edwardberman / shopt

Shear Optimization with ShOpt.jl, a julia package for empirical point spread function characterizations

Table of Contents

• About
  • Who Should Use ShOpt
  • Analytic Profile Fits
  • Pixel Grid Fits
  • Interpolation Across the Field of View
• Inputs and Outputs
  • Inputs
  • Outputs
• Running
  • Command
  • Dependencies
  • Set Up
  • Multi-Threading
  • Testing
• Program Architecture
• How Can I Contribute?
• Config / YAML Information
• Known Issues
• Contributors
• Further Acknowledgements

About

Shear Optimization with ShOpt.jl, a julia library for empirical point spread function characterizations. We aim to improve upon the current state of Point Spread Function Modeling by using Julia to leverage performance gains, use a different mathematical formulation than the literature to provide more robust analytic and pixel grid fits, improve the diagnostic plots, and add features such as wavelets and shapelets. At this projects conclusion we will compare to existing software such as PIFF and PSFex. Work done under McCleary's Group.

We release the related benchmarking code at https://github.com/mcclearyj/cweb_psf

Start by Cloning This Repository. Then see TutorialNotebook.ipynb or follow along the rest of this README.md to get started! Note that the commands in the tutorial notebook are meant to give a sense of procedure and can be executed with the Julia REPL itself.

Who Should Use

Users looking for empirical point spread function characterization software tailored for the data coming from the James Webb Space Telescope, or on a dataset with the similar characteristics. For example, the point spread function spans 100s of pixels because of the pixel scale of your camera, the point spread function is not well approximated by an analytic profile, or the point spread function varies alot across the field of view. For any of these reasons, you should consider using ShOpt.jl. ShOpt.jl is not a single function package, and we would encourage the user to explore the full functionality of ShOpt.jl in the sample config to tailor the software to their needs.

JWST Data is now publicly available at: https://cosmos.astro.caltech.edu/page/cosmosweb-dr. ShOpt was evaluated across all of the wavelengths in the 30mas pixel scale at this link.

Analytic Profile Fits

ShOpt.jl's analytic profile fitting takes inspiration from a number of algorithms outside of astronomy, notably SE-Sync, an algorithm that solves the robotic mapping problem by considering the manifold properties of the data. With sufficiently clean data, the SE-Sync algorithm will descend to a global minimum constrained to the manifold $SE(d)^n / SE(d)$. Following suit, we are able to put a constraint on the solutions we obtain to $[s, g_1, g_2]$ to a manifold. The solution space to $[s, g_1, g_2]$ is constrained to the manifold $$B_2(r) \times \mathbb{R}_{+}$$. The existence of the constraint on shear is well known; nevertheless, the parameter estimation task is usually framed as an unconstrained problem.

Path to $[s, g_1, g_2]$ in $B_2(r) \times \mathbb{R}_+$

Pixel Grid Fits

PCA Mode

We used the first n weights of a Principal Component Analysis and use that to construct our PSF in addition to a smoothing kernel to account for aliasing

Autoencoder Mode

For doing Pixel Grid Fits we use an autoencoder model to reconstruct the Star
PCA mode

function pca_image(image, ncomponents)    
  #Load img Matrix
  img_matrix = image

  # Perform PCA    
  M = fit(PCA, img_matrix; maxoutdim=ncomponents)    

  # Transform the image into the PCA space    
  transformed = MultivariateStats.transform(M, img_matrix)    

  # Reconstruct the image    
  reconstructed = reconstruct(M, transformed)    

  # Reshape the image back to its original shape    
  reconstructed_image = reshape(reconstructed, size(img_matrix)...)    
end    

Autoencoder mode

# Encoder    
encoder = Chain(    
                Dense(r*c, 128, leakyrelu),    
                Dense(128, 64, leakyrelu),    
                Dense(64, 32, leakyrelu),    
               )    
#Decoder
decoder = Chain(    
                Dense(32, 64, leakyrelu),    
                Dense(64, 128, leakyrelu),    
                Dense(128, r*c, tanh),    
               )    
#Full autoencoder
autoencoder = Chain(encoder, decoder)    

#x_hat = autoencoder(x)    
loss(x) = mse(autoencoder(x), x)    

# Define the optimizer    
optimizer = ADAM()    

Interpolation Across the Field of View

[s, g1, g2] are all interpolated across the field of view. Each Pixel is also given an interpolation across the field of view for an nth degree polynomial in (u,v), where n is supplied by the user

