FUNCTION dustem_brute_force_fit,sed $

							   							 ,table_name $
							   							 ,filters $
							   							 ,fixed_parameters_description=fixed_parameters_description $
							   							 ,fixed_parameters_values=fixed_parameters_values $
							   							 ,normalize=normalize $
							   							 ,params_hit=params_hit $
							   							 ,params_uncertainties=params_uncertainties $
							   							 ,params_min=params_min $
							   							 ,params_max=params_max $
							   							 ,fact=fact $
							   							 ,chi2=chi2 $
							   							 ,rchi2=rchi2 $
							   							 ,best_sed=best_sed $
							   							 ,best_grid_sed=best_grid_sed $
							   							 ,show_sed=show_sed $
							   							 ,nostop=nostop $

							   							 ,reset=reset $

							   							 ,weighted_params=weighted_params $

							   							 ,weighted_sed=weighted_sed $
							   							 ,weighted_chi2=weighted_chi2 $
							   							 ,rweighted_chi2=rweighted_chi2 $

							   							 ,all_chi2s=all_chi2s $

							   							 ,weightning=weightning $
							   							 ,marginalize_method=marginalize_method $

							   							 ,status=status

;+
; NAME:
;    dustem_brute_force_fit
; PURPOSE:
;    does brut force fit from a pre-computed model grid
; CATEGORY:
;    DustEM
; CALLING SEQUENCE:

;    res=dustem_brute_force_fit(sed,table_name,filters[,/normalize][,fact=][,chi2=][,rchi2=][rchi2][best_sed=][,/show_sed][,/nostop][,/reset])

; INPUTS:

;    sed  = sed to be fit. Should contain all filters provided to the routine trhough input variable filter

;    table_name = name of the fits file to be used as a grid (see dustem_make_sed_table.pro for creation)
;    filters    = filters to be used in the grid table (only used if !dustem_grid does not exist or if /reset)

; OPTIONAL INPUT PARAMETERS:

;    fixed_parameters_description = if set, describes parameters of the table that should be considered fixed.
;    fixed_parameters_values      = values of parameters of the table that should be considered fixed.

; OUTPUTS:

;    res = best values of the free parameters. Not that linear parameters to the sed should be multiplied by fact.

; OPTIONAL OUTPUT PARAMETERS:

;    fact = best scaling factor. Such that best_sed*fact is for NH=1.e20 H/cm2 (as in the grid). IF sed is already normalized to NH=1e20 H/cm2, no need to set /normalize

;    chi2 = chi2 of the best fit
;    rchi2 = reduced chi2 of the best fit

;    best_grid_sed = best SED found in the grid
;    best_sed = best fit SED matching input SED (such that best_sed=best_grid_sed*fact)
;    params_hit= 1 if parameter hits the top, -1 if hits the bottom, 0 otherwise
;    params_uncertainties = 90% confidence level uncertainty on parameters
;    params_min = 90% lower limit of parameters
;    params_max = 90% upper limit of parameters

;    weighted_params = likelyhood exp(-chi2/2) weighted best parameter values

;    weighted_sed = likelyhood exp(-chi2/2) weighted best sed

; ACCEPTED KEY-WORDS:
;    help                  = if set, print this help

;    normalize = if set, the best match is done based on the shape of the SEDs only, instead of their absolute values. fact will be computed.
;    show_sed = if set, plot the seds
;    nostop   = if set, do not stop

; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The DustEMWrap IDL code must be installed
; PROCEDURE:

;    Scans the grid SEDs to find the best fit (lowest chi2)

; EXAMPLES

;    

; MODIFICATION HISTORY:

;    Written by JPB Jan 2024

;    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_brute_force_fit'
  params=0.

  goto,the_end
ENDIF

;stop

;sed must be normalized to NH=1e20 H/cm2

;==== test for the existence of a loaded grid and initialize if needed.

;==== only filters contained in the filters variable are kept in the grid, and the sum of grid seds are recomputed accordingly

defsysv,'!dustem_grid',0,exist=exist

IF exist EQ 0 or keyword_set(reset)THEN BEGIN

    dustem_define_grid,table_name,filters=filters,status=status
    IF status NE 1 THEN BEGIN
    	params=-1
    	goto,the_end
    ENDIF

ENDIF

;seds=!dustem_grid.seds

