;
; 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*
;
docol=1      ; 0=black&white, 1=color
dops=1       ; 0=screen (though can do PS manually), 1=separate PS files
doper=[1922,1976]
;
; Prepare for plotting
;
if dops eq 0 then begin
  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
  if docol eq 1 then begin
    def_1color,51,color='red'
    def_1color,50,color='black'
  endif else begin
    def_1color,50,color='black'
    def_1color,51,color='white'
  endelse
  def_1color,52,color='lgrey'
endif
;
; 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 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'
;
if total(abs(obsyr-timey)) ne 0 then message,'Ooops!'
;
kper=where((timey ge doper(0)) and (timey le doper(1)),n1)
;oea=oea(kl,*)
;fillea=fillea(kl,*)
;
; 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(0:7,*)=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.
      ;
      nhalf1=fix(n1*0.4)
      nhalf2=n1-nhalf1
      nhalf2a=fix(nhalf2/2)
      nhalf2b=nhalf2-nhalf2a
      for iset = 0 , 2 do begin
        ;
        ; Select period to use
        ;
        case iset of
          0: begin
;            ktrain=indgen(nhalf2)
;            ktest=indgen(nhalf1)+nhalf2
            ktrain=[indgen(nhalf2a),indgen(nhalf2b)+nhalf2a+nhalf1]
            ktest=indgen(nhalf1)+nhalf2a
          end
          1: begin
;            ktest=indgen(nhalf1)
;            ktrain=indgen(nhalf2)+nhalf1
            ktest=indgen(nhalf1)
            ktrain=indgen(nhalf2)+nhalf1
          end
          2: begin
            ktrain=indgen(n1)
            ktest=!values.f_nan
          end
        endcase
        ;
        ; If first attempt, use screening regression to select predictors
        ;
        if iset eq 0 then begin
          dep1=depts(kper[ktrain])
          ind1=indts(*,kper[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(kper[ktrain])
        ind1=indts(*,kper[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(*)
            help,yfit,dep1
            r3=mkcorrelation(reform(yfit),dep1,filter=10.)
            r3str=string(round(r3*100.),format='(3I4)')
            print,r3str
            for iyr = 0 , nyr-1 do begin
              ind3=indts(ineq,*)
              mypred(iyr,imode)=pcrconst+total(pcrcoeff*ind3(*,iyr))
            endfor
            end
        endcase
        print,'Train r=',mulcorr
        ;
        ; Test it
        ;
        if iset ne 2 then begin
          dep2=depts(kper[ktest])
          ind2=indts(*,kper[ktest])
          ind2=ind2(ineq,*)
          estdep=fltarr(nhalf1)
          for ival = 0 , nhalf1-1 do begin
            estdep(ival)=pcrconst+total(pcrcoeff*ind2(*,ival))
          endfor
          testr=correlate(dep2,estdep)
          print,'Test  r=',testr
          if iset eq 0 then verr=testr
        endif
        ;
endif
        ;
      endfor
      ;
if finite(ineq(0)) ne 0 then begin
      ;
      if dops ne 0 then begin
        ploton,1
        loadct,39
        multi_plot,nrow=3,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=14
        endelse
        if docol eq 1 then begin
          def_1color,51,color='red'
          def_1color,50,color='black'
        endif else begin
          def_1color,50,color='black'
          def_1color,51,color='white'
        endelse
        def_1color,52,color='lgrey'
      endif
      ;
      ; Plot results
      ;
      p1=depts
      p2=reform(mypred(*,imode))
      mknormal,p1,timey,refper=doper,refsd=refsd,refmean=refmean
      p1=depts
      mknormal,p2,timey,refper=doper
      p2=(p2*refsd(0))+refmean(0)
      pause
      r1=mkcorrelation(p1,p2,timey,refperiod=doper)
      r2=mkcorrelation(p1,p2,timey,refperiod=[1977,1998])
      r3str=' ver_r='+strcompress(string(verr,format='(F5.2)'),/remove_all)+$
        ' fit_r='+strcompress(string(r1[0],format='(F5.2)'),/remove_all)+$
        ' recent_r='+strcompress(string(r2[0],format='(F5.2)'),/remove_all)
      plot,timey,p1,/nodata,$
        /xstyle,xtitle='Year',$
        /ystyle,yrange=[-24.,24.],ytitle='PC (obs & pred)',$
        title='Mode #'+string(imode+1,format='(I2)')+r3str
      polyfill,doper[[0,1,1,0,0]],!y.crange[[0,0,1,1,0]],color=52
      axis,xaxis=0,/xstyle
      axis,xaxis=1,/xstyle,xtickformat='nolabels'
      axis,yaxis=0,/ystyle
      oplot,!x.crange,[0.,0.]
      oplot,timey,p1,color=50,thick=1+3*(1-docol)
      oplot,timey,p1,color=51,thick=2
      oplot,timey,p2,color=50,thick=2
      filter_cru,10.,tsin=p1,tslow=p1low
      filter_cru,10.,tsin=p2,tslow=p2low
      oplot,timey,p1low,thick=5+3*(1-docol),color=50
      oplot,timey,p1low,thick=5,color=51
      oplot,timey,p2low,thick=5,color=50
      ;
      if dops ne 0 then begin
        plotoff
        if docol eq 0 then ctit='' else ctit='_col'
        spawn,'mv idl.ps1 holpt2_fig'+$
          strcompress(string(imode+11,format='(I2)'),/remove_all)+$
          'd'+ctit+'.ps'
      endif
endif
      ;
    endif
    ;
    print,'------------------------'
    ;
  endif
endfor
;
end