;+
;
; Project   : s3drs and Pixon reconstruction software
;                   
; Name      : p_pixonloop
;               
; 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
;               
;-            
; /best = hours=72, qcriteria=1.0, /saveeach=5, /endsave
; but any of the above can be overwritten by command line options.
;
; runs the step-by-step full Pixon solver (aka 'pseudo' plus Pixon,
; dubbed 'map' here) by default, or can use the built-in Full
; method ('/full') or the quick method ('/quick') or just the
; pseudo solver ('/pseudo'), or the older 'classic' method by Paul Reiser
; ('/classic'), or even the older 'anneal' method ('/anneal')
;
; Note that, currently, the 'quick' method in Pixon does not work
; for our sample problems.  We recommend avoiding it.
;
; Note also that /init runs the first loop after initing
;
; Can use optional parameter 'hours', where it will try to run
; loops for X hours (overridding 'nloops' but retaining 'istep' in
; case the user is plotting information).
;
; Or, can use the 'qcriteria=' limit, eg qcrit=1, qcrit=10, et cetera
; to set a limiting 'q' so that we end when Q goes below that limit.
;
; If you give both hours= and qcrit=, it will end when the first
; condition is met (i.e. either qcrit dips below the desired
; threshold or the clock runs out).
;
; Note both criteria hours= and qcrit= are in addition to the
; usual fit tolerance.  Runs end when the desired fit-to-noise
; ftol is reached, if that happens before the hours or qcrit is met.
;
; The parameter 'runtime' returns the time it took to run this program,
; and is usually used by 'hours' but can also be used by the user.
; Note 'hours' does not include the time it takes to init the system,
; nor run the handful of loops necessary for timing.  Therefore,
; you need to estimate that overhead when deciding whether to do
; a run.  Note that for long hours=X runs, this overhead is minor.
; On a 1.xGHz PC, a 64^3 run took to do an init/estimate of runtime >5min.
;
; /flush clears out any existing timage and puts in a neutral one
;
; /init is the same as using both /flush and /restart
;
; With 'saveeach=n', will save the results (pxwrapper, oPixon, plus text
; report) into an incrementing save file every 'saveeach' loops
; (remember that each loop lasts 'istep' mxit iterations.)
; With '/endsave', it saves the results at the end of the fit
; (recommended in case a run ends unattended.)
;

PRO p_pixonloop,pxwrapper,oPixon,nloops,istep,update=update,plot=plot,$
                init=init, flush=flush, restart=restart, hours=hours, $
                runtime=runtime, qcriteria=qcriteria, tag=tag, map=map, $
                pseudo=pseudo,quick=quick,classic=classic,full=full, $
                anneal=anneal,saveeach=saveeach,endsave=endsave,best=best,$
                resume=resume

   pixon_env,pxwrapper.object; force check to ensure pixon is loaded

   best=KEYWORD_SET(best)
   if (best) then begin
       if (n_elements(hours) eq 0) then hours=72.0
       if (n_elements(qcriteria) eq 0) then qcriteria=1.0
       if (n_elements(saveeach) eq 0) then saveeach=5
       if (n_elements(endsave) eq 0) then endsave=1
   end

   resume=KEYWORD_SET(resume)

   if (n_elements(nloops) eq 0) then nloops=1
   if (n_elements(istep) eq 0) then istep=pxwrapper.sp.mxit
   update=KEYWORD_SET(update)
   plot=KEYWORD_SET(plot)
   init=KEYWORD_SET(init); flush and restart
   restart=KEYWORD_SET(restart); Pixon's 'restart'
   flush=KEYWORD_SET(flush)
   pseudo=KEYWORD_SET(pseudo)
   quick=KEYWORD_SET(quick)
   classic=KEYWORD_SET(classic)
   anneal=KEYWORD_SET(anneal)
   full=KEYWORD_SET(full); the default, included as a key to allow explicit calls
   endsave=KEYWORD_SET(endsave)
   map=KEYWORD_SET(map)
   if (n_elements(tag) eq 0) then tag='pixonloop'
   if (n_elements(qcriteria) eq 0) then qcriteria=-1

   if (anneal eq 0 and pseudo eq 0 and classic eq 0 and $
       quick eq 0 and full eq 0) then begin
       map=1                   ; default method is map
   end else begin
       map=0; if another method is chosen
   end

   if (n_elements(saveeach) eq 0) then saveeach=0
   if (n_elements(hours) eq 0) then hours=0

   if (pxwrapper.object eq 't') then begin
       if (full or quick or map) then $
         print,"Warning, tetrahedral pixon must use 'pseudo'"
       quick=0
       full=0
       map=0
       pseudo=1
   end

   if (qcriteria gt 0) then print,'** Running until Q is less than ',qcriteria
   if (hours gt 0) then print,'** Running for ',hours,' hours'

   ; replace any earlier result with a blank image
   N=pxwrapper.sr.Ncube
   if (flush or init) then begin
       memtalk,'p_pixonloop pre-init'
       ; original was: oPixon->Replace,image=dblarr(N,N,N)*0+1/double(N)
       im=oPixon->Get(/im); as a pointer, we hope
       if (pxwrapper.sc.doubleprec) then const=1.0/double(N) else const=1.0/float(N)
       for il=0,N-1 do for jl=0,N-1 do for kl=0,N-1 do (*(im[0]))[il,jl,kl]=const
       oPixon->Replace,image=im
       memtalk,'p_pixonloop post-init'
   end

   ; put in a blank pixonmap
