;
; Calibrates the age-banded density data against observed temperature
; for a selected region (or quasi-hemispheric mean), using various
; options for focussing on low-frequency term only, and various
; options for quantifying uncertainty.
;
doplot=2      ; 0=none, 1=time series, 2=scatter of parameters,
              ; 3=mean parameters over timscales
;
; Select which version of the age-banded records to use.
fnband='_50-300'   ; age bands used.  Choices: '_50-300' ''=20-550
fnmint='_min2'     ; minimum trees in each band.  Choices: '_min2' ''=all used
fnepsw=''          ; weighting used for all.  Choices: '_eps' ''=no weighting
;
dosave=0     ; save results or not?  0=no, 1=IDLSAVE, 2=IDLSAVE & ASCII
;
; Define periods and time scales for calibration
;
calper=[1871,1958]        ; calibration period
;
caltscale=[-10,-9,-8,-7,-6,-5,-4,-3,-2,1,2,3,4,5,6,7,8,9,10]  ; smoothing to apply prior to calibration
;caltscale=[9]
caltscale=[1]
ncalts=n_elements(caltscale)
;
errtscale=[1]
nerrts=n_elements(errtscale)
;
; Select region(s) to calibrate (0-8=regions, 9=all, 10=nh)
;
listreg=[9]
ndoreg=n_elements(listreg)
;
; Select sequence of fixed grids to calibrate, ending in
; the "best coverage" grid.  Also list the times between
; which to interpolate the calibration error statistics,
; for generating a smoothly varying error through time.
;
; All available fixed grids
allfixedyr=[1400,1500,1600,1700,1800,1950,1988]
nfixed=n_elements(allfixedyr)
;
; Which to use and in what order
iord=[0,1,2,3,4,6,5]         ; order to process (keep 1950 last!)
iord=[5]
nord=n_elements(iord)
;
; When to apply the values, and which value to use
useyr= [1400,1500,1600,1700,1800,1900,1983,1994]
usefix=[1400,1500,1600,1700,1800,1950,1950,1988]
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=4,ncol=2,layout='large',/landscape
if !d.name eq 'X' then begin
  wintim,xsize=850,ysize=750
  !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]
;
; Get the hemispheric temperature series
;
fn='datastore/compbest_mxd_fixed1950.idlsave'
print,'Observed NH and ALL temperature'
print,fn
print
restore,filename=fn
nhtemp=reform(comptemp(*,3))
alltemp=reform(comptemp(*,2))
; Get rid of pre-1871 temperatures
knh=where(timey lt 1871)
nhtemp(knh)=!values.f_nan
alltemp[knh]=!values.f_nan
;
; Get all the age-banded MXD series (regions, hemi and back to 1400)
;
fn='/cru/u2/f055/treeharry/banding/bandregnh'+$
  fnband+fnmint+fnepsw+'_fixed.idlsave'
print,'Full tree time series from full coverage data'
print,fn
print
restore,filename=fn
; 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+2)*!values.f_nan
newfm(*,0:nreg-1)=transpose(fullmxd(*,*))
newfm(*,nreg)=nhts(*)     ; predictor is same for ALL and NH, so repeat
newfm(*,nreg+1)=nhts(*)
fullmxd=newfm
;
; Get regional temperature series
;
fn='datastore/regbest_mxd_fixed1950.idlsave'
print,'Observed regional temperature'
print,fn
print
restore,filename=fn
; 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)) ,ncomp )
;
; Combine hemispheric and regional temperatures
;
newrt=fltarr(nyr,nreg+2)*!values.f_nan
newrt(*,0:nreg-1)=regtemp(*,*)
newrt(*,nreg)=alltemp(*)
newrt(*,nreg+1)=nhtemp(*)
regtemp=newrt
tempyear=timey
tempnyr=nyr
;
; Need to read in all the series generated by various fixed grids,
; then cut down to just the period that overlaps with observations
;
for ifixed = 0 , nfixed-1 do begin
  ;
  fyr=allfixedyr(ifixed)
  fnadd=string(fyr,format='(I4)')
  if fnadd eq '1950' then fnadd=''
  ;
  ; Get all the age-banded MXD series (regions, hemi and back to 1400)
  ;
  fn='/cru/u2/f055/treeharry/banding/bandregnh'+$
    fnband+fnmint+fnepsw+'_fixed'+fnadd+'.idlsave'
  if ifixed eq 0 then print,'Tree time series from fixed grid coverage data'
  print,fn
  if ifixed eq nfixed-1 then print
  restore,filename=fn
  ; Gets: timey,nyr,nhts,allts,alleps,nreg,regname
  ;
  if ifixed eq 0 then begin
    newrm=fltarr(tempnyr,nreg+2,nfixed)*!values.f_nan
  endif
  newrm(0:ncomp-1,0:nreg-1,ifixed)=transpose(allts(*,kcomp))
  newrm(0:ncomp-1,nreg,ifixed)=nhts(kcomp)
  newrm(0:ncomp-1,nreg+1,ifixed)=nhts(kcomp)
  ;
