;pro ufit,setfname,debug,var,a
;fitting program to call model-generating programs that read
;setup-files.
;This version is build to call utm (universal transit modeler)
;debflag goes from 0 to 3
;var,a are output-params
;written by Hans J Deeg
;
;Version 9 Nov 2000   
;
;version history
; 3/5/99 first working version
; 5/5/99 includes wfit handling
; 3/8/00 dbrat and chitol are now optionals in parameter file
; 30/9/00 clarified screen-output of chi and stdev
; 9/10/00 inserted output params a,var; debflag 0 doesn't plot anymore
; 18/10/00 changed repmax from 10 to 3 and itmax from 20 to 10
; 8/11/00 repmax and itmax are now optional paramfile values without 
; warning on absence. Temporary fitsets have now a unique
; filename. This allows to run more then 1 ufit process at the same
; time in the same directory.
; 9/11/00 fixed bug introduced yesterday


@utm   ;compile utm  (compiles also setupfile and readlc procedures)

function frac, fnum
;returns the fractional part of a number
;example:
; print frac(-3.1)
;     0.9000000
;print frac(3.1)
;     0.1000000
;
;7 mar 1997
;hans joerg deeg
return, (fnum mod 1) +1.0 - ((fnum mod 1) ge 0)
end

function testf,x,a
;this must be a function callable by an array x, and has to return an array
;of the same dimension with computed y values.
;HERE CALL TO UTM. 
;THIS IS THE ONLY ROUTINE THAT NEEDS ADAPTION FOR OTHER FUNCTIONS
common parrs,pvarr,pnarr
common fitparnamecom,fitparname,fitsetmp

yval=dblarr(n_elements(x))
;update 'a' into parameterfile
adim=n_elements(a)
for i=0,adim-1 do begin
    paraput,fitparname[i],a[i]
endfor
savesetupfile,fitsetmp
;invoke utm
utm,fitsetmp,0,x,yval
return,yval
end 


function mqmfitf,x,a
common sharxy,alast,xarr,yarr,ydat,xind,pdrat,debflag

;preparaton routine for lmfit
;Calculates array function 'testf' and its partical derivatives once.
;Recalculates only if coefficients 'a' have changed.
;Return value is a vector of function-value at x and partial derivatives
;x must be a single value! (=scalar) and must be element of xarr

;print,'start mqm xind=',xind
if xind eq n_elements(xarr) then begin
    xind = 0                    ;new round of invocations
;print,'xind reset to 0'
endif
if total(a ne alast) gt 0.0 then begin ;new recalculation
    yarr[*,0]=testf(xarr,a)     ;values of function

    k= pdrat                    ;kappa in partial derivatives
    for i=0,n_elements(a)-1 do begin ;calc partial derivatives
        ap=a & ap[i] = a[i]*(1+k) ;param-vectors with 1 param changed
        yarr[*,i+1]=(testf(xarr,ap)-yarr[*,0])/(1*a[i]*k)
;    below, slower version to calc part derivs on both sides 
;       am=a & am[i] = a[i]*(1-k)
;       yarr[*,i+1]=(testf(xarr,ap)-testf(xarr,am))/(2*a[i]*k)
    endfor
    if debflag ge 1 then begin
        oplot,xarr,yarr[*,0],linestyle=2
        ydiff=ydat-yarr[*,0]
        oplot,xarr,ydiff,linestyle=1
    endif
    if debflag ge 2 then begin
        chi=total(ydiff*ydiff)
        sdev=sqrt(chi/(n_elements(ydat)-1))
        print,'        mqmfitf:params: ',a
        print,'        mqmfitf:chisq : ',chi
        print,'        mqmfitf:stdev : ',sdev
        print,'        -----'
    endif
    xind=0
endif 
if x ne xarr[xind] then print,'ERROR:mqmfitf out of sync: x,xarr:',x,xarr[xind]
;print,'x,xind,xarr[xind],yarr[xind,*]',x,xind,xarr[xind],yarr[xind,*]
xind=xind+1                     ;bookkeeping for next invocation
alast = a
;if xind eq 129 then stop
;print,'end mqm'

return,yarr[xind-1,*]           ;return vector
end




pro ufit,setfname,debug,var,a
common parrs,pvarr,pnarr
common sharxy,alast,xarr,yarr,ydat,xind,pdrat,debflag
common fitparnamecom,fitparname,fitsetmp

if n_elements(debug) ne 0 then debflag=debug else debflag =0 ;debug-flag
readsetupfile,setfname

infilename=parasign('infilename') ;get input lc
lcarr=readlc(infilename,2)
xdat=lcarr[*,0]                 ;x,y arrays with data to be fitted
ydat=lcarr[*,1]

wflag=fix(parasign('wflag',0))  ;determines if last fit is being saved
if wflag then lcfilename= parasign('lcfilename','lc.out')

;----
;make up vector of initial fit parameters from setup file
;first, get NAMES of fit parameters from setup
fitparname=strarr(8)            ;max of 8 fit parmaeters
adim=0
for i=0,n_elements(pnarr)-1 do begin
    if pnarr[i] eq 'fit' then begin
        if adim eq 8 then print,"Ufit: Too many fit parameters specified!"
        fitparname[adim] = pvarr[i]  
        adim = adim+1
    endif
