;
; Quantifies the time-dependence of the correlation and regression slope
; coefficients for regional MXD data.
; Low pass filtering at century and longer time scales never gets rid of the
; trend - so eventually I start to scale down the 120-yr low pass time series
; to mimic the effect of removing/adding longer time scales!
; EFFECTIVE SAMPLE SIZE is now computed by generating long time series and
; filtering them, to compute the lag-autocorrelations etc.
;
; Choose set of filters to use
;
fuse=1      ; 0 = lots of filters,  1 = fewer filters
;
; Define periods for calibration, filters etc.
;
calper=[1881,1960]        ; calibration period
trv=2   ; 0=MXD 1=TRW 2=MXD with TRW overlaid 3=ABD_MXD (gets obs from MXD)
case trv of
  0: begin
    fnvar='mxd'
    end
  1: begin
    fnvar='trw'
    end
  2: begin
    fnvar='mxd'
    end
  3: begin
    fnvar='mxd'
    end
endcase
titadd=strupcase(fnvar)
;
; We want lots of overlapping band-passed filters, whose width increases
; linearly with time scale (from a width of 10 for the 0-10 year timescale,
; to a width of 50 for the 50-100 year timescale, noting that only an 80-yr
; run of data is analysed).
;
if fuse eq 0 then begin
  allmin=[findgen(26),findgen(8)*2+26,replicate(40,12)]
  allmax=[findgen(26)*2+10,findgen(8)*5+65,findgen(12)*10+110]
  nband=n_elements(allmin)
  toohigh=where(allmax gt 120,n2hi)
  if n2hi gt 0 then allmax(toohigh)=(-findgen(n2hi)*0.1-0.1) > (-1.)
endif else begin
  allmin=[ 0.,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, $
           15, 17, 19, 22, 25, 28, 31, 34, 37,replicate(40,11)]
  allmax=[10., 12, 14, 17, 20, 23, 26, 29, 32, 36, 40, 44, 48, 52, 56, $
           60, 65, 70, 75, 80, 85, 90,100,110,120,replicate(130,10)]
  nband=n_elements(allmin)
  toohigh=where(allmax gt 120,n2hi)
  if n2hi gt 0 then allmax(toohigh)=(-findgen(n2hi)*0.1-0.1) > (-1.)
endelse
;
print,allmin
print
print,allmax
print
;
rband=fltarr(nband)
bband=fltarr(nband)
bseband=fltarr(nband)
;
; If we are to overlay the TRW results on the MXD figures, then restore
; the previously saved TRW values.  Similarly for the ABD values.
;
if trv eq 2 then begin
  restore,filename='quantify_tsdcal.idlsave'
;    allsfac,allbeta,allse,allr
  trwsfac=allsfac
  trwbeta=allbeta
  trwse=allse
  trwr=allr
  ;
  restore,filename='quantify_tsdcal_abd.idlsave'
;    allsfac,allbeta,allse,allr
  abdsfac=allsfac
  abdbeta=allbeta
  abdse=allse
  abdr=allr
endif
;
if trv eq 3 then begin
  ;
  ; 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+2)*!values.f_nan
  newfm(*,0:nreg-1)=transpose(fullmxd(*,*))
  newfm(*,nreg)=nhts(*)
  newfm(*,nreg+1)=nhts(*)
  fullmxd=newfm
endif
;
; Get the hemispheric series
;
restore,filename='../compbest_'+fnvar+'_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 regional series (temp, and truncated MXD)
;
restore,filename='../regbest_'+fnvar+'_fixed1950.idlsave'
; Gets: nreg,regname,regtemp,rawtemp,timey,bestmxd etc.
nyr=n_elements(timey)
;
; 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
;
if trv eq 3 then begin
  ;
  ; Replace MXD with age-banded MXD (cutting out an appropriate period!)
  ;
  kmxd=where((yrmxd ge timey(0)) and (yrmxd le timey(nyr-1)),nmxd)
  bestmxd=fullmxd(kmxd,*)
  print,yrmxd(kmxd)
  ;
