From 1868ca6cced04a5084bea6ab0f3fe23775f7091b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jo=C3=A3o=20Faria?= <joao.faria@astro.up.pt>
Date: Tue, 12 May 2020 17:23:31 +0100
Subject: [PATCH] a start to implement #62; requires changes to DNest4

---
 src/distributions/mixGaussianLogUniform.cpp | 48 +++++++++++++++++++++
 src/distributions/mixGaussianLogUniform.h   | 36 ++++++++++++++++
 2 files changed, 84 insertions(+)
 create mode 100644 src/distributions/mixGaussianLogUniform.cpp
 create mode 100644 src/distributions/mixGaussianLogUniform.h

diff --git a/src/distributions/mixGaussianLogUniform.cpp b/src/distributions/mixGaussianLogUniform.cpp
new file mode 100644
index 00000000..a4dde238
--- /dev/null
+++ b/src/distributions/mixGaussianLogUniform.cpp
@@ -0,0 +1,48 @@
+#include "mixGaussianLogUniform.h"
+
+namespace DNest4
+{
+
+/**
+ * Construct an mixGaussianLogUniform distribution
+*/
+mixGaussianLogUniform::mixGaussianLogUniform(double mean, double sigma, double lower, double upper)
+{
+    G = Gaussian(mean, sigma);
+    LU = LogUniform(lower, upper);
+}
+
+
+/* mixGaussianLogUniform distribution function */
+double mixGaussianLogUniform::cdf(double x) const
+{
+    return 0.5 * G.cdf(x) + 0.5 * LU.cdf(x);
+}
+
+/* Quantile estimate, midpoint interpolation */
+double mixGaussianLogUniform::cdf_inverse(double p) const
+{
+    // not implemented
+    throw std::runtime_error("cdf_inverse not implemented for `mixGaussianLogUniform`");
+}
+
+double mixGaussianLogUniform::log_pdf(double x) const
+{
+    return 0.5 * G.log_pdf(x) + 0.5 * LU.log_pdf(x);
+}
+
+double mixGaussianLogUniform::perturb(double& x, RNG& rng) const
+{
+    x = cdf(x);
+    x += rng.randh();
+    wrap(x, 0.0, 1.0);
+    if (rng.rand_int(2) == 0)
+        x = G.cdf_inverse(x);
+    else
+        x = LU.cdf_inverse(x);
+    return 0.0;
+}
+
+
+} // namespace DNest4
+
diff --git a/src/distributions/mixGaussianLogUniform.h b/src/distributions/mixGaussianLogUniform.h
new file mode 100644
index 00000000..90950e21
--- /dev/null
+++ b/src/distributions/mixGaussianLogUniform.h
@@ -0,0 +1,36 @@
+#pragma once
+
+// DNest4/code
+#include "DNest4.h"
+#include <limits>
+#include <stdlib.h>
+#include <stdint.h>
+#include <cmath>
+#include <iostream>
+#include <fstream>
+#include <vector> 
+#include <algorithm>
+
+using namespace std;
+
+namespace DNest4
+{
+
+class mixGaussianLogUniform:public ContinuousDistribution
+{
+    private:
+        DNest4::Gaussian G;
+        DNest4::LogUniform LU;
+
+    public:
+        mixGaussianLogUniform(double mean=0.0, double sigma=1.0, double lower=1.0, double upper=2.0);
+
+        double cdf(double x) const;
+        double cdf_inverse(double p) const;
+        double log_pdf(double x) const;
+        double perturb(double& x, RNG& rng) const override;
+};
+
+
+} // namespace DNest4
+