;
; Calibrates and plots regional MXD against regional temperature
; Computes standard errors on the reconstructions too.
; DOES IT FOR AGE-BANDED DENSITY REGIONS AND NH
;
; Define periods for calibration
;
calper=[1881,1960]        ; calibration period
verper=[1961,1994]        ; verification period
;
dosmooth=0         ; 0=plot (& do errors for) annual data, 1=smoothed data
thalf=10
donh=0                    ; 0=do NH 1=do ALL
;
if donh eq 0 then nhtit='NH' else nhtit='ALL'
;
; Get the hemispheric temperature series
;
restore,filename='datastore/compbest_mxd_fixed1950.idlsave'
if donh eq 0 then nhtemp=reform(comptemp(*,3)) $
             else nhtemp=reform(comptemp(*,2))
; Get rid of pre-1871 temperatures
knh=where(timey lt 1871)
nhtemp(knh)=!values.f_nan
;
; Get all the age-banded MXD series (regions, hemi and back to 1400)
;
restore,filename='/cru/u2/f055/treeharry/banding/bandregnh_fixed.idlsave'
; Gets: timey,nyr,nhts,allts,alleps,nreg,regname
yrmxd=timey
fullmxd=allts
; Combine hemispheric MXD with regional MXD
nfull=n_elements(yrmxd)
newfm=fltarr(nfull,nreg+1)*!values.f_nan
newfm(*,0:nreg-1)=transpose(fullmxd(*,*))
newfm(*,nreg)=nhts(*)
fullmxd=newfm
;
; Get regional temperature series
;
restore,filename='datastore/regbest_mxd_fixed1950.idlsave'
; Gets: nreg,regname,regtemp,rawtemp,timey,bestmxd etc.
;
; Identify overlap between instrumental and MXD series
;
nyr=n_elements(timey)
kcomp=where( (yrmxd ge timey(0)) and (yrmxd le timey(nyr-1)) )
bestmxd=fullmxd(kcomp,*)
;
; Combine hemispheric and regional temperatures
;
regname=[regname,nhtit]
newrt=fltarr(nyr,nreg+1)*!values.f_nan
newrt(*,0:nreg-1)=regtemp(*,*)
newrt(*,nreg)=nhtemp(*)
regtemp=newrt
nreg=nreg+1
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=3,layout='large',/landscape
if !d.name eq 'X' then begin
  window,xsize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=9
endelse
def_1color,10,color='lblue'
def_1color,11,color='red'
def_1color,12,color='vlblue'
;
range0=[-2.5,-2.0,-3.5,-2.5,-1.5,-0.5,-2.5,-2.0,-3.0,-1.5]
range1=[ 1.5, 2.0, 2.5, 1.5, 1.0, 1.0, 1.5, 1.5, 1.5, 1.0]
;
calregts=fltarr(nfull,nreg)*!values.f_nan
;
; Process each region separately
;
for ireg = 0 , nreg-1 do begin
  print,regname(ireg)
  ;
  ; Extract series
  ;
  tstemp=reform(regtemp(*,ireg))
  tsmxd=reform(bestmxd(*,ireg))
  tsfull=reform(fullmxd(*,ireg))
  ;
  ; Identify calibration and verification subsets
  ;
  calkl=where( finite(tstemp+tsmxd) and $
               (timey ge calper(0)) and (timey le calper(1)) , ncal )
  verkl=where( finite(tstemp+tsmxd) and $
               (timey ge verper(0)) and (timey le verper(1)) , nver )
  print,'Calibration period:'+string(calper,format='(2I5)')+$
    '  length:'+string(ncal,format='(I5)')
;  print,'Verification period:'+string(verper,format='(2I5)')+$
;    '  length:'+string(nver,format='(I5)')
  ;
  ; Calibrate the MXD
  ;
  mkcalibrate,tsmxd(calkl),tstemp(calkl),fitcoeff=fitcoeff,autocoeff=autocoeff
  print,'RAW'
  print,'     r   alpha    beta SEalpha  SEbeta    sig    rmse'
  print,fitcoeff,format='(F6.2,4F8.4,F7.2,F8.4)'
  print,' a_mxd a_tem a_resid   ess'
  print,autocoeff,format='(3F6.2,F8.1)'
  ; Scale up the standard error on the regression coefficient if there is
  ; autocorrelation in the residuals
  aresid=autocoeff(2)
  serialfac=(1.+abs(aresid))/(1.-abs(aresid))
  sebeta=fitcoeff(4)*serialfac
  print,'SEbeta scaled up=',sebeta
  ;
  ; Generate calibrated record and uncertainties etc.
  ; Use the extended MXD record (after checking it matches the short one!)
  ;
  xxx=tsfull(kcomp)
  yyy=tsmxd
  tse=total( (xxx-yyy)^2 , /nan )
