;
; Calibrates and plots regional MXD against regional temperature, testing
; for timescale-dependent calibration coefficients!
;
; Define periods for calibration
;
calper=[1881,1960]        ; calibration period
thalf=10.
bandlim=[0,999,0,5,15,30,999,30,80]
nband=n_elements(bandlim)-1
rband=fltarr(nband)
bband=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
multi_plot,nrow=4,ncol=3,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=9
endelse
def_1color,10,color='blue'
def_1color,11,color='red'
def_1color,12,color='mlgrey'
;
range0=[-2.5,-2.0,-3.5,-2.5,-1.5,-0.5,-2.5,-2.0,-3.0,-1.5,-1.5]
range1=[ 1.5, 2.0, 2.5, 1.5, 1.0, 1.0, 1.5, 1.5, 1.5, 1.0, 1.0]
;
calregts=fltarr(nfull,nreg)*!values.f_nan
;
; Process each region separately
;
listreg=[0,1,2,3,4,6,7,8,9,10]
nlist=n_elements(listreg)
for jreg = 0 , nlist-1 do begin
  ireg=listreg(jreg)
  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 )
  print,'Calibration period:'+string(calper,format='(2I5)')+$
    '  length:'+string(ncal,format='(I5)')
  ;
  ; Calibrate the MXD using a single band
  ;
  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)'
  ;
  ; 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 )
  if tse gt 0.03 then message,'Series do not match'
  calmxd=tsfull*fitcoeff(2)+fitcoeff(1)
  calregts(*,ireg)=calmxd(*)
  ;
  ; Calibrate the MXD, separately for each frequency band
  ;
  for iband = 0 , nband-1 do begin
    t1=reform(tstemp)
    m1=reform(tsmxd)
    mint=bandlim(iband)
    maxt=bandlim(iband+1)
    if mint ge maxt then begin
      rband(iband)=!values.f_nan
      bband(iband)=!values.f_nan
      bseband(iband)=!values.f_nan
    endif else begin
;      if mint ne 0 then begin
;        filter_cru,/nan,mint,tsin=t1,tslow=t2
;        t1=t2
;        filter_cru,/nan,mint,tsin=m1,tslow=m2
;        m1=m2
;      endif
;      if maxt ne 999 then begin
;        filter_cru,/nan,maxt,tsin=t1,tshigh=t2,tslow=t3
;        t1=t2
;        filter_cru,/nan,maxt,tsin=m1,tshigh=m2,tslow=m3
;        m1=m2
;      endif
      if mint ne 0 then begin
        filter_cru,/nan,mint,tsin=t1,tslow=t2
        filter_cru,/nan,mint,tsin=m1,tslow=m2
      endif else begin
        t2=t1
        m2=m1
      endelse
      if maxt ne 999 then begin
        filter_cru,/nan,maxt,tsin=t1,tslow=t3
        filter_cru,/nan,maxt,tsin=m1,tslow=m3
      endif else begin
        t3=t1 & t3(*)=0.
        m3=m1 & m3(*)=0.
      endelse
      t1=t2-t3
      m1=m2-m3
      if iband eq nband-1 then begin
        pause
        mknormal,t3,timey,refperiod=calper
        mknormal,m3,timey,refperiod=calper
        plot,timey,t3,/xstyle,xrange=calper,$
          /ystyle,yrange=[-2,2],ytitle='Normalised T or MXD',$
          title=regname(ireg)
        oplot,timey,m3,thick=3
        oplot,!x.crange,[0,0],linestyle=1
      endif
      mkcalibrate,m1(calkl),t1(calkl),fitcoeff=fitcoeff,autocoeff=autocoeff
      rband(iband)=fitcoeff(0)
      bband(iband)=fitcoeff(2)
      lag1r=autocoeff(2)
      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!
      ;
;      serfac=(1.+abs(lag1r))/(1.-abs(lag1r))     ; standard AR1 formula
      ;
;      lagnr=lagnr > 0.                           ; my old alternative
;      serfac=1.+2.*total(lagnr)
      ;
      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)
      bseband(iband)=sqrt( (fitcoeff(4)^2)*serfac )
      print,'PERIODICITY BAND=',mint,maxt
      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)'
      print,'SERIAL CORRELATION FACTOR=',serfac
    endelse
  endfor
  ;
  pause
  xx=findgen(nband)+1
  xtickname=string(bandlim(0:nband-1),format='(I2)')+'-'+$
    string(bandlim(1:nband),format='(I3)')
  ml=where(finite(rband) eq 0,nmiss)
  if nmiss gt 0 then xtickname(ml)=' '
  plot,xx,rband,psym=def_sym(10),$
    /xstyle,xrange=[0,nband+1],xticks=nband-1,$
    xtickv=xx,xtickname=xtickname,$
    /ystyle,yrange=[0,1],ytitle='Correlation',$
    title=regname(ireg)
