PRO pxn_test, PSOUT=psout

; Include pixon code in front of path:
; (do it one time only)
;
if STRPOS(!path,'pixon_IDL_source') EQ -1 then $
  !PATH = '/nfs/wiwaxia/d4/ASC/packages/pixon/pixon_IDL_source:'+!PATH

; Make simulated 1-D data...
;-------------------------------------
; This could be replace by a section to read in
; image_data and psf arrays from real data...
; (e.g., 
;    image_data from a section of a Chandra PHA file 
; and psf from a row of an rmf file...)
;-------------------------------------
;
; Tell us about the data:
n_dim = 1
n_array = 64

; Set up the true image features
; Hardcoded for 64 long array...
;
; Start with low-level background
image_true = fltarr(64) + 0.1
; A wide, low plateau
image_true[11:34]=1.5+image_true[11:34]
; Bright isolated source on the plateau
image_true[15]=9.0+image_true[15]
; A narrower, tall plateau
image_true[24]=3.0+image_true[24]
image_true[25:29]=5.0+image_true[25:29]
image_true[30]=3.0+image_true[30]
; A He-like tripplet
image_true[38]=7.0+image_true[38]
image_true[43]=2.0+image_true[43]
image_true[48]=5.0+image_true[48]

; Scale the image
image_true = 50.0 * image_true

; Setup a normalized psf (even number of points)
psf = [0.02,0.04,0.2,0.34,0.25,0.08,0.05,0.02]
psf = psf/TOTAL(psf)

; Apply the PSF to the image
image_nonoise = convolve(image_true, psf)

; Possion statistics added!
; or not...
;;if KEYWORD_SET(NONOISE) then begin
;;  image_data = image_nonoise
;;end else begin
  image_data = g_poiss(image_nonoise)
;;end

; Show the true image and the data...
if KEYWORD_SET(PSOUT) then begin
  pre_print_landscape
end else begin
  window,0
end

!p.multi = [0,1,2]
plot, image_true, PSYM=-4,xrange=[-1.,64.],/xstyle, $
	TITLE='True Image',linestyle=1
plot, image_data, PSYM=-2,xrange=[-1.,64.],/xstyle, $
	TITLE = 'Data: noisy(*), noise-free(diamond)'
oplot, image_nonoise, PSYM=-4, linestyle=1

if KEYWORD_SET(PSOUT) then begin
  device, /close
  SPAWN, 'cp idl.ps pxn_truedata.ps'
  set_plot, 'X'
end

;-------------------------------------
; At this point we have available:
;   n_dim, n_array, psf, image_data [, image_true]
;-------------------------------------


; wait for user...
print, ''
print, ' Here''s the input True and Noisy images...'
print, ''
print, ' Next, we''ll run pxn with a pixon map of all 1.0''s. '
print, '   This just checks that an unconstrained solution'
print, '   is possible, e.g., that PSF and noise estimates are sensible.'
print, ' --- hit any key to continue... ---'
user=get_KBRD(1)


; - - - First PIXON Step - - - 
; This step "performs a pure chi-square minimization with no
; [pixon] constraints on the structural scales in the image",
; good residuals should be obtained in this case.
;
; Basically there is a radius = 1.0 pixon at each pixel so
; the forward folding problem is underconstrained.
;

; keep track of how many iterations...
nit=0

; Create an error array
image_error = (SQRT(image_data) > 1.0)

; Put PSF in same size array as the image:
; centered too
n_psf = n_elements(psf)
psf_full = [ fltarr((n_array-n_psf)/2) , psf, $
	fltarr((n_array-n_psf)/2)]

; Flag valid parts of the data
; Set to all 1.0's
image_valid = 1.0+fltarr(n_elements(image_data))
image_valid(indgen(n_psf))=0.0
image_valid( indgen(n_psf) + (n_array - n_psf) )=0.0
n_valid = TOTAL(image_valid)

; initial pixon map...
; Each pixel has pixon of size 1.0
pxnmap = 1.0 + 0.0*image_data

; starting pseudo image:
ps_image = 0.0 * image_data


; first pxn running...
;
; PRO PXN, Data, PSF, PXN, Image, Residuals, Sig=Sig, Weights=Weights, $
;         DMask=DMask, IMask=IMask, Qtol=Qtol, npxns=npxns, PImage=PImage, $
;                 Rec=Rec, wsc=wsc, pow=pow
;
pxn, image_data, psf_full, pxnmap, $   ; <-- inputs
	npxns=1, $                     ; <-- input: number of pixon sizes
	Sig=image_error, $             ; <-- input: sigma of noise in data
        Dmask=image_valid, $           ; <-- input: mask valid data
        Imask=image_valid, $           ; <-- input: mask valid image
	PImage = ps_image, $           ; <-- in/out: pseudo image
	pxn_image, pxn_resid, $        ; <-- outputs: decon image, norm resids
	wsc=2.0                        ; <-- plot scale 


; Make some plots... 
; Show the true image and the data...
if KEYWORD_SET(PSOUT) then begin
  pre_print_portrait
end else begin
  window,0
end

!p.multi = [0,1,4]
csize = 1.8

  plot, image_true, PSYM=-4,xrange=[-1.,n_array],/xstyle, linestyle=1, $
	TITLE='Iteration = '+STRCOMPRESS(nit)+ $
		' : True (diamonds), Deconvolved(hist)', $
		CHARSIZE=csize
  oplot, pxn_image, PSYM=10

; Show the data
plot, image_data, PSYM=-2,xrange=[-1.,n_array],/xstyle, $
	TITLE='Data(--*--)', CHARSIZE=csize