The plot on the left shows the average cutout of all stars in a supplied catalog. The plot in the middle shows the average point spread function model for each star. The plot on the right shows the average normalized error between the observed star cutouts and the point spread function model.

Inputs and Outputs

Currently, the inputs are JWST Point Spread Functions source catalogs. The current outputs are images of these Point Spread Functions, Learned Analytic Fits, Learned Pixel Grid Fits, Residual Maps, Loss versus iteration charts, and p-value statisitcs. Not all functionality is working in its current state. Planned functionality for more Shear checkplots.

Inputs

Image	Description
	Star Taken From Input Catalog
shopt.yml	Config File for Tunable Parameters
* starcat.fits	Star Catalog to take vignets from

Outputs

Image	Description
summary.shopt	Fits File containing summary statistics and information to reconstruct the PSF
	Pixel Grid Fit for the Star Above
	Residual Map for Above Model and Fit
	s varying across the field of view
	g1 varying across the field of view
	g2 varying across the field of view
	Histogram for learned profiles for each star in an analytic fit with their residuals
	Same data recorded as a scatterplot with and without outliers removed and with error bars

NB: This is not a comprehensive list, only a few cechkplots are presented. See the shopt.yml to configure which plots you want to see and save!

Running

Command

To run shopt.jl

First use Source Extractor to create a catalog for ShOpt to accept and save this catalog in the appropriate directory

Run julia shopt.jl [configdir] [outdir] [catalog]

There is also a shell script that runs this command so that the user may call shopt from a larger program they are running

Dependencies

Not all of these will be strictly necessary depending on the checkplots you produce, but for full functionality of ShOpt the following are necessary. Source Extractor (or Source Extractor ++) is also not a strict dependency, but in practice one will inevitably install to generate a catalog.

Julia	Python	Binaries	Julia	Julia
Plots	matplotlib	SEx	ProgressBars	BenchmarkTools
ForwardDiff	astropy		UnicodePlots	Measures
LinearAlgebra	numpy		CSV	Dates
Random			FFTW	YAML
Distributions			Images	CairoMakie
SpecialFunctions			ImageFiltering	Flux
Optim			DataFrames	QuadGK
IterativeSolvers			PyCall	Statistics

Set Up

Start by cloning this repository. There are future plans to release ShOpt onto a julia package repository, but for now the user needs these files contents.

The dependencies can be installed in the Julia REPL. For example:

import Pkg; Pkg.add("PyCall")

We also provide dependencies.jl, which you can run to download all of the Julia libraries automatically by reading in the imports.txt file. Simply run julia dependencies.jl in the command line. For the three python requirements, you can similarly run python dependenciesPy.py. We also have a project enviornment for the Julia packages in the folder shopt_env. It would be straightforward to add the line activate("shopt_env") to the first line of shopt.jl to read in this environment. Note that Pycall would have to be built within this environment in some of the next steps.

For some functionality we need to use wrappers for Python code, such as reading in fits files or converting (x,y) -> (u,v). Thus, we need to use certain Python libraries. Thankfully, the setup for this is still pretty straightfoward. We use PyCall to run these snippets.

There are four different methods to get Julia and Python to interopt nicely. We provide all of them as some of these methods play better with different systems.

First install the required Python libraries via pip (which is what dependenciesPy.py does). Now, for method 1, invoke the following in the julia REPL:

using PyCall
ENV["PYTHON"] = "/path_desired_python_directory/python_executable"; import Pkg; Pkg.build("PyCall")
pyimport("astropy")

It is important that the command pyimport("astropy") should be run after re-loading the Julia REPL.

Method 2. If you have a Conda Enviornment setup, you may find it easier to run

using PyCall
pyimport_conda("astropy", "ap") #ap is my choice of name and astropy is what I am importing from my conda Enviornment

Method 3. julia also has a way of working with conda directly. Note that julia will create its own conda enviornment to read from.

using Conda
Conda.add("astropy", :my_env) #my conda enviornment is named my_env

You may also do using Conda; Conda.add("astropy", "/path/to/directory") or using Conda; Conda.add("astropy", "/path/to/directory"; channel="anaconda")

Method 4. On the off chance that none of these works, a final method may look like the following

using PyCall
run(`$(PyCall.python) -m pip install --upgrade cython`)
run(`$(PyCall.python) -m pip install astropy`) 

After the file contents are downloaded the user can run julia shopt.jl [configdir] [outdir] [catalog] as stated above. Alternatively, they can run the shellscript that calls shopt in whatever program they are working with to create their catalog. For example, in a julia program you may use run(`./runshopt.sh [configdir] [outdir] [catalog]`)

