;
; Attempts to reconstruct observed SLP summer modes of variability
; from multiple linear regression on MXD principle components.
; Only those modes selected as having a good control on land air temperatures
; are attempted, and the potential predictors are pre-determined to reduce
; the possibility of overfitting the regression.
; HAVE TO JUST USE THE OVERLAPPING PERIODS FOR PCR!!!!
;
; Prepare for plotting
;
!p.multi=[2,1,4,0,0]
def_1color,50,color='red'
;
; Restore full MXD gridded dataset
;
print,'Reading in MXD data'
restore,filename='mxd_infill.idlsave'
; ntime,maxbox,boxdens,g,keepi,keepj,timey,baselist,nbase
;
; First restore the summer SLP modes
;
restore,filename='obs_summer_modes.idlsave'
;  timey,nretain,useinfill,usecorr,userot,xlon,ylat,nyr,nx,ny,fillea,ev
onret=nretain
oea=fillea
obsyr=timey
;
; Next restore the MXD modes
;
restore,filename='mxd_modes.idlsave'
;  timey,nretain,usefilt,thalf,g,ntime,ea,ev,maxbox,keepi,keepj
;
; Try to estimate the amplitude timeseries associated with each MXD
; EOF when using the UNFILTERED MXD data
;
; First remove long-term mean (not 61-90 mean!)
mkanomaly,boxdens
fullea=fltarr(ntime,nretain)
for iretain = 0 , nretain-1 do begin
  for iyr = 0 , ntime-1 do begin
    fullea(iyr,iretain)=total(boxdens(*,iyr)*ev(*,iretain))
  endfor
endfor
;for iretain = 0 , nretain-1 do begin
;  pause
;  plot,timey,fullea(*,iretain)
;  oplot,!x.crange,[0.,0.]
;  filter_cru,10.,/nan,tsin=reform(fullea(*,iretain)),tslow=tslow
;  oplot,timey,tslow,thick=3
;endfor
;
; Specify which modes to try, and their potential predictors
;
trymode=[0,1,1,0,0,1,1,0,1]       ; 0=no 1=yes
usepred=intarr(nretain,onret)
; MODES2USE:  1 2 3 4 5 6 7 8 9101112
usepred(*,1)=[0,0,0,0,1,0,0,0,0,0,0,1]
usepred(*,2)=[0,0,1,0,0,0,0,1,0,0,0,0]
usepred(*,5)=[0,0,0,0,0,0,0,0,0,0,0,0]
usepred(*,6)=[0,0,0,0,0,0,1,0,0,0,0,0]
usepred(*,8)=[0,0,1,0,0,1,0,0,0,0,0,0]
;
;-------------------------------------- OLD PREDICTORS (REFINED & REDUCED NOW)
;usepred(*,1)=[1,0,0,0,1,1,0,1,0,0,0,1]
;usepred(*,2)=[0,0,1,0,1,1,0,1,1,0,0,0]
;usepred(*,5)=[0,1,1,0,0,0,0,0,0,0,0,0]
;usepred(*,6)=[0,0,1,1,0,0,1,0,1,0,0,0]
;usepred(*,8)=[0,0,1,1,0,1,1,1,0,0,0,0]
;--------------------------------------
;
; Locate overlap period
;
minyr=obsyr(0)
maxyr=timey(ntime-1)
kslp=where(obsyr le maxyr,n1)
kmxd=where(timey ge minyr,n2)
if n1 ne n2 then message,'oops!'
;
; Try each mode in turn
;
mypred=fltarr(ntime,onret)*!values.f_nan
for imode = 0 , onret-1 do begin
  if trymode(imode) eq 1 then begin
    print,'TRYING TO RECONSTRUCT MODE',imode
    ;
    ; Extract timeseries to use
    ;
    depraw=reform(oea(*,imode))
    ; *** FILTER THE SLP PC ***
    filter_cru,thalf,/nan,tsin=depraw,tshigh=depts
    iuse=reform(usepred(*,imode))
    uselist=where(iuse eq 1,nuse)
    if nuse gt 0 then begin
      indts=reform(ea(*,uselist))
      indts=transpose(indts)        ; make it (modes,time)
      fullts=reform(fullea(*,uselist))
      fullts=transpose(fullts)        ; make it (modes,time)
      ;
      ; Perform multiple linear regression over entire series
      ;
      wts=fltarr(n1)+1.
      aa=indts(*,kmxd)
      bb=depts(kslp)
      pcrcoeff=regress(aa,bb,wts,yfit,pcrconst,pcrsigma,pcrf,lincorr,$
            mulcorr,pcrchi,pcrstat,/relative_weight)
      pcrcoeff=reform(pcrcoeff)
      ;
