;+
; NAME:
;                  AUX_SFH
;
; RESTRICTIONS:   This routine is not approved part of the ULySS
;                 package. It is not tested and was designed for a
;                 particular purpose.
; 
;                 (Comments from PhP 2012/11/15) The determined SFHs
;                 are highly unstable and sensitive to degeneracies
;                 and the need for interations with various
;                 decompositions adds to the instability. The user can
;                 himself decide if it is reasonable to use this
;                 approach. 
;
; PURPOSE:
;                  Determine a star formation history
; USAGE:
;                  aux_sfh, <nami>, <model>, <namo>, 
;                  /CLEAN, DETECT_THRESHOLD=detect, STATUS=status, _EXTRA=_extra
;
; DESCRIPTION:
;   The spectrum <nami> is decomposed in a series of bursts described 
;   in <model>. If this decomposition is considered as inacceptable,
;   the model is modified (e. g. consecutive bursts are merged),
;   and the procedure is iteratively called (re-entrant call).
;
; ARGUMENTS:
;    model - model component
;    nami - input name
;    namo - output name
;
; KEYWORDS:
;   CLEAN - perform the fit twice - clean and then fit
;   _EXTRA - keywords to be passed to ulyss
;
; HISTORY: 
;   this program was send by Mina (2010/05/26)
;   It is analysing some spectra to find the SFH in 3 bursts
;   PhP: I modified this program to read the files from GMOS/Spolaor
;        Do the fit in two passes (clean and then global minimisation)
;   
;   Rewrite the program to work for any number and characteristics of bursts
;   in the SFH decomposition (PhP 2010/05/30 and 31st).
;
;   Generelize the program to work for arbirary components (not
;   only ssps). For the moment the component should consiste of
;   ssps+other definitions (MK 2012/05/26)
;-
; CATEGORY:
;                  AUX

pro aux_sfh, nami, model, namo, DETECT_THRESHOLD=detect, STATUS=status, _EXTRA=_extra, CLEAN=clean

  common uly_path, uly_root

  if n_elements(nami) eq 0 then begin
     print, 'Usage: AUX_SFH, nami, model, namo ...'
     return
  endif

  if n_elements(detect) eq 0 then detect = 0
  status = 0

  if file_test(namo+'.res') eq 1 then file_delete, namo+'.res'

  if n_elements(clean) eq 1 then begin
     ulyss, nami, model, _EXTRA=_extra, FILE_OUT=namo, CLEAN=clean
     file_delete, namo+'.res'
     undefine, clean
     ulyss, namo+'.fits', model, FILE_OUT=namo, _EXTRA=_extra
  endif else ulyss, nami, model, FILE_OUT=namo, _EXTRA=_extra

  result = uly_solut_tread(namo+'.res')
  if size(result, /TYPE) ne 8 then begin
     message, /INFO, 'No fit solution'
     status = 1
     return
  endif

;; number of SSPs in the composite model
  ssps = where(model.eval_fun eq 'SSP', nbursts, COMPLEMENT=others, NCOMPLE=nothers)
  print, 'Number of bursts:', nbursts, ' for ', nami 
  age_limits = dblarr(2, nbursts)
  feh_limits = dblarr(2, nbursts)
  fixed = intarr(3, nbursts)
  age_ng = intarr(nbursts)
  feh_ng = intarr(nbursts)
  fixw = intarr(nbursts)
  age = dblarr(nbursts)
  feh = dblarr(nbursts)
  for k=0,nbursts-1 do begin
     age_limits[*,k] = (*model[k].para)[0].limits
     feh_limits[*,k] = (*model[k].para)[1].limits
     age_ng[k] = n_elements(*(*model[k].para)[0].guess)
     feh_ng[k] = n_elements(*(*model[k].para)[1].guess)
     fixed[*,k] = (*model[k].para).fixed
     if model[k].lim_weig[0] eq 0 and model[k].lim_weig[1] eq 0 then fixw[k] = 1
     age[k] = (*(result.cmp[k]).para)[0].value
     feh[k] = (*(result.cmp[k]).para)[1].value
  endfor

  print, 'age limits1:', exp(age_limits)
  print, 'fixed1:', fixed
  print, 'fixw1:', fixw
  print, 'age_ng:', age_ng
  print, 'feh_ng:', feh_ng

  iter = 0

  for k=0,nbursts-1 do begin
     print, 'check burst ', k
     print, 'weight', result.cmp[k].weight
     if result.cmp[k].weight gt 0 then begin ; check only if positive contrib
        print, 'fixed', (*model[k].para)[0].fixed
        if (*model[k].para)[0].fixed eq 0 then begin
           d = (*model[k].para)[0].limits - (*(result[0].cmp[k]).para)[0].value
           if d[0] gt -0.001 and k ge 1 then begin
              if result.cmp[k-1].l_weight/result.cmp[k].l_weight gt 10. then begin
                 print, 'cancel this burst to zero ', k
                 fixw[k] = 1
                 fixed[*,k] = [1, 1, 1]
                 iter = 1
                 break
              endif 
              print, 'merge ',k, ' with ', k-1 
