large update feb 28, 2024

ECCCO_multiion_inversion.py

@@ -55,6 +55,7 @@
         #cube_file = response_dir + 'eccco_is_response_feldman_m_el_with_tables_lw_pm1230_'+str(psf)+'pix.fits'
         # cube_file = response_dir +'D16Feb2024_eccco_response_feldman_m_el_with_tables_s_i_slw_coopersun.fits'
         cube_file = response_dir + 'D16Feb2024_eccco_response_feldman_m_el_with_tables_s_i_lw_coopersun.fits'
+        #cube_file = response_dir + "response_img_normalized.fits"
         #cube_file = response_dir + 'D1Aug2023_eccco_response_feldman_m_el_with_tables_lw.fits'
         #cube_file = response_dir + 'D14Feb2024_eccco_response_feldman_m_el_with_tables_lw.fits'
 
@@ -64,7 +65,9 @@
         data_dir ='data/'
     #    summed_image  = data_dir + 'eccco_lw_forwardmodel_thermal_response_psf'+str(psf)+'pix_el_decon.fits'
         summed_image  = data_dir+'eccco_is_lw_forwardmodel_thermal_response_psf'+str(psf)+'pix_el.fits'
-        sample_weights_data  = data_dir +'eccco_is_lw_forwardmodel_sample_weights_psf'+str(psf)+'pix_el.fits'
+        #summed_img = data_dir+'forwardmodel_img_normalized.fits'
+        sample_weights_data = data_dir +'eccco_is_lw_forwardmodel_sample_weights_psf'+str(psf)+'pix_el.fits'
+        #sample_weights_data = data_dir + 'weights_img_normalized.fits'
         #summed_image  = data_dir+'eccco_lw_forwardmodel_thermal_response_psf'+str(psf)+'pix_el.fits'
         #sample_weights_data  = data_dir +'eccco_lw_forwardmodel_sample_weights_psf'+str(psf)+'pix_el.fits'
 
@@ -83,7 +86,7 @@
         # field_angle_range = None
 
         rsp_dep_name = 'logt'
-        rsp_dep_list = np.round((np.arange(57,78, 1) / 10.0), decimals=1)
+        rsp_dep_list = np.round((np.arange(57, 78, 1) / 10.0), decimals=1)
 
     #smooth_over = 'spatial'
         smooth_over = 'dependence'
@@ -92,8 +95,6 @@
                           rsp_dep_name=rsp_dep_name, rsp_dep_list=rsp_dep_list,
                           solution_fov_width=solution_fov_width,smooth_over=smooth_over,field_angle_range=field_angle_range)
 
-
-
         #inversion.initialize_input_data(summed_image)#,image_mask_file)
         #sample_weights_data = calculate_weights(summed_image, weight_file, 8., 1.)
         #print(sample_weights_data)

configs/full_aia_data.toml

@@ -0,0 +1,22 @@
+[paths]
+response = "data/ECCCO_speedtest_runs/D27Feb2024_eccco_response_feldman_m_el_with_tables_sw_lw_s_i_scaled.fits"
+weights = "data/ECCCO_speedtest_runs/combined_ECCCO_weights_sw_lw_s_i_scaled.fits"
+image = "data/ECCCO_speedtest_runs/combined_ECCCO_sw_lw_s_i_scaled.fits"
+output = "output/full_data_fewer_temps/"
+
+[settings]
+solution_fov_width = 2
+detector_row_range = [0, 792]
+field_angle_range = [-1227, 1227]
+response_dependency_name = "logt"
+response_dependency_list =[5.7, 5.8, 5.9,
+    6.0 , 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8]
+   #, 6.9, 7.0 , 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7]
+smooth_over = 'dependence'
+alphas = [0.01]
+rhos = [0.1]
+
+# alpha = 0.2 -> 162 seconds
+# alpha = 0.01 -> 760 seconds
+# alpha = 0.02 fewer temp bins -> 121 seconds
+# alpha = 0.01 fewer temp bins -> 587 seconds

configs/img_norm.toml

