From 7b2afc7ec8f0662ecc802973c19d3c644cd95a0e Mon Sep 17 00:00:00 2001
From: gully <igully@gmail.com>
Date: Tue, 4 Apr 2017 10:51:42 -0700
Subject: [PATCH] Save Flux_scalars, the average value of each model spectrum,
 required for absolute flux of models and accurate solid angles.

Note that this is backwards-INCOMPATIBLE.  We have to pass an extra positional argument to many of the classes and methods that are not expecting one.  So trying to use this code with previous versions will fail.  This change is not desirable but the old code was misleading in the solid angle support, so it can be considered a bug fix.
---
 Starfish/emulator.py | 60 +++++++++++++++++++-------------------------
 1 file changed, 26 insertions(+), 34 deletions(-)

diff --git a/Starfish/emulator.py b/Starfish/emulator.py
index 36dd4dae..ec02260e 100644
--- a/Starfish/emulator.py
+++ b/Starfish/emulator.py
@@ -74,7 +74,7 @@ class PCAGrid:
     Create and query eigenspectra.
     '''
 
-    def __init__(self, wl, dv, flux_mean, flux_std, eigenspectra, w, w_hat, gparams):
+    def __init__(self, wl, dv, flux_mean, flux_std, eigenspectra, w, w_hat, gparams, flux_scalars):
         '''
 
         :param wl: wavelength array
@@ -91,6 +91,8 @@ def __init__(self, wl, dv, flux_mean, flux_std, eigenspectra, w, w_hat, gparams)
         :type w: 2D np.array (m, M)
         :param gparams: The stellar parameters of the synthetic library
         :type gparams: 2D array of parameters (nspec, nparam)
+        :param flux_scalar: scalar average of each spectrum
+        :type flux_scalar: 1D np.array (M)
 
         '''
         self.wl = wl
@@ -98,6 +100,7 @@ def __init__(self, wl, dv, flux_mean, flux_std, eigenspectra, w, w_hat, gparams)
         self.flux_mean = flux_mean
         self.flux_std = flux_std
         self.eigenspectra = eigenspectra
+        self.flux_scalars = flux_scalars
         self.m = len(self.eigenspectra)
         self.w = w
         self.w_hat = w_hat
@@ -136,7 +139,8 @@ def create(cls, interface):
 
         # Normalize all of the fluxes to an average value of 1
         # In order to remove uninteresting correlations
-        fluxes = fluxes/np.average(fluxes, axis=1)[np.newaxis].T
+        flux_scalars = np.average(fluxes, axis=1)
+        fluxes = fluxes/flux_scalars[np.newaxis].T
 
         # Subtract the mean from all of the fluxes.
         flux_mean = np.average(fluxes, axis=0)
@@ -184,7 +188,7 @@ def create(cls, interface):
         # Calculate w_hat, Eqn 20 Habib
         w_hat = get_w_hat(eigenspectra, fluxes, M)
 
-        return cls(wl, dv, flux_mean, flux_std, eigenspectra, w, w_hat, gparams)
+        return cls(wl, dv, flux_mean, flux_std, eigenspectra, w, w_hat, gparams,flux_scalars)
 
     def write(self, filename=Starfish.PCA["path"]):
         '''
@@ -208,6 +212,9 @@ def write(self, filename=Starfish.PCA["path"]):
         pdset[2,:] = self.flux_std
         pdset[3:, :] = self.eigenspectra
 
+        scalarset = hdf5.create_dataset("flux_scalars", (self.M,), compression='gzip', dtype="f8", compression_opts=9)
+        scalarset[:] = self.flux_scalars
+
         wdset = hdf5.create_dataset("w", (self.m, self.M), compression='gzip',
             dtype="f8", compression_opts=9)
         wdset[:] = self.w
@@ -244,13 +251,16 @@ def open(cls, filename=Starfish.PCA["path"]):
         wdset = hdf5["w"]
         w = wdset[:]
 
+        scalarset = hdf5["flux_scalars"]
+        flux_scalars = scalarset[:]
+
         w_hatdset = hdf5["w_hat"]
         w_hat = w_hatdset[:]
 
         gdset = hdf5["gparams"]
         gparams = gdset[:]
 
-        pcagrid = cls(wl, dv, flux_mean, flux_std, eigenspectra, w, w_hat, gparams)
+        pcagrid = cls(wl, dv, flux_mean, flux_std, eigenspectra, w, w_hat, gparams, flux_scalars)
         hdf5.close()
 
         return pcagrid
@@ -327,36 +337,6 @@ def reconstruct_all(self):
 
         return recon_fluxes
 
-
-class F_bol_interp:
-    '''
-    Interpolate the F_bol from the grid for any arbitrary stellar parameter
-    '''
-    def __init__(self, GridInterface):
-        '''
-        Provide the emulation products.
-
-        :param pca: object storing the principal components, the eigenpsectra
-        :type pca: PCAGrid
-        :param eparams: Optimized GP hyperparameters.
-        :type eparams: 1D np.array
-        '''
-        
-        myHDF5 = HDF5Interface()
-        n_gps, n_dims = myHDF5.grid_points.shape
-        F_bol = np.zeros(n_gps)
-        for i in range(n_gps):
-            _, hdr = myHDF5.load_flux_hdr(myHDF5.grid_points[i])
-            F_bol[i] = hdr["F_bol"]
-        X_grid = myHDF5.grid_points
-        y_dat = F_bol
-        ld = LinearNDInterpolator(X_grid, y_dat)
-        self.ld = ld
-
-    def interp(self, grid_pars):
-        return self.ld(grid_pars)
-
-
 class Emulator:
     def __init__(self, pca, eparams):
         '''
@@ -384,6 +364,8 @@ def __init__(self, pca, eparams):
 
         self.V11 = self.iPhiPhi + Sigma(self.pca.gparams, self.h2params)
 
+        self.ld = LinearNDInterpolator(self.pca.gparams, self.pca.flux_scalars)
+
         self._params = None # Where we want to interpolate
 
         self.V12 = None
@@ -391,6 +373,9 @@ def __init__(self, pca, eparams):
         self.mu = None
         self.sig = None
 
+        # New in v0.3: absolute flux scalar
+        self.flux_scalar = None
+
     @classmethod
     def open(cls, filename=Starfish.PCA["path"]):
         '''
@@ -437,10 +422,17 @@ def params(self, pars):
         self.mu.shape = (-1)
         self.sig = self.V22 - self.V12.T.dot(np.linalg.solve(self.V11, self.V12))
 
+        # Calculate the flux_scalar
+        self.flux_scalar = self.ld(self._params)
+
     @property
     def matrix(self):
         return (self.mu, self.sig)
 
+    @property
+    def absolute_flux(self):
+        return self.flux_scalar
+
     def draw_many_weights(self, params):
         '''
         :param params: multiple parameters to produce weight draws at.