endif
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=4,ncol=2,layout='large'
if !d.name eq 'X' then begin
  wintim,ysize=850,xsize=700
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=14
endelse
def_1color,10,color='blue'
def_1color,11,color='red'
def_1color,12,color='mlgrey'
;
; Process each region separately
;
listreg=[0,1,2,3,4,6,7,8,9,10]
nlist=n_elements(listreg)
allr=fltarr(nband,nreg)
allr(*,*)=!values.f_nan
allbeta=fltarr(nband,nreg)
allbeta(*,*)=!values.f_nan
allse=fltarr(nband,nreg)
allse(*,*)=!values.f_nan
allsfac=fltarr(nreg)
allsfac(*)=!values.f_nan
for jreg = 0 , nlist-1 do begin
  ireg=listreg(jreg)
  print,regname(ireg)
  ;
  ; Extract series
  ;
  tstemp=reform(regtemp(*,ireg))
  tsmxd=reform(bestmxd(*,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)'
  sfac=fitcoeff(2)
  ;
  ; Calibrate the MXD, separately for each frequency band
  ;
  for iband = 0 , nband-1 do begin
    t1=reform(tstemp)
    m1=reform(tsmxd)
    mint=allmin(iband)
    maxt=allmax(iband)
    if maxt lt 0 then begin
      maxfac=maxt
      maxt=120
    endif else begin
      maxfac=1.
    endelse
    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
      if maxfac ne 1. then begin
        t3=t3*(1.+maxfac)
        m3=m3*(1.+maxfac)
      endif
    endif else begin
      t3=t1 & t3(*)=0.
      m3=m1 & m3(*)=0.
    endelse
    t1=t2-t3
    m1=m2-m3
    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 these are not well
    ; defined from such a short series, so we use 'compute_neff.pro' which
    ; generates a very long synthetic AR1 series with same AR1 coefficient
    ; as the residuals, then smooths this with the same filters and computes
    ; the autocorrelation function and degrees of freedom from that.  To do
    ; this we need to the the lag-1 autocorrelation of the residuals if they
    ; weren't filtered.  So we need to apply the calibration to unfiltered
    ; data and then compute the residuals from that.  That's complicated, so
    ; for now, let's assume residuals were uncorrelated prior to filtering.
    ; In order to keep ALL degrees of freedom, I compute high-pass and
    ; band-pass cases by subtracting the low-pass cases from unfiltered
    ; cases!
    ;
    repeat begin
      xresid=randomn(seed,ncal)
      ar1=a_correlate(xresid,1)
      ar1=ar1(0)
    endrep until (abs(ar1) lt 0.05)
;    tfi=[mint,maxt]
;    if maxfac eq 1 then tsc=[0,0] else tsc=[0,-maxfac]
;    neffdata=compute_neff(xresid,tfilter=tfi,tscale=tsc)
;    serfac=neffdata(1)
    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])
      if maxfac eq 1 then ne0=neffdata1(0)-neffdata2(0) $
                     else ne0=neffdata1(0)-(1.+maxfac)*neffdata2(0)
    endif
    serfac=float(ncal)/float(ne0)
    if iband eq nband-1 then begin
      lagnr=a_correlate(resid4,indgen(30)+1)
      sf1=1.+2.*total(abs(lagnr))
      zlist=where(abs(lagnr) lt 0.2,nzero)
      if nzero gt 0 then lagnr(zlist(0):29)=0.
      sf2=1.+2.*total(abs(lagnr))
      print,serfac,sf1,sf2
    endif
;      serfac=1.
    ;
      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
  endfor
  ;
  pause
  xx=findgen(nband)+1
  xn=fix(nband/(7-fuse*2))
  xgap=fix(nband/xn)
  ix=findgen(xn)*xgap
  trenlist=where(allmax(ix) lt 0.,ntren)
  if ntren gt 0 then ix(trenlist)=ix(trenlist)-1
  xtickname=string(allmax(ix),format='(I3)')
  xtickname=string(allmin(ix),format='(I2)')+'-'+xtickname
  trenlist=where(allmax(ix) lt 0.,ntren)
  if ntren gt 0 then xtickname(trenlist)=' '
  ml=where(finite(rband) eq 0,nmiss)
  if nmiss gt 0 then xtickname(ml)=' '
  xn=xn+1
  ix=[ix,nband-1]
  xtickname=[xtickname,'>'+string(allmin(nband-1),format='(I2)')]
  plot,xx,rband,psym=def_sym(10),symsize=0.75,$
    /xstyle,xrange=[0,nband+1],xticks=xn-1,$
    xtickv=xx(ix),xtickname=xtickname,$
    /ystyle,yrange=[0,1],ytitle='Correlation',$
    title=regname(ireg)
  if trv eq 2 then begin
    oplot,xx,trwr(*,ireg),psym=def_sym(15),symsize=0.75
    oplot,xx,abdr(*,ireg),psym=def_sym(7),symsize=0.75
  endif