;  oplot,!x.crange,[0,0],linestyle=1
  ;
  yr=[0,1]
  if ireg eq 2 then yr=[0,1.6]
  if ireg ge 9 then yr=[0,0.5]
  plot,xx,bband,psym=def_sym(10),$
    /xstyle,xrange=[0,nband+1],xticks=nband-1,$
    xtickv=xx,xtickname=xtickname,$
    /ystyle,yrange=yr,ytitle='Regression slope (!Uo!NC)',$
    title=regname(ireg)
  for iband = 0 , nband-1 do begin
    oplot,replicate(xx(iband),2),bband(iband)+bseband(iband)*[-2,2]
  endfor
  ;
  ; And low and high passed too!
  ;
;  filter_cru,/nan,thalf,tsin=tsmxd,tshigh=mxdhi,tslow=mxdlo
;  filter_cru,/nan,thalf,tsin=tstemp,tshigh=temphi,tslow=templo
;  mkcalibrate,mxdhi(calkl),temphi(calkl),fitcoeff=hicoeff,autocoeff=autohi
;  print,'HIGH-PASSED',thalf
;  print,'     r   alpha    beta SEalpha  SEbeta    sig    rmse'
;  print,hicoeff,format='(F6.2,4F8.4,F7.2,F8.4)'
;  print,' a_mxd a_tem a_resid   ess'
;  print,autohi,format='(3F6.2,F8.1)'
;  mkcalibrate,mxdlo(calkl),templo(calkl),fitcoeff=locoeff,autocoeff=autolo,$
;    nlag=5
;  print,'LOW-PASSED',thalf
;  print,'     r   alpha    beta SEalpha  SEbeta    sig    rmse'
;  print,locoeff,format='(F6.2,4F8.4,F7.2,F8.4)'
;  print,' a_mxd a_tem a_resid   ess'
;  print,autolo,format='(3F6.2,F8.1)'
;  filter_cru,/nan,thalf,tsin=tsfull,tshigh=fullhi,tslow=fulllo
;  sepmxd=(fullhi*hicoeff(2)+hicoeff(1))+(fulllo*locoeff(2)+locoeff(1))
  ;
  ; For overlap period with non-missing data, compute correlation between
  ; the temperature series and the 2-band reconstruction
  ;
;  kkk=calkl+(timey(0)-yrmxd(0))
;  print,'r(Two-band calibrated series,Temp)='+$
;    string(correlate(tstemp(calkl),sepmxd(kkk)),format='(F6.2)')
  ;
;  if regname(ireg) eq 'ALL' then !p.multi=[0,1,2,0,0]
;  pause
;  plot,timey,tstemp,$
;    /xstyle,xrange=[1400,1994],xtitle='Year',$
;    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
;    /ystyle,yrange=[range0(ireg),range1(ireg)],$
;    title=regname(ireg)
;  oplot,yrmxd,calmxd,color=10       ;,thick=3
;  oplot,yrmxd,sepmxd,color=11       ;,thick=3
;  oplot,!x.crange,[0.,0.],linestyle=1
;  filter_cru,/nan,thalf,tsin=tstemp,tslow=ylow
;  oplot,timey,ylow,thick=3
;  filter_cru,/nan,thalf,tsin=calmxd,tslow=ylow
;  oplot,yrmxd,ylow,thick=3,color=10
;  filter_cru,/nan,thalf,tsin=sepmxd,tslow=ylow
;  oplot,yrmxd,ylow,thick=3,color=11
  ;
;  if (regname(ireg) eq 'ALL') or (regname(ireg) eq 'NH') then begin
;  plot,timey,tstemp,thick=2,$
;    /xstyle,xrange=[1850,2000],xtitle='Year',$
;    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
;    title=regname(ireg)
;  oplot,yrmxd,calmxd,color=10,thick=3
;  oplot,yrmxd,sepmxd,color=11,thick=3
;  oplot,!x.crange,[0.,0.],linestyle=1
;  endif
  ;
endfor
;
; Now output the calibrated series
;
tempregts=regtemp
tempnyr=nyr
temptimey=timey
nyr=nfull
timey=yrmxd
save,filename='regtemp_tsd_calibrated.idlsave',$
  nyr,nreg,calregts,regname,timey,tempregts,tempnyr,temptimey
;
end