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 $
                       ,help=help,debug=debug

;+
; 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

; first decide what needs to be done
fpar_status=ptr_valid((*!dustem_fit).fixed_param_descs)
par_status=ptr_valid((*!dustem_fit).param_descs)
fpd_set = keyword_set(fpd)
pd_set = keyword_set(pd)

old_pd=[] & old_iv=[] & old_parinfo=[]
old_ulimed=[] & old_llimed=[] & old_ulims=[] & old_llims=[]
old_fixed=[] & old_tied=[]
nold=0.

new_pd=[] & new_iv=[] & new_parinfo=[]
new_ulimed=[] & new_llimed=[] & new_ulims=[] & new_llims=[]
new_fixed=[] & new_tied=[]
nnew=0.

old_fpd=[] & old_fiv=[]
nold_fixed=0.

new_fpd=[] & new_fiv=[]
nnew_fixed=0.

if fpar_status eq 1 then begin
   old_fpd=(*(*!dustem_fit).fixed_param_descs)
   old_fiv=(*(*!dustem_fit).fixed_param_init_values)
   nold_fixed=n_elements(old_fpd)
end

if par_status eq 1 then begin
   old_pd=(*(*!dustem_fit).param_descs)
   old_iv=(*(*!dustem_fit).param_init_values)
   nold=n_elements(old_pd)
   old_parinfo=*!dustem_parinfo
   old_fixed=old_parinfo.fixed
   old_tied=old_parinfo.tied
   noldpars=n_elements(old_parinfo)
   old_ulimed=intarr(nold)
   old_llimed=intarr(nold)
   old_ulims=fltarr(nold)
   old_llims=fltarr(nold)
   for i=0,nold-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
end

if pd_set eq 1 then begin
   new_pd=pd & new_iv=iv & nnew=n_elements(new_pd)
   new_ulimed=intarr(nnew)
   new_llimed=intarr(nnew)
   new_ulims=intarr(nnew)
   new_llims=intarr(nnew)
   new_fixed=intarr(nnew)
   new_tied=strarr(nnew)
   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
end

if fpd_set eq 1 then begin
   new_fpd=fpd & new_fiv=fiv & nnew_fixed=n_elements(new_fpd)
end

if keyword_set(force_reset) then goto, initialise

if pd_set eq 1 then begin
   if par_status eq 1 then begin
      match,strupcase(new_pd),strupcase(old_pd),a,b,count=ndups
      if ndups gt nold then begin
         message, 'Number of duplicates exceeds number of existing free parameters. This should not happen ?'/info
         stop
      end else if ndups eq 0 then begin
         message,'No duplicate free parameters found...',/info
      end else if ndups eq nold then begin
         message,'Updating all old values of free parameters...',/info
         old_pd=[]
         old_iv=[]
         old_parinfo=[]
         old_ulimed=[]
         old_llimed=[]
         old_ulims=[]
         old_llims=[]
         old_fixed=[]
         old_tied=[]
      end else if ndups gt 0 and ndups lt nold then begin
         message,string(ndups)+' matching parameters found. Updating values...',/info
         for i=0,ndups-1 do print,old_pd[b[i]]
         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
   endif

   if fpar_status eq 1 then begin
      match,strupcase(new_pd),strupcase(old_fpd),a,b,count=ndups
      if ndups gt nold_fixed then begin
         message, 'Number of duplicates exceeds number of existing fixed parameters. This should not happen ?'/info
         stop
      end else if ndups eq 0 then begin
         message,'No duplicate free parameters found in old fixed params...',/info
      end else if ndups eq nold_fixed then begin
         message,'Updating all old values of fixed parameters found in new free params...',/info
         old_fpd=[]
         old_fiv=[]
      end else if ndups gt 0 and ndups lt nold_fixed then begin
         message,string(ndups)+' matching fixed parameters found. Updating values...',/info
         for i=0,ndups-1 do print,old_fpd[b[i]]
         remove,b,old_fpd
         remove,b,old_fiv
      endif
   endif
endif

if fpd_set eq 1 then begin
   if par_status eq 1 then begin
      match,strupcase(new_fpd),strupcase(old_pd),a,b,count=ndups
      if ndups gt nold then begin
         message, 'Number of duplicates exceeds number of existing fixed parameters. This should not happen ?'/info
         stop
      end else if ndups eq 0 then begin
         message,'No duplicate fixed parameters found  in old free params...',/info
      end else if ndups eq nold then begin
         message,'Updating all old values of free parameters found in fixed params...',/info
         old_pd=[]
         old_iv=[]
         old_parinfo=[]
         old_ulimed=[]
         old_llimed=[]
         old_ulims=[]
         old_llims=[]
         old_fixed=[]
         old_tied=[]
      end else if ndups gt 0 and ndups lt nold then begin
         message,string(ndups)+' matching parameters found. Updating values...',/info
         for i=0,ndups-1 do print,old_pd[b[i]]
         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
   endif

   if fpar_status eq 1 then begin
      match,strupcase(new_fpd),strupcase(old_fpd),a,b,count=ndups
      if ndups gt nold_fixed then begin
         message, 'Number of duplicates exceeds number of existing fixed parameters. This should not happen ?'/info
         stop
      end else if ndups eq 0 then begin
         message,'No duplicate fixed parameters found in old fixed params...',/info
      end else if ndups eq nold_fixed then begin
         message,'Updating all old values of fixed parameters...',/info
         old_fpd=[]
         old_fiv=[]
      end else if ndups gt 0 and ndups lt nold_fixed then begin
         message,string(ndups)+' matching fixed parameters found. Updating values...',/info
         for i=0,ndups-1 do print,old_fpd[b[i]]
         remove,b,old_fpd
         remove,b,old_fiv
      endif
   endif
endif

;merge the pd 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]

fpd=[new_fpd,old_fpd]
fiv=[new_fiv,old_fiv]

initialise:
;== INITIALIZE DUST MODEL PARAMETERS and PLUGIN FREE PARAMETERS
; currently this works without checking if PD is non-NULL because we
; require there alway to be a PD vector... one day we should change this
npd=n_elements(pd) & nfpd=n_elements(fpd)

if npd ge 1 then dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims,fixed=fixed,tied=tied

if npd eq 0 and nold gt 0 then $
   dustem_init_parinfo,pd,iv,/clear

if nfpd ge 1 then dustem_init_fixed_params,fpd,fiv

if nfpd eq 0 and nold_fixed gt 0 then $
   dustem_init_fixed_params,fpd,fiv,/clear

;== INITIALIZE ANY PLUGINS
IF not keyword_set(no_reset_plugin_structure) THEN BEGIN
   dustem_init_plugins,pd,fpd=fpd
ENDIF

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

the_end:

END