; Pixon map values
plot, pxnmap, PSYM=4,xrange=[-1.0,n_array],/xstyle, $
	CHARSIZE=csize, YRANGE=[0.0,1.1*MAX(pxnmap)],/YSTYLE, $
	TITLE='Pixon Map'

; Pseudo image
;  oplot, (ps_image > 0.0), PSYM=2

; Residual values
; Residual values
chi2all = TOTAL((pxn_resid)^2*image_valid)/n_valid
plot, pxn_resid, PSYM=4,xrange=[-1.,n_array],/xstyle, $
	CHARSIZE=csize, YRANGE=[-5.0,5.0],/YSTYLE, $
	TITLE='Pixon Residuals, Chi per point, Chi^2/N ='+ $
		string(chi2all,format='(F7.2)')

if KEYWORD_SET(PSOUT) then begin
  device, /close
  SPAWN, 'cp idl.ps pxn_itr'+strcompress(nit,/remove)+'.ps'
  set_plot, 'X'
end

; - - - - - - - - - - -
; Will iterate pxncal and pxn below...
; Set pxnmap to have a large value instead of all 1's,
; this large value will set the initial maximum pixon size.
; Base the value on the size of the data array:
pxnmap(where(image_valid NE 0.0)) = n_array/4

; starting pseudo image:
ps_image = 0.0 * image_data

; starting pxn_image:
pxn_image = image_data

; Setting the noise larger than it is will cause larger
; sized pixon values - a good starting point.  On each iteration
; reduce the sigma_factor closer to 1.0 to approach
; the full-data significance...
; Do a square-root of current value at each iteration
; Starting value squared:
sigma_factor = 4.0^2

; wait for user...
print, ''
print, ' OK, now we''ll do real pixon iterations starting'
print, ' with very wide pixons and iterating to smaller sizes'
print, ' as the data allow.'
print, ' --- hit any key to continue... ---'
user=get_KBRD(1)

while sigma_factor GT 1.01 do begin
; = = = = = = iterate from here to bottom = = = = = =

; - - - Run pxncal - - - 
;
; Function pxncal, Data, Image, Psf, SIG, gdmsk=gdmsk
;
;   Now, instead of assigning all pixels to have a
; pixon value of 1.0, run this routine to estimate
; the pixon radii for each pixel in the image...
;

nit = nit + 1
sigma_factor = SQRT(sigma_factor)

print, ' - - - - -'
print, ' Pixon iteration number:', nit
print, '   sigma_factor = ',sigma_factor
print, ' - - - - -'
pxnmap = pxncal( image_data, pxn_image, psf_full, $
		sigma_factor*image_error, $
		gdmsk=image_valid, maxdels=MAX(pxnmap*image_valid) )

; Set border values to 0.0
pxnmap = pxnmap * image_valid

; Report the uniq values used
print, pxnmap(UNIQ(pxnmap,SORT(pxnmap)))

; wait for user...
print, ''
print, ' Just calculated a new pixon scale map for this iteration'
print, ' and will next run pxn to find best fit with this pixon map.'
print, ' --- hit any key to continue... ---'
user=get_KBRD(1)


; - - - Additional pxn iterations with non-trivial pixon map - - -
pxn, image_data, psf_full, pxnmap, $   ; <-- inputs
	npxns=12, $                    ; <-- input: number of pixon sizes
	Sig=sigma_factor*image_error, $ ; <-- input: sigma of noise in data
        Dmask=image_valid, $           ; <-- input: mask valid data
        Imask=image_valid, $           ; <-- input: mask valid image
	PImage = ps_image, $           ; <-- in/out: pseudo image
	pxn_image, pxn_resid, $        ; <-- outputs: decon image, norm resids
	wsc=2.0                        ; <-- plot scale 


; Make some plots... 
; Show the true image and the data...
if KEYWORD_SET(PSOUT) then begin
  pre_print_portrait
end else begin
  window,0
end

!p.multi = [0,1,4]
csize = 1.8

  plot, image_true, PSYM=-4,xrange=[-1.,n_array],/xstyle, linestyle=1, $
	TITLE='Iteration = '+STRCOMPRESS(nit)+ $
		' : True (diamonds), Deconvolved(hist)', $
		CHARSIZE=csize
  oplot, pxn_image, PSYM=10
; Pseudo image
;  oplot, (ps_image > 0.0), PSYM=2

; Show the data
plot, image_data, PSYM=-2,xrange=[-1.,n_array],/xstyle, $
	TITLE='Data(--*--)', CHARSIZE=csize

; Pixon map values
plot, pxnmap, PSYM=4,xrange=[-1.0,n_array],/xstyle, $
	CHARSIZE=csize, YRANGE=[0.0,1.1*MAX(pxnmap)],/YSTYLE, $
	TITLE='Pixon Map (sigma_factor = ' + $
		string(sigma_factor,FORMAT='(F7.2)')+')'

; Residual values
chi2all = TOTAL((pxn_resid)^2*image_valid)/n_valid
plot, pxn_resid, PSYM=4,xrange=[-1.,n_array],/xstyle, $
	CHARSIZE=csize, YRANGE=[-5.0,5.0],/YSTYLE, $
	TITLE='Pixon Residuals, Chi per point, Chi^2/N ='+ $
		string(chi2all,format='(F7.2)')

if KEYWORD_SET(PSOUT) then begin
  device, /close
  SPAWN, 'cp idl.ps pxn_itr'+strcompress(nit,/remove)+'.ps'
  set_plot, 'X'
end

print, ''
print, ' Iteration complete, see plot window...'
print, ' --- hit any key to continue... ---'
user=get_KBRD(1)

; = = = = = = end of iterations = = = = = =
end


RETURN
END