;      print,pcrcoeff
;      print,pcrsigma*1.5
      print,lincorr*100.
      print,pcrconst
      print,mulcorr*100.
      print,pcrf
      ;
      ; Also apply it to unfiltered predictors
      ;
      for ival = 0 , ntime-1 do begin
        mypred(ival,imode)=pcrconst+total(pcrcoeff*fullts(*,ival))
      endfor
      ;
      ; Now split into equal testing and training sets and redo regression
      ;
      nhalf1=fix(n1/2)
      nhalf2=n1-nhalf1
      for iset = 0 , 1 do begin
        ;
        ; Select period to use
        ;
        case iset of
          0: begin
            ktrain=indgen(nhalf2)
            ktest=indgen(nhalf1)+nhalf2
          end
          1: begin
            ktest=indgen(nhalf1)
            ktrain=indgen(nhalf2)+nhalf1
          end
        endcase
        ;
        ; Build model
        ;
        dep1=depts(kslp(ktrain))
        ind1=indts(*,kmxd(ktrain))
        wts=fltarr(nhalf2)+1.
        pcrcoeff=regress(ind1,dep1,wts,yfit,pcrconst,pcrsigma,pcrf,lincorr,$
              mulcorr,pcrchi,pcrstat,/relative_weight)
        pcrcoeff=reform(pcrcoeff)
if iset eq 0 then begin
  print,pcrcoeff
  print,pcrsigma*1.5
endif
        print,'Train r=',mulcorr
        ;
        ; Test it
        ;
        dep2=depts(kslp(ktest))
        ind2=indts(*,kmxd(ktest))
        estdep=fltarr(nhalf1)
        for ival = 0 , nhalf1-1 do begin
          estdep(ival)=pcrconst+total(pcrcoeff*ind2(*,ival))
        endfor
        print,'Test  r=',correlate(dep2,estdep)
        ;
      endfor
      ;
      ; Plot results (if required)
      ;
      p1=depraw
      p2=reform(mypred(*,imode))
      mknormal,p1,obsyr,refsd=refsd,refmean=refmean,refperiod=[minyr,maxyr]
      mknormal,p2,timey,refperiod=[minyr,maxyr]
      p2=(p2*refsd(0))+refmean(0)
    if imode eq 2 then begin
      plot,timey,p2,/nodata,$
        /xstyle,xrange=[1695,1995.],xtitle='Year',$
        /ystyle,yrange=[-15.,15.],ytitle='PC time series (obs & pred)',$
        title='(c)  Mode #'+string(imode+1,format='(I2)')
      oplot,!x.crange,[0,0]
      oplot,obsyr,depraw
      oplot,timey,p2,color=50
      filter_cru,20.,/nan,tsin=depraw,tslow=ppp
      oplot,obsyr,ppp,thick=3
      filter_cru,20.,/nan,tsin=p2,tslow=ppp
      oplot,timey,ppp,color=50,thick=3
    endif
      ; 
      ; Correlate full and low-pass over overlap period
      ;
      print,'RAW R=',correlate(p2(kmxd),depraw(kslp))
      filter_cru,10.,/nan,tsin=p2,tslow=lowmxd
      filter_cru,10.,/nan,tsin=depraw,tslow=lowslp
      print,'LOW R=',correlate(lowmxd(kmxd),lowslp(kslp))
      ;
    endif
    ;
    print,'------------------------'
    ;
  endif
endfor
;
end