Merge pull request #48 from dirac-institute/allow_y0_set

Allow y0 set

README.md

@@ -1,6 +1,6 @@
 # siggi 
-[![Build Status](https://travis-ci.org/jbkalmbach/siggi.svg?branch=master)](https://travis-ci.org/jbkalmbach/siggi)
-[![codecov.io](https://codecov.io/github/jbkalmbach/siggi/coverage.svg?branch=master)](https://codecov.io/github/jbkalmbach/siggi?branch=master)
+[![Build Status](https://travis-ci.org/dirac-institute/siggi.svg?branch=master)](https://travis-ci.org/dirac-institute/siggi)
+[![codecov.io](https://codecov.io/github/dirac-institute/siggi/coverage.svg?branch=master)](https://codecov.io/github/dirac-institute/siggi?branch=master)
 [![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](http://www.gnu.org/licenses/gpl-3.0)
 [![DOI](https://zenodo.org/badge/118051192.svg)](https://zenodo.org/badge/latestdoi/118051192)
 

examples/scripts/run_catsim_prior_lsst_plusOne.py

@@ -0,0 +1,105 @@
+## Example LSST+1 script
+# First import code
+import sys
+sys.path.append('..')
+import os
+import time
+from siggi import siggi, filters, spectra, Sed
+from siggi.lsst_utils import BandpassDict, Bandpass, PhotometricParameters
+import numpy as np
+
+def prior_z(z):
+    #normalization = 0.08333
+    a, b, z0 = np.array([ 0.7787937 ,  2.47523432,  1.47403178])
+    return (z**a)*np.exp(-(z/z0)**b)
+
+if __name__ == "__main__":
+
+    f = filters()
+    s = spectra()
+
+    sed_list = []
+    for sed_name in os.listdir('../data/cww_kin_lephare/'):
+        sed_obj = Sed()
+        sed_obj.readSED_flambda('../data/cww_kin_lephare/%s' % sed_name)
+        # Convert from angstroms to nm
+        sed_obj.wavelen /= 10.
+        sed_list.append(sed_obj)
+
+    num_filters = 1
+    total_non_zero = 0
+    num_trials = 510
+
+    bp_list = []
+    new_phot_params = {}
+    bp_dir = '/astro/store/epyc/users/brycek/siggi/siggi/data/lsst_baseline_throughputs'
+    for filter_name in ['u', 'g', 'r', 'i', 'z', 'y']:
+        current_bp = Bandpass()
+        print(os.path.join(bp_dir, 'filter_%s.dat' % filter_name))
+        current_bp.readThroughput(os.path.join(bp_dir, 'filter_%s.dat' % filter_name))
+        bp_list.append(current_bp)
+        new_phot_params[filter_name] = PhotometricParameters(bandpass=filter_name)
+
+    new_phot_params['filter_0'] = PhotometricParameters(nexp=160*2, bandpass='filter_6')
+    #bp_list.append(new_filt['filter_0'])                                                                                                                                                       
+    frozen_dict = BandpassDict(bp_list, ['u', 'g', 'r', 'i', 'z', 'y'])#, 'filter_0'])  
+
+    sed_weights = np.ones(len(sed_list))/len(sed_list)
+
+    sig_example = siggi(sed_list,
+                        sed_weights, prior_z, 
+                        z_min=0.00, z_max=2.3, z_steps=47, phot_params=new_phot_params)
+
+    x0 = None#[[300.,         414.99703239, 562.22744686, 630.98517889 ],
+          #[ 318.40555564,  408.17312354,  870.49778848, 1073.39207548] ]
+    y0 = None
+
+    start = time.time()
+    print('Starting at ', time.localtime())
+
+    rand_state = np.random.RandomState(1305)
+    ratio = None
+
+    res = sig_example.optimize_filters(num_filters=num_filters,
+                                       filt_min=300., filt_max=1100.,
+                                       set_ratio=ratio,
+                                       procs=2, n_opt_points=num_trials,
+                                       system_wavelen_min=300.,
+                                       system_wavelen_max=1200.,
+                                       starting_points=x0,
+                                       frozen_filt_dict=frozen_dict,
+                                       frozen_filt_eff_wavelen=[365., 477., 622., 765., 870., 1015],
+                                       #acq_func_kwargs_dict={'kappa':3.5},
+                                       optimizer_verbosity=10,
+                                       rand_state=rand_state,
+                                       save_optimizer='frozen_filter_opt.pkl')
+
+    trial_vals = np.array(res.yi)
+    trial_pts = np.array(res.Xi)
+    best_val = np.min(res.yi)
+    best_pt = trial_pts[np.argmin(res.yi)]
+    random_pts_used = res.random_pts_used
+
+    non_zero_pts = np.where(trial_vals != 0.)[0]
+    total_non_zero = len(non_zero_pts)
+
+    suffix = 'ff_10sed_catsim_long'
+
+    with open('results/run_results_%s.txt' % suffix, 'w') as f:
+        f.write('After %i trials: \n' % total_non_zero)
+        f.write('Random Points Used: %i \n' % random_pts_used)
+        f.write('Best point: ')
+        for pt_val in best_pt:
+            f.write('%.4f ' % pt_val)
+        f.write('\n')
+        f.write('Best Information Gain: %.4f' % (-1.*best_val))
+
+    np.savetxt('results/run_points_%s.txt' % suffix, trial_pts, fmt='%f')
+    np.savetxt('results/run_values_%s.txt' % suffix, -1.*np.array(trial_vals), fmt='%f')
+
+
+    finish = time.time()
+
+    print('Job finished in %.4f seconds.' % (finish-start))
+
+    print(best_pt, -1.*best_val)

examples/scripts/run_catsim_prior_ratio_0p1.py

@@ -0,0 +1,97 @@
+## Example 6 filter script
+# First import code
+import sys
+sys.path.append('..')
+import os
+import time
+from siggi import siggi, filters, spectra, Sed
+from siggi.lsst_utils import BandpassDict, Bandpass
+from scipy.stats import norm
+import numpy as np
+
+def prior_z(z):
+    a, b, z0 = np.array([ 0.7787937 ,  2.47523432,  1.47403178])
+    return (z**a)*np.exp(-(z/z0)**b)
+
+
+if __name__ == "__main__":
+
+    f = filters()
+    s = spectra()
+
+    sed_list = []
+    for sed_name in os.listdir('../data/cww_kin_lephare/'):
+        sed_obj = Sed()
+        sed_obj.readSED_flambda('../data/cww_kin_lephare/%s' % sed_name)
+        # Convert from angstroms to nm
+        sed_obj.wavelen /= 10.
+        sed_list.append(sed_obj)
+
+    sed_weights = np.ones(len(sed_list))/len(sed_list)
+
+    sig_example = siggi(sed_list,
+                        sed_weights, prior_z, 
+                        z_min=0.00, z_max=2.3, z_steps=47)
+
+    num_filters = 6
+    total_non_zero = 0
+    num_trials = 210
+
+    bp_list = []
+    bp_dir = '../siggi/data/lsst_baseline_throughputs'
+    for filter_name in ['u', 'g', 'r', 'i', 'z', 'y']:
+        current_bp = Bandpass()
+        print(os.path.join(bp_dir, 'filter_%s.dat' % filter_name))
+        current_bp.readThroughput(os.path.join(bp_dir, 'filter_%s.dat' % filter_name))
+        bp_list.append(current_bp)
+
+    x0 = None
+    y0 = None
+
+    start = time.time()
+    print('Starting at ', time.localtime())
+
+    rand_state = np.random.RandomState(1305)
+    ratio = 0.1
+
+    res = sig_example.optimize_filters(num_filters=num_filters,
+                                       filt_min=300., filt_max=1100.,
+                                       set_ratio=ratio,
+                                       procs=2, n_opt_points=num_trials,
+                                       system_wavelen_min=300.,
+                                       system_wavelen_max=1200.,
+                                       starting_points=x0,
+                                       optimizer_verbosity=10,
+                                       rand_state=rand_state,
+                                       max_search_factor=100,
+                                       save_optimizer='frozen_filter_opt_%i.pkl' % int(ratio*10))
+
+    trial_vals = np.array(res.yi)
+    trial_pts = np.array(res.Xi)
+    best_val = np.min(res.yi)
+    best_pt = trial_pts[np.argmin(res.yi)]
+    random_pts_used = res.random_pts_used
+
+    non_zero_pts = np.where(trial_vals != 0.)[0]
+    total_non_zero = len(non_zero_pts)
+
+    suffix = '6filter_catsim_%02i' % int(ratio*10)
+
+    with open('results/run_results_%s.txt' % suffix, 'w') as f:
+        f.write('After %i trials: \n' % total_non_zero)
+        f.write('Random Points Used: %i \n' % random_pts_used)
+        f.write('Best point: ')
+        for pt_val in best_pt:
+            f.write('%.4f ' % pt_val)
+        f.write('\n')
+        f.write('Best Information Gain: %.4f' % (-1.*best_val))
+
+    np.savetxt('results/run_points_%s.txt' % suffix, trial_pts, fmt='%f')
+    np.savetxt('results/run_values_%s.txt' % suffix, -1.*np.array(trial_vals), fmt='%f')
+
+
+    finish = time.time()
+
+    print('Job finished in %.4f seconds.' % (finish-start))
+
+    print(best_pt, -1.*best_val)

siggi/calcIG.py

@@ -1,5 +1,6 @@
 from __future__ import division
 
+import os
 import numpy as np
 from scipy import interpolate
 from scipy.spatial.distance import cdist
@@ -39,13 +40,18 @@ class calcIG(mathUtils):
     """
 
     def __init__(self, filter_dict, sed_list, y_probs, y_vals,
-                 n_pts=500000,
+                 n_pts=1000000,
                  sky_mag=20.47, ref_filter=None, phot_params=None,
                  fwhm_eff=0.8):
 
+        bp_dict_folder = os.path.join(os.path.dirname(__file__),
+                                      'data',
+                                      'lsst_baseline_throughputs')
+        bp_dict = BandpassDict.loadTotalBandpassesFromFiles(
+            bandpassDir=bp_dict_folder)
+
         if ref_filter is None:
-            ref_filter = Bandpass()
-            ref_filter.imsimBandpass()
+            ref_filter = bp_dict['i']
 
         self.sed_list = sed_list
         self.n_pts = n_pts

siggi/plotting.py

@@ -43,17 +43,17 @@ def __init__(self, sed_list, best_point,
             BandpassDict.addSystemBandpass(trap_dict)
 
         if frozen_filt_dict is None:
-            self.filter_dict = total_filt_dict
+            self.filter_dict = hardware_filt_dict#total_filt_dict
         else:
             if (type(frozen_filt_eff_wavelen) != list):
                 raise ValueError("If including frozen filters, " +
                                  "need list of eff. wavelengths.")
             filter_wavelengths = frozen_filt_eff_wavelen +\
                 self.find_filt_centers(filter_info)
             filter_names_unsort = frozen_filt_dict.keys() +\
-                total_filt_dict.keys()
+                hardware_filt_dict.keys()
             filter_list_unsort = frozen_filt_dict.values() +\
-                total_filt_dict.values()
+                hardware_filt_dict.values()
             sort_idx = np.argsort(filter_wavelengths)
             filter_names = [filter_names_unsort[idx] for idx in sort_idx]
             filter_list = [filter_list_unsort[idx] for idx in sort_idx]
@@ -295,7 +295,7 @@ def plot_ig_space(self, test_pts, test_vals, filter_idx,
 
         extent = [np.min(xi), np.max(xi), np.min(yi), np.max(yi)]
         plt.imshow(zi_lin, cmap=plt.cm.plasma, origin='lower',
-                   extent=extent, interpolation='bicubic')
+                   extent=extent, interpolation='bicubic', vmax=1.)
 
         if return_centers is True:
             return xi, yi

siggi/siggi.py

@@ -62,7 +62,7 @@ class siggi(_siggiBase):
 
     def __init__(self, spec_list, spec_weights, z_prior,
                  z_min=0.00, z_max=2.5, z_steps=51,
-                 calib_filter=None, calib_mag=25.,
+                 calib_filter=None, calib_mag=25.3,
                  phot_params=None):
 
         self.num_sed_types = len(spec_list)
@@ -104,6 +104,7 @@ def optimize_filters(self, filt_min=300., filt_max=1100.,
                          set_ratio=None,
                          width_min=30., width_max=120.,
                          starting_points=None,
+                         starting_vals=None,
                          system_wavelen_min=300., system_wavelen_max=1150.,
                          system_wavelen_step=0.1,
                          procs=1, n_opt_points=100, acq_func_kwargs_dict=None,
@@ -287,6 +288,21 @@ def optimize_filters(self, filt_min=300., filt_max=1100.,
         else:
             opt = load_optimizer
 
+        if starting_vals is not None:
+
+            len_y0 = len(starting_vals)
+            x_submit = []
+            for x_sample in x1[:len_y0]:
+                x_submit.append(list(x_sample))
+            opt.tell(x_submit, list(starting_vals))
+
+            if len_y0 < len(x1):
+                x0 = [list(x_on) for x_on in x1[len_y0:]]
+            else:
+                x0 = [list(x_on) for x_on in x0[len_y0:]]
+
+        print(x0)
+
         with Parallel(n_jobs=procs, batch_size=1, backend=parallel_backend,
                       verbose=self.verbosity) as parallel:
             while i < n_opt_points:

tests/test_calcIG.py

@@ -1,4 +1,5 @@
 import unittest
+import os
 import numpy as np
 from scipy.stats import norm, entropy, multivariate_normal
 from siggi import filters, spectra, calcIG
@@ -26,6 +27,12 @@ def setUpClass(cls):
         cls.imsimBand = Bandpass()
         cls.imsimBand.imsimBandpass()
 
+        cls.frozen_dict = BandpassDict.loadTotalBandpassesFromFiles(
+                        bandpassNames=['i'],
+                        bandpassDir=os.path.join(os.path.dirname(__file__),
+                                                 '../siggi/data',
+                                                 'lsst_baseline_throughputs'))
+
         phot_params = {}
         for idx in range(6):
             phot_params['filter_%i' % idx] = \
@@ -65,7 +72,7 @@ def test_calc_colors(self):
 
         np.testing.assert_equal(colors, np.zeros(np.shape(colors)))
 
-        sky_fn = self.sky_spec.calcFluxNorm(21.2, self.imsimBand)
+        sky_fn = self.sky_spec.calcFluxNorm(21.2, self.frozen_dict['i'])
         self.sky_spec.multiplyFluxNorm(sky_fn)
 
         test_error = calcMagError_sed(test_c.sed_list[0][0],

tests/test_siggi.py

@@ -71,8 +71,8 @@ def prior_z(z, z0=0.5):
         np.testing.assert_array_equal(t_1.Xi, t_2.Xi)
         np.testing.assert_array_equal(t_1.yi, t_2.yi)
         np.testing.assert_almost_equal(np.max(np.abs(t_1.yi[:10])),
-                                       1.63715332)
-        self.assertGreaterEqual(np.max(np.abs(t_1.yi)), 1.63715332)
+                                       1.2962127)
+        self.assertGreaterEqual(np.max(np.abs(t_1.yi)), 1.2962127)
 
         # Test pickling of optimization can replicate results