;   N=pxwrapper.sr.Ncube
;   pxnidx = FLTARR(N,N,N)+pxwrapper.sp.npxn-1
;   oPixon->Set, pxnidx=pxnidx
;   oPixon->Set, reltol=rsched[0]
;  pxnmap=oPixon->Get(/pxnidx, /nop)

   ; free all before start and create blank map and blank image
   if (init or restart) then begin
       oPixon->Restart
;       oPixon->Replace,mxit=1; just a quick hop to prep things and settle
;       if (pxwrapper.object eq 't') then begin
;           ; tetra lacks a full method
;           oPixon->pseudo,/restart
;       end else begin
;           oPixon->Full,/restart
;       end
   end

   if (hours ne 0) then begin

       if (pxwrapper.sc.looptime eq 0) then begin

           print,'starting time estimate'
       ; make a first time estimate thru quick invocation of this routine
       ; (sans init, hours, update, plot options)
           p_pixonloop,pxwrapper,oPixon,1,3, runtime=rtest, $
             pseudo=pseudo,quick=quick,classic=classic,full=full,map=map
           pxwrapper.sc.looptime = rtest/3.0; estimates seconds/loop
           if (pxwrapper.sc.looptime lt 1) then pxwrapper.sc.looptime=1
       end

       total_loops = hours * 3600.0 /pxwrapper.sc.looptime
       total_nloops = long(total_loops/istep)
       if (total_nloops lt 1) then begin
           ; if overzealous isteps request, readjustment needed...
           ; force into a 1-loop setup
           istep = long(total_loops/2)
           istep=min(1,istep) ; always ensure there's at least a few iterations
           total_nloops=1
       end
       print,'Each iteration takes ',pxwrapper.sc.looptime,' seconds.'
;       print,'Estimating ',total_nloops,' loops of ',istep,'
;       iterations will fit within ',$
       ts=strtrim(string(total_nloops),2)
       is=strtrim(string(istep),2)
       hs=strtrim(string(hours),2)
       print,'Estimating ',ts,' loops of ',is,' iterations will fit within ',$
         hs,' hours.'
       nloops=total_nloops  ; reset this and we are ready to run below
   end

   oPixon->Replace,mxit=istep

   starttime=systime(/seconds)

   cwin=!D.WINDOW               ; temporarily save current window

;;   tv_multi,stereo_project(pxwrapper,oPixon),title='FIT IN PROGRESS'
   if (update ne 1 and plot eq 1) then begin

       tv_multi,(stereo_project(pxwrapper,oPixon) < 5000 > 3000),/log,/border,title='FIT IN PROGRESS',res=196/pxwrapper.sr.Ndata[0],winid=winid

       ; need to create window for plot, too

;       plotx = [pxwrapper.sc.iterations]

;       ploty1 = [pxwrapper.sc.pcount]
;       ploty2 = [pxwrapper.sc.chisq]
;       ploty3 = [pxwrapper.sc.q]

       window,/free,xsize=600,ysize=400
;       !P.MULTI = [0, 0, 3]
;       plot,plotx,ploty1,title='Pixon Count',/ynozero,yticks=2
;       plot,plotx,ploty2,title='Chisq',/ynozero,yticks=2
;       plot,plotx,ploty3,title='Q (Chisq/DOF)',/ynozero,yticks=2
;       !P.MULTI= 0              ; reset back to normal
       wplot = !D.WINDOW

       clean_wset,cwin                ; now restore previous window

   end

   savecounter=0
   actualiterations=0
   for iloop=1,nloops do begin
       print,iloop,nloops,istep,format='("... running loop ",i," of ",i,",",i," iterations")'
       if (quick) then begin
           oPixon->Quick
       end else if (pseudo) then begin
           oPixon->Pseudo
       end else if (classic) then begin
           ; earlier Paul Reiser method, inside of his annealing
           oPixon->Map
           ; can smooth image /image, pseudomap /phi, or pixon map /pxnidx
           oPixon->PixonSmooth, /image
           oPixon->Pseudo
       end else if (anneal) then begin
           ; discarded 'anneal' method
           if (tag_exist(pxwrapper.sr,'N')) then N=pxwrapper.sr.N else $
               N=pxwrapper.sr.Ncube ; backward compat
           sp=pxwrapper.sp
           p_anneal,oPixon,sp.jsched,sp.msched,sp.rsched,sp.npxn,N,pxwrapper
       end else if (full) then begin
           oPixon->Full
       end else begin
           ; map method, recommended
           oPixon-> Pixon::NNLS
           oPixon-> Map
           oPixon-> pseudo
       end
       actualiterations=actualiterations+istep

       p_update_runstats,pxwrapper,oPixon,istep

       ++savecounter
       if (saveeach gt 0 and savecounter eq saveeach) then begin
           p_report,pxwrapper,oPixon,/savedata,runname=tag
           savecounter=0
       end

