;
; Plots calibration statistics for the MXD against temperature.
; Also infills the few grid boxes without calibration coefficients.
;
matchvar=1        ; 0=regression, 1=variance matching
if matchvar eq 0 then slopetit='Regression slope' else slopetit='Scaling factor'
;
if matchvar eq 0 then fnadd='_regress' else fnadd=''
restore,filename='calibmxd1'+fnadd+'.idlsave'
; Gets:  g,mxdyear,mxdnyr,fdcorr,fdalph,fdbeta,fdvexp,fdcalib,mxdfd2,$
;        fdrver,fdvver,timey
;
nrowcol=[3,2]
dopanel=[1,1,0,0,0,1,0,0]  ; rcal, rver, a, b, REcal, REver, offset, binfill
dopanel=[0,0,0,0,0,0,0,1]  ; rcal, rver, a, b, REcal, REver, offset, binfill
panlab='('+['a','b','c','d','e','f','g','h']+') '
ipan=0
jpan=0
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=nrowcol[0],ncol=nrowcol[1],layout='large'
if !d.name eq 'X' then begin
  window,ysize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=9
endelse
def_1color,20,color='red'
def_1color,21,color='orange'
def_1color,22,color='yellow'
def_1color,23,color='green'
def_1color,24,color='lblue'
def_1color,25,color='deepblue'
def_1color,26,color='vlpurple'
def_1color,28,color='black'
def_smearcolor,fromto=[26,28]
cc=28-findgen(9)
;
outfd=total(finite(mxdfd2),3)
list4=where(outfd eq 0)
list2=where(outfd gt 0)
outfd(list4)=4.
outfd(list2)=2.
;
map=def_map(/npolar) & map.limit(0)=25.
labels=def_labels(/off)
;
levs=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]
if dopanel[ipan] then begin
  ;
  inter_boxfd,fdcorr,g.x,g.y,$
    labels=labels,map=map,$
    c_colors=cc,$
    levels=levs,/scale,$
    title=panlab[jpan]+'Calibration correlation'
  ;
  fdin={ fd: outfd, x: g.x, y: g.y, nx: g.nx, ny: g.ny }
  whizz_fd,fdin=fdin,fdout=f,limit=map.limit
  boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
  ;
  jpan=jpan+1
endif
ipan=ipan+1
;
checkfd=fdcorr
print
print,'Calibration correlation'
for i = 0 , n_elements(levs)-2 do begin
  print,total(checkfd gt levs[i]),'>',levs[i]
endfor
;
levs=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]
if dopanel[ipan] then begin
  ;
  pause
  inter_boxfd,fdrver,g.x,g.y,$
    labels=labels,map=map,$
    c_colors=cc,$
    levels=levs,/scale,$
    title=panlab[jpan]+'Verification correlation'
  boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
  ;
  jpan=jpan+1
endif
ipan=ipan+1
;
checkfd=fdrver
print
print,'Verification correlation'
for i = 0 , n_elements(levs)-2 do begin
  print,total(checkfd gt levs[i]),'>',levs[i]
endfor
;
if dopanel[ipan] then begin
  ;
  pause
  inter_boxfd,fdalph,g.x,g.y,$
    labels=labels,map=map,$
    c_colors=cc,$
    levels=0.1*[-1.1,-0.6,-0.4,-0.2,-0.1,0.,0.1,0.2,0.4,1.1],/scale,$
    title=panlab[jpan]+'Regression intercept'
  boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
  ;
  jpan=jpan+1
endif
ipan=ipan+1
;
if dopanel[ipan] then begin
  ;
  pause
  inter_boxfd,fdbeta,g.x,g.y,$
    labels=labels,map=map,$
    c_colors=cc,$
    levels=[-1.,0.2,0.4,0.6,0.8,0.9,1.0,1.2,1.4,1.6],/scale,$
    title=panlab[jpan]+slopetit
  boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
  ;
  jpan=jpan+1
endif
ipan=ipan+1
;
levs=[-1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8]
if dopanel[ipan] then begin
  ;
  pause
  inter_boxfd,fdvexp,g.x,g.y,$
    labels=labels,map=map,$
    c_colors=cc,$
    levels=levs,/scale,$
    title=panlab[jpan]+'% variance explained (calibration)'
  boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
  ;
  jpan=jpan+1
