;
; Calibrates and plots regional MXD against regional temperature
; Computes standard errors on the reconstructions too.
; Errors are time dependent due to declining coverage early on.
;
dosave=1    ; 0=do not, 1=do save results (errors)
;
trv=0           ; selects tree-ring-variable: 0=MXD 1=TRW
case trv of
  0: begin
    fnvar='mxd'
    end
  1: begin
    fnvar='trw'
    end
endcase
titadd=strupcase(fnvar)

; Define periods for calibration
;
calper=[1881,1960]        ; calibration period
verper=[1961,1994]        ; verification period
;
dosmooth=1         ; 0=plot (& do errors for) annual data, 1=smoothed data
thalf=10
donh=1                    ; 0=do NH 1=do ALL
;
if donh eq 0 then nhtit='NH' else nhtit='ALL'
;
allfixedyr=[1400,1500,1600,1700,1800,1950,1988]
iord=[0,1,2,3,4,6,5]         ; order to process (keep 1950 last!)
nfixed=n_elements(allfixedyr)
useyr=[1400,1500,1600,1700,1800,1900,1983,1994]
;
; 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=14
endelse
def_1color,12,color='lsand'
def_1color,11,color='mblue'
def_1color,10,color='lyellow'
;
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]
;
; Get the extended hemispheric series
;
restore,filename='compts_'+fnvar+'_fixed1950.idlsave'
fullnhmxd=reform(compmxd(*,2))
;
; Get the extended MXD series
;
restore,filename='regts_'+fnvar+'_fixed1950.idlsave'
; Gets: timey,bestmxd etc.
yrmxd=timey
fullmxd=bestmxd
; Combine extended hemispheric MXD with these ones
nfull=n_elements(yrmxd)
newfm=fltarr(nfull,nreg+1)*!values.f_nan
newfm(*,0:nreg-1)=fullmxd(*,*)
newfm(*,nreg)=fullnhmxd(*)
fullmxd=newfm
;
; For the shorter calibration series, need to read in all the series
; generated by various fixed grids
;
for ifixed = 0 , nfixed-1 do begin
  ;
  fyr=allfixedyr(ifixed)
  fnadd=string(fyr,format='(I4)')
  ;
  ; Get the hemispheric series
  ;
  restore,filename='datastore/compbest_'+fnvar+'_fixed'+fnadd+'.idlsave'
  if donh eq 0 then nhtemp=reform(comptemp(*,3)) $
               else nhtemp=reform(comptemp(*,2))
  nhmxd=reform(compmxd(*,2,0))
  ; Get rid of pre-1871 temperatures
  knh=where(timey lt 1871)
  nhtemp(knh)=!values.f_nan
  ;
  ; Get regional series (temp, and truncated MXD)
  ;
  restore,filename='datastore/regbest_'+fnvar+'_fixed'+fnadd+'.idlsave'
  ; Gets: nreg,regname,regtemp,rawtemp,timey,bestmxd etc.
  ;
  ; Add in the hemispheric series
  ;
  if ifixed eq 0 then begin
    nyr=n_elements(timey)
    newrt=fltarr(nyr,nreg+1,nfixed)*!values.f_nan
    newrm=fltarr(nyr,nreg+1,nfixed)*!values.f_nan
  endif
  newrt(*,0:nreg-1,ifixed)=regtemp(*,*)
  newrt(*,nreg,ifixed)=nhtemp(*)
  newrm(*,0:nreg-1,ifixed)=bestmxd(*,*,1)
  newrm(*,nreg,ifixed)=nhmxd(*)
  ;
endfor
regtemp=newrt
bestmxd=newrm
regname=[regname,nhtit]
nreg=nreg+1
;
; Identify overlap between full and truncated MXD series
;
kcomp=where( (yrmxd ge timey(0)) and (yrmxd le timey(nyr-1)) )
;
calregts=fltarr(nfull,nreg)*!values.f_nan
calregse=fltarr(nfull,nreg)*!values.f_nan
;
; Repeat for each region separately
;
for ireg = 0 , nreg-1 do begin
  print,regname(ireg)+' '+titadd
  ;
  ; Repeat for each fixed grid calculation
  ;
  fixcoeff=fltarr(4,nfixed+1)*!values.f_nan
  for jfixed = 0 , nfixed-1 do begin
    ifixed=iord(jfixed)
    fyr=allfixedyr(ifixed)
    ;
    ; Extract series
    ;
    tstemp=reform(regtemp(*,ireg,ifixed))
    tsmxd=reform(bestmxd(*,ireg,ifixed))
    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 )
    if ifixed eq 0 then begin
      print,'Calibration period:'+string(calper,format='(2I5)')+$
        '  length:'+string(ncal,format='(I5)')
    endif
