<jupyter_start><jupyter_text>Import Statements, Preliminaries<jupyter_code>%autosave 180
%matplotlib inline
import numpy as np
from scipy.stats import multivariate_normal
from multiprocessing import Pool
from matplotlib import pyplot as plt
from astropy.utils.data import get_pkg_data_filename
import galsim
import os
import math
import time
import logging
import sys
from astropy.io import fits
simname = "/nfs/slac/des/fs1/g/sims/jderose/BCC/Chinchilla/Herd/Chinchilla-4/v1.9.2/addgalspostprocess/mag2/Chinchilla-4-y3wlpz.0.fits"
truthname = "/nfs/slac/des/fs1/g/sims/jderose/BCC/Chinchilla/Herd/Chinchilla-4/v1.9.2/addgalspostprocess/truth/truth/Chinchilla-4_lensed.4.fits"
simhdu = fits.open(simname)
truthhdu = fits.open(truthname)
sim_table = simhdu[1].data
truth_table = truthhdu[1].data<jupyter_output><empty_output><jupyter_text>Functions and List Creation<jupyter_code>cat_file_name = os.path.join('/nfs/slac/des/fs1/g/sims/jderose/BCC/Chinchilla/Herd/Chinchilla-4/v1.9.2/addgalspostprocess/truth/truth/','Chinchilla-4_lensed.4.fits')
cat = galsim.Catalog(cat_file_name)

noisy_file_name = os.path.join('/afs/slac.stanford.edu/u/ki/gwisk/CRNProject','noisy.fits')
noisefree_file_name = os.path.join('/afs/slac.stanford.edu/u/ki/gwisk/CRNProject','noisefree.fits')

rband = truth_table['LMAG'][:,1] #grizy
_ra = truth_table['RA']
_dec = truth_table['DEC']
_hlr = truth_table['SIZE']
te1 = truth_table['EPSILON'][:,0]
te2 = truth_table['EPSILON'][:,1]
gaussiannoise = galsim.GaussianNoise(sigma=0.005)
x, y = np.mgrid[-24.5:25.5:1, -24.5:25.5:1]
pos = np.zeros(x.shape + (2,))
pos[:, :, 0] = x; pos[:, :, 1] = y

def HSMfunc(integer):
    _nmoments, _gfluxes, _images, failures = ([] for i in range(4))
    for i in range(integer):
        # catalogs
        hlradius = _hlr[i]   
        elip1 = te1[i]
        elip2 = te2[i]
        lmag = rband[i]
        gflux = 10**((lmag-22.5)/(-2.5))
        _gfluxes.append(gflux)
        #make galaxy
        gal = galsim.Gaussian(flux=gflux, half_light_radius=hlradius)
        # shearing
        gal = gal.shear(galsim.Shear(e1=elip1, e2=elip2))
        # create image, convolve
        image = galsim.ImageF(50, 50)
        final = galsim.Convolve([gal])
        image = final.drawImage(image, scale=0.2)
        image.addNoise(gaussiannoise)
        #create image...again
        image2 = galsim.ImageF(50, 50)
        final2 = galsim.Convolve([gal])
        image2 = final2.drawImage(image2, scale=0.2)
        _images.append(image2)
        #HSM
        nmoment = image.FindAdaptiveMom(strict = False ,guess_sig=hlradius*5) #, precision=1e-6)
        _nmoments.append(nmoment)
        if nmoment.moments_n_iter  == 0:
            failures.append(i)
    return _nmoments, _gfluxes, _images, failures 


def cheatfunc(integer):
    _covariances = []
    _noisematrices = []
    for i in range(integer):
        noiseimage = galsim.ImageF(50, 50)
        # catalogs
        hlradius = _hlr[i] 
        sigma = np.sqrt(2*np.log(2)) * hlradius # Conversion from hlr to sigma necessary
        elip1 = te1[i]
        elip2 = te2[i]
        lmag = rband[i]
        gflux = 10**((lmag-22.5)/(-2.5))
        shearing2 = galsim.Shear(e1=-1*elip1, e2=elip2)
        matrix = (shearing2.getMatrix())*(sigma**2)*25 
        _covariances.append(matrix)
        noiseimage.addNoise(gaussiannoise)
        _noisematrices.append(noiseimage)
    return _covariances, _noisematrices<jupyter_output><empty_output><jupyter_text>Function Calls 1<jupyter_code>_nmoments1, _gfluxes1, _images1, failures  = HSMfunc(2000)
_covariances1, _noisematrices1 = cheatfunc(2000)

nmoments = [i for j, i in enumerate(_nmoments1) if j not in failures]
gfluxes = [i for j, i in enumerate(_gfluxes1) if j not in failures]
covariances = [i for j, i in enumerate(_covariances1) if j not in failures] # HSM doesn't always converge well, so the purpose of this is to remove galaxies for which HSM didn't converge.
images = [i for j, i in enumerate(_images1) if j not in failures]
noisematrices = [i for j, i in enumerate(_noisematrices1) if j not in failures]

cheatfluxes = np.zeros(len(nmoments))
hsmfluxes = np.zeros(len(nmoments))
weird_method = np.zeros(len(nmoments))

for j in range(len(nmoments)):
    rv = multivariate_normal([0, 0], covariances[j])
    r1 = np.sum(np.multiply(rv.pdf(pos), images[j].array + noisematrices[j].array)) # * (1/_hlr[j]**2)) noisematrices[j].array
    r2 = np.sum(np.multiply(rv.pdf(pos), rv.pdf(pos))) #/(_hlr[j])**2)
    cheatfluxes[j] = r1/r2
    
for i in range(len(nmoments)):
    hsmfluxes[i] = nmoments[i].moments_amp