Multithreading

Before running, we recommend that users run export JULIA_NUM_THREADS=4 on Unix machines set JULIA_NUM_THREADS=4 . By default Julia will run the program on a single thread, but the polynomial interpolation step is inherently parallelizable. The program is set to use all of threads available to it. You may do more than 4, just be cautious about what your system can provide and when you start getting dimminishing returns.

Testing

To test that everything works, running the dependencies.jl should test that everything is installed correctly in addition to downloading. Running

pyimport("astropy")
pyimport("matplotlib")
pyimport("numpy")

in the Julia REPL should ensure that Julia and Python are interopping correctly.

Additionally, in the Julia REPL, we may write

using Base.Threads
nthreads()

to make sure we are using 4 or more threads. The TutorialNotebook.ipynb will walk you through and test all of the functionality of the program.

Program Architecture

TutorialNotebook.ipynb

Run ShOpt inside of a Jupyter Notebook and learn both how to run the program and how to reconstruct the PSF

shopt.jl

A runner script for all functions in this software

dataPreprocessing.jl

A wrapper for python code to handle fits files and dedicated file to deal with data cleaning

dataOutprocessing.jl

Convert data into a summary.shopt file. Access this data with reader.jl. Produces some additional python plots.

reader.jl

Get Point Spread Functions at an arbitrary (u,v) by reading in a summary.shopt file

plot.jl

A dedicated file to handle all plotting

radialProfiles.jl

Contains analytic profiles such as a Gaussian Fit and a kolmogorov fit

analyticLBFGS.jl

Provides the necessary arguments (cost function and gradient) to the optimize function for analytic fits

pixelGridAutoencoder.jl

Houses the function defining the autoencoder and other machine learning functions supplied to Flux's training call

interpolate.jl

For Point Spread Functions that vary across the Field of View, interpolate.jl will fit a nth degree polynomial in u and v to show how each of the pixel grid parameters change across the (u,v) plane

outliers.jl

Contains functions for identifying and removing outliers from a list

powerSpectrum.jl

Computes the power spectra for a circle of radius k, called iteratively to plot P(k) / k

kaisserSquires.jl

Computes the Kaisser-Squires array to be plotted

runshopt.sh

A shell script for running Shopt. Available so that users can run a terminal command in whatever program they are writing to run shopt.

LICENSE

MIT LICENSE

README.md

User guide, Dependencies, etc.

index.md

For official website

_config.yml

Also for official website

imports.txt

List of Julia Libraries used in the program

packages.txt

List of Python Libraries used in the program

dependencies.jl

Download all of the imports from imports.txt automatically

dependenciesPy.py

Download all of the imports from packages.txt automatically

Config / YAML Information

saveYaml

• Set true if you want to save the YAML to the output directory for future reference, set to false otherwise

NNparams

• epochs
  • Set the Maximum Number of training epochs should the model never reach the minimum gradient of the loss function. Set to 1000 by default
• minGradientPixel
  • A stopping gradient of the loss function for a pixel grid fit. Set to 1e-6 by default

AnalyticFitParams

• minGradientAnalyticModel
  • A stopping gradient of the loss function for an analytic profile fit for input star vignets from a catalog. Set to 1e-6 by default
• minGradientAnalyticLearned
  • A stopping gradient of the loss function for an analytic profile fit for stars learned by a pixel grid fit. Set to 1e-6 by default
• analyticFitStampSize
  • The box size for the subset of your stamp (see stamp size) you wish to use for analytic profile fitting. Ensure to specify this to be smaller than the stamp size of the vignets themselves. Set to 64 by default, therefore fitting an analytic profile to the middle 64 x 64 pixels of the stamp.

dataProcessing

• SnRPercentile
  • Supply a float that represents the percentile below which stars will be filtered by on the basis of signal to noise ratio. Set to 0.33 by default
• sUpperBound
  • Stars fit with an analytic profile are filtered out if their s exceeds this upper bound. Set to 1 by default
• sLowerBound
  • Stars fit with an analytic profile are filtered out if their s falls beneath this lower bound. Set to 0.075 by default

plots

• Set true to plot and save a given figure, false otherwise

polynomialDegree

• The degree of the polynomial used to interpolate each pixel in the stamp across the field of view. Set to 3 by default

stampSize

