Added inverse emulators

Sometimes, plain old regression is required. So I added inverse
emulators (in goes reflectance, out comes parameters).

gp_emulator/__init__.py

@@ -1,7 +1,15 @@
+#!/usr/bin/env python
+
+__author__ = "J Gomez-Dans"
+__copyright__ = "Copyright 2017, 2018 J Gomez-Dans"
+__version__ = "1.6.5"
+__license__ = "GPLv3"
+__email__ = "[email protected]"
 
 from .GaussianProcess import GaussianProcess, k_fold_cross_validation
 from .multivariate_gp import MultivariateEmulator
 from .lhd import lhd
 from .emulation_helpers import create_training_set, create_validation_set
 from .emulation_helpers import create_emulator_validation
 from .emulation_helpers import create_single_band_emulators
+from .emulation_helpers import create_inverse_emulators

gp_emulator/emulation_helpers.py

@@ -7,6 +7,34 @@
 from .multivariate_gp import MultivariateEmulator
 
 
+def integrate_passbands(spectrum, band_pass):
+    """Integrate a spectrum over some band pass function.
+    Parameters
+    ----------
+    spectrum: iter
+        An array of spectra (n_spectra * n_wavelengths)
+    band_pass: iter
+        A 2D array of size `(n_bands x n_wavelengths)`. The sum over the
+        second dimension of the array should evaluate to 1.
+        
+    Returns
+    -------
+    An array of n_spectra * n_bands with the integrated spectra.
+    """
+
+    n_bands = band_pass.shape[0]
+    if band_pass.dtype == np.bool:
+        x_pband = [spectrum[:, band_pass[i, :]].mean(axis=1)
+                         for i in range(n_bands)]
+    else:
+        band_pass = band_pass / (band_pass.sum(axis=1)[:, None])
+        x_pband = [np.sum(spectrum[:, :] * band_pass[i, :],
+                                axis=1) for i in range(n_bands)]
+
+    x_train_pband = np.array(x_pband)
+    return x_train_pband
+
+
 def create_single_band_emulators(emulator, band_pass, n_tries=15):
     """This function creates per band emulators from the full-spectrum
     emulator. It requires passing an array of band pass functions of
@@ -32,25 +60,50 @@ def create_single_band_emulators(emulator, band_pass, n_tries=15):
     A list of `n_bands` `GaussianProcess` emulators.
     """
 
-    n_bands = band_pass.shape[0]
-    if band_pass.dtype == np.bool:
-        x_train_pband = [emulator.X_train[:, band_pass[i, :]].mean(axis=1)
-                         for i in range(n_bands)]
-    else:
-        band_pass = band_pass / (band_pass.sum(axis=1)[:, None])
-        x_train_pband = [np.sum(emulator.X_train[:, :] * band_pass[i, :],
-                                axis=1) for i in range(n_bands)]
-
-    x_train_pband = np.array(x_train_pband)
+    x_train_pband = integrate_passbands(emulator.X_train, band_pass)
     emus = []
-    for i in range(n_bands):
+    for i in range(band_pass.shape[0]):
         gp = GaussianProcess(emulator.y_train[:] * 1.,
                              x_train_pband[i, :])
         gp.learn_hyperparameters(n_tries=n_tries)
         emus.append(gp)
     return emus
 
 
+def create_inverse_emulators (original_emulator, band_pass, state_config,
+                              n_tries=10):
+    """
+    This function takes a multivariable output trained emulator
+    and "retrains it" to take input reflectances and report a
+    prediction of single input parameters (i.e., regression from
+    reflectance/radiance to state). This is a useful starting 
+    point for spatial problems.
+    
+    Parameters
+    ------------
+    original_emulator: emulator
+        An emulator (type gp_emulator.GaussianProcess)
+    srf: array
+        A 2d bandpass array (nbands, nfreq). Logical type
+    state_config: ordered dict
+        A state configuration ordered dictionary.
+    """
+
+    # For simplicity, let's get the training data out of the emulator
+    X = original_emulator.X_train*1.
+    y = original_emulator.y_train*1.
+    # Apply band pass functions here...
+    x_train_pband = integrate_passbands(original_emulator.X_train, band_pass)
+
+    # A container to store the emulators
+    gps = {}
+    for  i, param in enumerate (state_config) :
+        gp = GaussianProcess ( x_train_pband.T, y[:, i] )
+        gp.learn_hyperparameters( n_tries = n_tries )
+        gps[param] = gp 
+    return gps
+
+
 def create_training_set(parameters, minvals, maxvals, fix_params=None, n_train=200):
     """Creates a traning set for a set of parameters specified by
     ``parameters`` (not actually used, but useful for debugging

setup.py

@@ -1,14 +1,31 @@
 #!/usr/bin/env python
+import codecs
+import os
+import re
 from os import path
 from setuptools import setup
 
 
 this_directory = path.abspath(path.dirname(__file__))
+
 with open(path.join(this_directory, 'README.md'), "rb") as f:
     long_description = f.read().decode('utf-8')
 
+def read(*parts):
+    with codecs.open(os.path.join(this_directory, *parts), 'r') as fp:
+        return fp.read()
+
+def find_version(*file_paths):
+    version_file = read(*file_paths)
+    version_match = re.search(r"^__version__ = ['\"]([^'\"]*)['\"]",
+                              version_file, re.M)
+    if version_match:
+        return version_match.group(1)
+    raise RuntimeError("Unable to find version string.")
+
+
 setup(name='gp_emulator',
-      version='1.6.4',
+      version=find_version("gp_emulator", "__init__.py"),
       description='A Python Gaussian Process emulator software package',
       classifiers=[
 	'Development Status :: 4 - Beta',