PRO dustem_init_params,model,pd,iv, $
                       fpd=fpd,fiv=fiv $
                       ,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims $
                       ,tied=tied,fixed=fixed $
                       ,polarization=polarization $
                       ,no_reset_plugin_structure=no_reset_plugin_structure $
                       ,force_reset=force_reset
;+
; NAME:
;    dustem_init_params
;
; PURPOSE:
; This routine initialises free and fixed parameters of both the
; dust models and plugins. It does some basic sanity checking. By
; default, it preserves parameters and related information that are
; already stored in memory. To re-initialise, run with the
; /force_reset keyword
;  
; CATEGORY:
;    DustEMWrap, Distributed, Mid-Level, Initialization
;
; CALLING SEQUENCE:
;    dustem_init_params, model, pd, iv, $
;          [,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,lolims=llims,tied=tied,fixed=fixed $
;          ,pol=pol,no_reset_plugin_structure=no_reset_plugin_structure,force_reset=force_reset]
;  
; INPUTS:
;    model -- the dust model
;    pd -- the parameter description vector
;    iv -- the initial values of the free parameters 
;
; OPTIONAL INPUT PARAMETERS:
;    fpd -- the fixed parameter description vector
;    fiv -- the initial values of the fixed parameters 
;    ulimed -- flag 1/0 if upper limits ON/OFF
;    llimed -- flag 1/0 if lower limits ON/OFF
;    ulims -- upper limit values
;    llims -- lower limit values
;
; OUTPUTS:
;    None
;
; OPTIONAL OUTPUT PARAMETERS:
;    None
;
; ACCEPTED KEY-WORDS:
;    polarization = check that the model makes a polarisation prediction
;    no_reset_plugin_structure = if set, does not reset the plugin structure
;    force_reset = if set, resets existing parameters in memory to the
;                  default values of the model (removes all trace of
;                  limits, no-default values, free/fixed from previous
;                  instructions)
;    help      = if set, print this help
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None
;
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
;
; PROCEDURES AND SUBROUTINES USED:
;   dustem_init_parinfo
;   dustem_init_plugins
;   dustem_init_fixed_params
;  
; EXAMPLES
;
; MODIFICATION HISTORY:
;    Written by AH October-2022
;    Significant changes by AH March-2024:
;       addition of force_reset keyword and nicer management of
;       existing parinfo in memory
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_init_params'
  goto,the_end
END

Npar=n_elements(pd)
Niv=n_elements(iv)

if Npar ne Niv then begin
   message, 'Number of parameters does not equal number of initial values provided',/info
   stop
end

if keyword_set(ulimed) then begin
   Nulimed=n_elements(ulimed)
   if Npar ne Nulimed then begin
      message, 'Number of parameters does not equal number of upper limit flags',/info
      stop
   end
end

if keyword_set(ulims) then begin
   Nulims=n_elements(ulims)
   if Npar ne Nulims then begin
      message, 'Number of upper limits must equal number of parameters (use ulimed to specify ON/OFF for each parameter)',/info
      stop
   end
   ; compare ulims and ivs
   on=where(ulimed eq 1,ct)
   if ct gt 0 then begin
      test=total(ulims[on] lt iv[on])
      if test ne 0 then begin
         message, 'Upper limits are less than requested initial values',/info
         stop
      end
      
   end
end

if keyword_set(llimed) then begin
   Nllimed=n_elements(llimed)
   if Npar ne Nllimed then begin
      message, 'Number of parameters does not equal number of lower limit flags',/info
      stop
   end
end

if keyword_set(llims) then begin
   Nllims=n_elements(llims)
   if Npar ne Nllims then begin
      message, 'Number of lower limits must equal number of parameters (use llimed to specify ON/OFF for each parameter)',/info
      stop
   end
   ; compare llims and ivs
   on=where(llimed eq 1,ct)
   if ct gt 0 then begin
      test=total(llims[on] gt iv[on])
      if test ne 0 then begin
         message, 'Lower limits are greater than requested initial values',/info
         stop
      end
      
   end
end

