solver feature (#29)

* solver feature

* need scipy now

Author: mathematicalmichael

setup.cfg

@@ -30,9 +30,11 @@ package_dir =
 # DON'T CHANGE THE FOLLOWING LINE! IT WILL BE UPDATED BY PYSCAFFOLD!
 setup_requires = pyscaffold>=3.2a0,<3.3a0
 # Add here dependencies of your project (semicolon/line-separated), e.g.
-install_requires = numpy
-                   matplotlib
-                   pyerf
+install_requires =
+                numpy
+                scipy
+                matplotlib
+                pyerf
 
 # The usage of test_requires is discouraged, see `Dependency Management` docs
 # tests_require = pytest; pytest-cov

src/mud/__init__.py

@@ -9,3 +9,5 @@
     __version__ = 'unknown'
 finally:
     del get_distribution, DistributionNotFound
+
+from .base import DensityProblem  # noqa: F401

src/mud/base.py

@@ -0,0 +1,53 @@
+import numpy as np
+from scipy.stats import distributions as dist
+from scipy.stats import gaussian_kde as gkde
+
+
+class DensityProblem(object):
+    def __init__(self, X, y, domain=None):
+        self.X = X
+        self.y = y
+        self.domain = np.array(domain)
+        self._up = None
+        self._pr = None
+        self._in = None
+        self._ob = None
+
+    def set_observed(self, distribution=dist.norm()):
+        self._ob = distribution.pdf(self.y).prod(axis=1)
+
+    def set_initial(self, distribution=None):
+        if distribution is None:  # assume standard normal by default
+            if self.domain is not None:  # assume uniform if domain specified
+                mn = np.min(self.domain, axis=1)
+                mx = np.max(self.domain, axis=1)
+                distribution = dist.uniform(loc=mn, scale=mx - mn)
+            distribution = dist.norm()
+        initial_dist = distribution
+        self._in = initial_dist.pdf(self.X).prod(axis=1)
+
+    def set_predicted(self, distribution=None):
+        if distribution is None:
+            distribution = gkde(self.y.T)
+            pred_pdf = distribution.pdf(self.y.T).T
+        else:
+            pred_pdf = distribution.pdf(self.y)
+        self._pr = pred_pdf
+
+    def fit(self):
+        if not self._in:
+            self.set_initial()
+            self._pr = None
+        if not self._pr:
+            self.set_predicted()
+        if not self._ob:
+            self.set_observed()
+
+        up_pdf = np.divide(np.multiply(self._in, self._ob), self._pr)
+        self._up = up_pdf
+
+    def mud_point(self):
+        if self._up is None:
+            self.fit()
+        m = np.argmax(self._up)
+        return self.X[m, :]

src/mud/funs.py

@@ -10,6 +10,7 @@
 import numpy as np
 
 from mud import __version__
+from mud import DensityProblem
 
 __author__ = "Mathematical Michael"
 __copyright__ = "Mathematical Michael"
@@ -264,5 +265,34 @@ def iterate(A, b, y, initial_mean, initial_cov,
     return chain
 
 
+def mud_problem(lam, qoi, qoi_true, domain, sd=0.05, num_obs=None):
+    """
+    Wrapper around mud problem, takes in raw qoi + synthetic data and
+    performs WME transformation, instantiates solver object
+    """
+    if lam.ndim == 1:
+        lam = lam.reshape(-1, 1)
+
+    if qoi.ndim == 1:
+        qoi = qoi.reshape(-1, 1)
+    dim_output = qoi.shape[1]
+
+    if num_obs is None:
+        num_obs = dim_output
+    elif num_obs < 1:
+        raise ValueError("num_obs must be >= 1")
+    elif num_obs > dim_output:
+        raise ValueError("num_obs must be <= dim(qoi)")
+
+    # this is our data processing step.
+    data = qoi_true[0:num_obs] + np.random.randn(num_obs) * sd
+    q = wme(qoi[:, 0:num_obs], data, sd).reshape(-1, 1)
+
+    # this just implements density-based solutions + mud point method
+    d = DensityProblem(lam, q, domain)
+#     d.fit() # optional. will compute if invoked while empty.
+    return d
+
+
 if __name__ == "__main__":
     run()