;       print,pxwrapper.sc.iterations,pxwrapper.sc.pcount,pxwrapper.sc.chisq,pxwrapper.sc.q


       if (n_elements(plot_iterations) eq 0) then begin
                                ; initialize
           plot_iterations=pxwrapper.sc.iterations
           plot_pcount = pxwrapper.sc.pcount
           plot_chisq = pxwrapper.sc.chisq
           plot_q = pxwrapper.sc.q
       end else begin
                                ; accumulate
           plot_iterations = [plot_iterations, pxwrapper.sc.iterations]
           plot_pcount = [plot_pcount, pxwrapper.sc.pcount]
           plot_chisq = [plot_chisq, pxwrapper.sc.chisq]
           plot_q = [plot_q, pxwrapper.sc.q]
       end

       ; help,plot_iterations,plot_q
       if (plot eq 1) then begin

           if (n_elements(winid) ne 0) then $
;             tv_multi,(stereo_project(pxwrapper,oPixon) < 5000 >
;             3000),/log,/border,/update,res=196/pxwrapper.sr.Ndata[0],winid=winid
             tv_multi,stereo_project(pxwrapper,oPixon),/log,/border,/update,res=196/pxwrapper.sr.Ndata[0],winid=winid

           if (n_elements(wplot) ne 0) then begin

               clean_wset,wplot

               !P.MULTI = [0, 0, 3]
               ; note we skip the first data point as it is too often 0
               plot,plot_iterations,plot_pcount, $
                 title='Pixon Count',/ynozero,yticks=2
               plot,plot_iterations,plot_chisq, $
                 title='Chisq',/ynozero,yticks=2
               plot,plot_iterations,plot_q, $
                 title='Q (Chisq/DOF)',/ynozero,yticks=2
               !P.MULTI = 0; reset back to normal

           end

           clean_wset,cwin; now restore previous window
       end

       if (hours ne 0) then begin
           ; check on our predicted runtime to see how we are doing so we
           ; never run more than mxit iterations over our requested runtime.
           currenthours = (systime(/seconds)-starttime)/3600.0; runtime in hrs
           if (currenthours gt hours) then begin
               print,'Whoops, taking longer than expected, ending early after ',actualiterations,' iterations, in ',currenthours,' hours (instead of ',hours,').'

               break
           end
       end

       if (pxwrapper.sc.q lt qcriteria) then begin
           ; check on desired Q level
           print,'Q dropped below ',qcriteria,', ending as requested.'
           break
       end

   end

   runtime= systime(/seconds) - starttime; runtime in seconds

   ; newer, better estimate of sec/loop (as long as we do at least a few loops)
   ; note use of 'i-1' in case original loop breaked early for overtime
   if (actualiterations gt 8) then begin
       pxwrapper.sc.looptime = runtime/actualiterations
       print,'New update on seconds per loop: ',pxwrapper.sc.looptime
   end

   pxwrapper.sc.totaltime = pxwrapper.sc.totaltime + runtime

;;   tv_multi,stereo_project(pxwrapper,oPixon),title='fit to real data'

   minutes=strtrim(string(runtime/60),2)
   print,'Pixon fit ended, run took ',minutes,' minutes'

   ; now add in these run stats to any previous runs,
   ; order is always iterations, pcount, chisq, q
   if (ptr_valid(pxwrapper.sc.runstats[0])) then begin
       ; add to prior runs by storing them, then freeing the pointers
;       print,'storing then freeing earlier runs'
       plot_iterations= [*(pxwrapper.sc.runstats[0]),plot_iterations]
       plot_pcount= [*(pxwrapper.sc.runstats[1]),plot_pcount]
       plot_chisq= [*(pxwrapper.sc.runstats[2]),plot_chisq]
       plot_q= [*(pxwrapper.sc.runstats[3]),plot_q]
       for ip=0,3 do ptr_free,pxwrapper.sc.runstats[ip]
   end
   ; now store
   ; help,plot_iterations,plot_q
   pxwrapper.sc.runstats[0]=ptr_new(plot_iterations)
   pxwrapper.sc.runstats[1]=ptr_new(plot_pcount)
   pxwrapper.sc.runstats[2]=ptr_new(plot_chisq)
   pxwrapper.sc.runstats[3]=ptr_new(plot_q)
   if (endsave) then p_report,pxwrapper,oPixon,/savedata,runname=tag
   ;help,*pxwrapper.sc.runstats[0],*pxwrapper.sc.runstats[3]

END