;  oplot,!x.crange,[0,0],linestyle=1
  ;
  print,max(bband)
  yr=[0,1.2]
  if ireg eq 2 then yr=[0,1.7]
  if ireg ge 9 then yr=[0,0.6]
  plot,xx,bband,psym=def_sym(10),symsize=0.75,$
    /xstyle,xrange=[0,nband+1],xticks=xn-1,$
    xtickv=xx(ix),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)*[-1,1]
    oplot,replicate(xx(iband),2),bband(iband)+bseband(iband)*[-2,2],$
      psym=def_sym(10),symsize=0.3
  endfor
  if trv eq 2 then begin
    oplot,xx,trwbeta(*,ireg),psym=def_sym(15),symsize=0.75
    oplot,xx,abdbeta(*,ireg),psym=def_sym(7),symsize=0.75
  endif
  allr(*,ireg)=rband(*)
  allbeta(*,ireg)=bband(*)
  allse(*,ireg)=bseband(*)
  allsfac(ireg)=sfac
  ;
endfor
;
; Save the results for later if TRW (to overlay on MXD)
;
if trv eq 1 then begin
  save,filename='quantify_tsdcal.idlsave',$
    allsfac,allbeta,allse,allr
endif
if trv eq 3 then begin
  save,filename='quantify_tsdcal_abd.idlsave',$
    allsfac,allbeta,allse,allr
endif
;
; Plot the average regressions
;
for ireg = 0 , nreg-1 do begin
  allbeta(*,ireg)=allbeta(*,ireg)/allsfac(ireg)
  allse(*,ireg)=allse(*,ireg)/allsfac(ireg)
endfor
totbeta=total(allbeta,2,/nan)
numbeta=total(finite(allbeta),2)
bband=totbeta/float(numbeta)
totse=total(allse^2,2,/nan)
numse=total(finite(allse),2)
bseband=sqrt(totse)/float(numse)
if trv eq 2 then begin
  !p.multi(0)=0
  mtit='(a)'
endif else begin
  !p.multi(0)=!p.multi(0)-2
  mtit='Average'
endelse
pause
yr=[0,3]
plot,xx,bband,psym=def_sym(10),symsize=0.75,$
  /xstyle,xrange=[0,nband+1],xticks=xn-1,$
  xtickv=xx(ix),xtickname=xtickname,$
  /ystyle,yrange=yr,ytitle='Regression slope (TEMP,'+titadd+')',$
  title=mtit
for iband = 0 , nband-1 do begin
  oplot,replicate(xx(iband),2),bband(iband)+bseband(iband)*[-1,1]
  oplot,replicate(xx(iband),2),bband(iband)+bseband(iband)*[-2,2],$
    psym=def_sym(10),symsize=0.3
endfor
;
; Plot the average regressions for TRW if necessary & ABD too!
;
if trv eq 2 then begin
  for iii = 0 , 1 do begin
    if iii eq 1 then begin
      allbeta=trwbeta
      allse=trwse
      allsfac=trwsfac
      ttt='TRW'
      isym=15
      partit='(c)'
    endif else begin
      allbeta=abdbeta
      allse=abdse
      allsfac=abdsfac
      ttt='ABD_MXD'
      isym=7
      partit='(b)'
    endelse
    for ireg = 0 , nreg-1 do begin
      allbeta(*,ireg)=allbeta(*,ireg)/allsfac(ireg)
      allse(*,ireg)=allse(*,ireg)/allsfac(ireg)
    endfor
    totbeta=total(allbeta,2,/nan)
    numbeta=total(finite(allbeta),2)
    bband=totbeta/float(numbeta)
    totse=total(allse^2,2,/nan)
    numse=total(finite(allse),2)
    bseband=sqrt(totse)/float(numse)
    yr=[0,3]
    plot,xx,bband,psym=def_sym(isym),symsize=0.75,$
      /xstyle,xrange=[0,nband+1],xticks=xn-1,$
      xtickv=xx(ix),xtickname=xtickname,$
      /ystyle,yrange=yr,ytitle='Regression slope (TEMP,'+ttt+')',$
      title=partit
    for iband = 0 , nband-1 do begin
      oplot,replicate(xx(iband),2),bband(iband)+bseband(iband)*[-1,1]
      oplot,replicate(xx(iband),2),bband(iband)+bseband(iband)*[-2,2],$
        psym=def_sym(10),symsize=0.3
    endfor
  endfor
endif
;
end