PRO me_test, NONOISE=nonoise, ITALL=itall
;+
;  PRO me_test, NONOISE=nonoise, ITALL=itall
;
;PURPOSE:
; Demo the astro-lib Max_Entropy IDL routine...
; Currently assumes 1-D data but could do 2-D as well...
;
;USAGE:
; Compiling:
;  IDL> .run me_test  [will see error for g_poiss, do it again:]
;  IDL> .run me_test
; Running:
;  IDL> me_test
;
;OPTIONAL KEYWORDS:
;  /NONOISE = Don't add Poisson noise to the fake data
;  /ITALL   = Carry out all the iterations up to max_iters
;             (other wise there are stopping criteria in the code.)
; 
;EXTERNAL ROUTINES:
;  convolve  [in astro-lib ?]
;  g_poiss [for pseudo-data creation, appended below]
;  g_rand  [for g_poiss, appended below]
;
;REVISIONS:
; 3/19/2001 dd Orignal version
;-

; Parameters - - - 

max_iters = 50

; - - - - - - - - -

; 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...
image_true = fltarr(64) + 0.1
image_true[15]=9.0
image_true[25:29]=5.0
image_true[24]=3.0
image_true[30]=3.0
image_true[38]=7.0
image_true[43]=2.0
image_true[48]=5.0

; Scale the image
image_true = 20.0 * image_true

; Setup a psf
psf = [0.02,0.04,0.2,0.34,0.25,0.08,0.05,0.02]

; 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...
;;!p.multi = [0,1,2]
;;plot, image_true, PSYM=4,xrange=[-1.,64.],/xstyle
;;plot, image_data, PSYM=-2,xrange=[-1.,64.],/xstyle
;;oplot, image_nonoise, PSYM=-4, linestyle=1

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


; Before the first ME iteration signal we are starting over:
multipliers = 0  
nit = 0
; and set some things up:
lastchi2n = 1.E12
keep_going = 1


; Do the iterations...
while keep_going EQ 1 do begin
; ---- repeat these lines to iterate... -----


  Max_Entropy, image_data, psf, image_deconv, multipliers, FT_PSF=psf_ft
  nit = nit+1

  ; Do not allow 0.0's in image_deconv
  ; (image is nominally in "counts" so this limit is generous:)
  image_deconv = image_deconv > 1.0E-6

  ; For statistics, use only those points with multipliers GT 0
  ; (looks like multipliers EQ 0 at the edges due to psf length...)
  sel = where(multipliers NE 0.0)

  ; Check and print:
  ;  - Chi2 at each iteration... assuming counting statistics
  ;  (could be added to MAX_ENTROPY since Re_conv is the same as image_folded)
  ;  - "entropy" measure of model (image_deconv) at each iteration
  
  ; Chi^2
  image_folded = convolve(image_deconv, psf)
  chipp = (image_data - image_folded)/SQRT(image_data > 1.0)
  thischi2n = TOTAL(chipp(sel)^2/N_ELEMENTS(chipp(sel)))

  ; Entropy: different versions...
  ;
  ;  "- f ln(f)" with f in 0-1 range:
  frac = image_deconv(sel)/TOTAL(image_deconv(sel))
  ent1 = TOTAL( -1.0 * frac * alog(frac) )
  ;  "ln(I)" 
  ent2 = TOTAL( alog(image_deconv(sel)) )
  ;  "sqrt(I)" 
  ent3 = TOTAL( sqrt(image_deconv(sel)) )

  ; Print them
  print, 'Iter = '+STRING(nit, FORMAT='(I3)')+ $
	',    Chi2/N = '+STRING( thischi2n , FORMAT='(F7.3)') + $
	',    -f ln(f) = '+STRING( ent1 , FORMAT='(F7.3)') + $
	',    ln(I) = '+STRING( ent2 , FORMAT='(F7.1)') + $
	',    sqrt(I) = '+STRING( ent3 , FORMAT='(F7.1)')

  if NOT( KEYWORD_SET(ITALL) ) then begin
    ; Stop iterations based in Chi2
    if thischi2n LT 1.0 then keep_going = 0
    if thischi2n GT lastchi2n then keep_going = 0
    if (thischi2n/lastchi2n) GT 0.99 then keep_going = 0
    lastchi2n = thischi2n
  end

  ; otherwise stop when max_iters reached
  if nit GE max_iters then keep_going = 0

