;
; Calibrates and plots regional MXD against regional temperature, testing
; for timescale-dependent calibration coefficients!
;
docol=0         ; 0=black&white 1=color
;
; Define periods for calibration
;
calper=[1881,1960]        ; calibration period
bandlim=[0,15,999]
nband=n_elements(bandlim)
rband=fltarr(nband)
aband=fltarr(nband)
bband=fltarr(nband)
aseband=fltarr(nband)
bseband=fltarr(nband)
;
; Get the hemispheric series
;
restore,filename='../compbest_mxd_fixed1950.idlsave'
nhtemp=reform(comptemp(*,3))
alltemp=reform(comptemp(*,2))
nhmxd=reform(compmxd(*,2,0))
; Get rid of pre-1871 termperatures
knh=where(timey lt 1871)
nhtemp(knh)=!values.f_nan
alltemp(knh)=!values.f_nan
;
; Get the extended hemispheric series
;
restore,filename='../compts_mxd_fixed1950.idlsave'
fullnhmxd=reform(compmxd(*,2))
;
; Get the extended MXD series
;
restore,filename='../regts_mxd_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+2)*!values.f_nan
newfm(*,0:nreg-1)=fullmxd(*,*)
newfm(*,nreg)=fullnhmxd(*)
newfm(*,nreg+1)=fullnhmxd(*)
fullmxd=newfm
;
; Get regional series (temp, and truncated MXD)
;
restore,filename='../regbest_mxd_fixed1950.idlsave'
; Gets: nreg,regname,regtemp,rawtemp,timey,bestmxd etc.
;
; Identify overlap between full and truncated MXD series
;
nyr=n_elements(timey)
kcomp=where( (yrmxd ge timey(0)) and (yrmxd le timey(nyr-1)) )
;
; Add in the hemispheric series
;
regname=[regname,'ALL','NH']
newrt=fltarr(nyr,nreg+2)*!values.f_nan
newrt(*,0:nreg-1)=regtemp(*,*)
newrt(*,nreg)=alltemp(*)
newrt(*,nreg+1)=nhtemp(*)
regtemp=newrt
newrm=fltarr(nyr,nreg+2)*!values.f_nan
newrm(*,0:nreg-1)=bestmxd(*,*,1)
newrm(*,nreg)=nhmxd(*)
newrm(*,nreg+1)=nhmxd(*)
bestmxd=newrm
nreg=nreg+2
;
; Prepare for plotting
;
loadct,39
if docol eq 1 then multi_plot,nrow=5,ncol=2,layout='large' $
              else multi_plot,nrow=4,ncol=1,layout='large'
if !d.name eq 'X' then begin
  window,ysize=850,xsize=700
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=10+6*(1-docol)
endelse
if docol eq 0 then cname=replicate('black',3) else cname=['blue','red','mlgrey']
for i = 0 , n_elements(cname)-1 do def_1color,i+10,color=cname(i)
;
range0=[-2.5,-1.8,-3.5,-2.5,-1.5,-0.5,-2.5,-2.0,-3.0,-1.5,-1.5]
range1=[ 1.5, 1.8, 2.5, 1.5, 1.0, 1.0, 1.5, 1.5, 1.5, 1.0, 1.0]
;
calregts=fltarr(nfull,nreg,2)*!values.f_nan
calregse=fltarr(nfull,nreg,2)*!values.f_nan
;
; Process each region separately
;
ptit='('+['a','b','c','d','e','f','g','h','i','j','k','l','m','n','o']+') '
if docol eq 1 then listreg=[0,1,2,3,4,6,7,8,9,10] $
              else listreg=[1,0]
