;
; Plots calibration statistics for the MXD against temperature.
; calibmxd4 has calibration done on both decline-corrected and uncorrected
; MXD data.  Statistics are for the corrected dataset only, and are
; artificially high for the calibration period, but are fine for the
; verification period.
; Infilling of regression coefficients is done, and then calibration applied
; to those boxes previously without 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='calibmxd3'+fnadd+'.idlsave'
; Gets: fdcalib, fdcorrect (and other stuff!)
restore,filename='calibmxd4'+fnadd+'.idlsave'
; Gets:  g,mxdyear,mxdnyr,fdcorr,fdalph,fdbeta,fdvexp,mxdfd2,$
;        fdrver,fdvver,fdseas,timey,fdcalibu,fdcalibc
;
nrowcol=[3,2]
dopanel=[1,1,1,1,1,1,1,1]  ; rcal, rver, REcal, REver, a, b, offset, binfill
dopanel=[0,0,0,1,0,0,0,0]  ; rcal, rver, REcal, REver, a, b, offset, binfill
dopanel=[0,0,0,0,0,0,0,1]  ; rcal, rver, REcal, REver, a, b, offset, binfill
panlab='('+['a','b','c','d','e','f','g','h']+') '
panlab='('+['d','b','c','d','e','f','g','h']+') '
panlab='('+['b','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=750
  !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)
cc2=[25,24,23,21,20]
;
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 (CORRECTED)'
  ;
  if jpan eq 0 then begin
    fdin={ fd: outfd, x: g.x, y: g.y, nx: g.nx, ny: g.ny }
    whizz_fd,fdin=fdin,fdout=f,limit=map.limit
  endif
  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]+'Correlation over verification period (CORRECTED)'
  ;
  if jpan eq 0 then begin
    fdin={ fd: outfd, x: g.x, y: g.y, nx: g.nx, ny: g.ny }
    whizz_fd,fdin=fdin,fdout=f,limit=map.limit
  endif
  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
;
levs=[-1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.9]
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 CORRECTED)'
  ;
  if jpan eq 0 then begin
    fdin={ fd: outfd, x: g.x, y: g.y, nx: g.nx, ny: g.ny }
    whizz_fd,fdin=fdin,fdout=f,limit=map.limit
  endif
  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.9]
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 (decline-adjusted)'
  ;
  if jpan eq 0 then begin
    fdin={ fd: outfd, x: g.x, y: g.y, nx: g.nx, ny: g.ny }
    whizz_fd,fdin=fdin,fdout=f,limit=map.limit
  endif
  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
;
if dopanel[ipan] then begin
  ;
  pause
  inter_boxfd,fdalph,g.x,g.y,$
    labels=labels,map=map,$
    levels=[-1.,-0.5,-0.3,-0.2,-0.1,0.,0.2],/scale,$
    title=panlab[jpan]+'Regression intercept (CORRECTED)'
  ;
  if jpan eq 0 then begin
    fdin={ fd: outfd, x: g.x, y: g.y, nx: g.nx, ny: g.ny }
    whizz_fd,fdin=fdin,fdout=f,limit=map.limit
  endif
  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,$
    levels=[0.4,0.8,0.9,1.0,1.1,1.2,1.4],/scale,$
    title=panlab[jpan]+slopetit+' (CORRECTED)'
  ;
  if jpan eq 0 then begin
    fdin={ fd: outfd, x: g.x, y: g.y, nx: g.nx, ny: g.ny }
    whizz_fd,fdin=fdin,fdout=f,limit=map.limit
  endif
  boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
  ;
  jpan=jpan+1
endif
ipan=ipan+1
;
print
print,'Regression slope'
print,total(finite(fdbeta)),' values'
print,total((fdbeta ge 0.9) and (fdbeta le 1.1)),' are between 0.9 and 1.1'
print,total((fdbeta ge 0.8) and (fdbeta le 1.2)),' are between 0.8 and 1.2'
;
; Now need to infill the 6 boxes without regression coefficients.  This is
; done using a gaussian-weighted mean of nearby coefficients.
;
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
;
; Now plot the infilled coefficients
;
if dopanel[ipan] then begin
  ;
  pause
  inter_boxfd,fdalph,g.x,g.y,$
    labels=labels,map=map,$
    levels=[-1.,-0.5,-0.3,-0.2,-0.1,0.,0.2],/scale,$
    title=panlab[jpan]+'Regression intercept (CORRECTED & INFILLED)'
  ;
  if jpan eq 0 then begin
    fdin={ fd: outfd, x: g.x, y: g.y, nx: g.nx, ny: g.ny }
    whizz_fd,fdin=fdin,fdout=f,limit=map.limit
  endif
  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=cc2,$
    levels=[0.85,0.95,1.0,1.05,1.15,1.25],/scale,$
    title=panlab[jpan]+slopetit+' (decline-adjusted & infilled)'
  ;
  if jpan eq 0 then begin
    fdin={ fd: outfd, x: g.x, y: g.y, nx: g.nx, ny: g.ny }
    whizz_fd,fdin=fdin,fdout=f,limit=map.limit
  endif
  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
  fdcalibu(*,*,iyr)=fdalph(*,*)+fdbeta(*,*)*fdcalib(*,*,iyr)
  fdcalibc(*,*,iyr)=fdalph(*,*)+fdbeta(*,*)*fdcorrect(*,*,iyr)
endfor
;
; Now save the data for later analysis
;
save,filename='calibmxd5'+fnadd+'.idlsave',$
  g,mxdyear,mxdnyr,fdcalibu,fdcalibc,mxdfd2,timey,fdseas
;
end