;    verkl=where( finite(tstemp+tsmxd) and $
;               (timey ge verper(0)) and (timey le verper(1)) , nver )
;  print,'Verification period:'+string(verper,format='(2I5)')+$
;    '  length:'+string(nver,format='(I5)')
    ;
    if ncal gt 30 then begin
    ;
    ; Calibrate the MXD
    ;
    mkcalibrate,tsmxd(calkl),tstemp(calkl),fitcoeff=fitcoeff,autocoeff=autocoeff
    if ifixed eq 0 then print,' Fyr SEalpha  SEbeta    rmse a_resid'
    fixcoeff(*,ifixed)=[fitcoeff([3,4,6]),autocoeff(2)]
    print,fyr,fixcoeff(*,ifixed),format='(I4,3F8.4,F8.2)'
    ;
    ; Generate calibrated record and uncertainties etc.
    ; Use the extended MXD record (after checking it matches the short one!)
    ; And only for the best fixed grid
    ;
    if fyr eq 1950 then begin
      xxx=tsfull(kcomp)
      yyy=tsmxd
      tse=total( (xxx-yyy)^2 , /nan )
;      print,tse
      if (tse gt 0.4) or (finite(tse) eq 0) then message,'Series do not match'
      calmxd=tsfull*fitcoeff(2)+fitcoeff(1)
      calregts(*,ireg)=calmxd(*)
      ;
      ; To compute errors, 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
      ;
      ; Generate time dependent standard error
      ;
      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]
        prets=callow
      endif else begin
        prets=calmxd
      endelse
      fixcoeff(*,nfixed)=fixcoeff(*,nfixed-1)  ; move 1988 values to end
      fixcoeff(*,nfixed-1)=fixcoeff(*,nfixed-2)  ; repeat 1950 values for 1983
      for jjj = nfixed-3 , 0 , -1 do begin     ; fill in missing values
        if finite(fixcoeff(0,jjj)) eq 0 then begin
          fixcoeff(*,jjj)=fixcoeff(*,jjj+1)
        endif
      endfor
      if finite(fixcoeff(0,nfixed)) eq 0 then begin     ; and end value
        fixcoeff(0,nfixed)=fixcoeff(0,nfixed-1)
      endif
      for iyr = 0 , nfull-1 do begin
        onefc3=interpol(reform(fixcoeff(0,*)),useyr,yrmxd(iyr))
        onefc4=interpol(reform(fixcoeff(1,*)),useyr,yrmxd(iyr))
        onefc6=interpol(reform(fixcoeff(2,*)),useyr,yrmxd(iyr))
        oneac2=interpol(reform(fixcoeff(3,*)),useyr,yrmxd(iyr))
        onefc3=onefc3(0)
        onefc4=onefc4(0)
        onefc6=onefc6(0)
        oneac2=oneac2(0)
      ;
      ; Now compute the one standard error for each reconstructed value
      ; (optionally follow procedure for smoothed data), using error
      ; parameters interpolated from various fixed grid series.
      ;
      ; SE due to uncertainty in (a) slope and (b) intercept coefficients,
      ; and (c) due to unexplained variance
      ;
      ; Scale up the standard error on the regression coefficient if there is
      ; autocorrelation in the residuals
      ;
      ; PREVIOUSLY I USED:
      serialfac=(1.+abs(oneac2))/(1.-abs(oneac2))
      ; BUT NOW I SQUARE ROOT IT, BECAUSE IT IS MULTIPLYING SE NOT SE^2
      serialfac=sqrt(serialfac)
      ;
      sebeta=onefc4*serialfac
      ;
      if dosmooth ne 0 then begin
        se_a=lowpred(iyr)*sebeta
        varfac=serialfac*sqrt(total(allwt^2)) < 1. ; don't let errors get bigger
        se_c=onefc6*varfac
      endif else begin
        se_a=predanom(iyr)*sebeta
        se_c=onefc6
      endelse
      se_b=onefc3
      serr=sqrt( se_a^2 + se_b^2 + se_c^2 )
      calregse(iyr,ireg)=serr
      ;
    endfor
      serr=calregse(*,ireg)
      ;
      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=titadd+' '+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
        plot,timey,tstemp,thick=3,$
          /xstyle,xrange=[1850,2000],xtitle='Year',$
          ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
          title=titadd+' '+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,noclip=0
        ;
        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,noclip=0
        oplot,timey,tstemp,thick=3
        oplot,yrmxd,calmxd,color=11,thick=3
        oplot,!x.crange,[0.,0.],linestyle=1
      endif
      ;
    endif
    endif
  endfor
endfor
;
; Save standard errors (note they apply only to the time scale chosen!)
;
if dosave then begin
  timey=yrmxd
  if dosmooth eq 1 then smst='sm'+string(thalf,format='(I2.2)') else smst=''
  save,filename='datastore/regtemp'+smst+'_'+fnvar+'_errors.idlsave',$
    calregse,thalf,dosmooth,nreg,regname,timey
endif
;
end