nlist=n_elements(listreg)
for jreg = 0 , nlist-1 do begin
  ireg=listreg(jreg)
  print
  print,regname(ireg)
  ;
  ; Extract series
  ;
  tstemp=reform(regtemp(*,ireg))
  tsmxd=reform(bestmxd(*,ireg))
  tsfull=reform(fullmxd(*,ireg))
  calregts(*,ireg,*)=0.
  calregse(*,ireg,*)=0.
  ;
  ; Identify calibration and verification subsets
  ;
  calkl=where( finite(tstemp+tsmxd) and $
               (timey ge calper(0)) and (timey le calper(1)) , ncal )
  print,'Calibration period:'+string(calper,format='(2I5)')+$
    '  length:'+string(ncal,format='(I5)')
  ;
  ; Use the extended MXD record (after checking it matches the short one!)
  ; for later generation of the calibrated record.
  ;
  xxx=tsfull(kcomp)
  yyy=tsmxd
  tse=total( (xxx-yyy)^2 , /nan )
  if tse gt 0.03 then message,'Series do not match'
  ;
  ; Calibrate the MXD, separately for each frequency band
  ;
  for iband = 0 , nband-1 do begin
    t1=reform(tstemp)
    m1=reform(tsmxd)
    f1=reform(tsfull)
    ; Filter the series if necessary
    if iband eq nband-1 then begin
      mint=0
      maxt=999
    endif else begin
      mint=bandlim(iband)
      maxt=bandlim(iband+1)
      if mint ne 0 then begin
        filter_cru,/nan,mint,tsin=t1,tslow=t2
        filter_cru,/nan,mint,tsin=m1,tslow=m2
        filter_cru,/nan,mint,tsin=f1,tslow=f2
      endif else begin
        t2=t1
        m2=m1
        f2=f1
      endelse
      if maxt ne 999 then begin
        filter_cru,/nan,maxt,tsin=t1,tslow=t3
        filter_cru,/nan,maxt,tsin=m1,tslow=m3
        filter_cru,/nan,maxt,tsin=f1,tslow=f3
      endif else begin
        t3=t1 & t3(*)=0.
        m3=m1 & m3(*)=0.
        f3=f1 & f3(*)=0.
      endelse
      t1=t2-t3
      m1=m2-m3
      f1=f2-f3
    endelse
    ; Calibrate the filtered or unfiltered series
    mkcalibrate,m1(calkl),t1(calkl),fitcoeff=fitcoeff,autocoeff=autocoeff
    rband(iband)=fitcoeff(0)
    aband(iband)=fitcoeff(1)
    bband(iband)=fitcoeff(2)
    aseband(iband)=fitcoeff(3)
    ; For high-pass filtered series, we set the acoeff error to zero