;  print,tse
  if tse gt 0.03 then message,'Series do not match'
  calmxd=tsfull*fitcoeff(2)+fitcoeff(1)
  calregts(*,ireg)=calmxd(*)
  ;
  ; Now compute the one standard error for each reconstructed value
  ; (optionally follow procedure for smoothed data)
  ;
  ; To do so, need the predictor in terms of anomalies wrt to its mean
  ; over the calibration period
  predanom=tsfull
  mkanomaly,predanom,yrmxd,refperiod=calper,refmean=predmean
  print,'Predictor mean over cal. per. is',predmean
  ;
  ; SE due to uncertainty in (a) slope and (b) intercept coefficients,
  ; and (c) due to unexplained variance
  ;
  filter_cru,/nan,thalf,tsin=tstemp,tslow=templow
  filter_cru,/nan,thalf,tsin=calmxd,tslow=callow
  if dosmooth ne 0 then begin
    filter_cru,/nan,thalf,tsin=predanom,tslow=lowpred,weight=filwt
    nwt=n_elements(filwt)
    allwt=[reverse(filwt(1:nwt-1)),filwt]
    se_a=lowpred*sebeta
    varfac=serialfac*sqrt(total(allwt^2)) < 1. ; don't let errors get bigger
    print,'Unexplained variance scaled (due to smoothing) by',varfac
    se_c=fitcoeff(6)*varfac
    prets=callow
  endif else begin
    se_a=predanom*sebeta
    se_c=fitcoeff(6)
    prets=calmxd
  endelse
  se_b=fitcoeff(3)
  serr=sqrt( se_a^2 + se_b^2 + se_c^2 )
  ;
  if regname(ireg) eq nhtit then !p.multi=[0,1,2,0,0]
  pause
  ;
  plot,timey,tstemp,/nodata,$
    /xstyle,xrange=[1400,1994],xtitle='Year',$
    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
    /ystyle,yrange=[range0(ireg),range1(ireg)],$
    title=regname(ireg)
  ;
  xfill=[yrmxd,reverse(yrmxd)]
  yfill=[prets-serr*2.,reverse(prets+serr*2.)]
  kl=where(finite(yfill),nkeep)
  yfill=yfill(kl)
  xfill=xfill(kl)
  polyfill,xfill,yfill,color=12
  ;
  xfill=[yrmxd,reverse(yrmxd)]
  yfill=[prets-serr,reverse(prets+serr)]
  kl=where(finite(yfill),nkeep)
  yfill=yfill(kl)
  xfill=xfill(kl)
  polyfill,xfill,yfill,color=10
  ;
  axis,yaxis=0,/ystyle,ytickformat='nolabels'
  axis,yaxis=1,/ystyle,ytickformat='nolabels'
  ;
  if dosmooth eq 0 then begin
    oplot,timey,tstemp
    oplot,yrmxd,calmxd,color=11
  endif
  ;
  oplot,!x.crange,[0.,0.],linestyle=1
  oplot,timey,templow,thick=3
  oplot,yrmxd,callow,thick=3,color=11
  ;
  if (regname(ireg) eq nhtit) and (dosmooth eq 0) then begin
  !p.multi[0]=0
  plot,timey,tstemp,thick=2,$
    /xstyle,xrange=[1870,1965],xtitle='Year',$
    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
    title=regname(ireg)
  oplot,yrmxd,calmxd,color=11,thick=2
filter_cru,/nan,10.,tsin=tstemp,tslow=t1
;oplot,timey,t1,thick=6
filter_cru,/nan,10.,tsin=calmxd,tslow=t1
;oplot,yrmxd,t1,thick=6,color=11
;  oplot,yrmxd,calmxd+serr,color=10
;  oplot,yrmxd,calmxd-serr,color=10
;  oplot,yrmxd,calmxd+serr*2.,color=10
;  oplot,yrmxd,calmxd-serr*2.,color=10
  oplot,!x.crange,[0.,0.],linestyle=1
  endif
  ;
endfor
;
end