;+
;
; Project   : s3drs and Pixon reconstruction software
;                   
; Name      : p_run_test
;               
; Purpose   : part of the s3drs GUI
;               
; Explanation: (part of the testing and legacy code, not documented herein)
;               
; Use       : 
;    
; Inputs    : 
;               
; Outputs   : 

; Keywords  :
;               
; Common    : 
;               
; Restrictions: none
;               
; Side effects: 
;               
; Category    : 
;               
; Prev. Hist. : (see CVS logs)
;
; Written     : Sandy Antunes, NRC, 2007-2009
;               
;-            
; originally for fluxrope, now being modified to allow hshell2 (/shell)
;
; And, if /cor2a is given, saves data and sigma and sets normd akin to Cor2A
;
; If given '/halfmask', adds a mask on frame 2 for 1/2 the view, to
; test masking
;
; If given /detmask, adds a detector-derived mask

PRO p_run_test,Ncube=Ncube,$
               local=local,pregen=pregen,plot=plot,exp=exp,hours=hours,$
               debug=debug,dnorm=dnorm,iloops=iloops,puresaito=puresaito,$
               float=float,double=double,ps=ps,offset=offset,$
               stopearly=stopearly,forcesigma=forcesigma,forcedata=forcedata,$
               e0=e0,shell=shell,gamma=gamma,zhack=zhack,tag=tag,$
               theta=theta,z0=z0,LL=LL,cor2a=cor2a,dataonly=dataonly,$
               cube=cube,tetra=tetra,planar=planar,roco=roco,saito=saito,$
               halfmask=halfmask,detmask=detmask,fits=fits,longrun=longrun,$
               offthomson=offthomson,imagesave=imagesave,noiselevel=noiselevel

; /pregen=use pregenerated data if available

;; extremely simple test case for Pixon, create and render a fluxrope

;   .reset
   on_error,2 
;   plot=0;  1=plot 0=no plots

   memtalk,'startup'

   plot=KEYWORD_SET(plot)
   ps=KEYWORD_SET(ps)
   saito=KEYWORD_SET(saito)
   stopearly=KEYWORD_SET(stopearly)
   imagesave=KEYWORD_SET(imagesave)
   dataonly=KEYWORD_SET(dataonly)
   ; how many sigma of extra noise to add
   if (n_elements(noiselevel) eq 0) then noiselevel=0
   local=KEYWORD_SET(local)
   pregen=KEYWORD_SET(pregen)
   halfmask=KEYWORD_SET(halfmask)
   detmask=KEYWORD_SET(detmask)
   offset=KEYWORD_SET(offset)
   offthomson=KEYWORD_SET(offthomson)
   longrun=KEYWORD_SET(longrun)
   puresaito=KEYWORD_SET(puresaito)
   exp=KEYWORD_SET(exp)
   fits=KEYWORD_SET(fits)
   debug=KEYWORD_SET(debug)
   double=KEYWORD_SET(double)
   if (double eq 0) then float=1; uses float as default
   zhack=KEYWORD_SET(zhack)
   planar=KEYWORD_SET(planar); default is no, use classic Pixon instead
   tetra=KEYWORD_SET(tetra); default is no, use classic Pixon instead
   if (n_elements(dnorm) eq 0) then dnorm=1.0
   if (n_elements(e0) eq 0) then e0=1.0
   if (n_elements(iloops) eq 0) then iloops=1
   if (n_elements(forcesigma) eq 0) then forcesigma=-1.0
   if (n_elements(forcedata) eq 0) then forcedata=-1.0

   shell=KEYWORD_SET(shell)
   cube=KEYWORD_SET(cube)
   if (shell and cube) then cube=0; shell takes precedence over cube

   if (n_elements(gamma) eq 0) then gamma=[0.0,90.0]
   nviews=n_elements(gamma)

   if (n_elements(theta) eq 0) then theta=make_array(nviews,/double,value=0.0)
   if (n_elements(z0) eq 0) then z0=make_array(nviews,/double,value=0.0)
   if (n_elements(LL) eq 0) then LL=make_array(nviews,/double,value=1.0)

   if (n_elements(Ncube) eq 0) then Ncube=32;

   prep_3drv,pxwrapper,views=nviews,N=Ncube, $
     gamma=gamma,theta=theta,z0=z0,LL=LL,dnorm=dnorm

   cor2a=KEYWORD_SET(cor2a)

   if (cor2a) then p_modelcor,pxwrapper,instr='Cor2'

   ; turns off other toggles to prevent them adding non-Cor2A values
