PRO check_phangs_ssps_isrf_prediction,show_map=show_map,save=save

;check_phangs_ssps_isrf_prediction
;check_phangs_ssps_isrf_prediction,/from_restore
;check_phangs_ssps_isrf_prediction,/from_classes_restore,/save

;must have run make_phangs_ssps_isrf_prediction.pro before

;from_restore:

;make_phangs_isrf_classes,bidon,/save
;make_phangs_isrf_classes,bidon

dustem_init

data_dir=!phangs_data_dir+'/ISRF/WORK/'
product_dir=!phangs_data_dir+'/ISRF/PRODUCTS/'

;filename=product_dir+'ngc0628_isrf_classes_one_ratio.fits'
filename=data_dir+'isrf_classes_one_ratio_on_ngc0628.fits'
class0=mrdfits(filename,1,h1)
classes=mrdfits(filename,2,h2)

lambir=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths)

;cgplot,isrf_wavelengths,*class0[0].ISRF_TEMPLATE,/xlog,/ylog
cgplot,isrf_wavelengths,class0[0].ISRF_TEMPLATE,/xlog,/ylog
FOR i=0L,29 DO BEGIN
	;cgoplot,isrf_wavelengths,*classes[i].ISRF_TEMPLATE
	cgoplot,isrf_wavelengths,classes[i].ISRF_TEMPLATE
ENDFOR
;stop

;from_classes_restore:
;restore,data_dir+'isrf_classes_one_ratio.sav',/verbose

file=data_dir+'ngc0628_isrf_min_prediction.sav'
st_info=file_info(file)
IF st_info.exists NE 1 THEN BEGIN
	message,'Could not find '+file,/continue
	stop
ENDIF
restore,file,/verb
;% RESTORE: Restored variable: ISRFS.
;% RESTORE: Restored variable: G0S.
;% RESTORE: Restored variable: OBJECT_DISTANCE.
;% RESTORE: Restored variable: OBJECT_THICKNESS.
;% RESTORE: Restored variable: USE_SOURCE_NAME.
file=data_dir+'ngc0628_astrosat_voronoi_prediction_fast.sav'
st_info=file_info(file)
IF st_info.exists NE 1 THEN BEGIN
	message,'Could not find '+file,/continue
	stop
ENDIF
restore,file,/verb
;% RESTORE: Restored variable: ACTUAL_SEDS.
;% RESTORE: Restored variable: PREDICTED_SEDS.
;% RESTORE: Restored variable: HREF.
;% RESTORE: Restored variable: ALL_SEDS_INDICES.
;% RESTORE: Restored variable: DO_NORMALIZE.
;% RESTORE: Restored variable: ALL_FILTERS.
;% RESTORE: Restored variable: NHVOR.

Nclasses=n_elements(classes)
refs_wavs=[classes[0].FLUX_RATIO_WAV]
reference_wavelength=classes[0].reference_wavelength
ratios=classes.flux_ratio

lambir=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths)

Nlamb=(size(ISRFS))[1]
Nvor=(size(ISRFS))[2]

n_ISRFS=fltarr(Nlamb,Nvor)
n_wavs=fltarr(Nlamb,Nvor)
ind=where(isrf_wavelengths GT 1.)
indd=ind[0]
FOR vid=0L,Nvor-1 DO BEGIN
  n_ISRFS[*,vid]=ISRFS[*,vid]/ISRFS[indd,vid]
  n_wavs[*,vid]=isrf_wavelengths
ENDFOR

;compute flux ratio for each voronoi bin
vor_ratios=fltarr(Nvor)
FOR i=0L,Nvor-1 DO BEGIN
	vor_ratios[i]=interpol(n_ISRFs[*,i],isrf_wavelengths,refs_wavs[0])
ENDFOR
;fill in template ISRF of each class
;method='one_spectrum'
method='average_spectrum'
empty_classes=intarr(Nclasses)
CASE method OF
	'one_spectrum':BEGIN ;This is selection based on a single ISRF
		FOR i=0L,Nclasses-1 DO BEGIN
			ind=where(vor_ratios GT ratios[i],count)
			classes[i].ISRF_template=ptr_new(n_ISRFS[*,ind[0]])
		ENDFOR
	END
	'average_spectrum':BEGIN ;This is based on an average
		delta=1.0 ;This is a bit arbitrary. Chosen so that there is no empty classes
		FOR i=0L,Nclasses-1 DO BEGIN
			;ind=where(vor_ratios GE ratios[i]-ratios[i]*delta and vor_ratios LE ratios[i]+ratios[i]*delta,count)
			ratio_min=classes[i].flux_ratio_min
			ratio_max=classes[i].flux_ratio_max
			;ind=where(vor_ratios GT ratios_min[i] and vor_ratios LE ratios_max[i],count)
			ind=where(vor_ratios GT ratio_min and vor_ratios LE ratio_max,count)
			IF count NE 0 THEN BEGIN
				classes[i].ISRF_template=ptr_new(la_mean(n_ISRFS[*,ind],dim=-1))
			ENDIF ELSE BEGIN
				message,'Found no voronoi bin in class '+strtrim(i+1,2),/continue
				;stop
				empty_classes[i]=1
			ENDELSE
		ENDFOR
	END