if keyword_set(llims) and keyword_set(ulims) then begin
   on=where(llimed eq 1 and ulimed eq 1,ct)
   if ct gt 0 then begin
      test=total(llims[on] gt ulims[on])
      if test ne 0 then begin
         message, 'Requested lower limits are greater than upper limits',/info
         stop
      end
   end
end


if keyword_set(fpd) then begin
   Nfpar=n_elements(fpd)
   Nfiv=n_elements(fiv)
   if Nfpar ne Nfiv then begin
      message, 'Number of fixed parameters does not equal number of fixed initial values',/info
      stop
   end
;== CHECK IF USER HAS MADE SOMETHING FIXED AND FREE
   match,strupcase(pd),strupcase(fpd),a,b,count=ct
   if ct ne 0 then begin
      message, 'Parameter(s) that are both free and fixed?',/info
      print,'Matching free parameters: '+pd[a]
      print,'Matching fixed parameters: '+fpd[b]
      stop
   end
end


exists=dustem_test_model_exists(model,/silent)
if exists ne 1 then $
   message,'Unknown dust model '+model

polexists=dustem_test_model_exists(model,/pol,/silent)
if keyword_set(polarization) and polexists eq 0 then $
   message,'You have requested polarization mode, but the '+model+' model does not predict polarization'

; basic sanity checking of PLUGINS
; to be written

;== PLAY NICE -- check what already exists in the parameter structure
defsysv,'!dustem_parinfo',exists=par_status

if keyword_set(force_reset) and !dustem_verbose NE 0 then $
   message, 'Reinitialising all previous parameter options (FORCE_RESET=1)',/info

if par_status eq 0 and !dustem_verbose NE 0 then $
   message, 'No existing parameters found',/info

if par_status eq 1 and not keyword_set(force_reset) then begin
   if  !dustem_verbose NE 0 then message,'Updating/adding model parameters',/info
   new_pd=pd & new_iv=iv & nnewpars=n_elements(new_pd)
   new_ulimed=intarr(nnewpars)
   new_llimed=intarr(nnewpars)
   new_ulims=intarr(nnewpars)
   new_llims=intarr(nnewpars)
   new_fixed=intarr(nnewpars)
   new_tied=intarr(nnewpars)
   if keyword_set(ulimed) then new_ulimed=ulimed
   if keyword_set(llimed) then new_llimed=llimed
   if keyword_set(ulims) then new_ulims=ulims
   if keyword_set(llims) then new_llims=llims
   if keyword_set(fixed) then new_fixed=fixed
   if keyword_set(tied) then new_tied=tied

   old_pd=(*(*!dustem_fit).param_descs)
   old_iv=(*(*!dustem_fit).param_init_values)
   old_parinfo=*!dustem_parinfo
   old_fixed=old_parinfo.fixed
   old_tied=old_parinfo.tied
   noldpars=n_elements(old_parinfo)
   old_ulimed=intarr(noldpars)
   old_llimed=intarr(noldpars)
   old_ulims=fltarr(noldpars)
   old_llims=fltarr(noldpars)
   for i=0,noldpars-1 do begin
      old_ulimed[i]=old_parinfo[i].limited[1]
      old_llimed[i]=old_parinfo[i].limited[0]
      old_ulims[i]=(old_parinfo[i].limits[1])*old_iv[i]
      old_llims[i]=(old_parinfo[i].limits[0])*old_iv[i]
   endfor
   ; check for duplicates 
   noldpars=n_elements(old_pd)
   match,strupcase(new_pd),strupcase(old_pd),a,b,count=ndups
   if ndups gt noldpars then begin
      message, 'Number of duplicates exceeds number of existing parameters. This should not happen ?'/info
      stop
   endif
   if ndups gt 0 and ndups le noldpars then begin
      if  !dustem_verbose NE 0 then begin
         message,string(ndups)+' matching parameters found. Updating values...',/info
         for i=0,ndups-1 do print,old_pd[b[i]]
      end
      remove,b,old_pd
      remove,b,old_iv
      remove,b,old_parinfo
      remove,b,old_ulimed
      remove,b,old_llimed
      remove,b,old_ulims
      remove,b,old_llims
      remove,b,old_fixed
      remove,b,old_tied
   endif
   ;merge the structures
   pd=[new_pd,old_pd]
   iv=[new_iv,old_iv]
   ulimed=[new_ulimed,old_ulimed]
   llimed=[new_llimed,old_llimed]
   ulims=[new_ulims,old_ulims]
   llims=[new_llims,old_llims]
   fixed=[new_fixed,old_fixed]
   tied=[new_tied,old_tied]