@@ -0,0 +1,42 @@
+[paths]
+response = "data/response_img_normalized.fits"
+weights = "data/weights_img_normalized.fits"
+image = "data/forwardmodel_img_normalized.fits"
+output = "output/img_norm/"
+
+[settings]
+solution_fov_width = 2
+detector_row_range = [450, 900]
+field_angle_range = [-2160, 2160]
+response_dependency_name = "logt"
+response_dependency_list =[5.7, 5.8, 5.9,
+    6.0 , 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9,
+    7.0 , 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7]
+smooth_over = 'dependence'
+alphas = [6.02]
+rhos = [1.0]
+
+# Inversion Time = 123.58615684509277 # threaded execution
+# Inversion Time = 332.0509181022644  # default number of processes in process pool
+# Inversion Time = 440.9607219696045  # 4 processes in process pool (seems to go around 80% CPU with some fluctuation)
+#### LIMIT NUMPY TO 1 THREAD
+# Inversion Time = 263.8916232585907  # 4 processes with numpy limited to 1 thread
+# Inversion Time = 81.74327301979065  # threaded execution with numpy limited to 1 thread
+# Inversion Time = 148.79005193710327 # 11 processes with numpy limited to 1 thread (should match the threaded number of threads)
+# Inversion Time = 152.01935005187988 # same as above with a thread pool
+# Inversion Time = 147.7941930294037  # my thread pool executor with inits with 11 max workers
+# Inversion Time = 144.82306694984436 # my process pool executor with inits with 11 max workers
+# Inversion Time = 145.85570001602173 # same as above with max_iter decreased by factor of 10
+# Inversion Time = 75.56018900871277  # my no shared memory thread pool executor with 11 max workers
+# Inversion Time = 77.73930597305298  # my process pool executor with copy on write for global with 11 max workers
+# Inversion Time = 88.20389604568481  # same as above with numpy having 2 threads
+# Inversion Time = 42.56954908370972  # same as above wtih numpy having 1 thread and selection='cyclic'
+
+# changing alpha can dramatically slow things
+# set alpha = 0.1 then Inversion Time = 1208.5097889900208
+# set alpha = 0.1 and selection = 'cyclic' Inversion Time = 464.01495003700256
+
+# warm start test
+# 50 rows:
+# off = Inversion Time = 25.761871099472046
+# on  = Inversion Time = 22.579617023468018

original_config.toml → configs/original_config.toml

@@ -4,11 +4,11 @@
 response = "data/D16Feb2024_eccco_response_feldman_m_el_with_tables_s_i_lw_coopersun.fits"
 weights = "data/eccco_is_lw_forwardmodel_sample_weights_psf4pix_el.fits"
 image = "data/eccco_is_lw_forwardmodel_thermal_response_psf4pix_el.fits"
-output = "output/"
+output = "output/original/"
 
 [settings]
 solution_fov_width = 2