endfor
if adim eq 0 then print,"Ufit: At least 1 fit parameter needs to be specified!"
fitparname = fitparname[0:adim-1]
;now get VALUES of fit params from setup
a=dblarr(adim)
for i=0,adim-1 do begin
    a[i]=double(parasign(fitparname[i]))
endfor
aorg=a                          ;keep copy of original values
print,'ufit:  input parameters--------------'
for i=0,adim-1 do print,fitparname[i]," :",a[i]
if debflag ge 2 then print,'ufit:  -----------fitting------------'
;-------


;initialize remaining common-block vars
;adim=n_elements(atrue)
pdrat=double(parasign('pdrat',0.001)) ;ratio of xi in partial derivatives
xdim=n_elements(xdat)
xarr=xdat
yarr=dblarr(xdim,adim+1)        ;yarr are values of current invoke of test-fct
alast = dblarr(adim)            ;keep for next invoke
xind=0
;-------

;-------
;override some param values in the set-up arrays
paraput,'tflag','1'             ;use only time-points (x-vals in utm)
paraput,'onoise','0'
paraput,'wflag','0'             ;don't save temporary results

;create unique name for fitsetmp that is based on systime to 1/100s sec
stime=frac(systime(1)/10000)*10000
seed=stime  ;this, since randomu will modify seed
ranu=randomu(seed)*10000 ;random number

fitsetmp=strcompress(string(format='(f15.2,f7.1)',stime,ranu),/remove_all)+'ufit.tmp'
;print,fitsetmp
savesetupfile,fitsetmp


;prepare plots
if debflag ge 1 then begin
    oflag=fix(parasign('oflag',0))
    !x.title='time ('+parasign('tunit')+')'
    case oflag of
    0: begin& !y.title='sys.lum' & !y.range=[min(ydat),max(ydat)] &end
    1: begin&!y.title='rel.loss'& !y.range= [max(ydat),min(ydat)] &end
    2: begin&!y.title='rel.mag' & !y.range=[max(ydat),min(ydat)] &end
    endcase
    plot, xdat, ydat,psym=+1,title=infilename ;Plot the initial data
endif


;------------------------the fit--------------------------------
chi=0D & chilast=9999999D         
conv=0                          ;needs to be named var
iter=0 &iterct=iter
repct=0

repmax=fix(parasign('repmax',3,/nowarn))   ;max number of calls to LMFIT
itmax=fix(parasign('itmax',10,/nowarn))    ;max number of iterations by LMFIT
tol=float(parasign('chitol',1e-5)) ;chisq tolerance


repeat begin  
                                ;this repeat loop is a fix until
                                ;it is better checked out WHY lmfit often stops too early
;print,'start LM'
    yfit = LMFIT(Xdat, Ydat, a, FUNCTION_NAME = 'mqmfitf',chisq=chi,converge=conv,itmax=itmax,alpha=alpha,covar=covar,iter=iter)
;print,'end LM'
    iter=iter+1  ;iter appreas to be always 1 too low
    iterct=iterct+iter
    repct=repct+1

;print,'iter=',iter
;print,'iterct=',iterct
;print,'repct',repct
                                ;update 'a' into parameterfile
    for i=0,adim-1 do begin
        paraput,fitparname[i],a[i]
    endfor
    savesetupfile,fitsetmp
                                ;define break conditions
    chidiff=(chilast-chi)/chilast
    chilast=chi
    if debflag ge 2 then begin
        print,'LMFIT result:'
        for i=0,adim-1 do print,fitparname[i]," :",a[i]
        print,'chisq :',chi
        print,'stdev :',sqrt(chi/(n_elements(xarr)-1))
        print,'----------------------------------------'
    endif    
    if debflag ge 3 then stop
if repct ge repmax then print,'!!warning: No convergence of LMFIT in repmax=',str(repmax),' repeats'
endrep until (repct ge repmax) or (chidiff le tol) 
;iterct=iterct+repct
;---------------------------------------------------------------

;output stuff
;print,'stddev= ',stddev(yfit-ydat)
;print,'converge=',conv
ydiff=ydat-yfit
if debflag ge 1 then begin
    plot, xdat,ydat,psym=1,title=infilename+' and fit'
    oplot, xdat, yfit,linestyle=2 ;Overplot the fitted data.
    oplot, xdat, ydiff, linestyle=1
endif
savesetupfile,'fitset'
if debflag ge 2 then print,'ufit: saved parameters to fitset ----'
print,'------------ufit: fitted parameters--------------'
for i=0,adim-1 do print,fitparname[i]," :",a[i]
print,'number iterations: ',iterct
;    print,'chisq  : ',chi
var=total(ydiff*ydiff)/(n_elements(ydat)-1)
sdev=sqrt(var)
print,'stdev  : ',sdev
if debflag ge 2 then begin
    print,'----------------detailed output--------------'
    print,'convergence index: ',conv
;    print,'variance: ',var
    print,'alpha  :'
    print,alpha
    print,' '&print,'covar :'
    print,covar
endif
if debflag ge 3 then stop
if wflag then begin             ;save fitted values 
    paraput,'wflag','1'         ;put it back into fitsetmp
    savesetupfile,'fitset'
    print,'running utm once more with final parameters to save fitted'
    print,'values to file: ',lcfilename
    utm,fitsetmp,0,xdat
endif
;remove temp fitsetfile
exestr='\rm '+fitsetmp
spawn,exestr
end