Merge pull request #7 from jmbhughes/cleantime

First working version

.github/workflows/docs.yml

@@ -0,0 +1,36 @@
+name: deploy-docs
+
+# Only run this when the master branch changes
+on:
+  push:
+    branches:
+    - main
+
+# This job installs dependencies, builds the book, and pushes it to `gh-pages`
+jobs:
+  deploy-docs:
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v2
+
+    # Install dependencies
+    - name: Set up Python 3.10
+      uses: actions/setup-python@v2
+      with:
+        python-version: "3.10"
+
+    - name: Install dependencies
+      run: |
+        pip install ".[docs]"
+
+    # Build the book
+    - name: Sphinx build
+      run: |
+        cd ./docs; make html; cd ..
+
+    # Push the book's HTML to github-pages
+    - name: GitHub Pages action
+      uses: peaceiris/[email protected]
+      with:
+        github_token: ${{ secrets.GITHUB_TOKEN }}
+        publish_dir: ./docs/_build/html

.gitignore

@@ -163,3 +163,5 @@ cython_debug/
 # overlappogram data
 data/
 output/
+*.fits
+*/*.fits

README.md

@@ -1,3 +1,16 @@
 # overlappogram
 
 [![CI](https://github.com/jmbhughes/overlappogram/actions/workflows/CI.yml/badge.svg)](https://github.com/jmbhughes/overlappogram/actions/workflows/CI.yml)
+
+![overlappogram example](overlappogram.png)
+
+Overlappogram is a Python package for inverting overlappogram observations of the Sun, for examples MaGIXS observations
+or ECCCO observations.
+
+## How to Use
+
+## Getting Help
+
+## Cite
+
+## Code of Conduct

configs/full_aia_data.toml → configs/old/full_aia_data.toml

@@ -13,10 +13,12 @@ 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]
+alphas = [0.1, 0.005]
 rhos = [0.1]
 
 # alpha = 0.2 -> 162 seconds
 # alpha = 0.01 -> 760 seconds
-# alpha = 0.02 fewer temp bins -> 121 seconds
+# alpha = 0.2 fewer temp bins -> 121 seconds
+# alpha = 0.1 fewer temp bins -> 156 seconds
 # alpha = 0.01 fewer temp bins -> 587 seconds
+# alpha = 0.005 fewer temp bins ->  941 seconds

configs/old/precompute.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/photons/"
+
+[settings]
+solution_fov_width = 2
+detector_row_range = [300, 400]
+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]
+smooth_over = 'dependence'
+alphas = [0.2]
+rhos = [0.1]
+
+# alpha = 0.2 -> 46 seconds
+# alpha = 0.1 -> 63 seeconds
+# alpha = 0.01 -> 277 seconds
+# alpha = 0.005 -> 471 seconds

configs/trimmed_aia_data.toml

@@ -0,0 +1,24 @@
+[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/photons/"
+
+[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]
+smooth_over = 'dependence'
+alphas = [0.2, 0.1, 0.01, 0.005]
+rhos = [0.1]
+
+# alpha = 0.2 -> 46 seconds
+# alpha = 0.1 -> 63 seeconds
+# alpha = 0.01 -> 277 seconds ... with 30 threads = 229 seconds
+#               ... with 40 threads = 220 seconds ... with 50 threads = 217 seconds
+#               ... boost tol=1E-5 -> 333 seconds (no quality impact)
+#               ... drop tol=1E-3 -> 102 seconds (no quality impact)
+#               ... drop tol=1E-2 -> 58 seconds (starting to see deviation from earlier results... is it worse?)
+# alpha = 0.005 -> 471 seconds

example_config.toml

@@ -0,0 +1,28 @@
+[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"
+gnt = "data/ECCCO_speedtest_runs/master_gnt_eccco_inphotons_cm3persperpix_with_tables.fits"
+
+[output]
+prefix = ""
+make_spectral = true
+directory = "output/test/"
+
+[inversion]
+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]
+smooth_over = 'dependence'
+
+[model]
+alphas = [0.2] #[0.2, 0.1, 0.01, 0.005]
+rhos = [0.1]
+warm_start = true
+tol = 1E-2
+max_iter = 10_000
+selection = 'cyclic'
+
+[execution]
+num_threads = 50

experiment_scripts/compare_spectrally_pure.py

@@ -0,0 +1,67 @@
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+from astropy.io import fits
+from scipy.signal import convolve2d
+from scipy.ndimage import gaussian_filter
+
+output_dir = "output/photons/plots/"
+ground_truth_path = "output/photons/ground_truth_spectrally_pure_data_cube_reshaped.fits"
+dep_ref = "data/ECCCO_speedtest_runs/eccco_sp_puremaps_psf4pix_inelectrons_cm3perspersr_with_tables.fits"
+output_paths = ["output/photons/combined_ECCCO_trim_sw_lw_s_i_scaled_spectrally_pure_data_cube_x2_1.0_0.1_wpsf.fits",
+                "output/photons/combined_ECCCO_trim_sw_lw_s_i_scaled_spectrally_pure_data_cube_x2_1.0_0.01_wpsf.fits",
+                "output/photons/combined_ECCCO_trim_sw_lw_s_i_scaled_spectrally_pure_data_cube_x2_1.0_0.2_wpsf.fits",
+                "output/photons/combined_ECCCO_trim_sw_lw_s_i_scaled_spectrally_pure_data_cube_x2_1.0_0.005_wpsf.fits"]
+
+with fits.open(dep_ref) as hdul:
+    dep_table = hdul[2].data
+
+for output_path in output_paths:
+    short_path = os.path.splitext(os.path.basename(output_path))[0]
+    ground_truth = fits.getdata(ground_truth_path)
+    output_data = fits.getdata(output_path)
+
+    for i in range(len(dep_table)):
+        gt = np.flipud(ground_truth[i].T)
+        gt = gaussian_filter(gt, 4 / (2 * np.sqrt(2 * np.log(2))))
+
+        vmin, vmax = 0, np.nanpercentile(ground_truth[i], 99.9)
+        interval = vmax * 0.5
+
+        # initial image comparison
+        fig, axs = plt.subplots(ncols=3, figsize=(30, 10))
+        im1 = axs[0].imshow(gt, vmin=vmin, vmax=vmax, origin='lower')
+        axs[0].set_title("Ground truth")
+        axs[1].imshow(output_data[i], vmin=vmin, vmax=vmax, origin='lower')
+        axs[1].set_title("Output")
+        im2 = axs[2].imshow(gt - output_data[i],
+                            origin='lower', cmap='seismic', vmin=-interval, vmax=interval)
+        axs[2].set_title("Ground truth - output")
+        for ax in axs:
+            ax.set_aspect(0.5)
+        fig.colorbar(im1, ax=axs[:2])
+        fig.colorbar(im2, ax=axs[2])
+        fig.suptitle(dep_table[i][1])
+        fig.savefig(output_dir + f"spectral_pure_{dep_table[i][1]}_{short_path}.png")
+        plt.close()
+
+        # scatter plots
+        fig, ax = plt.subplots()
+        ax.plot(gt.flatten(), output_data[i].flatten(), '.', ms=1)
+        ax.plot([0, np.max(ground_truth[i])], [0, np.max(ground_truth[i])], 'r-')
+        ax.set_title(dep_table[i][1])
+        ax.set_xlabel("Ground truth")
+        ax.set_ylabel("Unfolded recreation")
+        fig.savefig(output_dir + f"scatter_{dep_table[i][1]}_{short_path}.png")
+        plt.close()
+
+        # histogram
+        Z, xedges, yedges = np.histogram2d(gt.flatten(), output_data[i].flatten(), 50)
+        fig, ax = plt.subplots()
+        im = ax.pcolormesh(xedges, yedges, np.log10(Z.T + 1E-3))
+        ax.plot([0, np.max(ground_truth[i])], [0, np.max(ground_truth[i])], 'r-')
+        ax.set_title(dep_table[i][1])
+        fig.colorbar(im, ax=ax)
+        fig.savefig(output_dir + f"histogram_{dep_table[i][1]}_{short_path}.png")
+        plt.close()

experiment_scripts/create_spectrally_pure_images.py

@@ -0,0 +1,18 @@
+import numpy as np
+
+from magixs_data_products import MaGIXSDataProducts
+
+dir_path = "output/test/"
+
+mdp = MaGIXSDataProducts()
+# image_list = ["output/photons/combined_ECCCO_trim_sw_lw_s_i_scaled_em_data_cube_x2_1.0_0.2_wpsf.fits",
+#               "output/photons/combined_ECCCO_trim_sw_lw_s_i_scaled_em_data_cube_x2_1.0_0.1_wpsf.fits",
+#               "output/photons/combined_ECCCO_trim_sw_lw_s_i_scaled_em_data_cube_x2_1.0_0.01_wpsf.fits",
+#               "output/photons/combined_ECCCO_trim_sw_lw_s_i_scaled_em_data_cube_x2_1.0_0.005_wpsf.fits"]
+#image_list = ["/Users/jhughes/Desktop/ECCCO_unfolding_share/output/combined_ECCCO_trim_sw_lw_s_i_scaled_em_data_cube_x2_1.0_0.2_wpsf.fits"]
+image_list = ["data/ECCCO_speedtest_runs/combined_ECCCO_sw_lw_s_i_scaled_em_data_cube_x2_1.0_0.005_wpsf.fits"]
+# gnt_file = "data/ECCCO_speedtest_runs/master_gnt_eccco_inelectrons_cm3perspersr_with_tables.fits"
+gnt_file = "data/ECCCO_speedtest_runs/master_gnt_eccco_inphotons_cm3persperpix_with_tables.fits"
+
+rsp_dep_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] #np.round((np.arange(57, 68+1, 1) / 10.0), decimals=1)
+mdp.create_level2_0_spectrally_pure_images(image_list, gnt_file, rsp_dep_list, dir_path)