-detector_row_range = [450, 455]
+detector_row_range = [450, 900]
 field_angle_range = [-2160, 2160]
 response_dependency_name = "logt"
 response_dependency_list =[5.7, 5.8, 5.9,

configs/spectra_only.toml

@@ -0,0 +1,20 @@
+# chop of the image to see the results of the spectra only
+[paths]
+response = "data/response_only_spectra.fits"
+weights = "data/weights_only_spectra.fits"
+image = "data/forward_model_only_spectra.fits"
+output = "output/only_spectra/"
+
+[settings]
+solution_fov_width = 2
+detector_row_range = [450, 900]
+field_angle_range = [-2160, 2160]
+response_dependency_name = "logt"
+response_dependency_list =[5.7, 5.8, 5.9,
+    6.0 , 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9,
+    7.0 , 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7]
+smooth_over = 'dependence'
+alphas = [1.01] #[5]
+rhos = [0.501] #[0.1]
+
+# 68 seconds

configs/trimmed_full_aia_data.toml

@@ -0,0 +1,20 @@
+[paths]
+response = "data/ECCCO_speedtest_runs/D27Feb2024_eccco_response_feldman_m_el_with_tables_trim_sw_lw_s_i_scaled.fits"
+weights = "data/ECCCO_speedtest_runs/combined_ECCCO_weights_trim_sw_lw_s_i_scaled.fits"
+image = "data/ECCCO_speedtest_runs/combined_ECCCO_trim_sw_lw_s_i_scaled.fits"
+output = "output/trimmed_full_data/"
+
+[settings]
+solution_fov_width = 2
+detector_row_range = [0, 792]
+field_angle_range = [-1227, 1227]
+response_dependency_name = "logt"
+response_dependency_list =[5.7, 5.8, 5.9,
+    6.0 , 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9,
+    7.0 , 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7]
+smooth_over = 'dependence'
+alphas = [0.01]
+rhos = [0.1]
+
+# alpha = 0.2 ->
+# alpha = 0.01 -> 690.062417268753

img_norm.toml → configs/uniform_weights.toml

@@ -1,10 +1,8 @@
-# version with the image normalized to look like the spectra
-
 [paths]
-response = "data/response_img_normalized.fits"
-weights = "data/forwardmodel_img_normalized.fits"
-image = "data/weights_img_normalized.fits"
-output = "output/"
+response = "data/D16Feb2024_eccco_response_feldman_m_el_with_tables_s_i_lw_coopersun.fits"
+weights = "data/weights_uniform.fits"
+image = "data/eccco_is_lw_forwardmodel_thermal_response_psf4pix_el.fits"
+output = "output/uniform_weights/"
 
 [settings]
 solution_fov_width = 2
@@ -18,4 +16,4 @@ smooth_over = 'dependence'
 alphas = [5]
 rhos = [0.1]
 
-# Inversion Time = 13.962096929550171
+# Inversion Time = 76.74999809265137

create_spectrally_pure_images.py

@@ -0,0 +1,12 @@
+import numpy as np
+
+from magixs_data_products import MaGIXSDataProducts
+
+dir_path = "output/full_data_fewer_temps/"
+
+mdp = MaGIXSDataProducts()
+image_list = ["output/full_data_fewer_temps/combined_ECCCO_sw_lw_s_i_scaled_em_data_cube_x2_1.0_0.01_wpsf.fits"]
+gnt_file = "data/ECCCO_speedtest_runs/master_gnt_eccco_inelectrons_cm3perspersr_with_tables.fits"
+
+rsp_dep_list = np.round((np.arange(56, 68, 1) / 10.0), decimals=1)
+mdp.create_level2_0_spectrally_pure_images(image_list, gnt_file, rsp_dep_list, dir_path)

experiment_scripts/restart_after_convergence_warm.py

@@ -0,0 +1,70 @@
+import time
+from multiprocessing import RawArray
+
+import numpy as np
+import toml
+from astropy.io import fits
+from sklearn.linear_model import ElasticNet
+
+from overlappogram.inversion_field_angles import Inversion
+
+ALPHA = 0.1
+
+if __name__ == '__main__':
+    with open("img_norm.toml") as f:
+        config = toml.load(f)
+
+    response_cube = fits.getdata(config['paths']['response'])
+    inversion = Inversion(rsp_func_cube_file=config['paths']['response'],
+                          rsp_dep_name=config['settings']['response_dependency_name'],
+                          rsp_dep_list=config['settings']['response_dependency_list'],
+                          solution_fov_width=config['settings']['solution_fov_width'],
+                          smooth_over=config['settings']['smooth_over'],
+                          field_angle_range=config['settings']['field_angle_range'])
+    inversion.initialize_input_data(config['paths']['image'],
+                                    None,
+                                    config['paths']['weights'])
+
+    # fits.writeto("response_matrix.fits", inversion.get_response_function())
+    X_d = inversion.get_response_function()
+    X_shape = inversion.get_response_function().shape
+    X = RawArray('d', X_shape[0] * X_shape[1])
+    X_np = np.frombuffer(X).reshape(X_shape)
+    np.copyto(X_np, inversion.get_response_function())
+
+    overlappogram_d = fits.getdata(config['paths']['image'])
+    overlappogram_shape = overlappogram_d.shape
+    overlappogram = RawArray('d', overlappogram_shape[0] * overlappogram_shape[1])
+    overlappogram_np = np.frombuffer(overlappogram).reshape(overlappogram_shape)
+    np.copyto(overlappogram_np, overlappogram_d)
+
+    weights_d = fits.getdata(config['paths']['weights'])
+    weights_shape = weights_d.shape
+    weights = RawArray('d', weights_shape[0] * weights_shape[1])
+    weights_np = np.frombuffer(weights).reshape(weights_shape)
+    np.copyto(weights_np, weights_d)
+
+    enet_model = ElasticNet(alpha=ALPHA,
+                            l1_ratio=0.1,
+                            max_iter=10_000,
+                            precompute=False,
+                            positive=True,
+                            copy_X=False,
+                            fit_intercept=False,
+                            warm_start=True,
+                            selection='cyclic')
+
+    i = 700
+    start = time.time()
+    enet_model.fit(X_d, overlappogram_d[i, :], sample_weight=weights_d[i, :])
+    data_out = enet_model.predict(X_d)
+    em = enet_model.coef_
+    end = time.time()
+    print("Inversion Time =", end - start)
+
+    # start = time.time()
+    # enet_model.fit(X_d, overlappogram_d[i, :], sample_weight=weights_d[i, :])
+    # data_out = enet_model.predict(X_d)
+    # em = enet_model.coef_
+    # end = time.time()
+    # print("Reconvergence Time =", end - start)

experiment_scripts/trim_imager.py

@@ -0,0 +1,20 @@
+# remove the imager from the response, weights, and overlappogram
+from astropy.io import fits
+
+response = "data/D16Feb2024_eccco_response_feldman_m_el_with_tables_s_i_lw_coopersun.fits"
+weights = "data/eccco_is_lw_forwardmodel_sample_weights_psf4pix_el.fits"
+image = "data/eccco_is_lw_forwardmodel_thermal_response_psf4pix_el.fits"
+
+
+response_hdul = fits.open(response)
+image_hdul = fits.open(image)
+weights_hdul = fits.open(weights)
+
+response_hdul_img_norm = response_hdul.copy()
+response_hdul_img_norm[0].data = response_hdul[0].data[:, :, :4096]
+response_hdul_img_norm.writeto("data/response_only_spectra.fits", overwrite=True)
+fits.writeto("data/forward_model_only_spectra.fits",
+             image_hdul[0].data[:, :4096],
+             header=image_hdul[0].header, overwrite=True)
+fits.writeto("data/weights_only_spectra.fits", weights_hdul[0].data[:, :4096],
+             header=weights_hdul[0].header, overwrite=True)

overlappogram/elasticnet_model.py

@@ -11,7 +11,7 @@ class ElasticNetModel(AbstractModel):
 
     def invert(self, response_function, data, sample_weights = None):
         #print(sample_weights)
-        self.model.fit(response_function, data, sample_weight=sample_weights)
+        self.model.fit(response_function, data, sample_weight=sample_weights, check_input=True)
         #self.model.fit(response_function, data)
         #score=(self.model.score(response_function, data, sample_weight=sample_weights))
         data_out = self.model.predict(response_function)
@@ -20,12 +20,12 @@ def invert(self, response_function, data, sample_weights = None):
         #return em, data_out, score
 
     def add_fits_keywords(self, header):
-        params = self.model.get_params()
-        #print(params)
-        header['INVMDL'] = ('Elastic Net', 'Inversion Model')
-        header['ALPHA'] = (params['alpha'], 'Inversion Model Alpha')
-        header['RHO'] = (params['l1_ratio'], 'Inversion Model Rho')
+        # params = self.model.get_params()
+        # #print(params)
+        # header['INVMDL'] = ('Elastic Net', 'Inversion Model')
+        # header['ALPHA'] = (params['alpha'], 'Inversion Model Alpha')
+        # header['RHO'] = (params['l1_ratio'], 'Inversion Model Rho')
+        pass
 
     def get_score(self, response_function, data):
-        score = self.model.score(response_function, data)
-        return score
+        return self.model.score(response_function, data)

overlappogram/ridge_model.py

@@ -27,5 +27,4 @@ def add_fits_keywords(self, header):
         # header['RHO'] = (params['l1_ratio'], 'Inversion Model Rho')
 
     def get_score(self, response_function, data):
-        score = self.model.score(response_function, data)
-        return score
+        return self.model.score(response_function, data)

overlappogram/sgd_model.py

@@ -21,10 +21,10 @@ def invert(self, response_function, data, sample_weights=None):
 
     def add_fits_keywords(self, header):
         params = self.model.get_params()
-        # print(params)
-        header['INVMDL'] = ('Elastic Net', 'Inversion Model')
-        header['ALPHA'] = (params['alpha'], 'Inversion Model Alpha')
-        header['RHO'] = (params['l1_ratio'], 'Inversion Model Rho')
+        # # print(params)
+        # header['INVMDL'] = ('Elastic Net', 'Inversion Model')
+        # header['ALPHA'] = (params['alpha'], 'Inversion Model Alpha')
+        # header['RHO'] = (params['l1_ratio'], 'Inversion Model Rho')
 
     def get_score(self, response_function, data):
         return self.model.score(response_function, data)

run_multiion_inversion.py

@@ -20,28 +20,34 @@ def run_inversion(config: dict):
     inversion.initialize_input_data(config['paths']['image'],
                                     None,
                                     config['paths']['weights'])
+    # inversion.sample_weights = None
 
     for alpha in config['settings']['alphas']:
         for rho in config['settings']['rhos']:
             enet_model = ElasticNet(alpha=alpha,
                                     l1_ratio=rho,
-                                    max_iter=100000,
-                                    precompute=True,
+                                    tol=1E-4,
+                                    max_iter=10_000,
+                                    precompute=False,
                                     positive=True,
+                                    copy_X=False,
                                     fit_intercept=False,
-                                    selection='cyclic')
+                                    selection='cyclic',
+                                    warm_start=False)
             inv_model = model(enet_model)
+
             basename = os.path.splitext(os.path.basename(config['paths']['image']))[0]
 
             start = time.time()
 
             postfix = 'x'+str(config['settings']['solution_fov_width'])+'_'+str(rho*10)+'_'+str(alpha)+'_wpsf'
             inversion.multiprocessing_invert(inv_model,
+            # inversion.invert(inv_model,
                                              config['paths']['output'],
                                              output_file_prefix=basename,
                                              output_file_postfix=postfix,
                                              detector_row_range=config['settings']['detector_row_range'],
-                                             score=False)
+                                             score=True)
 
             end = time.time()
             print("Inversion Time =", end - start)