;
; Attempts to reconstruct observed SLP summer modes of variability
; from stepwise screening multiple linear regression on calibrated, corrected
; 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.
; *THIS VERSION DOES NOT USE HIGH-PASS FILTERED DATA*
; Have to select just the overlapping period to perform the PCR on.
;
dopub=2        ; 0=all results etc.  1=publication graph only (calib period)
               ; 2=publication graph (full reconstruction period)
docol=1        ; 0=black&white, 1=colour
dops=1       ; 0=screen (though can do PS manually), 1=separate PS files
;
case dopub of
  0: begin
   nc=1
   xr=[1697,1998]
  end
  1: begin
   nc=2
   xr=[1922,1998]
  end
  2: begin
   nc=1
   xr=[1697,1998]
  end
endcase
;
; Prepare for plotting
;
if dops eq 0 then begin
  loadct,39
  multi_plot,nrow=4,ncol=nc,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 ne 0 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
;!y.margin=[0,0]
;!y.omargin=[4,2]
;
; 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='mxdmodes_16971976unrcoradjraw.idlsave'
;  mxdyear,nretain,usefilt,usecorr,usedecline,useper,userot,thalf
;  g,mxdnyr,ea,ev,lampct
ntime=mxdnyr
timey=mxdyear
fillea=ea
if nretain ne 12 then message,'Only want to provide 12 potential predictors'
;
; Specify which modes to try, and their potential predictors
;
trymode=[1,1,1,1,1,1,1,1]       ; 0=no 1=yes
;if dopub eq 1 then trymode=[0,1,1,0,1,0,1,0,1]
;trymode=[0,1,1,0,0,1,1,0,1]       ; 0=no 1=yes
;trymode=[0,1,0,0,0,0,0,0,0]       ; 0=no 1=yes
usepred=intarr(nretain,onret)
; MODES2USE:  1 2 3 4 5 6 7 8 9101112
;usepred(*,1)=[0,0,0,1,1,0,0,0,0,0,0,0]
;usepred(*,2)=[0,1,1,0,0,0,0,0,0,0,0,0]
;usepred(*,5)=[0,1,0,0,1,0,0,0,0,1,0,0]
;usepred(*,6)=[0,1,1,0,0,0,0,1,0,0,0,0]
;usepred(*,8)=[0,0,1,0,0,0,0,1,0,0,0,0]
usepred(*,*)=0
usepred(0:7,*)=1
;
; 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,'Oooops!'
;
;ptit='('+['a','b','c','d','e','f','g','h','i','j']+') '
;npart=0
;
; 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+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
;        print,ktrain,format='(20I4)'
;        print,ktest,format='(20I4)'
        ;
        ; If first attempt, use screening regression to select predictors
        ;
        if iset eq 0 then begin
          dep1=depts(kslp(ktrain))
          ind1=indts(*,kmxd(ktrain))
          ineq=!values.f_nan     ; NaN means do not force any into regression
          help,dep1
          stepwise,ind1,dep1,ineq,0.050,0.050,varnames=vn
        endif
        if finite(ineq(0)) then begin
          ;
          ; Build model using screened predictors
          ;
          dep1=depts(kslp(ktrain))
          ind1=indts(*,kmxd(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'
              r3=mkcorrelation(reform(yfit),dep1,filter=10.)
              fitr=r3[0]
              r3str=string(round(r3*100.),format='(3I4)')
              print,r3str
              for iyr = 0 , ntime-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(kslp(ktest))
            ind2=indts(*,kmxd(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
      ;
      ; Plot results
      ;
      if finite(ineq(0)) then begin
        ;
        if dops ne 0 then begin
          ploton,1
          loadct,39
          multi_plot,nrow=3,ncol=nc,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 ne 0 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
        ;
        p1=depts
        p2=reform(mypred(*,imode))
        mknormal,p1,obsyr,refsd=refsd,refmean=refmean,refperiod=[minyr,maxyr]
        p1=depts
        mknormal,p2,timey,refperiod=[minyr,maxyr]
        p2=(p2*refsd(0))+refmean(0)
        pause
        if dopub eq 1 then begin
          mtit='Mode #'+string(imode+1,format='(I2)')+$
            ' ver_r='+strcompress(string(verr,format='(F5.2)'),/remove_all)+$
            ' fit_r='+strcompress(string(fitr,format='(F5.2)'),/remove_all)
        endif else begin
          mtit=''
        endelse
        plot,obsyr,p1,/nodata,$
          /xstyle,xrange=xr,xtitle='Year',$
          /ystyle,yrange=[-24.,24.],ytitle='PC (obs & pred)',$
          title=mtit
        doper=[1922,1976]
        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
;        if dopub eq 0 then begin
;          xyouts,1703,12.,ptit(npart)+'Mode #'+$
;            string(imode+1,format='(I1)')+'  '+r3str
;        endif else begin
;          xyouts,1703,12.,ptit(npart)+'Mode #'+string(imode+1,format='(I1)')
;        endelse
;        npart=npart+1
;        if !p.multi(0) eq 0 then axis,xaxis=0,/xstyle,xtitle='Year'
        oplot,!x.crange,[0.,0.]
        if dopub ne 2 then begin
          if docol eq 0 then oplot,obsyr,p1,color=50,thick=1+3*dopub
          oplot,obsyr,p1,color=51,thick=2
        endif
        oplot,timey,p2,color=50,thick=2
        filter_cru,10.,tsin=p1,tslow=p1low
        filter_cru,10.,tsin=p2,tslow=p2low
        if docol eq 0 then oplot,obsyr,p1low,color=50,thick=5+3*dopub
        oplot,obsyr,p1low,thick=5,color=51
        oplot,timey,p2low,thick=5,color=50
        ;
        if dops ne 0 then begin
          plotoff
          if dopub eq 1 then ftit='e' else ftit='f'
          if docol ne 0 then ctit='_col' else ctit=''
          spawn,'mv idl.ps1 holpt2_fig'+$
            strcompress(string(imode+11,format='(I2)'),/remove_all)+$
            ftit+ctit+'.ps'
        endif
        ;
      endif
      ;
    endif
    ;
    print,'------------------------'
    ;
  endif
endfor
;
end