from scipy import integrate
rv = multivariate_normal([0, 0], covariances[j])
f, ((ax1, ax2, ax3), (ax4, ax5, ax6), (ax7, ax8, ax9)) = plt.subplots(3,3,figsize=(15,15))
im1 = ax1.imshow(rv.pdf(pos)*gfluxes[0])
im2 = ax2.imshow(images[0].array)
im3 = ax3.imshow(images[0].array-rv.pdf(pos)*gfluxes[0])
f.colorbar(im1,ax=ax1)
f.colorbar(im2,ax=ax2)
f.colorbar(im3,ax=ax3)

grid = np.array([[integrate.dblquad(lambda b, a : rv.pdf([a,b]), x, x+1, lambda a: y, lambda a: y+1)[0] for y in range(-25, 25)] for x in range(-25, 25)])

print(grid.shape)
im4 = ax4.imshow(grid*gfluxes[0])
im5 = ax5.imshow(images[0].array)
im6 = ax6.imshow(images[0].array-grid*gfluxes[0])
f.colorbar(im4,ax=ax4)
f.colorbar(im5,ax=ax5)
f.colorbar(im6,ax=ax6)

im7 = ax7.imshow(rv.pdf(pos)*gfluxes[0])
im8 = ax8.imshow(grid*gfluxes[0])
im9 = ax9.imshow(rv.pdf(pos)*gfluxes[0]-grid*gfluxes[0])
f.colorbar(im7,ax=ax7)
f.colorbar(im8,ax=ax8)
f.colorbar(im9,ax=ax9)

f.tight_layout()
print(np.sum(images[0].array), gfluxes[0])<jupyter_output>(50, 50)
0.11153627 0.11153710526316728<jupyter_text>Plot of Matrix Flux vs. HSM Flux and Plot of Matrix Flux vs. Galaxy Flux<jupyter_code>from sklearn.linear_model import LinearRegression
cheatfluxes = cheatfluxes.reshape((-1, 1))
model = LinearRegression().fit(cheatfluxes, hsmfluxes)
print('intercept:', model.intercept_)
print('slope:', model.coef_)
plt.figure()
z = np.linspace(0, 2, 100)
y = z*model.coef_ + model.intercept_

plt.plot(z, y, '-r')
plt.scatter(cheatfluxes, hsmfluxes)
plt.xlabel("matrix fluxes")
plt.ylabel("hsm fluxes")


plt.figure()
plt.scatter(cheatfluxes, gfluxes)
plt.xlabel("matrix fluxes")
plt.ylabel("galaxy fluxes")

cheatfluxes = cheatfluxes.reshape((-1, 1))
gfluxes = np.asarray(gfluxes)
model = LinearRegression().fit(cheatfluxes, gfluxes)
print('intercept:', model.intercept_)
print('slope:', model.coef_)
plt.scatter(gfluxes, hsmfluxes)
plt.xlabel("galaxy fluxes")
plt.ylabel("hsm fluxes")
plt.scatter(gfluxes, hsmfluxes - gfluxes)
plt.scatter(gfluxes, [cheatfluxes[i]-gfluxes[i] for i in range(len(cheatfluxes))])<jupyter_output><empty_output><jupyter_text>Delta Flux Function Calls<jupyter_code>_nmoments1, _gfluxes1, _images1, failures  = HSMfunc(2000)
_covariances1, _noisematrices1 = cheatfunc(2000)

nmoments = [i for j, i in enumerate(_nmoments1) if j not in failures]
gfluxes = [i for j, i in enumerate(_gfluxes1) if j not in failures]
covariances = [i for j, i in enumerate(_covariances1) if j not in failures]
images = [i for j, i in enumerate(_images1) if j not in failures]
noisematrices = [i for j, i in enumerate(_noisematrices1) if j not in failures]

cheatfluxes = np.zeros(len(nmoments))
hsmfluxes = np.zeros(len(nmoments))

for j in range(len(nmoments)):
    rv = multivariate_normal([0, 0], covariances[j])
    r1 = np.sum(np.multiply(rv.pdf(pos), noisematrices[j].array)) # * (1/_hlr[j]**2)) noisematrices[j].array
    r2 = np.sum(np.multiply(rv.pdf(pos), rv.pdf(pos))) #/(_hlr[j])**2)
    cheatfluxes[j] = r1/r2
    
for i in range(len(nmoments)):
    hsmfluxes[i] = nmoments[i].moments_amp - gfluxes[i]
print(np.amax(hsmfluxes), np.amin(np.abs(hsmfluxes)), np.sum(hsmfluxes)/len(hsmfluxes))

print(np.amax(cheatfluxes), np.amin(np.abs(cheatfluxes)), np.sum(cheatfluxes)/len(cheatfluxes))
plt.figure()
plt.scatter(cheatfluxes, gfluxes)
plt.xlabel("matrix fluxes")
plt.ylabel("galaxy fluxes")
plt.figure()
plt.scatter(cheatfluxes, hsmfluxes)
plt.xlabel("matrix fluxes")
plt.ylabel("hsm fluxes")
plt.figure()
plt.scatter(hsmfluxes, gfluxes)
plt.xlabel("hsm fluxes")
plt.ylabel("galaxy fluxes")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15296999387630128, 'fraction_numerical': 0.035762400489895894, 'mean_word_length': 4.686629526462395, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '283011', 'license': 'MIT', 'n_tokens_mistral': 3625, 'n_tokens_neox': 3344, 'n_words': 611, 'path': '.ipynb_checkpoints/Weight_Matrix_Test-checkpoint.ipynb', 'repo_name': 'cahebert/read-noise', 'size': '8165'}