• The size of the vignet for which use wish to fit. Used interpolation for oversampling and a simple crop for undersampling. Set to 131 by default to fit and interpolate 131 x 131 pixels

training_ratio

• Before doing a polynomial interpolation, the remaining stars will be divided into training and validation stars based off of this float. Set to 0.8 by default, indicating 80% training stars 20% validation stars

CommentsOnRun

• This is Where You Can Leave Comments or Notes To Self on the Run! Could be very useful if you save the yaml file with each run

How to Contribute

Do one of the following:

• Contact [email protected]
• Fork this repository or submit a pull request

If you would like to report an issue or problem please raise an issue on this repository. If you would like to seek support, again see [email protected] .

Known Issues and Enhancements

• We are working on more configurability in the YAML for a smoother user experience, however, everything in this repository is functional and ready to work out of the box
• We are working on a chisq pixel grid fit mode. Right now the gain value g is hardcoded for chi square pixel grid fits

Contributors

• Edward Berman
• Jacqueline McCleary

With help from collaborators at COSMOS-Web: The JWST Cosmic Origins Survey

Further Acknowledgements

• The Northeastern Cosmology Group for Their Continued Support and Guidance
• The Northeastern Physics Department and Northeastern Undergraduate Research and Fellowships, for making this project possible with funding from the Northeastern Physics Research Co-Op Fellowship and PEAK Ascent Award respectively
• David Rosen, who gave valuable input in the early stages of this project and during his course Math 7223, Riemannian Optimization
• The COSMOS Web Collaboration for providing data from the James Webb Space Telescope and internal feedback

Cite

Citation for the current Astronomical Journal (AJ) preprint:

@misc{berman2024efficient,
      title={Efficient PSF Modeling with ShOpt.jl: A PSF Benchmarking Study with JWST NIRCam Imaging}, 
      author={Edward Berman and Jacqueline McCleary and Anton M. Koekemoer and Maximilien Franco and Nicole E. Drakos and Daizhong Liu and James W. Nightingale and Marko Shuntov and Diana Scognamiglio and Richard Massey and Guillaume Mahler and Henry Joy McCracken and Brant E. Robertson and Andreas L. Faisst Caitlin M. Casey and Jeyhan S. Kartaltepe},
      year={2024},
      eprint={2401.11625},
      archivePrefix={arXiv},
      primaryClass={astro-ph.IM}
}

Citation for the current Journal of Open Source Software (JOSS) preprint:

@misc{berman2023shoptjl,
      title={ShOpt.jl: A Julia Package for Empirical Point Spread Function Characterization of JWST NIRCam Data}, 
      author={Edward Berman and Jacqueline McCleary},
      year={2023},
      eprint={2310.00071},
      archivePrefix={arXiv},
      primaryClass={astro-ph.IM}
}

shopt's Issues

[JOSS] Improving the documentation

This issue is part of the JOSS review process.

I have been following the installation instructions and running of the program on the test catalog. Here are some comments regarding the documentation.

• The part on how to compile Python libraries for use in Julia is not straightforward to follow. In particular, it is not clear that there are three, even four, ways to achieve the same thing. I advise to re-formulate this part, stating ahead that there are these X different ways (e.g., as a list), and guide the user to the right section, instead of having all the commands in a row.
• It is important that the command pyimport("astropy") should be run after re-loading the Julia REPL, so it would better to have this explicitely mentioned.
• Is the TutorialNotebook.ipynb made just for copy-pasting lines of codes into the Julia REPL, or can this be executed using Julia directly? This is not clear to me. Probably because this is the first time I run Julia code: it means that the documentation can be be improved.
• using Base.Threads is missing in TutorialNotebook.ipynb.

Relative Error Loss Function NaN

Mean Squared Error Loss Function can handle NaN out of the box with Flux but the Relative Error Loss Function is acting more detrimental. NaN's replace all model stars with I(x,y) < -1000 as a way of denoting other Stars in the image (Source Extractor sets these values to -e32)

Aligning Failed Stars

Shopt Currently produces robust Analytic and Pixel Grid fits most of the time. Additionally, on the chance that a fit fails, I have Shopt record such. I wanted to raise an issue as a self note to make sure not to include the failed instances in the plotting, because at the moment these failed fits are making good data seem indecipherable, at least, overshadowing what I would otherwise deem successful. I am assigning this issue to myself to fix.

Bugs/enhancements for ShOpt runs

I ran into the following issues when running ShOpt for the first time.

Bugs

