;
; 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=1        ; 0=all results etc.  1=publication graph only
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=5,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 dopub eq 0 then begin
  def_1color,50,color='red'
  def_1color,51,color='black'
endif else begin
  def_1color,50,color='black'
  def_1color,51,color='white'
endelse
!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,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(*,*)=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/2)
      nhalf2=n1-nhalf1
      for iset = 0 , 2 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
          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(kslp(ktrain))
          ind1=indts(*,kmxd(ktrain))
          ineq=!values.f_nan     ; NaN means do not force any into regression
          stepwise,ind1,dep1,ineq,0.075,0.075,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.)
              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
            print,'Test  r=',correlate(dep2,estdep)
          endif
          ;
        endif
        ;
      endfor
      ;
      ; Plot results
      ;
      if finite(ineq(0)) then begin
        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
        plot,obsyr,p1,/nodata,$
          /xstyle,xrange=[1697,1995],xtickformat='nolabels',$  ;xtitle='Year',$
          /ystyle,yrange=[-17.,17.],ytitle='PC (obs & pred)'   ;,$
;          title='Mode #'+string(imode+1,format='(I2)')+r3str
        xyouts,1703,12.,ptit(npart)+'Mode #'+string(imode+1,format='(I1)')
        npart=npart+1
        if !p.multi(0) eq 0 then axis,xaxis=0,/xstyle,xtitle='Year'
        oplot,!x.crange,[0.,0.]
        oplot,obsyr,p1,color=50,thick=1+3*dopub
        oplot,obsyr,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,obsyr,p1low,thick=5+3*dopub,color=50
        oplot,obsyr,p1low,thick=5,color=51
        oplot,timey,p2low,thick=5,color=50
      endif
      ;
    endif
    ;
    print,'------------------------'
    ;
  endif
endfor
;
end