ENDCASE
;===== remove empty classes
ind=where(empty_classes EQ 0,count)
classes=classes[ind]
Nclasses=count
;ratios_min=ratios_min[ind]
;ratios_max=ratios_max[ind]
;===== Fill in class number and name
FOR i=0L,Nclasses-1 DO BEGIN
	classes[i].number=i+1       ;because class 0 will be for Mathis ISRF
	str='I('+strtrim(refs_wavs[0],2)+'/I('+strtrim(reference_wavelength,2)+')='+strtrim(ratios[i],2)
	classes[i].name=str
ENDFOR

stop

;classify Voronoi bin according to flux ratio
;This could also be done on chi2, instead of just ratio ...
vor_class=lonarr(Nvor)
FOR i=0L,Nvor-1 DO BEGIN
  ;ratio=interpol(n_ISRFs[*,i],lambir,refs_wavs[0])
  this_ratio=vor_ratios[i]
  ind=where(classes.flux_ratio_min LT this_ratio and classes.flux_ratio_max GE this_ratio,count)
  ;ind=where(ratios_min LT this_ratio and ratios_max GE this_ratio,count)
  xv=ind[0]
  ;xv=round(interpol(xratios,ratios,this_ratio))
  IF ptr_valid(classes[xv].voronoi_members) THEN BEGIN
     classes[xv].voronoi_members=ptr_new([*classes[xv].voronoi_members,i])
  ENDIF ELSE BEGIN
  	  classes[xv].voronoi_members=ptr_new([i])
  ENDELSE
  classes[xv].Nvor=classes[xv].Nvor+1
  this_chi2=total((n_ISRFS[*,i])^2-(*classes[xv].ISRF_template)^2)
  IF ptr_valid(classes[xv].voronoi_chi2) THEN BEGIN
  	 classes[xv].voronoi_chi2=ptr_new([*classes[xv].voronoi_chi2,this_chi2])
  ENDIF ELSE BEGIN
  	 classes[xv].voronoi_chi2=ptr_new([this_chi2])
  ENDELSE
  ;This is a voronoi vector of the belonging class
  vor_class[i]=xv
ENDFOR
;classes[xv].voronoi_chi2=ptr_new((*classes[xv].voronoi_chi2)[1:*])
;classes[xv].voronoi_members=ptr_new((*classes[xv].voronoi_members)[1:*])

;==== plot histogram of voronoi classes
res=histogram(vor_class,min=0,max=Nclasses-1,Nbin=Nclasses,locations=x)

window,win & win=win+1
!p.multi=[0,1,2]
cgplot,x,res,psym=10,xrange=[-1,Nclasses+1],/xstyle,xtit='Class',ytit='Npix'
cgplot,ratios,res,psym=10,xrange=[1.e-5,10],/xstyle,xtit='Ratio',ytit='Npix',/xlog

;==== make a sky map of the classes
window,win & win=win+1
map_classes=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
FOR i=0L,Nvor-1 DO BEGIN
    map_classes[*all_seds_indices[i]]=vor_class[i]
ENDFOR
obp=[1.1,0.,1.15,1]
image_cont20,map_classes,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,imrange=[-1,Nclasses+1],off_bar=obp,tit='Class map'

;=== plot ISRF templates of all classes
window,win & win=win+1
xtit='Wavelength [mic]'
ytit='Normalized ISRF'
tit='ISRF classes'
loadct,13
!p.multi=0
med_isrf=la_median(n_ISRFS,dim=-1)
cgplot,isrf_wavelengths,med_isrf,color='black',/xlog,/ylog,xrange=[1.e-2,1.e2],/xsty,yrange=[1.e-5,10.],/nodata,xtit=xtit,ytit=ytit,tit=tit
colors=fix(range_gen(Nclasses,[10,250]))
FOR i=0L,Nclasses-1 DO cgoplot,isrf_wavelengths,*classes[i].ISRF_template,color=colors[i]
cgoplot,isrf_wavelengths,med_isrf,color='black',thick=2.
cgoplot,isrf_wavelengths,*class0.isrf_template,color='black',thick=2.
cgoplot,replicate(refs_wavs,2),10^!y.crange,linestyle=2
cgoplot,replicate(reference_wavelength,2),10^!y.crange,linestyle=2


;stop

IF keyword_set(save) THEN BEGIN
	file_save=data_dir+'ngc0682_isrf_classes_one_ratio.sav'
	message,'About to write '+file_save,/continue
	message,'Make sure you are allowed to do that and type .c to proceed.',/continue
	save,classes,vor_class,map_classes,file=file_save,/verb
	message,'Saved '+file_save,/continue
ENDIF

message,'========== Finished computing ISRF classes',/continue

stop


END