seds=(*!dustem_grid).st_seds
params=(*!dustem_grid).st_params
totals=(*!dustem_grid).sed_totals
Nseds=(*!dustem_grid).Nsed
Nfilters=(*!dustem_grid).Nfilters
Nparams=(*!dustem_grid).Nparams
table_params=(*!dustem_grid).params
table_pmins=(*!dustem_grid).pmin_values
table_pmaxs=(*!dustem_grid).pmax_values

use_weightning='rellikelyhood'
IF keyword_set(weightning) THEN BEGIN
	use_weightning=weightning
ENDIF

;=== marginalize table for fixed parameters, if any

;stop
use_marginalize_method='NEAREST'
IF keyword_set(marginalize_method) THEN use_marginalize_method=marginalize_method

IF keyword_set(fixed_parameters_description) THEN BEGIN

	dustem_marginalize_grid,fixed_parameters_description,fixed_parameters_values,seds,Nseds,params,Nparams,table_params,table_pmins,table_pmaxs,method=use_marginalize_method

ENDIF

;help,seds,Nseds,params,Nparams,table_params,table_pmins,table_pmaxs
;SEDS            STRUCT    = -> <Anonymous> Array[400]
;NSEDS           LONG      =          400
;PARAMS          STRUCT    = -> <Anonymous> Array[400]
;NPARAMS         LONG      =            3
;TABLE_PARAMS    STRING    = Array[3]
;TABLE_PMINS     DOUBLE    = Array[3]
;TABLE_PMAXS     DOUBLE    = Array[3]

;stop

;=== These are the full arrays needed for calculations

chi2s=dblarr(Nseds)

grid_seds=dblarr(Nseds,Nfilters)
total_seds=dblarr(Nseds)
grid_param_values=dblarr(Nseds,Nparams)

;=== compute sed grid Array
FOR i=0L,Nseds-1 DO BEGIN
	FOR j=0L,Nfilters-1 DO BEGIN

		;grid_seds[i,j]=seds[i].(j+Nparams)
		grid_seds[i,j]=seds[i].(j)

	ENDFOR

	;total_seds[i]=seds[i].total
	total_seds[i]=totals[i]

ENDFOR
;=== compute the parameter values Array
FOR i=0L,Nseds-1 DO BEGIN
	FOR j=0L,Nparams-1 DO BEGIN

		;grid_param_values[i,j]=seds[i].(j)
		grid_param_values[i,j]=params[i].(j)

	ENDFOR
ENDFOR
;stop

;=== loop over grid lines to compute chi2s

;we do actually store the total of grid sed for each grid sed into the grid fits file

;we could also store the total of observed seds in the observed sed save fits file (!)
facts=dblarr(Nseds)
facts[*]=1.d0

;stop

IF keyword_set(normalize) THEN BEGIN
	FOR i=0L,Nseds-1 DO BEGIN

		ref_sed=reform(grid_seds[i,*])
		;total_ref_sed1=total_seds[i]
		total_ref_sed2=total(ref_sed)
		;diff=abs(total_ref_sed1-total_ref_sed2)
		;print,diff
		;total_ref_sed=total_ref_sed1    ;This does not work. I'm not sure why, as it is supposed to be the same pre-computed ...
		total_ref_sed=total_ref_sed2    ;This does !!		
		;IF diff GT 0.1 THEN stop
		facts[i]=total(sed.stokesI)/total_ref_sed
    this_ref_sed=ref_sed*facts[i]
    this_sed=sed.stokesI
    this_var=sed.sigmaII
		;this_sed=sed.stokesI/facts[i]
		;this_var=sed.sigmaII/facts[i]^2

	  ;chi2s[i]=total((sed.stokesI-grid_seds[i,*])^2/sed.sigmaII)

	  chi2s[i]=total((this_sed-this_ref_sed)^2/this_var)

	ENDFOR
ENDIF ELSE BEGIN
	FOR i=0L,Nseds-1 DO BEGIN

	   chi2s[i]=total((sed.stokesI-grid_seds[i,*])^2/sed.sigmaII)

	ENDFOR
ENDELSE

all_chi2s=chi2s

chi2=min(chi2s)
ind=where(chi2s EQ chi2,count)

;stop
;print,facts[14406],facts[ind[0]]
;print,chi2s[14406],chi2s[ind[0]]