endif
ipan=ipan+1
;
checkfd=fdvexp
print
print,'Calibration RE'
for i = 0 , n_elements(levs)-2 do begin
  print,total(checkfd gt levs[i]),'>',levs[i]
endfor
;
levs=[-1.1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8]
if dopanel[ipan] then begin
  ;
  pause
  inter_boxfd,fdvver,g.x,g.y,$
    labels=labels,map=map,$
    c_colors=cc,$
    levels=levs,/scale,$
    title=panlab[jpan]+'Verification RE'
  boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
  ;
  jpan=jpan+1
endif
ipan=ipan+1
;
checkfd=fdvver
print
print,'Verification RE'
for i = 0 , n_elements(levs)-2 do begin
  print,total(checkfd gt levs[i]),'>',levs[i]
endfor
;
; Now need to infill the 6 boxes without regression coefficients.  This is
; done using a gaussian-weighted mean of nearby coefficients.
; Also infill the mean offsets (which are used to determine the mean level,
; rather than the regression intercept, anyway!).
;
print,'Infilling coefficients'
allx=fltarr(g.nx,g.ny)
for iy = 0 , g.ny-1 do allx(*,iy)=g.x(*)
ally=fltarr(g.nx,g.ny)
for ix = 0 , g.nx-1 do ally(ix,*)=g.y(*)
statx=allx(list2)
staty=ally(list2)
;
statval=fdalph(list2)
misslist=where(finite(statval) eq 0,nmiss)
if nmiss gt 0 then begin
  fd_extend,statval,statx,staty,search=1000.,wavelen=400.
  fdalph(list2(misslist))=statval(misslist)
endif
;
statval=fdbeta(list2)
misslist=where(finite(statval) eq 0,nmiss)
if nmiss gt 0 then begin
  fd_extend,statval,statx,staty,search=1000.,wavelen=400.
  fdbeta(list2(misslist))=statval(misslist)
endif
;
statval=fdcaloffset(list2)
misslist=where(finite(statval) eq 0,nmiss)
if nmiss gt 0 then begin
  fd_extend,statval,statx,staty,search=1000.,wavelen=400.
  fdcaloffset(list2(misslist))=statval(misslist)
endif
;
; Now plot the infilled coefficients
;
if dopanel[ipan] then begin
  ;
  pause
  inter_boxfd,fdcaloffset,g.x,g.y,$
    labels=labels,map=map,$
    c_colors=cc,$
    levels=0.5*[-1.1,-0.6,-0.4,-0.2,-0.1,0.,0.1,0.2,0.4,1.1],/scale,$
    title=panlab[jpan]+'Calibration mean offset (infilled)'
  boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
  ;
  jpan=jpan+1
endif
ipan=ipan+1
;
if dopanel[ipan] then begin
  ;
  pause
  inter_boxfd,fdbeta,g.x,g.y,$
    labels=labels,map=map,$
    c_colors=cc,$
    levels=[0.,0.2,0.4,0.6,0.7,0.8,1.0,1.2,1.4,1.6],/scale,$
    title=panlab[jpan]+slopetit+' (infilled)'
  boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
  ;
  jpan=jpan+1
endif
ipan=ipan+1
;
; Now compute the calibrated values now that all boxes have coefficients
;
;for iyr = 0 , mxdnyr-1 do begin
;  fdcalib(*,*,iyr)=fdalph(*,*)+fdbeta(*,*)*mxdfd2(*,*,iyr)
;endfor
print,'Making reconstruction'
i=0
for ix = 0 , g.nx-1 do begin
  for iy = 0 , g.ny-1 do begin
    if finite(fdbeta[ix,iy]) then begin
      ts1=fdbeta[ix,iy]*reform(mxdfd2[ix,iy,*])
      mkanomaly,ts1,mxdyear,refperiod=calper
      fdcalib[ix,iy,*]=ts1+fdcaloffset[ix,iy]
    endif
  endfor
endfor
print
;
; Now save the data for later analysis
;
save,filename='calibmxd2'+fnadd+'.idlsave',$
  g,mxdyear,mxdnyr,fdcorr,fdalph,fdbeta,fdvexp,fdcalib,mxdfd2,$
  fdrver,fdvver,timey,fdseas,fdcaloffset,calper
;
end