The atlantics

.gitignore

@@ -0,0 +1,5 @@
+data
+data/*
+plots/*
+wl_*.npz
+*.pkl

ComplexSOM.yaml

@@ -0,0 +1,57 @@
+metrics:  [SNR_3x2]
+bands: riz
+training_file: data_buzzard/training.hdf5
+validation_file: data_buzzard/validation.hdf5
+output_file: example/ComplexSOM_output.txt
+# Backend implementing the metrics, either: "firecrown" (default), "jax-cosmo"
+metrics_impl: firecrown
+
+run:
+  # This is a class name which will be looked up
+  ComplexSOM:
+    run_5:
+      # These settings are sent to the classifier
+      testing: False
+      bins: 5
+      som_dim: [51,51]
+      num_groups: 200
+      num_threads: 32
+      group_type: 'redshift'
+      whiten: False
+      redshift_propset: 2
+      data_threshold: [0,30]
+      sparse_frac: 1.0
+      plots: True 
+      plot_dir: 'plots' 
+      metric: "SNR_3x2"
+      use_inbin: False
+      use_outbin: False
+      # These special settings decide whether the
+      # color and error colums are passed to the classifier
+      # as well as the magnitudes
+      colors: False
+      errors: False
+      sizes: True
+    run_10:
+      # These settings are sent to the classifier
+      testing: False
+      bins: 10
+      som_dim: [51,51]
+      num_groups: 200
+      num_threads: 32
+      group_type: 'redshift'
+      redshift_propset: 2
+      data_threshold: [0,30]
+      sparse_frac: 1.0
+      plots: True 
+      plot_dir: 'plots' 
+      metric: "SNR_3x2"
+      use_inbin: False
+      use_outbin: False
+      # These special settings decide whether the
+      # color and error colums are passed to the classifier
+      # as well as the magnitudes
+      colors: False
+      errors: False
+      sizes: True
+

bin/challenge.py

@@ -27,12 +27,15 @@ def main(config_yaml):
     # Decide if anyone needs the colors calculating and/or errors loading
     anyone_wants_colors = False
     anyone_wants_errors = False
+    anyone_wants_sizes = False
     for run in config['run'].values():
         for version in run.values():
             if version.get('errors'):
                 anyone_wants_errors = True
             if version.get('colors'):
                 anyone_wants_colors = True
+            if version.get('sizes'):
+                anyone_wants_sizes = True
 
 
     bands = config['bands']
@@ -41,14 +44,16 @@ def main(config_yaml):
         config['training_file'],
         bands,
         errors=anyone_wants_errors,
-        colors=anyone_wants_colors
+        colors=anyone_wants_colors,
+        size=anyone_wants_sizes
     )
 
     validation_data = tc.load_data(
         config['validation_file'],
         bands,
         errors=anyone_wants_errors,
-        colors=anyone_wants_colors
+        colors=anyone_wants_colors,
+        size=anyone_wants_sizes
     )
 
     training_z = tc.load_redshift(config['training_file'])
@@ -59,7 +64,7 @@ def main(config_yaml):
     else:
         metrics_fn = tc.compute_scores
 
-    with open(config['output_file'],'w') as output_file:
+    with open(config['output_file'],'a') as output_file:
         for classifier_name, runs in config['run'].items():
             for run, settings in runs.items():
                 scores = run_one(classifier_name, bands, settings,
@@ -77,15 +82,16 @@ def run_one(classifier_name, bands, settings, train_data, train_z, valid_data,
     if classifier.wants_arrays:
         errors = settings.get('errors')
         colors = settings.get('colors')
-        train_data = tc.dict_to_array(train_data, bands, errors=errors, colors=colors)
-        valid_data = tc.dict_to_array(valid_data, bands, errors=errors, colors=colors)
+        sizes = settings.get('sizes')
+        train_data = tc.dict_to_array(train_data, bands, errors=errors, colors=colors, sizes=sizes)
+        valid_data = tc.dict_to_array(valid_data, bands, errors=errors, colors=colors, sizes=sizes)
 
     print ("Executing: ", classifier_name, bands, settings)
 
     ## first check if options are valid
     print (settings, classifier.valid_options)
     for key in settings.keys():
-        if key not in classifier.valid_options and key not in ['errors', 'colors']:
+        if key not in classifier.valid_options and key not in ['errors', 'colors', 'sizes']:
             raise ValueError(f"Key {key} is not recognized by classifier {classifier_name}")
 
     print ("Initializing classifier...")
@@ -102,7 +108,7 @@ def run_one(classifier_name, bands, settings, train_data, train_z, valid_data,
 
     print ("Making some pretty plots...")
     name = str(classifier.__name__)
-    tc.metrics.plot_distributions(valid_z, results, f"plots/{name}_{settings}_{bands}.png")
+    tc.metrics.plot_distributions(valid_z, results, f"plots/{name}_{bands}.png")
 
     return scores
 

example/SimpleSOM.yaml

@@ -0,0 +1,28 @@
+metrics:  [SNR_3x2]
+bands: riz
+training_file: data/training.hdf5
+validation_file: data/validation_small.hdf5
+output_file: example/example_output.txt
+# Backend implementing the metrics, either: "firecrown" (default), "jax-cosmo"
+metrics_impl: firecrown
+
+run:
+  # This is a class name which will be looked up
+  SimpleSOM:
+    run_3:
+      # These settings are sent to the classifier
+      bins: 10
+      som_dim: [151,151]
+      num_groups: 200
+      num_threads: 32
+      group_type: 'redshift'
+      data_threshold: [0,30]
+      sparse_frac: 0.05
+      plots: True 
+      plot_dir: 'plots' 
+      # These special settings decide whether the
+      # color and error colums are passed to the classifier
+      # as well as the magnitudes
+      colors: False
+      errors: False
+

the_atlantics/example.py

@@ -0,0 +1,74 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from snc import SnCalc
+
+
+# Let's start with a "large" set of 12 samples
+n_bins = 12
+z_cent = np.linspace(0.2, 1.4, n_bins+1)+0.5/n_bins
+sigma_z = 0.05 * (1+z_cent)
+z_arr = np.linspace(0, 2, 1024)
+# Total number of galaxies in each sample
+n_gals = (z_cent/0.5)**2*np.exp(-(z_cent/0.5)**1.5) * 10. * 1.8E4 * 60**2
+
+
+def pz_gau(z, zc, sz):
+    return np.exp(-0.5*((z-zc)/sz)**2)/np.sqrt(2*np.pi*sz**2)
+
+
+# Create the N(z) arrays, which should integrate to the total number
+# of galaxies in each sample
+nzs = [ng * pz_gau(z_arr, zc, sz)
+       for ng, zc, sz in zip(n_gals, z_cent, sigma_z)]
+
+# Let's plot them
+plt.figure()
+for nz in nzs:
+    plt.plot(z_arr, nz)
+plt.xlabel(r'$z$', fontsize=16)
+plt.ylabel(r'$N(z)$', fontsize=16)
+
+
+c = SnCalc(z_arr, nzs)
+# Now let's do a dumb thing and bunch them up into
+# 4, 3, 2 and 1 subsamples
+assign4 = [np.array([0, 1, 2]),
+           np.array([3, 4, 5]),
+           np.array([6, 7, 8]),
+           np.array([9, 10, 11])]
+assign3 = [np.array([0, 1, 2, 3]),
+           np.array([4, 5, 6, 7]),
+           np.array([8, 9, 10, 11])]
+assign2 = [np.array([0, 1, 2, 3, 4, 5]),
+           np.array([6, 7, 8, 9, 10, 11])]
+assign1 = [np.arange(n_bins)]
+
+
+# OK, caclulate S/N for all of these groupings
+snr_d4 = c.get_sn(assign4, full_output=True)
+print('4 bins, S/N = %.1lf' % (snr_d4['snr']))
+snr_d3 = c.get_sn(assign3, full_output=True)
+print('3 bins, S/N = %.1lf' % (snr_d3['snr']))
+snr_d2 = c.get_sn(assign2, full_output=True)
+print('2 bins, S/N = %.1lf' % (snr_d2['snr']))
+snr_d1 = c.get_sn(assign1, full_output=True)
+print('1 bins, S/N = %.1lf' % (snr_d1['snr']))
+
+
+# Note that the calculator also returns the power
+# spectra, if you wanna look at them.
+plt.figure()
+for i in range(4):
+    for j in range(i, 4):
+        if i == j:
+            plt.plot(snr_d4['ls'],
+                     snr_d4['cl'][:, i, j],
+                     'k-')
+            plt.plot(snr_d4['ls'],
+                     snr_d4['nl'][0, i, j]*np.ones_like(snr_d4['ls']), 'k--')
+        else:
+            plt.plot(snr_d4['ls'],
+                     snr_d4['cl'][:, i, j],
+                     '-', c='#AAAAAA')
+plt.loglog()
+plt.show()

the_atlantics/example_minimizer.py

@@ -0,0 +1,86 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from snc import SnCalc
+from scipy.optimize import minimize
+
+
+# Let's start with a "large" set of 12 samples
+n_bins = 100
+z_cent = np.linspace(0.2, 1.4, n_bins+1)+0.5/n_bins
+sigma_z = 0.05 * (1+z_cent)
+z_arr = np.linspace(0, 2, 1024)
+# Total number of galaxies in each sample
+n_gals = (z_cent/0.5)**2*np.exp(-(z_cent/0.5)**1.5) * 10. * 1.8E4 * 60**2
+
+
+def pz_gau(z, zc, sz):
+    return np.exp(-0.5*((z-zc)/sz)**2)/np.sqrt(2*np.pi*sz**2)
+
+
+# Create the N(z) arrays, which should integrate to the total number
+# of galaxies in each sample.
+nzs = [ng * pz_gau(z_arr, zc, sz)
+       for ng, zc, sz in zip(n_gals, z_cent, sigma_z)]
+
+# Let's plot them
+plt.figure()
+for nz in nzs:
+    plt.plot(z_arr, nz)
+plt.xlabel(r'$z$', fontsize=16)
+plt.ylabel(r'$N(z)$', fontsize=16)
+
+# Now initialize a S/N calculator for these initial groups.
+c_wl = SnCalc(z_arr, nzs, use_clustering=False)
+c_gc = SnCalc(z_arr, nzs, use_clustering=True)
+print("Initializing WL")
+c_wl.get_cl_matrix(fname_save='wl_nb%d.npz' % n_bins)
+print("Initializing GC")
+c_gc.get_cl_matrix(fname_save='gc_nb%d.npz' % n_bins)
+
+# Let's focus on WL for now.
+# OK, let's brute-force explore the total S/N as a function
+# of the bin edge for a 2-bin scenario.
+print("Exploring 2-bins")
+zs = np.linspace(0.01, 1.5, 128)
+sns_wl = np.array([c_wl.get_sn_from_edges(np.array([z])) for z in zs])
+sns_gc = np.array([c_gc.get_sn_from_edges(np.array([z])) for z in zs])
+plt.figure()
+plt.plot(zs, sns_wl/np.amax(sns_wl), 'r-', label='Lensing')
+plt.plot(zs, sns_gc/np.amax(sns_gc), 'k-', label='Clustering')
+plt.xlabel('$z$')
+plt.ylabel(r'$(S/N)/(S/N)_{\rm max}$')
+plt.legend(loc='upper left', frameon=False)
+
+# Now let's do a 3-bin case (more coarsely-grained so we can do it quickly)
+print("Exploring 3-bins")
+zs = np.linspace(0.01, 1.5, 32)
+sns = np.array([[c_wl.get_sn_from_edges(np.array([z1, z2])) for z1 in zs]
+                for z2 in zs])
+
+plt.figure()
+plt.imshow(sns, vmin=1440, extent=[0.01, 1.5, 0.01, 1.5],
+           origin='lower')
+plt.xlabel('$z_1$')
+plt.ylabel('$z_2$')
+cb = plt.colorbar()
+cb.set_label('$S/N$')
+
+
+# Finally, let's write the function to minimize and optimize for a 4-bin case.
+def minus_sn(edges, calc):
+    return -calc.get_sn_from_edges(edges)
+
+
+print("Optimizing WL")
+edges_0 = np.array([0.3, 0.6, 0.9])
+res = minimize(minus_sn, edges_0, method='Powell', args=(c_wl,))
+print("WL final edges: ", res.x)
+print("Maximum S/N: ", c_wl.get_sn_from_edges(res.x))
+print(" ")
+
+# That was for weak lensing. Let's do the same thing for clustering.
+print("Optimizing GC")
+res = minimize(minus_sn, edges_0, method='Powell', args=(c_gc,))
+print("GC final edges: ", res.x)
+print("Maximum S/N: ", c_gc.get_sn_from_edges(res.x))
+plt.show()