;   if (cor2a) then dnorm=1.0
;   if (cor2a) then forcesigma=-1.0
;   if (cor2a) then for i=0,pxwrapper.sr.Nd-1 do $
;     pxwrapper.sr.detector[i]='COR2'

;   help,pxwrapper.sr.L0tau,pxwrapper.sr.tau
;   print,pxwrapper.sr.L0tau,pxwrapper.sr.tau
;   if (cor2a) then pxwrapper.sr.L0tau=15.0/float(0.5*pxwrapper.sr.Ncube)
;   if (cor2a) then pxwrapper.sr.tau=pxwrapper.sr.L0tau/(2.0*pxwrapper.sr.L0)
;   help,pxwrapper.sr.L0tau,pxwrapper.sr.tau
;   print,pxwrapper.sr.L0tau,pxwrapper.sr.tau

   if (exp) then print,'Using experimental memory model' else $
     print,'Using original memory model'

;   N=256
;   N=32
;   N=512

   modelname='fluxrope'
   pregenfile='fluxrope.dat'

   if (shell) then modelname='hshell2'
   if (shell) then stuff=[0,0]
   if (offset) then stuff=[45.0,45.0]
   if (shell) then pregenfile='shell.dat'

   if (puresaito) then modelname='saito'

   if (cube) then modelname='ncube'
   if (cube) then stuff=[Ncube/3]
   if (cube) then pregenfile='cube.dat'

   if (pregen eq 1 and file_exist(pregenfile) eq 0 ) then pregen=0
   if (pregen eq 1) then restore,/verb,pregenfile
   if (n_elements(datum) eq 0) then pregen=0; restore failed
   if (pregen eq 1) then if (n_elements(datum[0,0,*]) ne nviews) then pregen=0

   ; for offset, shrink and embed
   if (offset) then Nm=Ncube/2 else Nm=Ncube
   if (pregen eq 0) then $
     timage=pixonmodels(Nm,modelname,ignore,stuff,stuff1)
;   help,datum

   if (offset) then timage2=fltarr(Ncube,Ncube,Ncube);
   if (offset) then timage2[Ncube/2:(Ncube/2)+Nm-1,0:Nm-1,0:Nm-1]=timage
   if (offset) then timage=timage2

   if (saito) then saitoimg=pixonmodels(Ncube,'saito',pxwrapper)
   if (shell and saito) then saitoimg=saitoimg/9E9
   if (modelname eq 'fluxrope' and saito) then saitoimg=saitoimg/30.0

   if (saito) then timage=timage + saitoimg

   if (imagesave) then save,file='timage.sav',timage

;   threeview,timage

;;   if (offset) then $
;;     timage2=timage[Nm/4:(3*Nm/4)-1,Nm/4:(3*Nm/4)-1,Nm/4:(3*Nm/4)-1]
;;   if (offset) then timage3=timage*0.0
;;   if (offset) then timage3[Nm/2:Nm-1,0:(Nm/2)-1,0:(Nm/2)-1]=timage2
;;   if (offset) then timage=fltarr(Ncube,Ncube,Ncube)

   ; cludgy way to move image to an upper quarter for off-center rendering
;   timage2=timage
;   timage2[Ncube/2+1:Ncube-1,*,*]=$
;     timage[1+Ncube/4:Ncube/2+Ncube/4-1,*,*]
;   timage2[0:Ncube/2,*,*]=min(timage)
;
;   timage3=timage2
;   timage3[*,Ncube/2+1:Ncube-1,*]=$
;     timage2[*,1+Ncube/4:Ncube/2+Ncube/4-1,*]
;   timage3[*,0:Ncube/2,*]=min(timage2)
;   if (offset) then timage=timage3

;   e0=max(timage)
;  prep_3drv,pxwrapper,views=nviews,N=Ncube,enorm=e0

   if (offthomson) then for i=0,pxwrapper.sr.Nd-1 do pxwrapper.sr.polar[i]='n'

   ; ability to test Paul's planar z-viewpoint, only
   if (n_elements(roco) ne 0) then pxwrapper.sr.roco=roco

  pxwrapper.sc.note='p_run_test'
  pxwrapper.sc.gname=modelname

  pxwrapper.sr.n0=e0

;  pxwrapper.sr.interp=3

  ; ZHACKS STILL DO NOT WORK!
  if (zhack) then pxwrapper.sr.theta[2]=90.0
  if (zhack) then pxwrapper.sr.phi[2]=90.0
  if (zhack) then pxwrapper.sr.z0[2]=pxwrapper.sr.lla0[2]
;  if (zhack) then pxwrapper.sr.z0[2]=pxwrapper.sr.L0