;                fehl = feh_limits[1,k-1:k]
;                feh_limits[1,k-1] = [min(fehl), max(fehl)]
              if fixw[k-1] eq 1 then begin
                 print, 'free the weight of cmp ', k-1
                 fixw[k-1] = 0
              endif
              age_limits[1,k-1] = age_limits[1,k]
              if (*model[k-1].para)[0].fixed eq 1 then begin
                 age_limits[0,k-1] = (*(*model[k-1].para)[0].guess)[0]
                 fixed[0,k-1] = 0
              endif
              print, 'age_ng 1:', exp(age_ng[k-1:k])
              n = 1 & if age_ng[k-1] eq 1 or age_ng[k] eq 1 then n = 0
              age_ng[k-1] += age_ng[k] - n

              ngd = feh_ng[k-1:k]/(feh_limits[1,k-1:k]-feh_limits[0,k-1:k])
              fehl = feh_limits[*,k-1:k]
              feh_limits[*,k-1] = [min(fehl), max(fehl)]
              feh_ng[k] = fix(min(ngd) * (feh_limits[1,k]-feh_limits[0,k]))

              print, 'ng 2:',age_ng[k-1], feh_ng[k-1]
              if k+1 le nbursts-1 then begin
                 age_limits = [[age_limits[*,0:k-1]], [age_limits[*,k+1:*]]]
                 feh_limits = [[feh_limits[*,0:k-1]], [feh_limits[*,k+1:*]]]
                 age_ng = [age_ng[0:k-1], age_ng[k+1:*]]
                 feh_ng = [feh_ng[0:k-1], feh_ng[k+1:*]]
                 age = [age[0:k-1], age[k+1:*]]
                 feh = [feh[0:k-1], feh[k+1:*]]
                 fixed = [[fixed[*,0:k-1]], [fixed[*,k+1:*]]]
                 fixw = [fixw[0:k-1], fixw[k+1:*]]
              endif else fixed = fixed[*,0:k-1]
              iter = 1
              break
           endif else if d[1] lt 0.001 and k lt nbursts-1 then begin
              if result.cmp[k+1].l_weight/result.cmp[k].l_weight gt 10. then begin
                 print, 'cancel this burst to zero ', k
                 fixw[k] = 1
                 fixed[*,k] = [1, 1, 1]
                 iter = 1
                 break
              endif 
              print, 'merge ',k, ' with ', k+1 
              print, age_limits
              if (*model[k+1].para)[0].fixed eq 0 then $
                 age_limits[1,k] = age_limits[1,k+1] $
              else $
                 age_limits[1,k] = (*(*model[k+1].para)[0].guess)[0]
              n = 1 & if age_ng[k] eq 1 or age_ng[k+1] eq 1 then n = 0
              age_ng[k] += age_ng[k+1] - n

              ngd = feh_ng[k:k+1]/(feh_limits[1,k:k+1]-feh_limits[0,k:k+1])
              fehl = feh_limits[*,k:k+1]
              feh_limits[*,k] = [min(fehl), max(fehl)]
              feh_ng[k] = fix(min(ngd) * (feh_limits[1,k]-feh_limits[0,k])+.1)
              print, 'ngd, feh_limits', k, ngd, feh_limits[*,k], feh_ng[k]

              if k+2 le nbursts-1 then begin
                 age_limits = [[age_limits[*,0:k]], [age_limits[*,k+2:*]]]
                 feh_limits = [[feh_limits[*,0:k]], [feh_limits[*,k+2:*]]]
                 age_ng = [age_ng[0:k], age_ng[k+2:*]]
                 feh_ng = [feh_ng[0:k], feh_ng[k+2:*]]
                 age = [age[0:k], age[k+2:*]]
                 feh = [feh[0:k], feh[k+2:*]]
                 fixed = [[fixed[*,0:k]], [fixed[*,k+2:*]]]
                 fixw = [fixw[0:k], fixw[k+2:*]]
              endif else fixed = fixed[*,0:k]
              iter = 1
              break
           endif
        endif else print, 'cmp', k, ' is fixed'
     endif else print, 'cmp',k,' has no contribution'
  endfor
  if size(fixed,/N_DIM) eq 1 then nbursts = 1 else $
     nbursts = (size(fixed,/DIM))[1]
  print, 'nbursts2: ', nbursts
  print, 'age limits2:', exp(age_limits)
  print, 'age ng2:', exp(age_ng)
  print, 'fixed2:', fixed
  print, 'feh_ng2:', feh_ng

  if iter eq 0 then begin
     print, 'check for low weights '
     for k=0,nbursts-1 do begin
        print, 'check burst ', k
        print, 'weight', result.cmp[k].weight, result.cmp[k].e_weight
        if result.cmp[k].weight gt 0 and $
           result.cmp[k].weight lt 3 * result.cmp[k].e_weight then begin 
           print, 'set this burst to zero ', k
           fixw[k] = 1
           fixed[*,k] = [1, 1, 1]
           iter = 1
           break
        endif
     endfor
  endif


  if iter eq 0 then begin
     print, 'check for low weights2 '
     tw = total(result.cmp.l_weight)
     for k=0,nbursts-1 do begin
        print, 'check burst ', k, result.cmp[k].l_weight/tw
        if result.cmp[k].weight gt 0 and $
           result.cmp[k].l_weight lt detect*tw then begin 
           print, 'constraint this burst to zero ', k
           fixw[k] = 1
           fixed[*,k] = [1, 1, 1]
           iter = 1
           break
        endif
     endfor
  endif

  if iter eq 0 and nbursts gt 1 then begin
     print, 'check the metallicity of the old burst'
     for k=nbursts-2,0,-1 do begin
        if exp(age[k]) gt 2000 and feh[k] lt feh[k+1] then begin
           print, 'alter the FeH limit ', k, feh_limits[*,k]
           feh_limits[0,k] = min([feh_limits[0,k]+0.3, feh[k+1]])
           iter = 1
           break
        endif
     endfor
  endif

  if iter eq 1 then begin
     for k=0,nbursts-1 do begin
        al = exp(age_limits[*, k])
        if fixed[0,k] eq 1 then begin                                                  ; if fixed
           if exp(age[k]) gt al[0] and exp(age[k]) lt al[1] then ag = exp(age[k]) $    ; the previous solution
           else begin 
              if exp(age[k]) lt al[0] then ag = al[0]+1
              if exp(age[k]) gt al[1] then ag = al[1]-1
           endelse
        endif else begin 
           ag = exp(range(alog(al[0]),alog(al[1]), age_ng[k], /OPEN))
        endelse
        zl = feh_limits[*, k]
        if fixed[1,k] eq 1 then begin
           zg = feh[k]
        endif else begin
           zg = range(zl[0],zl[1], feh_ng[k], /OPEN)
        endelse
        if fixw[k] eq 1 then wl=[0.,0.] else undefine, wl
        print, 'model', k, ag, zg, fixed[*,k]
        print, 'limits', al, zl
        modlk = uly_ssp(MODEL=strtrim(strmid(model[0].descr, 6),2), AG=ag, AL=al, ZL=zl, ZG=zg, $
                        FIX=fixed[*,k], WL=wl)
        if k eq 0 then modl = modlk else modl = [modl, modlk]
        
     endfor
;; add other components (for example emission lines)
     if nothers gt 0 then modl = [modl,model[others]] 
     sfh, namo+'.fits', modl, namo, DETECT=detect, STATUS=status ; re-entrant call to sfh
     heap_free, modl
  endif else print, "Solution accepted for ", nami

  return
end
;----------------------- end ---------------------------