• dataPreprocessing.jl assumes VIGNET exists in HDU 2 of a FITS file. This is a very common issue; it would be nice if ShOpt handled this gracefully.
• reader.jl looks for an extension that doesn't exist in the default truncated mode (3)
• Star cutouts are off-center in "comfort plot" cutouts produced during run. Is this systematic? It could be minor but if that off-centered star is fed to the pixel grid fitter, that could propagate into the model and cause real problems down the line.

Suggested enhancements

• It would be nice if one could supply other config files to ShOpt. Half the point of using configuration files is not having to change the same file for different use cases.
• Not a bug, but the polynomial degree of the analytic (Gaussian) profile fit should be a tunable parameter. At the very least, it should match the pixel grid fit polynomial degree.

[JOSS] Reproducibility: Figure 2

Hi,

I am trying to reproduce Figure 2 of the paper, and while it appears visually appealing, it would be helpful to have some additional clarifications:

1. It would be beneficial if you could specify in the repository (e.g. README) or on the GitHub site the precise JWST data file you downloaded and whether any modifications were made to the default SExtractor configuration files.
2. Is it possible for users to generate Figure 2 by simply providing the exact catalog and using TutorialNotebook.ipynb, or would additional code be required?

Thanks a lot for your time!

[JOSS] Default output directory naming on macos

The default naming of the output directory follows something like 2024-01-15T11:20:40.848, but : characters are forbidden in macos path names (and are replaced by the / when displayed in a file explorer). Hence I would suggest not to use these characters but - or _ instead.

[JOSS] Paper Comments

The tool and paper look great! Here are some comments (mostly minor):

Summary

• 16: You say "...essentially point sources before their light passes through the atmosphere and telescope..." Obviously this is normally correct from the ground, but since this is specifically a JWST paper you may want to remove the atmosphere part, or change it to "...through the atmosphere (when observing from the ground) and telescope..."
• 23: Since you've restricted this to NIRCam for now, I'd add that in the last paragraph: "...tailored to JWST imaging" --> "...tailored to JWST NIRCam imaging"

Statement of Need

• 32: Don't need to redefine JWST here
• 33: Please put an "e.g.," in front of the COSMOS citation, NIRCam offers a lot more scientific opportunities than just that!
• 67: less --> fewer
• While I agree that figure 2 highlights the need for non-Gaussian PSF models, would it be possible to make this a 4-panel figure showing what happens if you use a 2D Gaussian? I know that people are doing this, and it may be nice to show how much larger the residuals are
• 131 x 131 pixels is very large...I would be surprised if your fitter is actually able to pull anything useful out at that scale. Even your Figure 2 is just showing 75 pixels, and already it's basically noise at the edge even on a log scale. Even at 20 pixels or so (radius) you should have almost all the flux, so I'm not sure if/why you would need to go all the way to 4-5 arc seconds. Contamination from other objects would be pretty common at that level as well.

State of the Field

• 87: spelling out PIFF (and PSFex for that matter) should go to the first mention in the previous section
• 95: Best to change "twenty-five years ago" to a year, so that in future years the line still makes sense
• 125: "takes up" --> "covers"

[JOSS] Missing references

In the current state of the JOSS paper, I suggest to mention STARRED as another efficient PSF reconstruction technique, as it has been shown to outperform other programs like PSFex (mentioned in the paper) and PSFr (which would need to be mentioned as well), and uses high the just-in-time compiled and auto-diff library JAX.

kolmogorov fit

Fitting a kolmogorov radial fit is taking unfeasibly long. I know that there is an added degree of difficulty in numerical integration but perhaps something else is causing a bottleneck? I have tried printing loss with every iteration and it seems to be working, but I am working it may never even converge enough to satisfy optim.jl's requirements.

[JOSS] Paper comments

Hi,

Your JOSS paper is pretty well written and clear! I have some minor suggestions, including fixing a few typos:

• I would remove the Introduction title, leaving simply Summary;
• Typo: 'unprecendented' (line 35);
• Typo: ‘biasses’ (line 67);
• Typo: ‘emblamatic’ (line 76);
• For the State of the Field section: I wonder how ShOpt.jl would compare with PSFr (https://github.com/sibirrer/psfr), which is quite recent and has been used for JWST data;
• Beichman et al., 2012a and Beichman et al., 2012b are exactly the same paper.

Code Cleanup is overdue. Actionable things include putting in strict types for functions and adding keyword arguments, fixing naming conventions, adding explanatory code comments, and adding a config file. Unsure if the config should be a YAML or if I can get away with a TXT