;    if maxt ne 999 then aseband(iband)=0.
    m4=m1(calkl)
    t4=t1(calkl)
    pred4=fitcoeff(1)+fitcoeff(2)*m4
    resid4=t4-pred4
    ;
    ; We need to inflate the uncertainty in the regression slope due to
    ; autocorrelation in the residuals.  We cannot use the standard
    ; formula since the smoothed cases are not AR1 models.  We have to
    ; sum over the leading N lag autocorrelations.  But what is N?  Given
    ; that the sample size is 80 years, then a lag1 autocorrelation of
    ; +-0.2 is significant at the 95% level for a one-tailed test (since
    ; we know that it will be negative for high-pass and positive for band
    ; or low pass filters.  So, we keep going until lag(i)r is >-0.2 and
    ; <+0.2, and then stop!
    ;
;    nlag=30
;    lagnr=a_correlate(resid4,indgen(nlag)+1)
;    if abs(lagnr(0)-autocoeff(2)) gt 0.001 then message,'Should be the same!'
;    lagabs=abs(lagnr)
;    notsiglist=where(lagabs lt 0.2,nnot)
;    if nnot gt 0 then lagabs(notsiglist(0):nlag-1)=0.
;    serfac=1.+2.*total(lagabs)
    ;
    ; Now use an alternative method of computing effective independent
    ; sample size
    ;
    if iband eq nband-1 then begin
      neffdata=compute_neff(resid4)
      serfac=neffdata(1)
    endif else begin
      repeat begin
        xresid=randomn(seed,ncal)
        ar1=a_correlate(xresid,1)
        ar1=ar1(0)
      endrep until (abs(ar1) lt 0.05)
      neffdata=compute_neff(xresid)
      neall=neffdata(0)
      if (mint eq 0) and (maxt eq 999) then ne0=neall
      if (mint eq 0) and (maxt ne 999) then begin
        neffdata=compute_neff(xresid,tfilter=[maxt,999])
        ne0=neall-neffdata(0)
      endif
      if (mint gt 0) and (maxt eq 999) then begin
        neffdata=compute_neff(xresid,tfilter=[mint,999])
        ne0=neffdata(0)
      endif
      if (mint gt 0) and (maxt ne 999) then begin
        neffdata1=compute_neff(xresid,tfilter=[mint,999])
        neffdata2=compute_neff(xresid,tfilter=[maxt,999])
        ne0=neffdata1(0)-neffdata2(0)
      endif
      serfac=float(ncal)/float(ne0)
    endelse
    ;
    bseband(iband)=sqrt( (fitcoeff(4)^2)*serfac )
    print,'PERIODICITY BAND=',mint,maxt
    print,'     r   alpha    beta SEalpha  SEbeta    sig    rmse'
    print,rband(iband),aband(iband),bband(iband),aseband(iband),$
      bseband(iband),format='(F6.2,4F8.4,F7.2,F8.4)'
;    print,' a_mxd a_tem a_resid   ess'
;    print,autocoeff,format='(3F6.2,F8.1)'
    print,'SERIAL CORRELATION FACTOR=',serfac
    ;
    ; Now generate the calibrated time series for this frequency band
    ;
    calmxd=f1*bband(iband)+aband(iband)
    if iband eq nband-1 then begin     ; single band
      calregts(*,ireg,0)=calmxd(*)
    endif else begin
      calregts(*,ireg,1)=calregts(*,ireg,1)+calmxd(*)
    endelse
    ;
    ; Finally compute the one standard error for each reconstructed value
    ; (not suitable for computing errors on smoothed data!).  To do so, we
    ; need the predictor series in terms of anomalies wrt to its mean
    ; over the calibration period.
    ;
    f1anom=f1
    mkanomaly,f1anom,yrmxd,refperiod=calper
    ;
    ; SE^2 = SE_A^2 + (SE_B*F1ANOM)^2 + SE_RESID^2      ; single band
    ; SE^2 = SE_A^2 + sum(SE_B*FILANOM)^2 + SE_RESID^2  ; multi band
    ; We can actually sum the SE_A^2 too, since they're set to zero for all
    ; except the low-pass case.
    ; But we cannot combine the separate RMSE of the residuals from the
    ; separate frequency bands, since the separate bands have a low, but not
    ; quite zero, cross-correlation.  We can, however, combine the residuals
    ; and then compute an RMSE of them.
    ;
    if iband eq nband-1 then begin
      sea2=aseband(iband)^2
      seb2=(f1anom*bseband(iband))^2
      dummy=moment(resid4,sdev=onesd)
      print,'SINGLE-BAND reconstruction residual rmse',onesd
      seresid2=onesd^2
      calregse(*,ireg,0)=sqrt( sea2 + seb2 + seresid2 )
      r1single=autocoeff(2)
    endif else begin
      ; Do residuals later, and square root later
      sea2=aseband(iband)^2
      seb2=(f1anom*bseband(iband))^2
      calregse(*,ireg,1)=calregse(*,ireg,1)+sea2+seb2
    endelse
    ;
  endfor
  ;
  ; For the multi-band reconstruction, add on the standard error^2 due to
  ; the overall residual variability, and then square root to get SE
  ;
  kkk=calkl+(timey(0)-yrmxd(0))
  resid4=tstemp(calkl)-calregts(kkk,ireg,1)
  dummy=moment(resid4,sdev=onesd)
  print,'MULTI-BAND reconstruction residual rmse',onesd
  calregse(*,ireg,1)=sqrt(calregse(*,ireg,1)+onesd^2)
  r1multi=a_correlate(resid4,[1])
  print,'Lag-1 autocorrelation in residuals: single,multi=',$
    r1single,r1multi,format='(A,2F6.2)'
  ;
  ; For overlap period with non-missing data, compute correlation between
  ; the temperature series and the multi-band reconstruction
  ;
  print,'r(multi-band calibrated series,temp)='+$
    string(correlate(tstemp(calkl),calregts(kkk,ireg,1)),format='(F6.2)')
  ;
  pause
  plot,timey,tstemp,linestyle=1-docol,thick=1+2*(1-docol),$
    /xstyle,xrange=[1801,1960],xtitle='Year',$
    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
    /ystyle,yrange=[range0(ireg),range1(ireg)],$
    title=ptit(jreg*2)+regname(ireg)
  calmxd=reform(calregts(*,ireg,0))
  sepmxd=reform(calregts(*,ireg,1))
  oplot,yrmxd,calmxd,color=10
  oplot,yrmxd,sepmxd,color=11,thick=(1-docol)*2+1
  oplot,!x.crange,[0.,0.],linestyle=1
  ;
  filter_cru,/nan,10.,tsin=tstemp,tslow=ylow
  oplot,timey,ylow,thick=2+docol
  if docol eq 0 then oplot,timey,ylow,psym=def_sym(10),symsize=0.4
  filter_cru,/nan,10.,tsin=calmxd,tslow=ylow
  oplot,yrmxd,ylow,thick=2+docol,color=10
  filter_cru,/nan,10.,tsin=sepmxd,tslow=ylow
  oplot,yrmxd,ylow,thick=3+2*(1-docol),color=11
  ;
  plot,yrmxd,calregse(*,ireg,0),$
    /xstyle,xrange=[1801,1960],xtitle='Year',$
    ytitle='Temperature reconstruction standard error  (!Uo!NC)',$
    title=ptit(jreg*2+1)+regname(ireg)
  oplot,yrmxd,calregse(*,ireg,1),thick=5
  ;
endfor
;
end