params=reform(grid_param_values[ind[0],*])  ;these are the best fit parameters

fact=facts[ind[0]]   ; This is the corresponding scaling factor

Nfreedom=Nfilters-Nparams-1   ;This is the degree of freedom

rchi2=chi2/Nfreedom        ;This is the reduced chi2
best_grid_sed=reform(grid_seds[ind[0],*])   ;This is the best SED fround in the grid.
best_sed=best_grid_sed*fact  ;This is the best scaled grid SED matching the input SED.

;=== compute Likelyhood weighted params values

epsilon=1.d-20

weighted_params=dblarr(Nparams)

CASE use_weightning OF
	'1/chi2'      :weights=1.d0/(chi2s+epsilon)
	'1/sqrt(chi2)':weights=1.d0/sqrt(chi2s+epsilon)
	'likelyhood'  :weights=exp(-1.d0*chi2s/2.)
	'rellikelyhood':weights=exp(-1.d0*(chi2s/chi2)/2.)
ENDCASE

;weights=1.d0/(chi2s+epsilon)

;weights=1.d0/sqrt(chi2s+epsilon)

;IF use_weightning EQ 'likelyhood' THEN weights=exp(-1.d0*chi2s/2.)    ;This is the likelyhood

tot_weights=total(weights)

;stop

FOR i=0L,Nparams-1 DO BEGIN

	;weighted_params[i]=total(reform(1.d0*grid_param_values[*,i])*weights)/tot_weights
	;Trying a log approach
  logweighted_param=total(reform(1.d0*alog10(grid_param_values[*,i])*weights))/tot_weights
	weighted_params[i]=10^logweighted_param

ENDFOR

;===compute best weighted SED
weighted_sed=fltarr(Nfilters)
FOR i=0L,Nfilters-1 DO BEGIN
	weighted_sed[i]=total(grid_seds[*,i]*weights)/tot_weights
ENDFOR
;=== computed weighted chi2
weighted_chi2=total((sed.stokesI-weighted_sed)^2/sed.sigmaII)
rweighted_chi2=weighted_chi2/Nfreedom
;stop

;stop

;==== Check bumping to table edges.

params_hit=intarr(Nparams)
FOR i=0L,Nparams-1 DO BEGIN

	diff_min=abs(params[i]-table_pmins[i])/table_pmins[i]*100.
	diff_max=abs(params[i]-table_pmaxs[i])/table_pmaxs[i]*100.
	IF diff_min LE 1. THEN params_hit[i]=-1
	IF diff_max LE 1. THEN params_hit[i]=1
  ;IF params[i] LE table_pmins[i] THEN params_hit[i]=-1
  ;IF params[i] GE table_pmaxs[i] THEN params_hit[i]=1

ENDFOR

;==== Compute parameter uncertainties

;stop

conf_level=90./100.

;print,Nfreedom,conf_level
IF Nfreedom NE 0. THEN BEGIN
	dchi2=delta_chi2(Nfreedom,conf_level)
ENDIF ELSE BEGIN
	message,'Nfreedom = 0. Beware of meaningless fit',/continue
	stop
	message,'seting Nfreedom = 1. Beware of meaningless fit',/continue
  dchi2=delta_chi2(1.,conf_level)
ENDELSE

ind=where(chi2s LE chi2+dchi2,count)

params_min=fltarr(Nparams)
params_max=fltarr(Nparams)

FOR i=0L,Nparams-1 DO BEGIN

	params_min[i]=min(grid_param_values[ind[0],i])
	params_max[i]=max(grid_param_values[ind[0],i])

ENDFOR
params_uncertainties=params_max-params_min

IF keyword_set(show_sed) THEN BEGIN
	xtit='Wavelength [mic]'
	ytit='SED (model=red,data=blue)'

	;should also plot uncertainties

	minv=la_min([sed.stokesI,best_sed])
	maxv=la_max([sed.stokesI,best_sed])
	cgplot,sed.wave,sed.stokesI,psym='Filled Circle',color='blue',/ylog,/ysty,xtit=xtit,ytit=ytit,/xlog,yrange=[minv,maxv]
	cgoplot,sed.wave,best_sed,psym='Filled Square',color='red'

	if not keyword_set(nostop) THEN stop

ENDIF

the_end:

RETURN,params
END