PRO dustem_marginalize_grid,fixed_parameters_description $
							   					 ,fixed_parameters_values $
												   ,seds $
												   ,Nseds $
												   ,params $
												   ,Nparams $
												   ,table_params $
												   ,table_pmins $
												   ,table_pmaxs $
												   ,method=method $
												   ,help=help

;+
; NAME:
;    dustem_marginalize_grid
; PURPOSE:
;    marginalize a grid on given fixed parameter
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_marginalize_grid,fixed_parameters_description,fixed_parameters_values,seds,Nseds,Nfilters,Nparams,table_params,table_pmins,table_pmaxs
; INPUTS:
;    fixed_parameters_description : fixed parameters description
;    fixed_parameters_values      : fixed parameters values
;    seds                         : sed grid to be marginalized (also an output)
;    Nseds                        : Number of seds (also an output)
;    params                       : sed grid parameter values to be marginalized (also an output)
;    Nparams                      : Number of parameters in sed grid (also an output)
;    table_params                 : description of parameters in sed grid (also an output)
;    table_pmins                  : min values of parameters in sed grid (also an output)
;    table_pmaxs                  : max values of parameters in sed grid (also an output)
; OPTIONAL INPUT PARAMETERS:
;    NONE
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    seds                         : sed grid to be marginalized (also an output)
;    Nseds                        : Number of seds (also an output)
;    Nparams                      : Number of parameters in sed grid (also an output)
;    table_params                 : description of parameters in sed grid (also an output)
;    table_pmins                  : min values of parameters in sed grid (also an output)
;    table_pmaxs                  : max values of parameters in sed grid (also an output)
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    ;Caution: grids should have regularly spaced parameter values for this to work
; PROCEDURE:
;    None
; EXAMPLES
;    
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard (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_marginalize_grid'
  goto,the_end
ENDIF


N_fixed_params=n_elements(fixed_parameters_description)

;Nparams=!dustem_grid.Nparams

use_method='NEAREST'
IF keyword_set(method) THEN BEGIN
	use_method=method
ENDIF

Nfilters=(*!dustem_grid).Nfilters
filters=(tag_names(((*!dustem_grid).st_seds)[0]))[0:Nfilters-1]
FOR i=0L,Nfilters-1 DO BEGIN
	filters[i]=strmid(filters[i],1,1000)
ENDFOR
;Nfilters=n_elements(filters)

use_dist_min_perc=1.   ;parameter values which are a this percent different will be considered equal in marginalization ...

CASE use_method OF
	'NEAREST': BEGIN
		FOR i=0L,N_fixed_params-1 DO BEGIN
			ind=where(table_params EQ fixed_parameters_description[i],count,complement=indc,ncomplement=nc)
			IF count NE 0 THEN BEGIN
				;fixed_parameter_value=fixed_parameters_values[ind[0]]
				fixed_parameter_value=fixed_parameters_values[i]
				;compute absolute distance between table parameter and the requested fixed value
				this_param_values=params[*].(ind[0])
				;=== This is a trick to make equal parameter values which are different below a percent level.
				;=== Needed to allow for slightly different values of a given parameter in merged tables (at numerical accuracy)
				this_param_values=1.*round(this_param_values/this_param_values[0]*use_dist_min_perc*100.)/100.*this_param_values[0]
				;order=sort(this_param_values)
				;this_param_indiv_values=this_param_values[order]
				;un=uniq(this_param_indiv_values)
				;this_param_indiv_values=this_param_indiv_values[un]
				;print,this_param_indiv_values
				;compute distance between provided parameter and parameter values
				distance=abs(fixed_parameter_value-this_param_values)
				;pick the smalest distance
				indd=where(distance EQ min(distance),ccount)
				;stop
				IF ccount EQ 0 THEN BEGIN
					message,'This should not happen',/continue
					stop
				ENDIF
				seds=seds[indd]
				params=params[indd]
				Nseds=ccount
			ENDIF
		ENDFOR
	END
	'RELLIKELYHOOD':BEGIN
		Nfilters=n_tags(seds[0])-1
		FOR i=0L,N_fixed_params-1 DO BEGIN
				ind=where(table_params EQ fixed_parameters_description[i],count,complement=indc,ncomplement=nc)
				params_v=params[*].(ind[0])
				sigma=(params_v*0.1)^2
				IF count NE 0 THEN BEGIN
					fixed_parameter_value=fixed_parameters_values[i]
					print,min(params_v),fixed_parameter_value,max(params_v)
					;compute chi2 between table parameter and the requested fixed value
					chi2=(fixed_parameter_value-params_v)^2/sigma
					min_chi2=min(chi2)
					weights=exp(-1.d0*chi2/min_chi2/2.)   ;relative Likelyhood
					;cgplot,params_v,weights,psym=-4,/ylog
					;distance=abs(fixed_parameter_value-params[*].(ind[0]))
				  ;indd=where(distance EQ min(distance),ccount)
					indd=where(chi2 EQ min(chi2),ccount)
				  ;print,ccount
					;stop
					IF ccount EQ 0 THEN BEGIN
						message,'This should not happen',/continue
						stop
					ENDIF
					in_seds=seds[indd]
					out_seds=in_seds
					tot_weights=total(weights)
					;stop
					;FOR i=0L,count-1 DO BEGIN
					ind_not_zero=where(weights NE 0.,count_not_zero)
		      FOR k=0L,n_elements(out_seds)-1 DO BEGIN
						FOR j=0L,Nfilters-1 DO BEGIN
							out_seds[k].(j)=total(seds[ind_not_zero].(j)*weights[ind_not_zero])/tot_weights
						ENDFOR
					ENDFOR
					;stop
		      ;==== recompute normalization
		      FOR k=0L,n_elements(out_seds)-1 DO BEGIN
		      	total_value=0.D0
		      	FOR j=0L,Nfilters-1 DO BEGIN
		      		total_value=total_value+out_seds[k].(j)
		      	ENDFOR
		      	out_seds[k].(Nfilters)=total_value
		      ENDFOR
					seds=out_seds
					params=params[indd]
					params.(ind[0])=fixed_parameter_value
					Nseds=ccount
					;stop
				ENDIF
		ENDFOR
	END
	'PREDICTOR':BEGIN
		;stop
		FOR i=0L,N_fixed_params-1 DO BEGIN
			ind=where(table_params EQ fixed_parameters_description[i],count,complement=indc,ncomplement=nc)
			IF count NE 0 THEN BEGIN
				fixed_parameter_value=fixed_parameters_values[i]
				;compute distance between table parameter and the requested fixed value
				distance=abs(fixed_parameter_value-params[*].(ind[0]))
				indd=where(distance EQ min(distance),ccount)             ;This will be ne thew number of parameter combinations
				;print,ccount
				out_seds=seds[indd]
        FOR k=0L,ccount-1 DO BEGIN
        	print,k,ccount
        	params_v=[params[indd[k]].(0)]
        	FOR kk=1L,Nparams-1 DO params_v=[params_v,params[indd[k]].(kk)]          
					params_v[ind[0]]=fixed_parameters_values[i]
					parameters_description=(*!dustem_grid).params
					parameters_values=params_v
					;stop
					out_sed=dustem_brute_force_fit_prediction_jpb(parameters_description $
											  ,parameters_values $
											  	;,filters=filters $
							   				  ,best_sed =best_sed $
							   				  ,best_parameters =best_parameters $
							   				  ,chi2_p=chi2_p  $
							   				  ,rchi2_p=rchi2_p  $
							   				  ,total_chi2s_p=total_chi2s_p $
							   				  ,status=status $
							   				  ,show_sed=show_sed $
							   				  ,nostop=nostop)
						;stop
					  FOR j=0L,Nfilters-1 DO BEGIN
							out_seds[k].(j)=out_sed[j]
						ENDFOR
						;stop
        	ENDFOR
					seds=out_seds
					params=params[indd]
					params.(ind[0])=fixed_parameter_value
					Nseds=ccount
					;stop
			ENDIF
		ENDFOR
  END
ENDCASE

;stop

the_end:

END