endif


;; stop
;; message,'line 262',/info
;; if ptr_valid((*!dustem_fit).param_descs) then print,(*(*!dustem_fit).param_descs)
;; if ptr_valid((*!dustem_fit).fixed_param_descs) then print,(*(*!dustem_fit).fixed_param_descs)
;; stop

;== INITIALIZE DUST MODEL PARAMETERS and PLUGIN FREE PARAMETERS
dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims,fixed=fixed,tied=tied

;; stop
;; message,'line 271',/info
;; if ptr_valid((*!dustem_fit).param_descs) then print,(*(*!dustem_fit).param_descs)
;; if ptr_valid((*!dustem_fit).fixed_param_descs) then print,(*(*!dustem_fit).fixed_param_descs)
;; stop

;== Initialize ANY FIXED PARAMETERS FOR DUST MODEL AND PLUGINS
IF keyword_set(fpd) THEN BEGIN
;== PLAY NICE -- check what already exists in the fixed parameter structure
   fpar_status=ptr_valid((*!dustem_fit).fixed_param_descs)
   
   if keyword_set(force_reset) and  !dustem_verbose NE 0 then $
      message, 'Reinitialising all previous fixed parameter options (FORCE_RESET=1)',/info
   
   if fpar_status eq 0  and  !dustem_verbose NE 0 then $
      message, 'No existing fixed parameters found',/info

   if fpar_status eq 1 and not keyword_set(force_reset) then begin
      if  !dustem_verbose NE 0 then message,'Updating/adding fixed parameters',/info
      new_fpd=fpd & new_fiv=fiv & nnewfixed=n_elements(new_fpd)
      
      old_fpd=(*(*!dustem_fit).fixed_param_descs)
      old_fiv=(*(*!dustem_fit).fixed_param_init_values)
      noldfixed=n_elements(old_fpd)
      
      match,strupcase(new_fpd),strupcase(old_fpd),a,b,count=ndups
      if ndups gt noldfixed then begin
         message, 'Number of duplicates exceeds number of fixed parameters. This should not happen ?'/info
         stop
      endif
      if ndups gt 0 and ndups le noldfixed then begin
         if  !dustem_verbose NE 0 then message,string(ndups)+' matching fixed parameters found. Updating values for '+old_fpd[b],/info
         remove,b,old_fpd
         remove,b,old_fiv
      endif
      fpd=[new_fpd,old_fpd]
      fiv=[new_fiv,old_fiv]
   endif
   dustem_init_fixed_params,fpd,fiv
ENDIF
;; stop
;; message,'line 312',/info
;; if ptr_valid((*!dustem_fit).param_descs) then print,(*(*!dustem_fit).param_descs)
;; if ptr_valid((*!dustem_fit).fixed_param_descs) then print,(*(*!dustem_fit).fixed_param_descs)
;; stop

;== INITIALIZE ANY PLUGINS
IF not keyword_set(no_reset_plugin_structure) THEN BEGIN
   dustem_init_plugins,pd,fpd=fpd
ENDIF
;; stop
;; message,'line 322',/info
;; if ptr_valid((*!dustem_fit).param_descs) then print,(*(*!dustem_fit).param_descs)
;; if ptr_valid((*!dustem_fit).fixed_param_descs) then print,(*(*!dustem_fit).fixed_param_descs)
;; stop

if !dustem_verbose NE 0 then message,'Finished initializing free and fixed parameters of dust model and plugins',/info
;;stop

the_end:

END