;  pxwrapper.sr.interp=5; 5 = use C version
;  pxwrapper.sr.interp=0; 0 = force native IDL for big runs

  pxwrapper.sp.noshm=1; 1 = do not use shm
  pxwrapper.sp.noshm=0; 0 = use shm
;  pxwrapper.sp.noshm=-1; -1 = sometimes use shm

  ; sets and zeroes common variables for experimental mem use
;;  pxn_memset,1,/debug
;  pxn_memset,exp,/debug
  if (exp) then pxn_memset,exp,debug=debug

  noct=2; half as many, check memory usage for large N
  noct=4 ; default, will be set to 4 in Pixon
  if (Ncube gt 128) then noct=2;

  npo=1; half as many, check memory usage for large N
  npo=2 ; default, will be set to 2 in Pixon
  if (Ncube gt 128) then npo=1;
  pxwrapper.sp=p_tweak_kernel(pxwrapper.sp,noct,npo,/debug)

  print,'npxn: ',pxwrapper.sp.npxn

   ; optional rescaling test
;   pxwrapper.sr.n0=max(timage)   fails.  n0 = 1 fails.  Why?

; offset tests, showed that p_center_data is not needed for real data.
;;  offset=6.0; test of p_center_data's effects
;;  for m=0,n_elements(pxwrapper.sr.k0)-1 do $
;;    pxwrapper.sr.k0[m] = ptr_new( p_getk0(pxwrapper.sr,m,shift=offset) )
;;  for m=0,n_elements(pxwrapper.sr.k0)-1 do $
;;    pxwrapper.sr.k0[m] = ptr_new( p_getk0(pxwrapper.sr,m,shift=offset,sun) )

  memtalk,'pre-p_generate'

  ; just in case I want to toggle to 'tetra' 
  if (tetra) then classic='tetra' else if (planar) then classic='planar' else $
    classic = 'classic'
  if (pregen eq 1) then $
    p_generate_pixon,classic,pxwrapper,oPixon,local=local else $
    p_generate_pixon,classic,pxwrapper,oPixon,image=timage,local=local

  memtalk,'post-p_generate'

;   img=oPixon->Get(/im,/nop)
;   threeview,image,res=6
;   p_generate_pixon,'tetra',pxwrapper,oPixon,image=timage

   if (pregen eq 1) then datum=rebin(datum,Ncube,Ncube,nviews) else $
     datum=stereo_project(pxwrapper,oPixon)

   print,'datum max(datum), min(datum): ',max(datum),min(datum)


   save,file=pregenfile,datum,pxwrapper
   print,'saved ',pregenfile
   if (dataonly) then return

   sdnorm=strtrim(string(floor(pxwrapper.sr.normd)),2)
   if (pxwrapper.sr.normd lt 1.) then $
     sdnorm=$
     'E-'+strtrim(string((strsplit(string(pxwrapper.sr.normd),'.'))[1]),2)

   if (n_elements(tag) eq 0) then tag=''
   tag=modelname+'_dnorm'+sdnorm+'_'+tag

   pxwrapper.sc.gname=tag

   if (float) then datum=float(datum)

   print,'orig max(datum), min(datum): ',max(datum),min(datum)
   if (forcedata lt 0) then $
     datum=datum*pxwrapper.sr.normd else $
     datum=(forcedata/max(datum))*datum

   if (forcesigma lt 0) then $
     sigma=0.05*pxwrapper.sr.normd else $
     sigma=forcesigma

   sigscalartest=0
   if (size(sigma,/n_dimensions) ne 3) then sigscalartest=1
   ; sigma is likely a scalar
   if (sigscalartest and n_elements(sigma) eq 1) then $
     sigma=replicate(sigma,pxwrapper.sr.Nd)
   if (sigscalartest) then sigscalar=sigma
   if (sigscalartest) then $
     sigma=fltarr(pxwrapper.sr.Ncube,pxwrapper.sr.Ncube,pxwrapper.sr.Nd)
   if (sigscalartest) then $
     for i=0,pxwrapper.sr.Nd-1 do sigma[*,*,i]=sigscalar[i]

   ; special Cor2A options
;;   dbest=4.0E4; our ideal Pixon max data value
;;   maxcor2a=4E-7
;;   dn2photon=15
;;   biasmean=688.0
;;   exptime=6.0
;;   MSBperDN=2.7E-12; Bsun/DN
   ;
;;   datamsb=datum * maxcor2a / max(datum)
;;   d5 = datamsb * exptime/MSBperDN
;;   sigmad5 = sqrt(d5*dn2photon)
;;   sigmad5plusbias = sqrt(sigmad5^2 + biasmean)
;;   izero = where (d5 le 0)
;;   inotzero = where (d5 gt 0)
;;   mind5=sqrt(min(d5[inotzero])); assigns an arbitrary min
;;   d5[izero]=mind5
;;   fractional_err = sigmad5plusbias/d5
;;   fractional_noise = datamsb * fractional_err
   ;
