;
; Tries to reconstruct Apr-Sep temperatures, on a box-by-box basis, from the
; EOFs of the MXD data set.  This is PCR, although PCs are used as predictors
; but not as predictands.  This PCR-infilling must be done for a number of
; periods, with different EOFs for each period (due to different spatial
; coverage).  *BUT* don't do special PCR for the modern period (post-1976),
; since they won't be used due to the decline/correction problem.
; Certain boxes that appear to reconstruct well are "manually" removed because
; they are isolated and away from any trees.
;
minr=0.40     ; minimum validated correlation acceptable for accepting
              ; 50%=0.71 40%=0.63 30%=0.55 20%=0.45 10%=0.32
;;;;;npred=8       ; no. of EOFs to use as predictors
dotrim=1      ; 0=don't 1=do trim away isolated grid boxes
;
loadct,39
multi_plot,nrow=2
if !d.name eq 'X' then window,ysize=800
def_1color,10,color='white'
def_1color,12,color='mgrey'
;
for iper = 0 , 6 do begin
;
case iper of
  0: begin & perst='14001976' & npred=2 & end
  1: begin & perst='14531976' & npred=3 & end
  2: begin & perst='15831976' & npred=5 & end
  3: begin & perst='16601976' & npred=6 & end
  4: begin & perst='16971976' & npred=8 & end
  5: begin & perst='17431976' & npred=10 & end
  6: begin & perst='18221976' & npred=12 & end
