;
; Attempts to reconstruct observed SLP summer modes of variability
; from stepwise screening multiple linear regression on observed land air
; temperature 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.
; *THIS VERSION DOES NOT USE HIGH-PASS FILTERED DATA*
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=4,ncol=2,layout='large'
if !d.name eq 'X' then begin
  window,ysize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=9
endelse
def_1color,50,color='red'
;
; 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
;
; Next restore the summer land air temperature modes
;
restore,filename='olat_modes.idlsave'
;  timey,nretain,usefilt,usecoloc,thalf,xlon,ylat,nyr,nx,ny,fillea,ev
if nretain ne 12 then message,'I think we only want to offer 12 predictors'
;
; Specify which modes to try, and their potential predictors
;
;trymode=[0,1,1,0,0,1,1,0,1]       ; 0=no 1=yes
trymode=[1,1,1,1,1,1,1,1,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,1,0,0,0,0]
usepred(*,2)=[0,0,1,0,0,1,0,0,0,0,0,0]
usepred(*,5)=[0,0,0,0,1,1,0,0,0,0,0,0]
usepred(*,6)=[0,0,1,1,0,0,1,1,0,0,0,0]
usepred(*,8)=[0,0,1,1,0,0,0,0,0,0,0,0]
usepred(*,*)=1
;
; Try each mode in turn
;
mypred=fltarr(nyr,onret)*!values.f_nan
for imode = 0 , onret-1 do begin
  if trymode(imode) eq 1 then begin
    print,'TRYING TO RECONSTRUCT MODE',imode+1
    ;
    ; Extract timeseries to use (simple cos both predictors and predictands
    ; already cover exactly the same period)
    ;
    depts=reform(oea(*,imode))
    iuse=reform(usepred(*,imode))
    uselist=where(iuse eq 1,nuse)
    if nuse gt 0 then begin
      indts=reform(fillea(*,uselist))
      indts=transpose(indts)        ; make it (modes,time)
      vn='PC#'+string(uselist+1,format='(I2.2)')
      ;
      ; Perform multiple linear regression over first half, second half,
      ; and full record.  Use stepwise screening on the first half, then
      ; subsequently use the same screening predictors.
      ;
      nhalf=nyr/2
      for iset = 0 , 2 do begin
        ;
        ; Select period to use
        ;
        case iset of
          0: begin
            ktrain=indgen(nhalf)
            ktest=ktrain+nhalf
          end
          1: begin
            ktest=indgen(nhalf)
            ktrain=ktest+nhalf
          end
          2: begin
            ktrain=indgen(nyr)
            ktest=!values.f_nan
          end
        endcase
        ;
        ; If first attempt, use screening regression to select predictors
        ;
        if iset eq 0 then begin
          dep1=depts(ktrain)
          ind1=indts(*,ktrain)
          ineq=!values.f_nan     ; NaN means do not force any into regression
          stepwise,ind1,dep1,ineq,0.05,0.05,varnames=vn
          print,ineq
        endif
        ;
        ; Build model using screened predictors
        ;
if finite(ineq(0)) ne 0 then begin
        ;
        dep1=depts(ktrain)
        ind1=indts(*,ktrain)
        ind1=ind1(ineq,*)
        wts=fltarr(n_elements(dep1))+1.
        pcrcoeff=regress(ind1,dep1,wts,yfit,pcrconst,pcrsigma,pcrf,lincorr,$
              mulcorr,pcrchi,pcrstat,/relative_weight)
        pcrcoeff=reform(pcrcoeff)
        case iset of
          0: begin
            print,'First half'
            end
          1: begin
            print,'Second half'
            end
          2: begin
            print,'Full period'
            mypred(*,imode)=yfit(*)
            r3=mkcorrelation(reform(yfit),dep1,filter=10.)
            r3str=string(round(r3*100.),format='(3I4)')
            print,r3str
            end
        endcase
        print,'Train r=',mulcorr
        ;
        ; Test it
        ;
        if iset ne 2 then begin
          dep2=depts(ktest)
          ind2=indts(*,ktest)
          ind2=ind2(ineq,*)
          help,ind2
          estdep=fltarr(nhalf)
          for ival = 0 , nhalf-1 do begin
            estdep(ival)=pcrconst+total(pcrcoeff*ind2(*,ival))
          endfor
          print,'Test  r=',correlate(dep2,estdep)
        endif
        ;
endif
        ;
      endfor
      ;
      ; Plot results
      ;
if finite(ineq(0)) ne 0 then begin
      p1=depts
      p2=reform(mypred(*,imode))
      mknormal,p1,refsd=refsd,refmean=refmean
      p1=depts
      mknormal,p2
      p2=(p2*refsd(0))+refmean(0)
      pause
      plot,timey,p1,$
        /xstyle,xtitle='Year',$
        /ystyle,yrange=[-15.,15.],ytitle='PC (obs & pred)',$
        title='Mode #'+string(imode+1,format='(I2)')+r3str
      oplot,!x.crange,[0.,0.]
      oplot,timey,p2,color=50
      filter_cru,10.,tsin=p1,tslow=p1low
      filter_cru,10.,tsin=p2,tslow=p2low
      oplot,timey,p1low,thick=5
      oplot,timey,p2low,thick=5,color=50
endif
      ;
    endif
    ;
    print,'------------------------'
    ;
  endif
endfor
;
end