; --------------------------------------------
; could end the iteration loop here to avoid a
; plot per iteration...
;;end


; Make some plots... 

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

if n_elements(image_true) GT 0 then begin
  plot, image_true, PSYM=4,xrange=[-1.,n_array],/xstyle, $
	TITLE='Iteration = '+STRCOMPRESS(nit)+ $
		' : Deconv(hist) and True(diamonds)', CHARSIZE=csize
  oplot, image_deconv, PSYM=10
end else begin
  plot, image_deconv, PSYM=10,xrange=[-1.,n_array],/xstyle, $
	TITLE='Iteration = '+STRCOMPRESS(nit)+ $
		' : Deconvolution', CHARSIZE=csize
end

; Show the data and the deconvolved model folded through the PSF
plot, image_data, PSYM=-2,xrange=[-1.,n_array],/xstyle, $
	TITLE='Folded-Deconv(hist) and Data(--*--)', CHARSIZE=csize
oplot, image_folded, PSYM=10

; Residuals
plot, chipp, PSYM=-2,xrange=[-1.,n_array],/xstyle, $
	TITLE='Chi per point:  (Data - Folded-Deconv) / SQRT(Data)' + $
		',  Chi2/N = '+STRING( thischi2n, FORMAT='(F7.3)' ), $
	CHARSIZE=csize, YRANGE=[-5.0,5.0],/YSTYLE


; --------------------------------------------
; or end the iteration loop here to plot each iteration
end


RETURN
END
;=====================================================================

;=====================================================================
FUNCTION g_poiss, mu
; Subroutine to return samples from poisson distribution
; If mu < 11. the correct poisson values are calculated and a random
; value chosen, otherwise a gaussian approximation is used.
; 5/27/92 dd

COMMON rand, seed

; Comment this out to use systime for seed
; if n_elements(seed) eq 0 then seed = 17
if n_elements(mu) eq 0 then mu = 1.0

; pre-calculate the factorials, 0-20
fact = [1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880.,3628800.,$
	39916800.,479001600.,6227020800.,8.71783E10,1.30767E12,$
	2.09228E13,3.55687E14,6.40237E15,1.21645E17,2.43290E18]
probx = fact

x = float(mu)

FOR i=0,n_elements(mu)-1 DO BEGIN

  IF(mu(i) GE 11.0) THEN BEGIN
; Gaussian approximation
    x(i) = FLOAT(mu(i) + sqrt(mu(i))*g_rand())
    x(i) = FLOAT(LONG(x(i)))
  ENDIF ELSE BEGIN
; Poisson calculation
    ; Fill array of cumulative probabilities for 0 to 20 events
    probx(0) = 1.0
    FOR ip=1,20 DO BEGIN
      probx(ip) = probx(ip-1) + (mu(i)^ip)/fact(ip)
    END
    p_norm = EXP(-1.*mu(i))
    probx = p_norm*probx
    r_val = randomu(seed)
    out_val = FLOAT(indgen(23))
    x(i) = INTERPOL(out_val,[0.,probx,1.0],[r_val])
    x(i) = FLOAT(LONG(x(i)))
  ENDELSE

END

RETURN, x
END
;=====================================================================

;=====================================================================
FUNCTION g_rand, dum
; Subroutine to return a normally distributed deviate with zero mean
; and unit variance (a la Abramowitz/Stegun p953, Acceptance-rejection(4).)
;
; 6/7/91 dd
;

COMMON rand, seed

; Comment this out to use systime for seed
; if n_elements(seed) eq 0 then seed = 17

REPEAT BEGIN
  u1 = randomu(seed)
  u2 = randomu(seed)
  x = -1.0*ALOG(u1)
ENDREP UNTIL ( (x-1.0)*(x-1.0) LT -2.0*ALOG(u2) )

IF randomu(seed) LT 0.5 THEN x = -1.0*x

RETURN, x
END
;=====================================================================