UNDER DEVELOPMENT

The software package is currently being dramatically improved and updated.
This is a (very) old version, potentially with bugs.
  Shapelets web page     Shapelets IDL code     Installation     Help pages     Simulated images  
List of routines
Instrutions for general code
Instructions for image simulation code

Source code for simage_resample_pdf.pro:

You can also view the help page for this routine.

;$Id: simage_resample_pdf.pro, v1.0$ ; ; Copyright © 2004 Richard Massey and Alexandre Refregier. ; ; This file is a part of the Shapelets analysis code. ; www.ast.cam.ac.uk/~rjm/shapelets/ ; ; The Shapelets code is free software; you can redistribute it and/or ; modify it under the terms of the GNU General Public License as published ; by the Free Software Foundation; either version 2 of the License, or ; (at your option) any later version. ; ; The Shapelets code is distributed in the hope that it will be useful, ; but WITHOUT ANY WARRANTY; without even the implied warranty of ; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ; GNU General Public License for more details. ; ; You should have received a copy of the GNU General Public License ; along with the Shapelets code; if not, write to the Free Software ; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ; ; ; ************************************************************************ ; ************************************************************************ ; ; NAME: ; SIMAGE_FLAG_PDF_TOPOLOGY ; ; CATEGORY: ; A component of simage_resample_pdf. ; ; PURPOSE: ; Flag which dimensions of a PDF (the phase dimensions if using polar ; shapelets) wrap around at 2pi, and which extend to infinity. ; ; INPUTS: ; n_max - Maximum n_max used in the shapelet catalogue. ; ; OPTIONAL INPUTS: ; none. ; ; KEYWORD PARAMETERS: ; /POLAR - Assume the PDF contains polar shapelet coefficients. ; ; OUTPUTS: ; topology - Vector containing a series of flags, one for each dimension ; of a PDF of shapelet coefficients: ; (0) Dimension is empty and should be skipped ; (1) Standard, open-ended topology ; (2) Wraps around at 2pi ; (3) Once a coefficient is sampled, collapse the dimension ; reorder - Order in which the dimensions should be done. ; (useful in old create.pro but now possibly redundant) ; ; MODIFICATION HISTORY: ; Mar 04 - Tidied up by RM. ; May 02 - Written by Richard Massey pro simage_flag_pdf_topology, n_max, topology, reorder, POLAR=polar COMPILE_OPT idl2, HIDDEN shapelets_make_nvec,n_max,n1,n2,n_coeffs topology = intarr(n_coeffs) if keyword_set(polar) then begin ; Polar shapelet are used in PDF - will need to ensure that the dimensions ; corresponding to phases wrap around, but the dimensions corresponding to ; the magnitudes extend to infinity. ; Modulus-type dimensions are open (topology #1) start_mods = 0 end_mods = n_elements(where(n1 ge n2 and (n1+n2) le n_max))-1 topology[start_mods:end_mods] = 1 ; Phase-type dimensions are closed (topology #2) start_phases = n_elements(where(n1 ge n2)) end_phases = start_phases+n_elements(where(n1 gt n2 and (n1+n2) le n_max))-1 topology[start_phases:end_phases] = 2 ; Re-order, so cut slice through phase dimensions in between mag dimensions reorder = intarr(n_coeffs) mods = indgen(end_mods-start_mods+1)+start_mods phases = indgen(end_phases-start_phases+1)+start_phases index = indgen(2*n_max+1)/4+1 either = fix(sin(indgen(2*n_max+1)*!pi/2)^2) counter = 0L & mcount = 0L & pcount = 0L for i=0,2*n_max do begin for j=0,index[i]-1 do begin if either[i] eq 0 then begin reorder[counter]=mods[mcount] mcount=mcount+1 endif else begin reorder[counter]=phases[pcount] pcount=pcount+1 endelse counter=counter+1 endfor endfor if counter lt n_coeffs then reorder[counter:*]=n_coeffs-2 ; Cope with additional dimensions of size and magnitude topology = [1,topology] reorder = [0,reorder+1] endif else begin ; Cartesian shapelets are used, so trivially reproduce the notation that will ; ensure all of the coefficients will be used, but the code can stay ; compatible with the polar shapelet PDF option. topology = [1,intarr(n_coeffs-1)] topology[0:(n_max+1)*(n_max+2)/2-1] = 1 reorder=indgen(n_coeffs+1) endelse end ; ; ************************************************************************ ; ************************************************************************ ; ; NAME: ; SIMAGE_PDF_KERNEL_WIDTH ; ; CATEGORY: ; A component of simage_resample_pdf. ; ; PURPOSE: ; Decides how much to smooth a discretely sampled PDF (ie the kernel width ; "h") for the Kernel method of Silverman (1984). Will be different in ; differnet directions and may even be adaptive (ie different for each ; sample point, as a function of the local density). ; ; INPUTS: ; pdf - an array in the form [#sample points,#dimensions] ; ; OPTIONAL INPUTS: ; topology - array containing a flag describing the topology of each ; dimension. Default: everything 1s. ; (0) dimension ignored (skip) ; (1) open - regular Cartesian space, [-inf,inf] ; (2) closed - eg phases wrap around at 2pi (but range need not ; be specified ; ; OPTIONAL INPUTS: ; none. ; ; KEYWORD PARAMETERS: ; /PAPER - Use the smoothing kernel verified by Massey et al. (2004) ; MNRAS 348, 214 to produce realistic new galaxy morphologies. ; /POLAR - Assume the PDF contains polar shapelet coefficients. ; /ADAPTATIVE: adaptive smoothing. Returned value will then be an array ; rather than just a vector. NOT CURRENTLY IMPLEMENTED! ; ; OUTPUTS: ; kernel_width - "h" in every dimension (non-adaptive) or for each sample ; point in every direction (adaptive). ; NOTES: ; A constant smoothing amount is calculated for each dimension. There are several ; options here. ; * Leave as is. Simple, but may not get the wings right and blurrs out cuspy peaks/ ; zero to give anomalously high coefficients when slightly underpopulated. ; Note that the maxwidth attempts to counter this effect. ; * Could ignore all zero values. Not sure what effect this would have. ; * I could put it in the loop and use only nearby (current goodobj) objects. I don't ; like this poor man's adaptative smoothing because then the order in which the ; coefs are chosen alters the PDF. ; * Better would be adaptative smoothing using a choice of whetever "h" is and storing ; a 2D array of smoothing widths: one for every particle in every direction. ; ; MODIFICATION HISTORY: ; Mar 04 - Tidied up by RM. ; May 02 - Written by Richard Massey ; function simage_pdf_kernel_width,pdf,topology,$ PAPER=paper,ADAPTATIVE=adaptative,SILENT=silent COMPILE_OPT idl2,HIDDEN if keyword_set(paper) then begin if (size(pdf))[2] ne 137 then message,"Can only use smoothing kernel from MNRAS 348, 214 with n_max=15" if max(topology) ne 2 then message,"Can only use smoothing kernel from MNRAS 348, 214 with smoothing in polar shapelet space" kernel_width=[0.0704679,0.224751,0.00290480,0.0114870,0.00802608,0.0185946,0.0211331,0.0475693,0.0133025,$ 0.0177635,0.0265865,0.0399162,0.0327700,0.0268335,0.0983199,0.0378721,0.0284180,0.0369072,$ 0.0431431,0.0403410,0.0319501,0.0427166,0.0389435,0.108788,0.0379268,0.0295396,0.0203575,$ 0.0480356,0.0339141,0.0229076,0.0213328,0.00682178,0.0427869,0.0302036,0.0726666,0.0231790,$ 0.0212116,0.0108982,0.0190023,0.0316582,0.00495220,0.00585339,0.0114293,0.00450972,0.00403608,$ 0.0215383,0.0126624,0.0295790,0.00772182,0.00831089,0.00136459,0.00740194,0.00918036,0.0112132,$ 0.00603026,0.00364463,0.00595399,0.000155495,0.00152507,0.000988727,0.00493227,0.00240407,0.00574509,$ 0.00149079,0.00194335,0.000124179,0.000941016,0.00181656,0.00199292,0.00188102,0.00158170,0.000772038,$ 0.00128351,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,$ 0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,$ 0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,$ 0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,$ 0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,$ 0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,$ 0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327,0.251327] message,"Using smoothing kernel verified in Massey et al. (2004) MNRAS 348, 214",/info;,/noprint=silent endif else begin ; Set defaults minwidth = 1e-9 ; Error catching (don`t want any divisions by zero) maxwidth = 0.25 ; Don`t want huge spread in PDF either. ARBITRARY. ; Find size of the probability space that we're dealing with message,"Calculating smoothing lengths in the PDF",/info,noprint=silent arrsize = size(pdf) n_obj = arrsize[1] n_coeffs = arrsize[2]-1 if not keyword_set(topology) then topology=intarr(n_coeffs)+1 ; Create empty variable to contain answer kernel_width = fltarr(n_coeffs) ; Manually set smoothing size in object size and magnitude directions kernel_width[0] = mean(simage_pair_separation(pdf[where(finite(pdf[*,0])),0]))*50.;15. ; Size (beta) smoothing lenth kernel_width[1] = mean(simage_pair_separation(pdf[where(finite(pdf[*,1])),1]))*150. ; Magnitude smoothing ; Remaining shapelet coefficient directions for i=2,n_coeffs-1 do begin message,strtrim(i,1)+'/'+strtrim(n_coeffs,1),/info,noprint=silent centre = pdf[*,i] if topology[i] eq 1 then begin ; Magnitude coef smoothing length kernel_width[i] = mean(simage_pair_separation(centre[where(finite(centre) eq 1)]))/15.;*0.3 endif else kernel_width[i] = !pi/12.5 ; Phase coef smoothing length endfor endelse ; Return kernel size that has been decided upon return,kernel_width end ; ; ************************************************************************ ; ************************************************************************ ; ;+ ; NAME: ; SIMAGE_PAIR_SEPARATION ; ; CATEGORY: ; A component of simage_resample_pdf. ; ; PURPOSE: ; Finds the nearest neighbours of a list of points - so that create.pro can use ; the median or max value as a smoothing length. Works for 1D numbers. ; ; INPUTS: ; points - A catalogue of objects read in by, eg: shapelets_read_shapeshapecat.pro or ; shapelets_read_shapecat_hdf.pro ; ; OUTPUTS: ; separations - Distance of every point to its nearest neighbour. ; ; MODIFICATION HISTORY: ; Jan 02 - Written by Richard Massey function simage_pair_separation,points compile_opt idl2, HIDDEN n=n_elements(points) separations=points-points ; to get varible type correct (long, double etc) separations[0] = min( abs( [points[1:n-1]] - points[0] ) ) for i=1,n-2 do separations[i] = min( abs( [points[0:i-1],points[i+1:n-1]] - points[i] ) ) separations[n-1] = min( abs( [points[0:n-2]] - points[n-1] ) ) return,separations end ; ; ************************************************************************ ; ************************************************************************ ; ;+ ; NAME: ; SIMAGE_RESAMPLE_PDF ; ; CATEGORY: ; Shapelets image simulation. ; ; PURPOSE: ; Generates a brand new object from a shapelet catalogue by drawing a ; bootstrapped sample from an (n_a+1)-dimensional probability space. This is ; smoothed using the Kernel method, where the widths of the kernels is decided ; by simage_pdf_kernel_width.pro and may even be adaptive to the local density. The new ; object is essentially the same as one of the originals, but perturbed by ; multivariate elliptical Gaussians in shapelet space. ; ; INPUTS: ; shapecat- A catalogue of objects read in by, eg: shapelets_read_shapeshapecat.pro or ; shapelets_read_shapecat_hdf.pro ; ; OPTIONAL INPUTS: ; width - size of kernel smoothing in each direction. This is slow to calculate, ; so is returned to speed up a second run. ; n_max - n_max to create new object up to (if catalogue goes that far) ; It is recommended that you go as far as the decompositions, because ; the size & magnitude of the objects are maintained. If, for example, ; an n_max=12 object is modelled by this code as an n_max=4 object, the ; output will be a large, bright version of the object`s nucleus. ; Default: whatever the catalogue goes up to. ; seed - Seed for random number generator. ; bs_obj - The catalogue number of an object to morph (bootstrap). ; Default: one picked at random. ; ; KEYWORD PARAMETERS: ; /PLOT - Plot newly generated object to screen at end. ; /SILENT - Operate silently. ; ; OUTPUTS: ; new_object - Cartesian decomp structure, representing a new, morphed galaxy. ; bs_obj - The catalogue number of the morphed object. ; ; NOTE: ; Covariance matrix and even errors of coefficients are ignored, but ; their existence does hint that some smoothing in shapelet space is ; philosophically justified. ; ; MODIFICATION HISTORY: ; Mar 04 - Tidied up by RM. ; May 02 - Complete rewrite by RM: renamed & now uses fast bootstrap method ; from Silverman (1984) rather than integrating the entire PDF. ; Apr 02 - RM - faithful size-magnitude plane ensured by fitting it by hand. ; Feb 02 - RM - now copes with polar shapelet decompositions. ; Dec 01 - Create.pro written by Richard Massey ;- pro simage_resample_pdf, shapecat, new_object, width=width, seed=seed, $ n_max_input=n_max_input, plot=plot, silent=silent, bs_obj=bs_obj COMPILE_OPT idl2 ;on_error,2 ; ; ; SET UP PARAMETERS ; ; message,'Initialising...',/info,noprint=silent if not keyword_set(seed) then seed=1009L ; Random seed if not keyword_set(n_max_input) then n_max_input=shapecat.maxn_max ; Nmax for the created object. if n_max_input eq 999 then n_max_input=0 ; Enable us to sample just size/mag plane n_max_cat = n_max_input < shapecat.maxn_max ; Make sure the coefs we want to sample exist! shapelets_make_nvec,n_max_cat,n1_cat,n2_cat,n_coeffs_cat ; Indices of each coefficient. n_objects = shapecat.n ; Number of sample points in PDF. ; ; ; LOCATE KNOWN POINTS IN PROBABILITY SPACE ; (corresponding to galaxies in the input catalogue) ; ; pdf = shapecat.coeffs ;for i=0,shapecat.n-1 do pdf[i,*] = pdf[i,*] / pdf[i,0] ; THIS HAS ALREADY BEEN DONE! ; Insert size, R, as the zeroth dimension pdf = transpose([transpose(sqrt(abs(shapecat.rsquared))),transpose(pdf)]) n_coeffs_cat = n_coeffs_cat+1 ; Change first dimension to be the magnitude (rather than just f_00) pdf[*,1] = shapecat.mag ; This should ensure that size-mag diagram will be automatically correct ; (albeit with some smoothing in that plane too) ; ; ; FLAG THE TOPOLOGY OF DIFFERENT DIMENSIONS OF THE PROBABILITY SPACE ; ; Needs optimising: ; Topology need only be calculated within simage_pdf_kernel_width.pro ; Reorder can be in any order - just needed to make "good" variable below. ; ; simage_flag_pdf_topology,shapecat.maxn_max,topology,reorder,polar=shapecat.polar ; ; ; FIND OPTIMUM KERNEL WITH WHICH TO SMOOTH THE PDF ; ; ; If possible, use the kernel size verified by Massey et al. (2004) MNRAS 348, ; 214 to produce realistic new galaxy morphologies. This is also a lot faster! paper=shapecat.polar ; Slow, so allow process to be done once only by returning kernel size as a vector if not keyword_set(width) then $ width=simage_pdf_kernel_width(pdf,topology,/SILENT,PAPER=paper) ; ; ; BOOTSTRAP FROM PDF ; ; message,'Generating a brand new shapelet object...',/info,noprint=silent coeffs = fltarr(n_coeffs_cat) ; Final answer will go here! ; Follows fast method from Silverman (1986), page 143. ; Step (i) - Choose I uniformly with replacement from {1,...,n} if not keyword_set(bs_obj) then begin baddecomp=1 while baddecomp do begin bs_obj=fix(randomu(seed)*n_objects) ;& bs_obj=213 & bs_obj=2113 & bs_obj=212 if shapecat.flag[bs_obj,0] le 3 then begin if ( shapecat.chisq[bs_obj] lt 2.5 ) or ( shapecat.flag[bs_obj,1] ge 3 and shapecat.chisq[bs_obj] lt 5. ) then baddecomp=0 endif endwhile endif ; Step (ii) - Generate epsilon to have probability density function K n_max_new=shapecat.n_max[bs_obj] < n_max_cat shapelets_make_nvec,n_max_new,n1_new,n2_new,n_coeffs_new good=reorder[0:n_coeffs_new] ; Step (iii) - Set Y = X_I + h*epsilon epsilon=randome(seed,n_coeffs_new+1) coeffs[good]=reform(pdf[bs_obj,good])+( width[good] * epsilon ) ; Strictly speaking, the phases are no longer in [0,2pi] but the trig ; functions needed to rotate them will automatically wrap the phases ; around to their proper range. ; ; ; EXTRACT BETA & a_00 INFORMATION ; ; message,'Renormalising size and magnitude of new object...',/info,noprint=silent rsq_aim=coeffs[0]^2 ; This is the size... mag_aim=coeffs[1] ; ...and magnitude we want the object to have. flux_aim=10.^(-.4*mag_aim) ; Express magnitude back in units of flux. coeffs=coeffs[1:*] & n_coeffs_cat=n_coeffs_cat-1 ; Remove beta from coefficient lists. coeffs[0]=1. ; Remove magnitude from coefficient lists. Set default f_00. Will renormalise later. ; ; ; ROTATE RANDOMLY ; ; ; Transform only the relevant coefficients to polars, for speed good=(good[sort(good)])[1:*]-1 ; WHAT THE HELL IS THIS?!! if shapecat.polar then coeffs=shapelets_polar_expand(coeffs[good]) else coeffs=shapelets_c2p(coeffs[0:n_coeffs_new-1]) ; Also works wih polars: ;if shapecat.polar then coeffs=polar_expand(coeffs) else coeffs=shapelets_c2p(coeffs[0:n_coeffs_new-1]) ;good = indgen(n_coeffs_cat-1)+1 phase = atan(imaginary(coeffs),double(coeffs)) modulus = abs(coeffs) ; Parity flip? yesno = 2*(randomu(seed) gt 0.5)-1 phase = phase*yesno ; Random rotation angle = randomu(seed)*2.*!dpi m = n1_new-n2_new phase = phase+m*angle coeffs = modulus*complex(cos(phase),sin(phase)) ; Transform back to Cartesian shapelets (should remember matrix for speed!!!) coeffs = shapelets_p2c(coeffs) ; Pad with zeros ;if n_coeffs_cat gt n_coeffs_new then coeffs = [coeffs,fltarr(n_coeffs_cat-n_coeffs_new)] ; ; ; STORE NEWLY GENERATED OBJECT (in a decomp structure) ; new_object={name:"Brand new object", $ type:"decomp", $ history:"Morphed from object #"+strtrim(string(bs_obj),1)+" in "+shapecat.name+" catalogue.",$ x:[15,15], $ beta:1., $ n_max:n_max_new, $ n_coeffs:n_coeffs_new, $ coeffs:coeffs[0:n_coeffs_new-1], $ error:fltarr(n_coeffs_new), $ n1:n1_new, $ n2:n2_new, $ n_pixels:[30,30], $ over:1B, $ integrate:1B, $ chisq:shapecat.chisq[bs_obj], $ skyfit:0., $ sex:0B, $ moments:{flux:flux_aim,mag:mag_aim,rsquared:sqrt(abs(rsq_aim))}} ; ; ; STORE NEWLY GENERATED OBJECT (in a one-element shapecat structure: obfuscated) ; ; ;new_object={name:"Brand new object", $ ; type:"shapecat", $ ; n:1, $ ; maxn_max:n_max_new, $ ; polar:0B, $ ; seeing:shapecat.seeing, $ ; x:[0,0], $ ; beta:beta, $ ; n_max:n_max_new, $ ; n_coeffs:n_coeffs_cat, $ ; coeffs:coeffs, $ ; error:fltarr(n_coeffs_cat), $ ; flux:flux_aim, $ ; flux_error:0, $ ; snr:0, $ ; mag:mag_aim, $ ; centroid:[0,0], $ ; centroid_error:[0,0], $ ; rsquared:sqrt(abs(rsq_aim)), $ ; rsquared_error:0, $ ; ellipticity:[0,0], $ ; ellipticity_error:[0,0], $ ; quadrupole:[[0,0],[0,0]], $ ; quadrupole_error:[[0,0],[0,0]], $ ; flag:[0,0], $ ; chisq:shapecat.chisq[bs_obj], $ ; sexid:shapecat.sexid[bs_obj], $ ; sexfwhm:shapecat.sexfwhm[bs_obj], $ ; sexflux:shapecat.sexflux[bs_obj], $ ; sexmag:shapecat.sexmag[bs_obj], $ ; sexclass:shapecat.sexclass[bs_obj], $ ; sexa:shapecat.sexa[bs_obj], $ ; sexb:shapecat.sexb[bs_obj], $ ; sextheta:shapecat.sextheta[bs_obj], $ ; sexe1:shapecat.sexe1[bs_obj], $ ; sexe2:shapecat.sexe2[bs_obj], $ ; sexx:shapecat.sexx[bs_obj], $ ; sexy:shapecat.sexy[bs_obj], $ ; sexarea:shapecat.sexarea[bs_obj], $ ; objx_max:shapecat.objx_max[bs_obj], $ ; objx_min:shapecat.objx_min[bs_obj], $ ; objy_min:shapecat.objy_min[bs_obj], $ ; objy_max:shapecat.objy_max[bs_obj]} ; ; ; SET FINAL COEFFICIENTS: RENORMALISE ; ; ; Check actual properties of generated object ; (coefs are still in form a_ij/a_00 but this won't matter since temporarily a_00=1) rsq_true = shapelets_rsquared(new_object,flux=flux_true) new_object.beta = float(sqrt(abs(rsq_aim/rsq_true))) ; so can renormalise to desired size... new_object.coeffs = new_object.coeffs*flux_aim/flux_true/new_object.beta ; ...and magnitude! ; ; ; DISPLAY NEW OBJECT ; ; if keyword_set(plot) then begin old_object=shapelets_shapecat2decomp(shapecat,bs_obj) new_object.x=old_object.x new_object.n_pixels=old_object.n_pixels print,' beta rsquared mag flux a(0)' print,'Original was',shapecat.beta[bs_obj],shapecat.rsquared[bs_obj],shapecat.mag[bs_obj],shapecat.flux[bs_obj],shapecat.coeffs[bs_obj,0] print,'Aiming for ', ' -------',rsq_aim,mag_aim,double(flux_aim),' ---------'; Used to be different because of disretely sampled PDF, but now same as above print,'Intermediate', 1.,float(rsq_true),' -------',flux_true,1. print,'End up with ',new_object.beta,float(shapelets_rsquared(new_object,flux=new_flux)),$ -2.5*alog10(float(new_flux)),new_flux,new_object.coeffs[0] window,0 shapelets_recomp,new_object,new_recomp plt_image,new_recomp,/frame,/col, title="!6Newly generated object", $ xtitle="!6Size: "+strtrim(new_object.moments.rsquared*.1,1)+'" '+$ "Magnitude: "+strtrim(mag_aim,1)+' '+$ "n!dmax!n: "+strtrim(n_max_new,1)+' '+$ "!4e!6: "+strtrim(sqrt(total(epsilon^2)),1) ;plt_decomp,new_object,/coef window,2 shapelets_recomp,old_object,old_recomp plt_image,old_recomp,/frame,/col, title="!6Original object" ;plt_decomp,old_object,/coef read,junk endif end

View the help page for this routine, return to the shapelets web page or return to the code help menu.


Last modified on 02nd May 2008 by Richard Massey.

Valid HTML 4.01!