diff --git a/.gitignore b/.gitignore index eaa66aa..d2cd122 100644 --- a/.gitignore +++ b/.gitignore @@ -4,7 +4,7 @@ test_arctic libarctic.so *.cpython*.so lib_test -arcticpy/src/wrapper.cpp +python/arcticpy/wrapper.cpp __pycache__/ build/ gsl/ diff --git a/makefile b/makefile index 8236998..a18d762 100644 --- a/makefile +++ b/makefile @@ -67,11 +67,11 @@ DIR_INC := $(DIR_ROOT)/include DIR_TEST := $(DIR_ROOT)/test #DIR_GSL ?= /cosma/local/gsl/2.5/lib #DIR_OMP ?= /cosma/local/openmpi/gnu_11.1.0/4.1.4/lib -#DIR_GSL ?= $(DIR_HOMEBREW) -#DIR_OMP ?= $(DIR_HOMEBREW) +DIR_GSL ?= $(DIR_HOMEBREW) +DIR_OMP ?= $(DIR_HOMEBREW) #DIR_OMP ?= $(DIR_MACPORTS)/libomp -DIR_OMP ?= $(DIR_MACPORTS) -DIR_GSL ?= $(DIR_MACPORTS) +#DIR_OMP ?= $(DIR_MACPORTS) +#DIR_GSL ?= $(DIR_MACPORTS) # Fallback self-installing GSL #DIR_GSL ?= $(DIR_ROOT)/gsl DIR_WRAPPER := $(DIR_ROOT)/python/arcticpy @@ -99,9 +99,9 @@ LIBARCTIC := -L $(DIR_ROOT) -Wl,-rpath,$(DIR_ROOT) -l$(TARGET) # Add multithreading to reduce runtime (requires OpenMP to have been installed) CXXFLAGS += -Xpreprocessor -fopenmp # Use this on a mac -# LIBS += -L $(DIR_OMP)/lib -lomp + LIBS += -L $(DIR_OMP)/lib -lomp # Use the following on cosma (can also use with macports) -LIBS += -L $(DIR_OMP)/lib -lgomp +#LIBS += -L $(DIR_OMP)/lib -lgomp diff --git a/python/arcticpy/read_noise.c b/python/arcticpy/read_noise.c new file mode 100644 index 0000000..eb0d5e4 --- /dev/null +++ b/python/arcticpy/read_noise.c @@ -0,0 +1,149 @@ +#include +#include +#include + +#define SQUARE(x) ((x) * (x)) + +void determine_noise_model(const double* imageIn, const double* imageOut, const int rows, const int cols, const double readNoiseAmp, const double readNoiseAmpFraction, const int smoothCol, double* output) { + + double* dval0 = output; // using already allocated output buffer to temporarily store the results to avoid allocating more memory (safe as long as we do not modify output before storing dval0 there and not use dval0 after modifying output for the first time) + double* dval9 = (double*)malloc(rows * cols * sizeof(double)); + + double mod_clip = readNoiseAmp * readNoiseAmpFraction; + double readNoiseAmp2 = SQUARE(readNoiseAmp); + + int idx; + + for (int i = 0; i < rows * cols; i++) { + dval0[i] = imageIn[i] - imageOut[i]; + dval9[i] = dval0[i]; + } + + // get the dval9 average + // inner part + for (int i = 1; i < rows - 1; i++) { +#pragma omp parallel for private(idx) + for (int j = 1; j < cols - 1; j++) { + idx = i * cols + j; + dval9[idx] += dval0[idx + cols + 1]; // comparison with bottom-left neighbour + dval9[idx] += dval0[idx + cols]; // comparison with bottom-central neighbour + dval9[idx] += dval0[idx + cols - 1]; // comparison with bottom-right neighbour + dval9[idx] += dval0[idx + 1]; // comparison with middle-left neighbour + dval9[idx] += dval0[idx - 1]; // comparison with middle-right neighbour + dval9[idx] += dval0[idx - cols + 1]; // comparison with top-left neighbour + dval9[idx] += dval0[idx - cols]; // comparison with top-central neighbour + dval9[idx] += dval0[idx - cols - 1]; // comparison with top-right neighbour + dval9[idx] /= 9; + } + } + // edges (excl. corners) +#pragma omp parallel for private(idx) + for (int i = 1; i < rows - 1; i++) { + // left edge (j=0) + idx = i * cols; + dval9[idx] += dval0[idx + cols + 1]; // comparison with bottom-left neighbour + dval9[idx] += dval0[idx + cols]; // comparison with bottom-central neighbour + dval9[idx] += dval0[idx + 1]; // comparison with middle-left neighbour + dval9[idx] += dval0[idx - cols + 1]; // comparison with top-left neighbour + dval9[idx] += dval0[idx - cols]; // comparison with top-central neighbour + dval9[idx] /= 6.0; + // right edge (j=cols-1) + idx = (i + 1) * cols - 1; + dval9[idx] += dval0[idx + cols]; // comparison with bottom-central neighbour + dval9[idx] += dval0[idx + cols - 1]; // comparison with bottom-right neighbour + dval9[idx] += dval0[idx - 1]; // comparison with middle-right neighbour + dval9[idx] += dval0[idx - cols]; // comparison with top-central neighbour + dval9[idx] += dval0[idx - cols - 1]; // comparison with top-right neighbour + dval9[idx] /= 6.0; + } + +#pragma omp parallel for private(idx), shared(dval0, dval9) + for (int j = 1; j < cols - 1; j++) { + // top edge (i=0) + idx = j; + dval9[idx] += dval0[idx + cols + 1]; // comparison with bottom-left neighbour + dval9[idx] += dval0[idx + cols]; // comparison with bottom-central neighbour + dval9[idx] += dval0[idx + cols - 1]; // comparison with bottom-right neighbour + dval9[idx] += dval0[idx + 1]; // comparison with middle-left neighbour + dval9[idx] += dval0[idx - 1]; // comparison with middle-right neighbour + dval9[idx] /= 6.0; + // bottom edge (i=rows-1) + idx = (rows - 1) * cols + j; + dval9[idx] += dval0[idx + 1]; // comparison with middle-left neighbour + dval9[idx] += dval0[idx - 1]; // comparison with middle-right neighbour + dval9[idx] += dval0[idx - cols + 1]; // comparison with top-left neighbour + dval9[idx] += dval0[idx - cols]; // comparison with top-central neighbour + dval9[idx] += dval0[idx - cols - 1]; // comparison with top-right neighbour + dval9[(rows - 1) * cols + j] /= 6.0; + } + // corners + dval9[0] += dval0[cols + 1]; // comparison with bottom-left neighbour + dval9[0] += dval0[cols]; // comparison with bottom-central neighbour + dval9[0] += dval0[1]; // comparison with middle-left neighbour + dval9[0] /= 4.0; + + idx = cols - 1; + dval9[idx] += dval0[idx + cols]; // comparison with bottom-central neighbour + dval9[idx] += dval0[idx + cols - 1]; // comparison with bottom-right neighbour + dval9[idx] += dval0[idx - 1]; // comparison with middle-right neighbour + dval9[idx] /= 4.0; + + idx = (rows - 1) * cols; + dval9[idx] += dval0[idx + 1]; // comparison with middle-left neighbour + dval9[idx] += dval0[idx - cols + 1]; // comparison with top-left neighbour + dval9[idx] += dval0[idx - cols]; // comparison with top-central neighbour + dval9[idx] /= 4.0; + + idx = rows * cols - 1; + dval9[idx] += dval0[idx - 1]; // comparison with middle-right neighbour + dval9[idx] += dval0[idx - cols]; // comparison with top-central neighbour + dval9[idx] += dval0[idx - cols - 1]; // comparison with top-right neighbour + dval9[idx] /= 4.0; + +#pragma omp parallel for shared(output, dval0, dval9) + for (int i = 0; i < cols * rows; i++) { + + // setting dval0u to output + output[i] = fmin(1.0, fmax(-1.0, dval0[i])) * SQUARE(dval0[i]) / (SQUARE(dval0[i]) + 4.0 * readNoiseAmp2); + // adding dval9u to output + output[i] += fmax(fmin(dval9[i], mod_clip), -mod_clip) * SQUARE(dval9[i]) / (SQUARE(dval9[i]) + 18.0 * readNoiseAmp2); + + int rcol = i % cols; + double dmod1, dmod2; + if (i < cols) { // first row + dmod1 = 0.0; + } else { + dmod1 = imageOut[i - cols] - imageOut[i]; + } + output[i] += fmax(fmin(dmod1, mod_clip), -mod_clip) * 4 * readNoiseAmp2 / (SQUARE(dmod1) + 4.0 * readNoiseAmp2); + + if (i >= (rows - 1) * cols) { // last row + dmod2 = 0.0; + } else { + dmod2 = imageOut[i + cols] - imageOut[i]; + } + output[i] += fmax(fmin(dmod2, mod_clip), -mod_clip) * 4 * readNoiseAmp2 / (SQUARE(dmod2) + 4.0 * readNoiseAmp2); + + if (smoothCol) { + double cmod1, cmod2; + if (rcol == 0) { // first column + cmod1 = 0.0; + } else { + cmod1 = imageOut[i - 1] - imageOut[i]; + } + output[i] += fmax(fmin(cmod1, mod_clip), -mod_clip) * 4 * readNoiseAmp2 / (SQUARE(cmod1) + 4.0 * readNoiseAmp2); + if (rcol == cols - 1) { // last column + cmod2 = 0.0; + } else { + cmod2 = imageOut[i + 1] - imageOut[i]; + } + output[i] += fmax(fmin(cmod2, mod_clip), -mod_clip) * 4 * readNoiseAmp2 / (SQUARE(cmod2) + 4.0 * readNoiseAmp2); + + output[i] /= 6.0; + } else { + output[i] /= 4.0; + } + } + + free(dval9); +} diff --git a/python/arcticpy/read_noise.py b/python/arcticpy/read_noise.py index 34baa03..027326c 100644 --- a/python/arcticpy/read_noise.py +++ b/python/arcticpy/read_noise.py @@ -1,4 +1,16 @@ import numpy as np +#import arcticpy as ac +from read_noise_c import determine_noise_model_c + +try: + from arcticpy import wrapper as w +except ImportError: + import wrapper as w + +from arcticpy.cti import add_cti,remove_cti +import matplotlib as mpl +import numpy as np +from scipy.optimize import curve_fit """ CTI correction moves trailed electrons back to their proper location, but also @@ -7,6 +19,70 @@ adjacent pixels. This class provides a set of routines to predict this effect, and counteract it. + +################# +SAMPLE USAGE CASE +################# + +If you have an image as a 2D numpy array and an arctic model (with parameters such as parallel_roe, +parallel_ccd, parallel_traps, etc...) run the following: + +1.) readnoise = arctic.ReadNoise(4.5) + The argument sets the level of readnoise (in units of electrons) + +2.) [optional] readnoise.set_arctic_parameters(**kwargs) + Input each arctic parameter as a keyword argument. + These will then be passed to subsequent function calls automatically + (Note: these choices can always be overridden in functions, if desired) + +3.) [optional] readnoise.optimise_SR_fraction_from_image(image) + Input your data image. The function will determine the optimum + fraction of S+R separation to apply to minimise CTI covariance + +4.) generate_SR_frames_from_image(image) + Do the actual S+R splitting on your science image. If the SR + optimisation routine is not run, assume a correction level of 100% + +5.) [optional] readnoise.measure_simulated_covariance_corners() + estimate the resulting covariance matrices in each ccd corner, using the optimised simulation + +######POSSIBLE ALTERNATIVE OPTION?##### +5.) [optional] readnoise.calculate_covariance_corners_from image(image) + measure covariance matrices using the real image data + (not yet implemented) +####################################### + +6.) [optional] readnoise.plot_covariance_matrices() + plot the covariance matrices of all four corners + (not fully implemented) + + +################# +EXAMPLE COMMANDS +################# +skyImage, readNoiseImage = readnoise.generate_SR_frames_from_image(image, sr_fraction=None, + parallel_roe, parallel_ccd, parallel_traps) + +image_corrected = arctic.cti_correct(skyImage, + parallel_roe, parallel_ccd, parallel_traps) + +image_corrected += readNoiseImage + + +#to be updated + +readnoise.covariance + +covariance = readnoise.etimate_residual_covariance_from_image(image, + matrix_size=5, fprSize=5, + parallel_roe, parallel_ccd, parallel_traps) + +covariance = readnoise.etimate_residual_covariance(sky_level, sky_sigma, + matrix_size=5, fprSize=5, + parallel_roe, parallel_ccd, parallel_traps) + + + """ class ReadNoise: @@ -16,8 +92,11 @@ def __init__( adjacency=0.3, noise_model_scaling=1.0, # originally 0.75 amplitude_scale=0.2, # originally 0.33 - n_iter=200 - serial=True + n_iter=200, + sr_fraction=None, + serial=True, + covariance_matrices=None, + **kwargs ): self.sigmaRN = sigma_readnoise @@ -26,6 +105,42 @@ def __init__( self.outScale = noise_model_scaling self.n_iter = n_iter self.smoothCol = serial + self.sigmaSky = None + self.skyFrameSim = None + self.readnoiseFrameSim = None + self.SRfrac_optimised = None + self.arcKwargs = {} + self.covarianceMatrices = {} + self.dataFrames = {} + +# self._sr_fraction = sr_fraction +# self._sr_fraction_optimised = None + +# @property +# def sr_fraction(self): +# """ +# RMS value of read noise, in units of electrons +# """ +# if self._sr_fraction is not None: return self._sr_fraction +# #if self._sr_fraction_optimised is None: self.optimise_sr_fraction() +# return self._sr_fraction_optimised +# +# @sr_fraction.setter +# def sr_fraction(self, value): +# self._sr_fraction = value +# +# @property +# def sr_fraction_optimised(self, **kwargs): +# """ +# Value of S+R splitting fraction that will minimise pixel-to-pixel covariance +# after CTI correction +# """ +# if self._sr_fraction_optimised is None: self.optimise_sr_fraction(**kwargs) +# return self._sr_fraction_optimised +# +# @sr_fraction_optimised.setter +# def sr_fraction_optimised(self, value): +# self._sr_fraction_optimised = value @property def sigma(self): @@ -39,122 +154,63 @@ def sigma(self, value): if value < 0: value = 0. self._sigma = value + ############### + # + #USER-FACING FUNCTIONS + # ############### - def determine_noise_model(self, imageIn, imageOut): - readNoiseAmp = self.sigmaRN - - dval0 = imageIn - imageOut - dval0u = dval0.copy() - - dval0u[dval0u>1] = 1 - dval0u[dval0u<-1] = -1 - - dval9 = imageIn*0 - dcount = imageIn*0 - - dval9[:-1,:-1]+=(imageIn[1: ,1: ]-imageOut[1: ,1: ]); dcount[:-1,:-1]+=1 - dval9[:-1,: ]+=(imageIn[1: ,: ]-imageOut[1: ,: ]); dcount[:-1,: ]+=1 - dval9[:-1,1: ]+=(imageIn[1: ,:-1]-imageOut[1: ,:-1]); dcount[:-1,1: ]+=1 - dval9[: ,:-1]+=(imageIn[: ,1: ]-imageOut[: ,1: ]); dcount[: ,:-1]+=1 - dval9[: ,: ]+=(imageIn[: ,: ]-imageOut[: ,: ]); dcount[: ,: ]+=1 - dval9[: ,1: ]+=(imageIn[: ,:-1]-imageOut[: ,:-1]); dcount[: ,1: ]+=1 - dval9[1: ,:-1]+=(imageIn[:-1,1: ]-imageOut[:-1,1: ]); dcount[1: ,:-1]+=1 - dval9[1: ,: ]+=(imageIn[:-1,: ]-imageOut[:-1,: ]); dcount[1: ,: ]+=1 - dval9[1: ,1: ]+=(imageIn[:-1,:-1]-imageOut[:-1,:-1]); dcount[1: ,1: ]+=1 - - dval9/=dcount - - readNoiseAmpFraction = self.ampScale - dval9u = dval9.copy() - if type(readNoiseAmp)==np.ndarray: - dval9u = dval9u.flatten() - readNoiseAmp = readNoiseAmp.flatten() - dval9u[dval9u > readNoiseAmp*readNoiseAmpFraction] = (readNoiseAmp*readNoiseAmpFraction)[dval9u > readNoiseAmp*readNoiseAmpFraction] - dval9u[dval9u < readNoiseAmp*-readNoiseAmpFraction] = (readNoiseAmp*-readNoiseAmpFraction)[dval9u < readNoiseAmp*-readNoiseAmpFraction] - dval9u = dval9u.reshape(imageIn.shape) - readNoiseAmp = readNoiseAmp.reshape(imageIn.shape) - else: - dval9u[dval9u > readNoiseAmp*readNoiseAmpFraction] = readNoiseAmp*readNoiseAmpFraction - dval9u[dval9u < readNoiseAmp*-readNoiseAmpFraction] = readNoiseAmp*-readNoiseAmpFraction - - dmod1 = imageIn*0 - dmod1[1:,:] = imageOut[:-1,:] - imageOut[1:,:] - dmod2 = imageIn*0 - dmod2[:-1,:] = imageOut[1:,:] - imageOut[:-1,:] - - if self.smoothCol: - cmod1 = imageIn*0 - cmod1[:,1:] = imageOut[:,:-1] - imageOut[:,1:] - cmod2 = imageIn*0 - cmod2[:,:-1] = imageOut[:,1:] - imageOut[:,:-1] - - dmod1u = dmod1.copy() - if type(readNoiseAmp)==np.ndarray: - dmod1u = dmod1u.flatten() - readNoiseAmp = readNoiseAmp.flatten() - dmod1u[dmod1u>readNoiseAmp*readNoiseAmpFraction] = (readNoiseAmp*readNoiseAmpFraction)[dmod1u>readNoiseAmp*readNoiseAmpFraction] - dmod1u[dmod1ureadNoiseAmp*readNoiseAmpFraction] = readNoiseAmp*readNoiseAmpFraction - dmod1u[dmod1ureadNoiseAmp*readNoiseAmpFraction] = (readNoiseAmp*readNoiseAmpFraction)[dmod2u>readNoiseAmp*readNoiseAmpFraction] - dmod2u[dmod2ureadNoiseAmp*readNoiseAmpFraction] = readNoiseAmp*readNoiseAmpFraction - dmod2u[dmod2ureadNoiseAmp*readNoiseAmpFraction] = (readNoiseAmp*readNoiseAmpFraction)[cmod1u>readNoiseAmp*readNoiseAmpFraction] - cmod1u[cmod1ureadNoiseAmp*readNoiseAmpFraction] = readNoiseAmp*readNoiseAmpFraction - cmod1u[cmod1ureadNoiseAmp*readNoiseAmpFraction] = (readNoiseAmp*readNoiseAmpFraction)[cmod2u>readNoiseAmp*readNoiseAmpFraction] - cmod2u[cmod2ureadNoiseAmp*readNoiseAmpFraction] = readNoiseAmp*readNoiseAmpFraction - cmod2u[cmod2u + + sr_fraction: float + scale fraction specifying what percentage of read noise should be separated into the R + frame. Default is set to the optimised value measured in , or 100% + if has not been previously run. + + + ''' ampReadNoise = self.sigmaRN + #if sr_fraction is unset, use system value + if sr_fraction == None: + sr_fraction = self.SRfrac_optimised + # if optimisation is not set, default to 100% + if sr_fraction == None: + sr_fraction = 1.0 + imageIn = image imageAdj = np.zeros_like(imageIn) #adjustment image (for read noise modeling) imageOut = imageIn.copy() #output ("smoothed") image @@ -166,9 +222,10 @@ def estimate_read_noise_model_from_image(self, image): nrmsGlobal = 0. smoother = 1 + print('iter','model_rms','target_rms','residual') for s in range(self.n_iter): oldchk = ampReadNoise-rmsGlobal - imageAdj= self.determine_noise_model(imageIn,imageOut) + imageAdj= self._determine_noise_model(imageIn,imageOut) if (ampReadNoise-rmsGlobal) > 0: imageOut += imageAdj * self.outScale / smoother @@ -182,7 +239,8 @@ def estimate_read_noise_model_from_image(self, image): nrmsGlobal = noiseImage[cond].size rmsGlobal = np.sqrt(rmsGlobal/nrmsGlobal) - print(s,rmsGlobal, (ampReadNoise-rmsGlobal)) + print("\033[K",end='\r') + print('%4d %f %5.2f %f'%(s,rmsGlobal, self.sigmaRN, (ampReadNoise-rmsGlobal)),end='') chk = ampReadNoise-rmsGlobal if chk*oldchk < 0: @@ -194,121 +252,547 @@ def estimate_read_noise_model_from_image(self, image): else: if abs(ampReadNoise - rmsGlobal) < 0.0001: break - return imageOut,noiseImage - #raise NotImplementedError - #return image + + #scale the smooth and noise frames if != 1 + if sr_fraction !=1: + imageOut+=noiseImage + noiseImage*=sr_fraction + imageOut-=noiseImage + + return imageOut,noiseImage #These are the "S" and "R" frames, respectively + ############### + ############### + def optimise_SR_fraction( + self, + background_level, + background_sigma=None, + n_pixels=2048, + subchip_size=500, + matrix_size=5, + fom_method='box', + **kwargs + ): + ''' + Adjust parameters of the S+R algorithm that estimates the read noise in an image. + The image itself is not required; merely its dimensions and sky level parameters. + + It works by simulating a subchip_size x subchip_size regions of the image, one + at each corner of the CCD, then tries to minimise the pixel-to-pixel covariance + after CTI correction in the most distant corner. + + This merely adjusts the self parameters contained within the ReadNoise instance. + The final optimised S+R value is stored as the class variable + + Parameters + ---------- + background_level : Float + The mean (sky) background level across the image. + + background_sigma : Float + The rms variation per pixel around the background level. + If not specified, it will assume sqrt(background_level) for shot noise. + + n_pixels : Float or (Float, Float) + The number of pixels in the CCD. If a single number is supplied, it assumes + the CCD is square. If two numbers are supplied, the assumed order is (n_y,n_x). + + subchip_size : Float + size of the cutout region, to test CTI behaviour in different regimes (e.g. far from readout, close to readout, etc.) + + matrix_size : Int + Size of covariance matrix used to estimate covariance. By default the program creates a 5x5 matrix + + fom_method : String + Method used to estimate the mean value of a covariance matrix, for figure-of-merit purposes. Can be: + box : take a square (3x3) region around the central pixel [DEFAULT] + row : take all cells along the central row + column : take all cells along the central column + + **kwargs : parameters that characterise CTI features in arCTIc (trap density, CCD, read out electronics (ROE) descriptions) + ''' + + # Parse inputs + if type(n_pixels) == int: + n_pixels=(n_pixels,n_pixels) + chip_size = n_pixels + if background_sigma==None: + background_sigma = np.sqrt(background_level) + self.sigmaSky = background_sigma + + #set prescan regions (this simulates the maximum possible CTI trailing, averaged at the far edge of the real CCD) + if 'parallel_roe' in kwargs: + kwargs['parallel_roe'].prescan_offset = n_pixels[0]+(subchip_size//2) + if 'serial_roe' in kwargs: + kwargs['serial_roe'].prescan_offset = n_pixels[1]+(subchip_size//2) + + #container for optimised values + self.optVals = {} + + #if they don't exist, make the sky background and readnoise images for the simulations + if (self.skyFrameSim is None) & (self.readnoiseFrameSim is None): + self._create_initial_sim_images(background_level,background_sigma,image_size=subchip_size) + + #create "simulation noise" covariance matrix (estimated from the sky+readnoise image) + #in this case, the fpr decrement is set to zero, because there is no CTI trailing + sim_matrix = self.covariance_matrix_from_image(self.skyFrameSim+self.readnoiseFrameSim,matrix_size=matrix_size,fprSize=5) + + #ideally, the matrix should be all zeros except for the centre square (== sigmaSky**2 + sigmaReadnoise**2) + #therefore, subtracting this from sim_matrix will provide the simulation noise covariance + sim_matrix[matrix_size//2,matrix_size//2] -= (self.sigmaRN**2+self.sigmaSky**2) + + self.optVals['pre-benchmark'] = sim_matrix + + #run S+R routine at 0% and 100% S+R fraction + result_array = np.array([]) + matrix_array = np.zeros((2,matrix_size,matrix_size)) + sr_frac = np.array([0.0,1.0]) + for frac in sr_frac: + print('\nestimating S+R covariance at %d percent level'%(100*frac)) + raw_matrix,correction_matrix = self._estimate_residual_covariance(self.skyFrameSim,self.readnoiseFrameSim,frac,matrix_size=matrix_size,**self.arcKwargs,**kwargs) + matrix_array[int(frac)] = correction_matrix + result = self.figure_of_merit(correction_matrix,fom_method=fom_method) #TODO: check to see if changing subgrid_size actually matters + result_array = np.append(result_array,result) + #interpolate between the results to determine the point of minimum residual covariance + a_fit,cov=curve_fit(self._fitter_function,sr_frac,result_array,absolute_sigma=True) + m = a_fit[1] + b = a_fit[0] + xint = -b/m #this intercept will be the optimised S+R fraction (where mean covariance=0) + + #store final optimised fraction as a class variable + self.SRfrac_optimised = xint + #self.optVals['optSRfrac'] = xint ############### ############### - def estimate_residual_covariance( - self, - background_level, - # Parallel - parallel_ccd=None, #can we do this with **kwargs ? I'm not sure in python - parallel_roe=None, - parallel_traps=None, - parallel_express=0, - parallel_window_offset=0, - parallel_window_start=0, - parallel_window_stop=-1, - parallel_time_start=0, - parallel_time_stop=-1, - parallel_prune_n_electrons=1e-10, - parallel_prune_frequency=20, - # Serial - if self.smoothCol: #should be more agnostic about direction or rows/cols (for variable names...) - serial_ccd=None, - serial_roe=None, - serial_traps=None, - serial_express=0, - serial_window_offset=0, - serial_window_start=0, - serial_window_stop=-1, - serial_time_start=0, - serial_time_stop=-1, - serial_prune_n_electrons=1e-10, - serial_prune_frequency=20, - # Pixel bounce - pixel_bounce=None, - # Output - verbosity=1, + def optimise_SR_fraction_from_image( + self, + image, + subchip_size=500, #size of the cutout region, to test CTI behaviour (a subset of pixels far form the readout) + matrix_size=5, #covariance matrix size (NxN array) + fom_method='box', + **kwargs ): - - raise NotImplementedError + ''' + A wrapper function for running the S+R optimisation routine using a real astronomical image rather than theoretical values. + This routine still calls but the sky paramerters and image size are determined from the data itself. + + Parameters + ---------- + image : numpy array + 2D numpy array containing the pixel values of the image + + subchip_size : Float + size of the cutout region, to test CTI behaviour in different regimes (e.g. far from readout, close to readout, etc.) + + matrix_size : Int + Size of covariance matrix used to estimate covariance. By default the program creates a 5x5 matrix + fom_method : String + Method used to estimate the mean value of a covariance matrix, for figure-of-merit purposes. Can be: + box : take a square (3x3) region around the central pixel [DEFAULT] + row : take all cells along the central row + column : take all cells along the central column + + **kwargs : parameters that characterise CTI features in arCTIc (trap density, CCD, read out electronics (ROE) descriptions) + + ''' + + image_shape = image.shape + sky_level = np.median(image) + sky_back = np.sqrt(sky_level) #assume pure shot noise for now + + self.optimise_SR_fraction(sky_level,sky_back,image_shape,subchip_size,matrix_size,fom_method,**kwargs) + ############### + ############### + def calculate_covariance_corners_from_image( + self, + image, + ): + raise NotImplementedError ############### - ############### - def optimise_parameters( - self, - image, - background_level, - #figure_of_merit=figure_of_merit(), # should probably pass a function here - **kwargs - ## Parallel - #parallel_ccd=None, #again, should do this with **kwargs - #parallel_roe=None, - #parallel_traps=None, - #parallel_express=0, - #parallel_window_offset=0, - #parallel_window_start=0, - #parallel_window_stop=-1, - #parallel_time_start=0, - #parallel_time_stop=-1, - #parallel_prune_n_electrons=1e-10, - #parallel_prune_frequency=20, - ## Serial - #serial_ccd=None, - #serial_roe=None, - #serial_traps=None, - #serial_express=0, - #serial_window_offset=0, - #serial_window_start=0, - #serial_window_stop=-1, - #serial_time_start=0, - #serial_time_stop=-1, - #serial_prune_n_electrons=1e-10, - #serial_prune_frequency=20, - ## Pixel bounce - #pixel_bounce=None, - ## Output - #verbosity=1, + ############### + def plot_matrix( + self, + covariance_matrix, + supressCentralPix=True, + title='MyPlotMatrix', + clearFrame=True, + **kwargs ): - #things to compare over iterations - outputFrame,noiseFrame = self.estimate_read_noise_model_from_image(image) - covar = self.covariance_matrix_from_image(image) - raise NotImplementedError + if supressCentralPix == True: + matrix = covariance_matrix.copy() + matrix_size = matrix.shape[0] + matrix[matrix_size//2,matrix_size//2] = 0 # set central pixel to zero, to better see correlation dynamic range + else: + matrix = covariance_matrix.copy() + if clearFrame: + mpl.pyplot.clf() + mpl.pyplot.imshow(matrix,cmap='bwr',**kwargs) + mpl.pyplot.colorbar() + mpl.pyplot.title(title) + mpl.pyplot.show() + ############### ############### + def plot_optimised_covariance_matrices( + self + ): + mpl.pyplot.figure() + mpl.pyplot.subplot(221) + self.plot_matrix(self.optVals['noCTIQuadMatrix_diff'],title='Readout Corner',clearFrame=False,vmin=-0.01,vmax=0.01) + mpl.pyplot.subplot(222) + self.plot_matrix(self.optVals['serialQuadMatrix_diff'],title='Serial Corner',clearFrame=False,vmin=-0.01,vmax=0.01) + mpl.pyplot.subplot(223) + self.plot_matrix(self.optVals['parallelQuadMatrix_diff'],title='Parallel Corner',clearFrame=False,vmin=-0.01,vmax=0.01) + mpl.pyplot.subplot(224) + self.plot_matrix(self.optVals['combinedQuadMatrix_diff'],title='Combined Corner',clearFrame=False,vmin=-0.01,vmax=0.01) ############### - def figure_of_merit(self,covariance_matrix): - xlen = covariance_matrix.shape[1]//2 - ylen = covariance matrix.shape[0]//2 - fullsum = np.sum(covariance_matrix) #sum over all cells - rowsum = np.sum(covariance_matrix[ylen,:]) #sum over central row (for serial CTI) - colsum = np.sum(covariance_matrix[:,xlen]) #sum over central column (for parallel CTI) - raise NotImplementedError - return fullsum,rowsum,colsum + ############### + # + #BACKGROUND FUNCTIONS + # ############### - def covariance_matrix_from_image(self, image, n_pixels=5): - raise NotImplementedError - covariance_matrix = np.zeros(n_pixels,n_pixels) + def _determine_noise_model(self, imageIn: np.ndarray, imageOut: np.ndarray): + ''' + Method for estimating readnoise on image (in S+R fashion) + assumes parallel+serial CTI trailing by default + ''' + + # the following line returns the result from the C implementation instead (comment out to switch back to native python code) + return determine_noise_model_c(imageIn, imageOut, self.sigmaRN, self.ampScale, self.smoothCol) + + #set target read noise amplitude to be equal to the value specified in the input function + readNoiseAmp = self.sigmaRN + + ''' + Define up the "comparison" variables for each S+R realisation. In every run, the current value of each pixel in the array is compared + to its neighbours; these comparisons will determine how much the pixel value will change in the next interation of the loop. There are + six such comparisons, some of which want the pixel to keep its current value, and others that want it to change. The result is a tug-of-war + between different modifiers, preventing the pixel value from varying by too much too quickly. + + The comparisons are: + + dval0 -- compare the pixel to itself. This comparison tries to keep the pixel value unchanged + + dval9 -- compare the pixel to its "local" average (the block of pixels surrounding it). In most + cases this is a (3x3) block, but pixels at the edge of the array are appropriately truncated. + Once again, this comparison tries to keep the pixel value unchanged. + + dmod1 and dmod2 -- compare pixel to its upper (dmod1) and lower (dmod2) row reighbours. + These comparisons try to bring the pixel closer to the average of its neighbours + + cmod1 and cmod2 -- compare pixel to its right (cmod1) and left (cmod2) column reighbours. + These comparisons try to bring the pixel closer to the average of its neighbours + + Finally, each comparison has a "clipped" version (e.g., dval0u) which truncates the comparison difference + to be within some upper/lower limits. These clipped values will determine the actual pixel modification. + + Note: all of these comparisons *should* be done with python broadcasting tricks...but I'm not sure this is actually the case + ''' + # initialize dval0 comparison + dval0 = imageIn - imageOut + + # set the clipped dval0 comparison + # note, this clipping is more stringent than the others + dval0u = np.clip(dval0, -1, 1) + + # initialize dval9 value and count (average will be (summed value)/(summed count)) + dval9 = dval0.copy() + + # do the dval9 calculation + dval9[:-1, :-1] += dval0[1:, 1:] # comparison with bottom-left neighbour + dval9[:-1, :] += dval0[1: ,:] # comparison with bottom-central neighbour + dval9[:-1, 1:] += dval0[1:, :-1] # comparison with bottom-right neighbour + dval9[:, :-1] += dval0[:, 1:] # comparison with middle-left neighbour + dval9[:, 1:] += dval0[:, :-1] + dval9[1:, :-1] += dval0[:-1, 1:] + dval9[1:, :] += dval0[:-1, :] + dval9[1:, 1:] += dval0[:-1, :-1] + + # get the dval9 average + dval9[1:-1,1:-1] /= 9 + # edges (excl. corners) + dval9[1:-1,0] /= 6 + dval9[1:-1,-1] /= 6 + dval9[0,1:-1] /= 6 + dval9[-1,1:-1] /= 6 + # corners + dval9[0, 0] /= 4 + dval9[-1, 0] /= 4 + dval9[0, -1] /= 4 + dval9[-1, -1] /= 4 + + # set the clipping modifier + # limit any difference to be at most the 1-sigma readnoise value (readNoiseAmp) scaled by a preset modifier (readNoiseAmpFraction) -- currently hard-coded to 0.2 + readNoiseAmpFraction = self.ampScale + mod_clip = readNoiseAmp * readNoiseAmpFraction + + # set the clipped dval9 comparison + dval9u = np.clip(dval9, -mod_clip, mod_clip) + + # initialise and set the dmod comparisons + dmod1 = np.zeros(imageIn.shape) + dmod1[1:,:] = imageOut[:-1,:] - imageOut[1:,:] + dmod2 = np.zeros(imageIn.shape) + dmod2[:-1,:] = imageOut[1:,:] - imageOut[:-1,:] + + # if specified, also initialise and set the cmod comparions + if self.smoothCol: + cmod1 = np.zeros(imageIn.shape) + cmod1[:,1:] = imageOut[:,:-1] - imageOut[:,1:] + cmod2 = np.zeros(imageIn.shape) + cmod2[:,:-1] = imageOut[:,1:] - imageOut[:,:-1] + + # set the clipped dmod values + dmod1u = np.clip(dmod1, -mod_clip, mod_clip) + dmod2u = np.clip(dmod2, -mod_clip, mod_clip) + + # if specified, also set the clipped cmod values + if self.smoothCol: + cmod1u = np.clip(cmod1, -mod_clip, mod_clip) + cmod2u = np.clip(cmod2, -mod_clip, mod_clip) + + # calulate weight parameters for each comparison (these were taken from the STScI code) + readNoiseAmp2 = readNoiseAmp**2 + dval0u *= np.square(dval0) / (np.square(dval0) + 4.0 * readNoiseAmp2) + dval9u *= np.square(dval9) / (np.square(dval9) + 18.0 * readNoiseAmp2) + dmod1u *= 4 * readNoiseAmp2 / (np.square(dmod1) + 4.0 * readNoiseAmp2) + dmod2u *= 4 * readNoiseAmp2 / (np.square(dmod2) + 4.0 * readNoiseAmp2) + if self.smoothCol: + cmod1u *= 4 * readNoiseAmp2 / (np.square(cmod1) + 4.0 * readNoiseAmp2) + cmod2u *= 4 * readNoiseAmp2 / (np.square(cmod2) + 4.0 * readNoiseAmp2) + + # return the appropriately-modified array (which can be sent to the next S+R iteration) + if self.smoothCol: + return (dval0u + dval9u + dmod1u + dmod2u + cmod1u + cmod2u) / 6 + else: + return (dval0u + dval9u + dmod1u + dmod2u) / 4 + + ############### + ############### + def _create_initial_sim_images( + self, + background_level, + background_sigma, + image_size=500, + overwrite=False #figure out how to implement this (it fails currently) + ): + ''' + A subroutine to generate sky_background and readnoise frames for S+R simulations + + These components will create a baseline "clean" image (sky + readnoise) that we + would expect to observe in the absence of CTI smearing effects + in each simulated corner. Keeping this fixed will control random noise fluctuations + ("simulation noise") when estimating covariances throughout the simulation. + + Note: This routine will be normally be called in + rather than through direct user input + ''' + + #create sky frame + skyFrame = np.random.normal(background_level, background_sigma, (image_size,image_size)) + self.skyFrameSim = skyFrame + + #create readnoise_frame + noiseFrame = np.random.normal(0, self.sigmaRN, (image_size,image_size)) + self.readnoiseFrameSim = noiseFrame + + #generate covariance matrix from clean image + + + return + ############### + ############### + def _estimate_residual_covariance( + self, + skyFrame, + noiseFrame, + sr_fraction, + matrix_size=5, #size of covariance grid, in pixels (also assumes square) + **kwargs + ): + ''' + Estimate a covariance matrix for an optimising simulation, assuming the sky/noise frames have already been created + ''' + #add CTI effects to skyFrame + ctiTrailedFrame = add_cti(skyFrame,**kwargs) + #add read noise to CTI-trailed image + readNoiseAddedFrame = ctiTrailedFrame + noiseFrame + #do S+R routine + sFrame, rFrame = self.generate_SR_frames_from_image(readNoiseAddedFrame,sr_fraction) + ###rescale S+R correction to a fractional value (for optimization) + ##sFrame+=rFrame + ##rFrame*=sr_fraction + ##sFrame-=rFrame + #clean CTI from S frame + ctiCorrectedFrame = remove_cti(sFrame,1,**kwargs) + #re-add R frame to correction + outputFrame = ctiCorrectedFrame + rFrame + + #determine covariance FOMs + fomSR = self.covariance_matrix_from_image(outputFrame,matrix_size=matrix_size) + fomBenchmark = self.covariance_matrix_from_image(skyFrame+noiseFrame,matrix_size=matrix_size) + fomDiff = fomSR - fomBenchmark + + return fomSR,fomDiff + ############### + ############### + def _fitter_function(self,x,intercept,slope): + ''' + simple functional form to optimize a (linear) S+R fit using covariance matrix data + passed to python optiser in optimise_SR_fraction, not needed by the user + ''' + y = intercept + slope * x + return y + ############### + ############### + def figure_of_merit(self, covariance_matrix, fom_method='box',subgrid_size=None): + ''' + function for estimating the covariance figure of merit: the average value over some section of a covariance matrix + if method == 'box' take a square region around the central pixel; 'row' is along the central row; column is along the central column + if subgrid size is selected, use only an NxN subset of pixels in the matrix, starting from the center outward + ''' + xlen = covariance_matrix.shape[1] + ylen = covariance_matrix.shape[0] + xcenter = xlen//2 + ycenter = ylen//2 + xstart,xend = 0,xlen + ystart,yend = 0,ylen + if subgrid_size is not None: + xdiff = (xlen-subgrid_size)//2 + ydiff = (ylen-subgrid_size)//2 + xstart+=xdiff + xend-=xdiff + ystart+=ydiff + yend-=ydiff + + covariance_matrix[ycenter,xcenter] = 0 #supress autocorrelation pixel + + if fom_method == 'box': + result_mean = (np.sum(covariance_matrix[ystart:yend,xstart:xend]) - covariance_matrix[ycenter,xcenter])/8 #average over central box + if fom_method == 'row': + result_mean = (np.sum(covariance_matrix[ycenter,xstart:xend]) - covariance_matrix[ycenter,xcenter])/4 #average over central row (for serial CTI) + if fom_method == 'column': + result_mean = (np.sum(covariance_matrix[ystart:yend,xcenter]) - covariance_matrix[ycenter,xcenter])/4 #sum over central column (for parallel CTI) + #raise NotImplementedError + return result_mean + ############### + ############### + def covariance_matrix_from_image(self, image, matrix_size=5,fprSize=5): + + matRange= matrix_size//2 #set positions of correlation pixels and set limits to remove "roll-over" region + + covariance_matrix = np.zeros((matrix_size,matrix_size)) #calcluate mean stats on image - x = image.flatten() + image2 = image[fprSize:,fprSize:] #remove FPR decrement + x = image2[matRange:-matRange,matRange:-matRange].flatten() #remove possible roll-over region xbar = np.mean(x) #roll image to get cross correlation - matRange= n_pixels//2 for i in range(-matRange,matRange+1,1): for j in range(-matRange,matRange+1,1): - y = np.roll(img,(-j,-i),axis=(0,1)).flatten() + y = np.roll(image2,(-j,-i),axis=(0,1))[2:-2,2:-2].flatten() ybar = np.mean(y) #calculate covariance - covar = np.sum((x-xbar)*(y-ybar))/(pl.std(x-xbar)*pl.std(y-ybar))/(x.size) - + covar = np.sum((x-xbar)*(y-ybar))/x.size #/(np.std(x-xbar)*np.std(y-ybar))/(x.size) #removing noise-level normalization + #covar = np.sum((x-xbar)*(y-ybar))/(np.std(x-xbar)*np.std(y-ybar))/(x.size) #populate matrix cells - covariance_matrix[i+2,j+2] = covar - + covariance_matrix[i+matRange,j+matRange] = covar + #covariance_matrix[matRange,matRange] = 0 #switch this normalization to plotting + return covariance_matrix + ############### + + + + + ############### + # + # FUNCTIONS TO BE REFORMATTED / RETHOUGHT? (WORK IN PROGRESS) + # + ############### + def measure_simulated_covariance_corners( + self, + n_pixels=None, #the size of the entire CCD over which the procedure will be run + subchip_size=500, #size of the cutout region, to test CTI behaviour (a subset of pixels far form the readout) + matrix_size=5, #covariance matrix size (NxN array) + fom_method='box', + **kwargs + ): + ''' + Apply optimised S+R correction to the four corners of the CCD (close to readout, far from parallel register, far from serial register, far from both) + and estimate a covariance matrix for each corner. These matrices can be used to correct any residual covariance in weak lensing shapes. + ''' + + #update any keyword arguments + allKwargs = kwargs|self.arcKwargs + + frac_opt = self.SRfrac_optimised #self.optVals['optSRfrac'] #set the optimized S+R fraction + if frac_opt == None: + frac_opt = 1.0 + + if type(n_pixels) == int: + n_pixels=(n_pixels,n_pixels) + chip_size = n_pixels + + ##set prescan regions (this simulates the maximum possible CTI trailing, averaged at the far edge of the real CCD) + #if 'parallel_roe' in allKwargs: + # kwargs['parallel_roe'].prescan_offset = n_pixels[0]+(subchip_size//2) + #if 'serial_roe' in allKwargs: + # kwargs['serial_roe'].prescan_offset = n_pixels[1]+(subchip_size//2) + + #re-run the S+R routine with the optimised level + raw_matrix_opt,correction_matrix_opt = self._estimate_residual_covariance(self.skyFrameSim,self.readnoiseFrameSim,frac_opt,matrix_size=matrix_size,**self.arcKwargs,**kwargs) + raw_matrix_opt_noSR,correction_matrix_opt_noSR = self._estimate_residual_covariance(self.skyFrameSim,self.readnoiseFrameSim,0,matrix_size=matrix_size,**self.arcKwargs,**kwargs) #this is dumb...do it better + + #run a second S+R routine in the region close to the readout registers, to identify/remove stochastic simulation noise + kwargs_noCTI = kwargs|self.arcKwargs + kwargs_noCTI['parallel_traps']=None + kwargs_noCTI['serial_traps']=None + #effectively, we are removing the effects of CTI, but still calculating a covariance matrix + raw_matrix_opt_noCTI,correction_matrix_opt_noCTI = self._estimate_residual_covariance(self.skyFrameSim,self.readnoiseFrameSim,frac_opt,matrix_size=matrix_size,**kwargs_noCTI) + raw_matrix_opt_noCTI_noSR,correction_matrix_opt_noCTI_noSR = self._estimate_residual_covariance(self.skyFrameSim,self.readnoiseFrameSim,0,matrix_size=matrix_size,**kwargs_noCTI) + + #create the "benchmark" covariance matrix: a purely non-correlated array with the central pixel equal to (readnoise_sigma**2 + background_sigma**2) + #Any difference between this and "raw_matrix_opt_noCTI" is the simulation noise, and can be eliminated + benchmark_matrix = raw_matrix_opt_noCTI*0 + centrePix = matrix_size//2 + benchmark_matrix[centrePix,centrePix] = self.sigmaRN**2 + self.sigmaSky**2 + #calculate the simulation-noise residual + residual_benchmark_matrix = raw_matrix_opt_noCTI - benchmark_matrix + + #Finally, if both parallel and serial CTI are present, calculate covariances of the "cross-region" chip corners (far from serial readout but close to parallel readout, and vice-versa) + if ('parallel_traps' in allKwargs)&('serial_traps' in allKwargs): + if (allKwargs['parallel_traps'] is not None)&(allKwargs['serial_traps'] is not None): + kwargs_parallelOnly = kwargs|self.arcKwargs + kwargs_parallelOnly['serial_traps']=None #make a parallel-only corner + kwargs_serialOnly = kwargs|self.arcKwargs + kwargs_serialOnly['parallel_traps']=None #make a serial-only corner + + raw_matrix_opt_parallel,correction_matrix_opt_parallel = self._estimate_residual_covariance(self.skyFrameSim,self.readnoiseFrameSim,frac_opt,matrix_size=matrix_size,**kwargs_parallelOnly) + raw_matrix_opt_serial,correction_matrix_opt_serial = self._estimate_residual_covariance(self.skyFrameSim,self.readnoiseFrameSim,frac_opt,matrix_size=matrix_size,**kwargs_serialOnly) + raw_matrix_opt_parallel_noSR,correction_matrix_opt_parallel_noSR = self._estimate_residual_covariance(self.skyFrameSim,self.readnoiseFrameSim,0,matrix_size=matrix_size,**kwargs_parallelOnly) + raw_matrix_opt_serial_noSR,correction_matrix_opt_serial_noSR = self._estimate_residual_covariance(self.skyFrameSim,self.readnoiseFrameSim,0,matrix_size=matrix_size,**kwargs_serialOnly) + #package all relevant quantities into optVals and return + self.optVals['combinedQuadMatrix'] = raw_matrix_opt + self.optVals['combinedQuadMatrix_diff'] = correction_matrix_opt + self.optVals['combinedQuadMatrix_noSR'] = raw_matrix_opt_noSR + self.optVals['combinedQuadMatrix_noSR_diff'] = correction_matrix_opt_noSR + self.optVals['parallelQuadMatrix'] = raw_matrix_opt_parallel + self.optVals['parallelQuadMatrix_diff'] = correction_matrix_opt_parallel + self.optVals['parallelQuadMatrix_noSR'] = raw_matrix_opt_parallel_noSR + self.optVals['parallelQuadMatrix_noSR_diff'] = correction_matrix_opt_parallel_noSR + self.optVals['serialQuadMatrix'] = raw_matrix_opt_serial + self.optVals['serialQuadMatrix_diff'] = correction_matrix_opt_serial + self.optVals['serialQuadMatrix_noSR'] = raw_matrix_opt_serial_noSR + self.optVals['serialQuadMatrix_noSR_diff'] = correction_matrix_opt_serial_noSR + self.optVals['noCTIQuadMatrix'] = raw_matrix_opt_noCTI + self.optVals['noCTIQuadMatrix_diff'] = correction_matrix_opt_noCTI + self.optVals['noCTIQuadMatrix_noSR'] = raw_matrix_opt_noCTI + self.optVals['noCTIQuadMatrix_diff_noSR'] = correction_matrix_opt_noCTI + self.optVals['benchmark'] = benchmark_matrix + self.optVals['benchmark_resid'] = residual_benchmark_matrix + ############### diff --git a/python/arcticpy/read_noise_wrapper.pyx b/python/arcticpy/read_noise_wrapper.pyx new file mode 100644 index 0000000..a87e914 --- /dev/null +++ b/python/arcticpy/read_noise_wrapper.pyx @@ -0,0 +1,12 @@ +import numpy as np +cimport numpy as np +cimport cython + +cdef extern void determine_noise_model(const double* imageIn, const double* imageOut, const int rows, const int cols, const double readNoiseAmp, const double readNoiseAmpFraction, const int smoothCol, double* output); + +def determine_noise_model_c(np.ndarray[np.float64_t, ndim=2] arr1, np.ndarray[np.float64_t, ndim=2] arr2, readNoiseAmp, readNoiseAmpFraction, smoothCol): + cdef int rows = arr1.shape[0] + cdef int cols = arr1.shape[1] + cdef np.ndarray[np.float64_t, ndim=2] out = np.zeros((rows, cols), dtype=np.float64) + determine_noise_model(&arr1[0,0], &arr2[0,0], rows, cols, readNoiseAmp, readNoiseAmpFraction, smoothCol, &out[0,0]) + return out diff --git a/setup.py b/setup.py index bbbdf5c..ba5e088 100644 --- a/setup.py +++ b/setup.py @@ -32,8 +32,8 @@ os.remove(file) # Build -if "CC" not in os.environ: - os.environ["CC"] = "g++" +#if "CC" not in os.environ: +# os.environ["CC"] = "g++" ext_headers = [os.path.join(dir_include, header) for header in os.listdir(dir_include)] ext_sources = [os.path.join(dir_src, src) for src in os.listdir(dir_src)] @@ -46,9 +46,17 @@ libraries=["gsl"], runtime_library_dirs=[dir_gsl_lib], include_dirs=[dir_wrapper_include, dir_include, np.get_include(), dir_gsl_include], - extra_compile_args=["-std=c++17", "-O3"], + extra_compile_args=["-std=c++17", "-O3", "-fopenmp"], + extra_link_args=["-fopenmp"], define_macros=[('NPY_NO_DEPRECATED_API', 0)], ), + Extension( + 'read_noise_c', + sources=[os.path.join(dir_wrapper, "read_noise_wrapper.pyx"), os.path.join(dir_wrapper, "read_noise.c")], + include_dirs=[np.get_include()], + extra_compile_args=["-O3", "-fopenmp"], + extra_link_args=["-fopenmp"] + ), ] setup( @@ -56,5 +64,7 @@ extensions, compiler_directives={"language_level": "3"} ), + packages=['arcticpy'], + package_dir={'arcticpy': dir_wrapper}, headers=ext_headers, # currently ignored (?) ) diff --git a/test/test_arcticpy.py b/test/test_arcticpy.py index bd817d0..aee21d5 100644 --- a/test/test_arcticpy.py +++ b/test/test_arcticpy.py @@ -1091,3 +1091,31 @@ def print_help(): print_help() except IndexError: print_help() + +class TestReadNoise: + + def test_noiseEstimation(): + #make simple "sky" and "noise" frames and combine them + sky = np.array([[104.62689576, 96.26344068, 102.05736548, 76.16124587, 107.55689551], + [ 96.87789779, 87.37248667, 110.70696472, 112.3185982 , 84.7109744 ], + [104.06312714, 99.2960313 , 105.16551541, 97.52700805, 104.66187477], + [100.90670311, 108.95078214, 100.96009345, 100.06043671, 92.82319682], + [ 87.42418643, 109.87867837, 108.59791761, 102.64102994, 77.73118275]]) + noise = np.array([[ 8.21416447, -4.40936237, -1.18755795, 0.06791093, -2.33996983], + [ 4.93207717, 6.17463103, -3.60709428, 7.50694942, -2.27221275], + [ 0.63142708, 4.19195127, 4.99175746, 6.83683726, -6.91249208], + [-1.50630185, 1.91698306, 6.03575168, 11.84120206, 1.78779769], + [ 0.21060601, 4.81669569, 2.89002324, -2.14453015, -1.87141227]]) + combo = sky+noise + + #create a read_noise object with sigma_noise == 4.5 electrons (matches input noise model) + testRN = cti.read_noise.ReadNoise(4.5) + + #separate combination into S ("smooth") + R ("readnoise") components + S,R = testRN.generate_SR_frames_from_image(combo) + + #check that S+R components are close to input images (within typical simulation noise limits) + assert abs(np.mean(R) - np.mean(noise)) < 2 + assert abs(np.std(R) - np.std(noise)) <= 0.001 + assert abs(np.mean(S) - np.mean(sky)) < 2 + assert abs(np.std(S) - np.std(sky)) < 1