endfor
;
bestmxd=newrm
regname=[regname,'ALL','NH']
nreg=nreg+2
nyr=tempnyr
timey=tempyear
;
; Create arrays for storing calibrations and uncertainties
;
calregts=fltarr(nfull,ndoreg)*!values.f_nan
calregse=fltarr(nfull,ndoreg)*!values.f_nan
;
; Repeat for each region separately
;
for jreg = 0 , ndoreg-1 do begin
  ireg=listreg[jreg]
  print,regname[ireg]
  ;
  meanfitcoeff=fltarr(nord,7,ncalts)*!values.f_nan
  meanautocoeff=fltarr(nord,4,ncalts)*!values.f_nan
  ;
  ; Repeat for each fixed grid calculation
  ;
  for jfixed = 0 , nord-1 do begin
    ifixed=iord(jfixed)
    fyr=allfixedyr(ifixed)
    print,'FIXED for',fyr
    ;
    ; Extract series
    ;
    tstemp=reform(regtemp(*,ireg))
    tsmxd=reform(bestmxd(*,ireg,ifixed))
    tsfull=reform(fullmxd(*,ireg))
    ;
    for icalts = 0 , ncalts-1 do begin
      ;
      ; If requested, filter the series prior to calibration
      ; Two possible options are (1) running smoothing (e.g. gaussian weighting),
      ; and (2) discrete means (not running means).  The latter has the disadvantage
      ; that some data is being thrown away (by taking e.g. non-overlapping decadal
      ; mean values), but the advantage that the filtering does not impose any
      ; extra autocorrelation on the data.  The solution used here is to do (2), but
      ; repeated for all options of the starting point for computing discrete means,
      ; and then the regression parameters are averaged over all outcomes.
      ;
      if abs(caltscale[icalts]) le 1 then begin
        ;
        ; No smoothing required
        ;
        nser=1
        tempser=reform(tstemp,tempnyr,nser)
        mxdcalser=reform(tsmxd,tempnyr,nser)
        mxdallser=reform(tsfull,nfull,nser)
        timeyser=timey
        ;
      endif
      ;
      if caltscale[icalts] lt -1 then begin
        ;
        ; Simple Gaussian-weighted filtering is done
        ;
        filter_cru,/nan,-caltscale[icalts],tsin=tstemp,tslow=tempser
        filter_cru,/nan,-caltscale[icalts],tsin=tsmxd,tslow=mxdcalser
        filter_cru,/nan,-caltscale[icalts],tsin=tsfull,tslow=mxdallser
        nser=1
        tempser=reform(tempser,tempnyr,nser)
        mxdcalser=reform(mxdcalser,tempnyr,nser)
        mxdallser=reform(mxdallser,nfull,nser)
        timeyser=timey
        ;
      endif
      ;
      if caltscale[icalts] gt 1 then begin
        ;
        ; Non-overlapping running means, for all registration positions
        ;
        nser=caltscale[icalts]
        tempser=mkallmeans(tstemp,caltscale[icalts])
        mxdcalser=mkallmeans(tsmxd,caltscale[icalts])
        mxdallser=mkallmeans(tsfull,caltscale[icalts])
        ;
        dummy=reform(tempser[*,0])
        ndummy=n_elements(dummy)
        timeyser=timey[0]+findgen(ndummy)*caltscale[icalts]+0.5*(caltscale[icalts]-1)
        ;
        if (fyr eq '1950') and (doplot eq 1) then begin
          for i = 0 , caltscale[icalts]-1 do begin
            pause
            plot,timey,tstemp,/xstyle
            oplot,timeyser+i,tempser[*,i],psym=10
            oplot,timey,tsmxd*0.4,color=240
            oplot,timeyser+i,mxdcalser[*,i]*0.4,psym=10,color=240
          endfor
        endif
        ;
      endif
      ;
      ; Now calibrate each of the generated series
      ;
      allfitcoeff=reform(fltarr(7,nser),7,nser)*!values.f_nan
      allautocoeff=reform(fltarr(4,nser),4,nser)*!values.f_nan
      help,allfitcoeff,allautocoeff
      for iser = 0 , nser-1 do begin
        ;
        ; Identify calibration and verification subsets
        ;
        x=reform(tempser[*,iser])
        y=reform(mxdcalser[*,iser])
        ty=timeyser+iser
        calkl=where(finite(x+y) and (ty ge calper(0)) and (ty le calper(1)),ncal)
        if jfixed eq 0 then begin
          print,'Calibration period:'+string(calper,format='(2I5)')+$
            '  length:'+string(ncal,format='(I5)')
        endif
        ;
        ; Calibrate the MXD
        ;
        mkcalibrate,y(calkl),x(calkl),fitcoeff=fitcoeff,autocoeff=autocoeff
        ;
        allfitcoeff[*,iser]=fitcoeff
        allautocoeff[*,iser]=autocoeff
        ;
        if iser+jfixed eq 0 then print,' Fyr SEalpha  SEbeta    rmse a_resid       r   alpha    beta'
        fixcoeff=[fitcoeff([3,4,6]),autocoeff(2),fitcoeff([0,1,2])]
        print,fyr,fixcoeff,format='(I4,3F8.4,2F8.2,2F8.4)'
        ;
        ; Generate the calibrated record and the residuals
        ;
        xpred=fitcoeff[1]+fitcoeff[2]*y
        xresid=xpred-x
        if doplot eq 1 then begin
          pause
          plot,ty,x,/xstyle,psym=10
          oplot,ty,xpred,color=240,psym=10
          pause
          plot,ty,xresid,/xstyle,psym=10
        endif
        ;
        ; Analyse the autoregressive and autocorrelation structure of residuals
        ;
        maxlag=10+20*(caltscale[icalts] lt 3)
        maxlag=maxlag < (ncal-1)
        xlag=findgen(maxlag+1)
        maxord=7
        tjo_tsauto,xresid[calkl],maxlag=maxlag,maxorder=maxord,acf=acf,$
          akaike=akaike,printakaike=(nser eq 1),wantorder=[1],arwant=ar1,$
          bestorder=bestorder,arbest=arbest
        print,'BEST ORDER =',bestorder
        if iser eq 0 then allakaike=fltarr(nser,2,maxord+1)*!values.f_nan
        allakaike[iser,*,*]=akaike[*,*]
        ;
        ; Plot autocorrelation function
        ;
        if doplot eq 2 then begin
          pause
          noer=((iser gt 0) and (doplot eq 0))
          if noer gt 0 then !p.multi[0]=!p.multi[0]+1
          plot,xlag,acf,psym=def_sym(10),yrange=[-0.4,1],/ystyle,$
            /xstyle,noclip=1,color=20+iser*20
          oplot,!x.crange,[0,0],linestyle=1
          ;
          ; Compute and plot the AR1 model
          ;
          acf_ar1=tjo_ar2acf(ar1[0,2],maxlag=maxlag)
          oplot,xlag,acf_ar1,color=20+iser*20
          ;
          ; Compute and plot the best order AR model
          ;
          if bestorder gt 0 then begin
            acf_arbest=tjo_ar2acf(arbest[2:*],maxlag=maxlag)
            oplot,xlag,acf_arbest,color=20+iser*20,thick=3
          endif
          ;
        endif
        ;
      endfor
      ;
      meanakaike=total(allakaike,1)/float(nser)
      print
      print,'Mean Akaike results'
      print
      print,'Autoregressive model order selection'
      print
      print,'  Lag: m    r(m)  se2(m)  AIC(m)'
      for ilag = 0 , maxord do begin
        print,ilag,acf[ilag],meanakaike[*,ilag],format='(I8,F8.3,F8.4,F8.3)'
      endfor
      print
      ;
      meanfitcoeff[jfixed,*,icalts]=total(allfitcoeff,2)/float(nser)
      meanautocoeff[jfixed,*,icalts]=total(allautocoeff,2)/float(nser)
      if doplot eq 2 then begin
        pause
        plot,allfitcoeff[0,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='r',/ynozero
        oplot,!x.crange,replicate(meanfitcoeff[jfixed,0,icalts],2)
        pause
        plot,allfitcoeff[1,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='intercept',/ynozero
        oplot,!x.crange,replicate(meanfitcoeff[jfixed,1,icalts],2)
        pause
        plot,allfitcoeff[2,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='slope',/ynozero
        oplot,!x.crange,replicate(meanfitcoeff[jfixed,2,icalts],2)
        pause
        plot,allfitcoeff[3,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='SE intercept',/ynozero
        oplot,!x.crange,replicate(meanfitcoeff[jfixed,3,icalts],2)
        pause
        plot,allfitcoeff[4,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='SE slope',/ynozero
        oplot,!x.crange,replicate(meanfitcoeff[jfixed,4,icalts],2)
        pause
        plot,allfitcoeff[6,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='rmse',/ynozero
        oplot,!x.crange,replicate(meanfitcoeff[jfixed,6,icalts],2)
        pause
        plot,allautocoeff[2,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],$
          title='residuals autocorrelation',/ynozero
        oplot,!x.crange,replicate(meanautocoeff[jfixed,2,icalts],2)
      endif
      ;
    endfor
    ;
  endfor
    ;
    if doplot eq 3 then begin
      !x.range=[min(caltscale)-1,max(caltscale)+1]
      !x.style=1
      pause
      plot,caltscale,meanfitcoeff[0,0,*],psym=def_sym(10),title='r',/ystyle,yrange=[0.3,0.9]
      for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,0,*],psym=def_sym(10),color=50+30*jfixed
      pause
      plot,caltscale,meanfitcoeff[0,1,*],psym=def_sym(10),title='intercept',/ystyle,yrange=[-0.4,-0.14]
      for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,1,*],psym=def_sym(10),color=50+30*jfixed
      pause
      plot,caltscale,meanfitcoeff[0,2,*],psym=def_sym(10),title='slope',/ystyle,yrange=[0.1,0.5]
      for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,2,*],psym=def_sym(10),color=50+30*jfixed
      oplot,!x.crange,[0.42,0.42],linestyle=1
      pause
      plot,caltscale,meanfitcoeff[0,3,*],psym=def_sym(10),title='SE intercept',/ystyle,yrange=[0.02,0.09]
      for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,3,*],psym=def_sym(10),color=50+30*jfixed
      pause
      plot,caltscale,meanfitcoeff[0,4,*],psym=def_sym(10),title='SE slope',/ystyle,yrange=[0.02,0.16]
      for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,4,*],psym=def_sym(10),color=50+30*jfixed
      pause
      ar1=reform(meanautocoeff[0,2,*])
      plot,caltscale,meanfitcoeff[0,4,*]*sqrt((1.+abs(ar1))/(1.-abs(ar1))),$
        psym=def_sym(10),title='SE slope (badly adjusted for AR1)',/ystyle,yrange=[0.04,0.3]
      for jfixed = 0 , nord-1 do begin
        ar1=reform(meanautocoeff[jfixed,2,*])
        plots,caltscale,meanfitcoeff[jfixed,4,*]*sqrt((1.+abs(ar1))/(1.-abs(ar1))),$
          psym=def_sym(10),color=50+30*jfixed
      endfor
      pause
      plot,caltscale,meanfitcoeff[0,6,*],psym=def_sym(10),title='rmse',/ystyle,yrange=[0.11,0.32]
      for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,6,*],psym=def_sym(10),color=50+30*jfixed
      pause
      plot,caltscale,meanautocoeff[0,2,*],psym=def_sym(10),$
        title='residuals autocorrelation',/ystyle,yrange=[-0.2,1]
      for jfixed = 0 , nord-1 do plots,caltscale,meanautocoeff[jfixed,2,*],psym=def_sym(10),color=50+30*jfixed
      for jfixed = 0 , nord-1 do xyouts,-10,-0.1+0.15*jfixed,string(allfixedyr[iord[jfixed]],format='(I4)'),$
        color=50+30*jfixed
      oplot,!x.crange,[0.1,0.1],linestyle=1
      !x.range=[0,0]
      !x.style=0
    endif
  ;
endfor
  ;
if 1 eq 0 then begin
    ;
    if ncal gt 30 then begin
    ;
    ; 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.03) 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
      ;
    serialfac=(1.+abs(oneac2))/(1.-abs(oneac2))
    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=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=regname(ireg)
        oplot,yrmxd,calmxd,color=11,thick=3
        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
      ;
    endif
    endif
endif
;
if dosave gt 0 then begin
  ;
  ; Save standard errors (note they apply only to the time scale chosen!)
  ;
  timey=yrmxd
  save,filename='bandtemp_mxd_errors.idlsave',$
    calregse,thalf,dosmooth,nreg,regname,timey
  ;
endif
;
end