Skip to content

Commit

Permalink
Merge pull request #36 from LSSTDESC/no_analysis
Browse files Browse the repository at this point in the history 
No analysis (#35)
  • Loading branch information
@fjaviersanchez
fjaviersanchez authored Mar 28, 2019
2 parents ef05d96 + af726e1 commit 73cf852
Show file tree
Hide file tree
Showing 4 changed files with 111 additions and 94 deletions.
9 changes: 8 additions & 1 deletion descwl/analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -784,7 +784,7 @@ def residuals(p):
bestfit_values[i,5] = parameters['dg2_%d'%i].value
return bestfit_values

def finalize(self,verbose,trace,calculate_bias):
def finalize(self,verbose,trace,calculate_bias,no_analysis):
"""Finalize analysis of all added stars and galaxies.
Args:
Expand Down Expand Up @@ -867,6 +867,13 @@ def finalize(self,verbose,trace,calculate_bias):
])

data = np.empty(num_galaxies, dtype=dtype)
if no_analysis:
# Return empty data table
table = astropy.table.Table(data, copy=False)
num_slices, h, w = self.stamps[0].shape
results = OverlapResults(self.survey, table, self.stamps,
self.bounds, num_slices)
return results

trace('allocated table of %ld bytes for %d galaxies' % (data.nbytes,num_galaxies))

Expand Down
129 changes: 67 additions & 62 deletions descwl/render.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -198,7 +198,8 @@ def description(self):
('PSF dilution factor is %.6f.' % self.psf_dilution)
])

def render_galaxy(self,galaxy,no_partials = False, calculate_bias = False):
def render_galaxy(self,galaxy,no_partials = False, calculate_bias = False,
no_analysis=False):
"""Render a galaxy model for a simulated survey.
Args:
Expand Down Expand Up @@ -281,69 +282,73 @@ def render_galaxy(self,galaxy,no_partials = False, calculate_bias = False):
if survey_overlap.area() == 0:
raise SourceNotVisible
self.survey.image[survey_overlap] += cropped_stamp[survey_overlap]
if not no_analysis:
# Give this Galaxy its own GalaxyRenderer.
galaxy.renderer = GalaxyRenderer(galaxy,cropped_stamp,self.survey)

# Define the parameter variations we consider for building Fisher matrices.
# The names appearing below are args of Galaxy.get_transformed_model().
# We do not include 'flux' below since the nominal image is already the
# partial derivative wrt flux (after dividing by flux).
variations = [
('dx',self.survey.pixel_scale/3.), # arcsecs
('dy',self.survey.pixel_scale/3.), # arcsecs
('ds',0.05), # relative dilation (flux preserving)
('dg1',0.03), # + shear using |g| = (a-b)/(a+b) convention
('dg2',0.03), # x shear using |g| = (a-b)/(a+b) convention
]

# Prepare the datacube that we will return.
ncube = 1

# Give this Galaxy its own GalaxyRenderer.
galaxy.renderer = GalaxyRenderer(galaxy,cropped_stamp,self.survey)

# Define the parameter variations we consider for building Fisher matrices.
# The names appearing below are args of Galaxy.get_transformed_model().
# We do not include 'flux' below since the nominal image is already the
# partial derivative wrt flux (after dividing by flux).
variations = [
('dx',self.survey.pixel_scale/3.), # arcsecs
('dy',self.survey.pixel_scale/3.), # arcsecs
('ds',0.05), # relative dilation (flux preserving)
('dg1',0.03), # + shear using |g| = (a-b)/(a+b) convention
('dg2',0.03), # x shear using |g| = (a-b)/(a+b) convention
]

# Prepare the datacube that we will return.
ncube = 1

positions = descwl.analysis.make_positions()
if not no_partials:
ncube = 1+len(variations)
if calculate_bias:
#15 is number of second partials for 5 parameters (no flux).
ncube = 1 + len(variations) + 15

height,width = cropped_stamp.array.shape
datacube = np.empty((ncube,height,width))
datacube[0] = cropped_stamp.array #flux partial is the same.

#calculate partials, if requested.
if not no_partials:
#The nominal image doubles as the flux partial derivative.
for i,(pname_i,delta_i) in enumerate(variations):
variation_stamp = (galaxy.renderer.draw(**{pname_i: +delta_i}).copy() -
galaxy.renderer.draw(**{pname_i: -delta_i}))
datacube[positions[pname_i]] = variation_stamp.array/(2*delta_i)

#calculate second partials, if requested.
positions = descwl.analysis.make_positions()
if not no_partials:
ncube = 1+len(variations)
if calculate_bias:
for j,(pname_j,delta_j) in enumerate(variations):
if(i==j):
galaxy_2iup = galaxy.renderer.draw(**{pname_i: +2*delta_i}).copy()
galaxy_2idown = galaxy.renderer.draw(**{pname_i: -2*delta_i}).copy()
variation_i_i = galaxy_2iup + galaxy_2idown - 2*galaxy.renderer.draw()
datacube[positions[pname_i,pname_i]] = ((variation_i_i).array/
(2*delta_i)**2)

elif(j>i):
galaxy_iup_jup = galaxy.renderer.draw(**{pname_i: +delta_i,
pname_j: +delta_j}).copy()
galaxy_iup_jdown = galaxy.renderer.draw(**{pname_i: +delta_i,
pname_j: -delta_j}).copy()
galaxy_idown_jup = galaxy.renderer.draw(**{pname_i: -delta_i,
pname_j: +delta_j}).copy()
galaxy_idown_jdown = galaxy.renderer.draw(**{pname_i: -delta_i,
pname_j: -delta_j})

variation_i_j = (galaxy_iup_jup + galaxy_idown_jdown -
galaxy_idown_jup - galaxy_iup_jdown)
datacube[positions[pname_i,pname_j]] = (variation_i_j.array /
(4*delta_i*delta_j))

#15 is number of second partials for 5 parameters (no flux).
ncube = 1 + len(variations) + 15

height,width = cropped_stamp.array.shape
datacube = np.empty((ncube,height,width))
datacube[0] = cropped_stamp.array #flux partial is the same.

#calculate partials, if requested.
if not no_partials:
#The nominal image doubles as the flux partial derivative.
for i,(pname_i,delta_i) in enumerate(variations):
variation_stamp = (galaxy.renderer.draw(**{pname_i: +delta_i}).copy() -
galaxy.renderer.draw(**{pname_i: -delta_i}))
datacube[positions[pname_i]] = variation_stamp.array/(2*delta_i)

#calculate second partials, if requested.
if calculate_bias:
for j,(pname_j,delta_j) in enumerate(variations):
if(i==j):
galaxy_2iup = galaxy.renderer.draw(**{pname_i: +2*delta_i}).copy()
galaxy_2idown = galaxy.renderer.draw(**{pname_i: -2*delta_i}).copy()
variation_i_i = galaxy_2iup + galaxy_2idown - 2*galaxy.renderer.draw()
datacube[positions[pname_i,pname_i]] = ((variation_i_i).array/
(2*delta_i)**2)

elif(j>i):
galaxy_iup_jup = galaxy.renderer.draw(**{pname_i: +delta_i,
pname_j: +delta_j}).copy()
galaxy_iup_jdown = galaxy.renderer.draw(**{pname_i: +delta_i,
pname_j: -delta_j}).copy()
galaxy_idown_jup = galaxy.renderer.draw(**{pname_i: -delta_i,
pname_j: +delta_j}).copy()
galaxy_idown_jdown = galaxy.renderer.draw(**{pname_i: -delta_i,
pname_j: -delta_j})

variation_i_j = (galaxy_iup_jup + galaxy_idown_jdown -
galaxy_idown_jup - galaxy_iup_jdown)
datacube[positions[pname_i,pname_j]] = (variation_i_j.array /
(4*delta_i*delta_j))
else:
if self.verbose_render:
print("No analysis performed")
height, width = cropped_stamp.array.shape
datacube = np.empty((1, height, width))
if self.verbose_render:
print('Rendered galaxy model for id = %d with z = %.3f' % (
galaxy.identifier,galaxy.redshift))
Expand Down
58 changes: 30 additions & 28 deletions descwl/survey.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ class Survey(object):
Raises:
RuntimeError: Missing or extra arguments provided or unable to calculate PSF size.
"""
def __init__(self,**args):
def __init__(self, no_analysis=False, **args):
if set(args.keys()) != set(Survey._parameter_names):
raise RuntimeError('Missing or extra arguments provided to Survey constructor.')
self.args = args
Expand Down Expand Up @@ -79,32 +79,33 @@ def __init__(self,**args):
psf_size_pixels = 2*int(math.ceil(10*atmospheric_psf_fwhm/self.pixel_scale))
self.psf_image = galsim.Image(psf_size_pixels,psf_size_pixels,scale = self.pixel_scale)
self.psf_model.drawImage(image = self.psf_image)
# Draw a (temporary) high-resolution (10x) image covering the same area.
zoom = 10
hires_psf_image = galsim.Image(zoom*psf_size_pixels,zoom*psf_size_pixels,scale = self.pixel_scale/zoom)
self.psf_model.drawImage(image = hires_psf_image)
# Calculate the unweighted second moments in arcsecs**2 of the hi-res PSF image.
hires_sum = np.sum(hires_psf_image.array)
hires_grid = (self.pixel_scale/zoom)*(np.arange(zoom*psf_size_pixels) - 0.5*zoom*psf_size_pixels + 0.5)
hires_x,hires_y = np.meshgrid(hires_grid,hires_grid)
psf_x = np.sum(hires_psf_image.array*hires_x)/hires_sum
psf_y = np.sum(hires_psf_image.array*hires_x)/hires_sum
hires_x -= psf_x
hires_y -= psf_y
psf_xx = np.sum(hires_psf_image.array*hires_x**2)/hires_sum
psf_xy = np.sum(hires_psf_image.array*hires_x*hires_y)/hires_sum
psf_yy = np.sum(hires_psf_image.array*hires_y**2)/hires_sum
self.psf_second_moments = np.array(((psf_xx,psf_xy),(psf_xy,psf_yy)))
# Calculate the corresponding PSF sizes |Q|**0.25 and (0.5*trQ)**0.5
self.psf_sigma_m = np.power(np.linalg.det(self.psf_second_moments),0.25)
self.psf_sigma_p = np.sqrt(0.5*np.trace(self.psf_second_moments))
# Also calculate the PSF size as |Q|**0.25 using adaptive weighted second moments
# of the non-hires PSF image.
try:
hsm_results = galsim.hsm.FindAdaptiveMom(self.psf_image)
self.psf_size_hsm = hsm_results.moments_sigma*self.pixel_scale
except RuntimeError as e:
raise RuntimeError('Unable to calculate adaptive moments of PSF image.')
if not no_analysis:
# Draw a (temporary) high-resolution (10x) image covering the same area.
zoom = 10
hires_psf_image = galsim.Image(zoom*psf_size_pixels,zoom*psf_size_pixels,scale = self.pixel_scale/zoom)
self.psf_model.drawImage(image = hires_psf_image)
# Calculate the unweighted second moments in arcsecs**2 of the hi-res PSF image.
hires_sum = np.sum(hires_psf_image.array)
hires_grid = (self.pixel_scale/zoom)*(np.arange(zoom*psf_size_pixels) - 0.5*zoom*psf_size_pixels + 0.5)
hires_x,hires_y = np.meshgrid(hires_grid,hires_grid)
psf_x = np.sum(hires_psf_image.array*hires_x)/hires_sum
psf_y = np.sum(hires_psf_image.array*hires_x)/hires_sum
hires_x -= psf_x
hires_y -= psf_y
psf_xx = np.sum(hires_psf_image.array*hires_x**2)/hires_sum
psf_xy = np.sum(hires_psf_image.array*hires_x*hires_y)/hires_sum
psf_yy = np.sum(hires_psf_image.array*hires_y**2)/hires_sum
self.psf_second_moments = np.array(((psf_xx,psf_xy),(psf_xy,psf_yy)))
# Calculate the corresponding PSF sizes |Q|**0.25 and (0.5*trQ)**0.5
self.psf_sigma_m = np.power(np.linalg.det(self.psf_second_moments),0.25)
self.psf_sigma_p = np.sqrt(0.5*np.trace(self.psf_second_moments))
# Also calculate the PSF size as |Q|**0.25 using adaptive weighted second moments
# of the non-hires PSF image.
try:
hsm_results = galsim.hsm.FindAdaptiveMom(self.psf_image)
self.psf_size_hsm = hsm_results.moments_sigma*self.pixel_scale
except RuntimeError as e:
raise RuntimeError('Unable to calculate adaptive moments of PSF image.')
# Calculate the mean sky background level in detected electrons per pixel.
self.mean_sky_level = self.get_flux(self.sky_brightness)*self.pixel_scale**2
# Create an empty image using (0,0) to index the lower-left corner pixel.
Expand Down Expand Up @@ -478,4 +479,5 @@ def from_args(cls,args):
for parameter_name in ctor_params:
if parameter_name in args_dict and args_dict[parameter_name] is not None:
ctor_params[parameter_name] = args_dict[parameter_name]
return Survey(survey_name=survey_name,filter_band=filter_band,**ctor_params)
return Survey(no_analysis=args['no_analysis'], survey_name=survey_name,
filter_band=filter_band,**ctor_params)
9 changes: 6 additions & 3 deletions simulate.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@ def main():
parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('--verbose', action = 'store_true',
help = 'Provide verbose output.')
parser.add_argument('--no-analysis', action = 'store_true',
help = 'Don\'t run analysis.')
parser.add_argument('--survey-defaults', action = 'store_true',
help = 'Print survey camera and observing parameter defaults and exit.')
parser.add_argument('--memory-trace', action = 'store_true',
Expand Down Expand Up @@ -83,7 +85,9 @@ def main():

try:
galaxy = galaxy_builder.from_catalog(entry,dx,dy,survey.filter_band)
stamps,bounds = render_engine.render_galaxy(galaxy,args.no_partials,args.calculate_bias)
stamps, bounds = render_engine.render_galaxy(
galaxy, args.no_partials, args.calculate_bias,
args.no_analysis)
analyzer.add_galaxy(galaxy,stamps,bounds)
trace('render')

Expand All @@ -100,8 +104,7 @@ def main():

except (descwl.model.SourceNotVisible,descwl.render.SourceNotVisible):
pass

results = analyzer.finalize(args.verbose,trace,args.calculate_bias)
results = analyzer.finalize(args.verbose,trace,args.calculate_bias,args.no_analysis)
output.finalize(results,trace)

except RuntimeError as e:
Expand Down

0 comments on commit 73cf852

Please sign in to comment.