endcase
;
; Restore Apr-Sep temperature gridded dataset
;
restore,filename='../summer_modes/obs_temp_as.idlsave'
;  timey,fdseas,xlon,ylat,nyr,nx,ny,fdltm,fdsd,missfrac,landonly
;  seaonly,mainlysea
;
; Now get the EOFs of the MXD data that are to be used as predictors
;
restore,filename='mxdmodes_abdlow_'+perst+'unrcoradjraw.idlsave'
;  mxdyear,nretain,usefilt,thalf,g,mxdnyr,ea,ev,lampct,usecorr,$
;  usedecline,useper,userot
nretain=npred
ea=ea(*,0:nretain-1)
;
; Identify the overlapping period, used for calibration & verification
;
kmxd=where(mxdyear ge timey(0),nmxd)
ktem=where(timey le mxdyear(mxdnyr-1),ntem)
if nmxd ne ntem then message,'Oooops!'
;
fdseas=fdseas(*,*,ktem)
timey=timey(ktem)
easub=ea(kmxd,*)
;
; Make a field to contain status of each box (0=can't reconstruct 1=can)
; and one to contain the reconstruction
;
fdstatus=fltarr(nx,ny)
recontemp=fltarr(nx,ny,mxdnyr)
recontemp(*,*,*)=!values.f_nan
;
; For each box, try to reconstruct with PCR
;
for ix = 0 , nx-1 do begin
  for iy = 0 , ny-1 do begin
    print,xlon(ix),ylat(iy),format='($,2F8.1)'
    ;
    ; Extract temp series and check for sufficient non-missing data
    ;
    onets=reform(fdseas(ix,iy,*))
    kgot=where(finite(onets),ngot)
    ;
    ; Because of the need for separate testing/training sets, need at least 40
    ; years of data
    ;
    if ngot ge 40 then begin
      ;
      ; Split data into two equal segments
      ;
      setlen1=fix(ngot/2)
      setlen2=ngot-setlen1
      kset1=kgot(0:setlen1-1)
      kset2=kgot(setlen1:ngot-1)
      ;
      for iter = 0 , 2 do begin  ; repeat for test,training & vice versa & full
        case iter of
          0: begin
            ktrain=kset1
            ktest=kset2
          end
          1: begin
            rmean=rskill
            ktrain=kset2
            ktest=kset1
          end
          2: begin
            rmean=(rmean+rskill)*0.5
            print,rmean,format='($,F7.2)'
            ktrain=kgot
            ktest=!values.f_nan ; when using the full period, cannot test it
          end
        endcase
        ;
        ; Build PCR (PC-based multiple linear regression)
        ;
        depts=onets(ktrain)
        indts=transpose(easub(ktrain,*))
        wts=fltarr(n_elements(ktrain))+1.
        pcrcoeff=regress(indts,depts,wts,yfit,pcrconst,pcrsigma,pcrf,lincorr,$
          mulcorr,pcrchi,pcrstat,/relative_weight)
        ;
        ; If regression failed, set skill to zero, else evaluate skill
        ;
  if iter lt 2 then begin
        if pcrstat ne 0 then begin
          rskill=-9.99
        endif else begin
          depts=onets(ktest)
          indts=transpose(easub(ktest,*))
          nval=n_elements(ktest)
          prets=fltarr(nval)
          for ival = 0 , nval-1 do begin
            prets(ival)=pcrconst+total(pcrcoeff*indts(*,ival))
          endfor
          rskill=correlate(depts,prets)
        endelse
        print,mulcorr,rskill,format='($,2F7.2)'
        ;
        ; Determine whether skill is high enough, and store reconstruction
        ; if it is
        ;
  endif else begin
        if rmean ge minr then begin
          fdstatus(ix,iy)=1.
          print,'  SUCCESS'
          indts=transpose(ea)
          for ival = 0 , mxdnyr-1 do begin
            recontemp(ix,iy,ival)=pcrconst+total(pcrcoeff*indts(*,ival))
          endfor
        endif else begin
          print,'  FAILED'
        endelse
  endelse
        ;
      endfor
      ;
    endif else begin
      print,'  TOO SHORT:',ngot
    endelse
    ;
  endfor
endfor
;
; Manually remove some unwanted isolated reconstructions
;
case iper of
  0: begin
    xremove=[ 177.5,  97.5]
    yremove=[  57.5,  37.5]
    end
  1: begin
    xremove=[ 177.5, -37.5, -62.5, -87.5]
    yremove=[  57.5,  57.5,  67.5,  57.5]
    end
  2: begin
    xremove=[ 137.5, -37.5,-157.5]
    yremove=[  17.5,  57.5,  42.5]
    end
  3: begin
    xremove=[ 142.5, 137.5, 132.5,  12.5,  -2.5, -37.5,-127.5,-152.5,-157.5]
    yremove=[  17.5,  17.5,   7.5,  27.5,  27.5,  57.5,   2.5,  42.5,  42.5]
    end
  4: begin
    xremove=[ 132.5,  67.5,  -2.5,  -7.5,-177.5]
    yremove=[   7.5,  17.5,  27.5,  27.5,  12.5]
    end
  5: begin
    xremove=[ 132.5,  82.5,  67.5,  67.5,  62.5,  -2.5]
    yremove=[   7.5,   2.5,   2.5,  17.5,   2.5,  27.5]
    end
  6: begin
    xremove=[ 177.5, 152.5, 127.5,  82.5,  67.5,  -2.5,-157.5]
    yremove=[  27.5,  22.5,   2.5,   2.5,   7.5,  27.5,  47.5]
    end
endcase
nrem=n_elements(xremove)
if dotrim ne 0 then begin
  for irem = 0 , nrem-1 do begin
    iix=where(xlon eq xremove(irem))
    iiy=where(ylat eq yremove(irem))
    fdstatus(iix,iiy)=0.
    recontemp(iix,iiy,*)=!values.f_nan
  endfor
endif else begin
  nrem=0
endelse
;
; Now save this PCR-based reconstruction
;
save,filename='calibmxd5_abdlow_pcr'+perst+'.idlsave',$
  mxdyear,mxdnyr,nx,ny,xlon,ylat,fdstatus,recontemp
;
dummy=where(fdstatus eq 1,ninfill)
print,'Total reconstructed = ',ninfill,'  out of total = ',nx*ny
;
; Plot layout of missing & infilled and complete data boxes
;
pause
map=def_map(/npolar) & map.limit(0)=0.
coast=def_coast(/get_device) & coast.double=0
sm=def_sm(/off) & sm.cyclic=0
labels=def_labels(/off)
inter_boxfd,fdstatus,xlon,ylat,map=map,$
  labels=labels,coast=coast,sm=sm,$
  levels=[-0.5,0.5,1.5],c_colors=[10,12],$
  title='PCR-reconstruction for: '+perst+$
    '  boxes:'+string(ninfill,format='(I4)')+' + '+string(nrem,format='(I2)')
;if dotrim ne 0 then map_plots,xremove,yremove,psym=def_sym(10)
map_plots,xremove,yremove,psym=def_sym(10)
;
endfor
;
end