;;   msbnorm = dbest/max(datamsb)
;   pxwrapper.sr.msbperdn=1.0
;   pxwrapper.sr.exptime=1.0

    if (cor2a) then sigma=p_quicknoise(datum,pxwrapper=pxwrapper)
    print,'sigma is max ',max(sigma),' vs datum max ',max(datum)
;   if (cor2a) then pxwrapper.sr.normd=normd

;;   if (cor2a) then datum=datamsb
;;   if (cor2a) then sigma = fractional_noise * msbnorm
;;   if (cor2a) then pxwrapper.sr.normd[*]=msbnorm
;;   if (cor2a) then pxwrapper.sr.exptime[*]=exptime   
;;   if (cor2a) then pxwrapper.sr.dn2photon[*]=dn2photon
;;   if (cor2a) then pxwrapper.sr.MSBperDN[*]=MSBperDN
   ; Note that we are keeping units for cor2A as MSB and then
   ; settind sr.normd to transform to a best pixon range, which is
   ; then automatically handled in p_insert_data

   print,'max(datum), max(sigma) and sr.normd: ',$
     max(datum),max(sigma),pxwrapper.sr.normd
;   help,datum,sigma

   ; note if noiselevel is 0, we get unperturbed data  
   datum = datum + sigma*noiselevel*randomn(seed,Ncube,Ncube,pxwrapper.sr.Nd)
   print,'postnoise max dat,sig: ',max(datum),max(sigma)
   ; recalc sigma if we added noise
   if (noiselevel ne 0) then sigma=p_quicknoise(datum,pxwrapper=pxwrapper);
   print,'postquicknoise max dat,sig: ',max(datum),max(sigma)

   if (plot) then tv_multi,datum,res=6,/axes,/border

;   hack=0
;   if (hack) then pxwrapper.sr.roco=[0,0]; converts S to P

   masks=sigma*0.0 + 1.0
   if (halfmask) then masks[Ncube/2:*,*,1]=0.0; put a half-mask on this puppy
   if (detmask) then for i=0,pxwrapper.sr.Nd-1 do masks[*,*,i]=$
     p_getmask(sr=pxwrapper.sr,iele=i)

   print,'p_run_test recap: insert sigma, datum'
;   help,sigma,datum
   print,'sigma: ',max(sigma),', datum: ',max(datum)

   if (cor2a) then pxwrapper.sr.normd=1.0

   if (halfmask or detmask) then $
     p_insert_data,pxwrapper,datum,oPixon,sigma=sigma,masks=masks else $
     p_insert_data,pxwrapper,datum,oPixon,sigma=sigma
   
   datacheck=oPixon->get(/data,/nop)
   sigcheck=oPixon->get(/sigma,/nop)
   wtcheck=oPixon->get(/wt,/nop)
   print,'check: oPixon sigma, data, wt'
   print,max(sigcheck),max(datacheck),max(wtcheck)

;   oPixon->Full,/restart ; works

   p_report,pxwrapper,oPixon,/save,runname=tag

   ; extra saves here
   if (fits) then pixon2fits,pxwrapper,oPixon,/data; extra save
   if (fits) then pixon2fits,pxwrapper,oPixon,/sigma; extra save

   memtalk,'pre-p_pixonloop'

   if (stopearly) then save,pxwrapper,oPixon,datum,sigma,file='stopearly.sav'

   if (stopearly) then return

   if (longrun) then oPixon->set,mxit=4000
     
   if (longrun) then opixon->full,/restart else $
     if (n_elements(hours) ne 0) then $
     p_pixonloop,pxwrapper,oPixon,/init,/full,hours=hours,$
     tag=tag,saveeach=5,/endsave,qcrit=0.5 else $
     p_pixonloop,pxwrapper,oPixon,iloops,20,/init,/full,saveeach=5,tag=tag,$
     /endsave,qcrit=0.5

   memtalk,'post-p_pixonloop'

   print,'post-p_pixonloop memory: ',memory(/current),memory(/highwater)

   p_report,pxwrapper,oPixon,/save,/echo,runname=tag

   if (ps) then p_standardplots,pxwrapper,oPixon,/no

; global constants to set for physically real problems:
;sr.n0, sr.chl, sr.L0, sr.lla0, sr.d, sr.interp, sr.psp
;sp.reltol, sp.abstol, sp.mxit, sp.npo

print,'End of p_run_test'

END