;
; ***MODIFIED TO ALSO DO CLASSICAL REGRESSION!!!!***
;
; Calibrates and plots regional MXD against regional temperature
; DOES IT FOR AGE-BANDED DENSITY REGIONS AND NH
;
; Must first select which version of the age-banded records to use.
; All options affect the "all" or "nh" series, but some also affect
; the individual regions.
fnband='_20-550'   ; age bands used.  Choices: '_50-300' ''=20-550
fnmint='_min2'     ; minimum trees in each band.  Choices: '_min2' ''=all used
fnepsw='_eps'          ; weighting used for all.  Choices: '_eps' ''=no weighting
;
doclassic=0        ; 0=normal inverse regression, 1=classical regression
;
dosave=2     ; save results or not?  0=no, 1=IDLSAVE, 2=IDLSAVE & ASCII
;
; Define periods for calibration
;
calper=[1881,1960]        ; calibration period
verper=[1961,1994]        ; verification period
thalf=10.
donh=1                    ; 0=do NH 1=do ALL
;
if donh eq 0 then nhtit='NH' else nhtit='ALL'
;
; Get the hemispheric temperature series
;
restore,filename='datastore/compbest_mxd_fixed1950.idlsave'
if donh eq 0 then nhtemp=reform(comptemp(*,3)) $
             else nhtemp=reform(comptemp(*,2))
; Get rid of pre-1871 temperatures
knh=where(timey lt 1871)
nhtemp(knh)=!values.f_nan
;
; Get all the age-banded MXD series (regions, hemi and back to 1400)
;
restore,filename='/cru/u2/f055/treeharry/banding/bandregnh'+$
  fnband+fnmint+fnepsw+'_fixed.idlsave'
; Gets: timey,nyr,nhts,allts,alleps,nreg,regname
yrmxd=timey
fullmxd=allts
; Combine hemispheric MXD with regional MXD
nfull=n_elements(yrmxd)
newfm=fltarr(nfull,nreg+1)*!values.f_nan
newfm(*,0:nreg-1)=transpose(fullmxd(*,*))
newfm(*,nreg)=nhts(*)
fullmxd=newfm
;
; Get regional temperature series
;
restore,filename='datastore/regbest_mxd_fixed1950.idlsave'
; Gets: nreg,regname,regtemp,rawtemp,timey,bestmxd etc.
;
; Identify overlap between instrumental and MXD series
;
nyr=n_elements(timey)
kcomp=where( (yrmxd ge timey(0)) and (yrmxd le timey(nyr-1)) )
bestmxd=fullmxd(kcomp,*)
;
; Combine hemispheric and regional temperatures
;
regname=[regname,nhtit]
newrt=fltarr(nyr,nreg+1)*!values.f_nan
newrt(*,0:nreg-1)=regtemp(*,*)
newrt(*,nreg)=nhtemp(*)
regtemp=newrt
nreg=nreg+1
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=3,layout='large',/landscape
if !d.name eq 'X' then begin
  window,xsize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=9
endelse
def_1color,10,color='blue'
def_1color,11,color='red'
;
range0=[-2.5,-2.0,-3.5,-2.5,-1.5,-0.5,-2.5,-2.0,-3.0,-1.5]
range1=[ 1.5, 2.0, 2.5, 1.5, 1.0, 1.0, 1.5, 1.5, 1.5, 1.0]
;
calregts=fltarr(nfull,nreg)*!values.f_nan
;
; Process each region separately
;
for ireg = 0 , nreg-1 do begin
  print,regname(ireg)
  ;
  ; Extract series
  ;
  tstemp=reform(regtemp(*,ireg))
  tsmxd=reform(bestmxd(*,ireg))
  tsfull=reform(fullmxd(*,ireg))
  ;
  ; Identify calibration and verification subsets
  ;
  calkl=where( finite(tstemp+tsmxd) and $
               (timey ge calper(0)) and (timey le calper(1)) , ncal )
  verkl=where( finite(tstemp+tsmxd) and $
               (timey ge verper(0)) and (timey le verper(1)) , nver )
  print,'Calibration period:'+string(calper,format='(2I5)')+$
    '  length:'+string(ncal,format='(I5)')
  print,'Verification period:'+string(verper,format='(2I5)')+$
    '  length:'+string(nver,format='(I5)')
  ;
  ; Calibrate the MXD
  ;
dummy=where(finite(tsmxd),nnnnnn)
if nnnnnn gt 0 then begin
  if doclassic then print,'CLASSICAL REGRESSION!!!!!!!' else print,'INVERSE REGRESSION!!!!!'
  mkcalibrate,tsmxd(calkl),tstemp(calkl),fitcoeff=fitcoeff,autocoeff=autocoeff,classic=doclassic
  print,'RAW'
  print,'     r   alpha    beta SEalpha  SEbeta    sig    rmse'
  print,fitcoeff,format='(F6.2,4F8.4,F7.2,F8.4)'
  print,' a_mxd a_tem a_resid   ess'
  print,autocoeff,format='(3F6.2,F8.1)'
  ;
  ; And low and high passed too!
  ;
  filter_cru,/nan,thalf,tsin=tsmxd,tshigh=mxdhi,tslow=mxdlo
  filter_cru,/nan,thalf,tsin=tstemp,tshigh=temphi,tslow=templo
  mkcalibrate,mxdhi(calkl),temphi(calkl),fitcoeff=hicoeff,autocoeff=autohi,classic=doclasssic
  print,'HIGH-PASSED',thalf
  print,'     r   alpha    beta SEalpha  SEbeta    sig    rmse'
  print,hicoeff,format='(F6.2,4F8.4,F7.2,F8.4)'
  print,' a_mxd a_tem a_resid   ess'
  print,autohi,format='(3F6.2,F8.1)'
  mkcalibrate,mxdlo(calkl),templo(calkl),fitcoeff=locoeff,autocoeff=autolo,$
    nlag=5,classic=doclassic
  print,'LOW-PASSED',thalf
  print,'     r   alpha    beta SEalpha  SEbeta    sig    rmse'
  print,locoeff,format='(F6.2,4F8.4,F7.2,F8.4)'
  print,' a_mxd a_tem a_resid   ess'
  print,autolo,format='(3F6.2,F8.1)'
  ;
  ; Generate calibrated record and uncertainties etc.
  ; Use the extended MXD record (after checking it matches the short one!)
  ;
  xxx=tsfull(kcomp)
  yyy=tsmxd
  tse=total( (xxx-yyy)^2 , /nan )
;  print,tse
  if tse gt 0.03 then message,'Series do not match'
  if doclassic eq 0 then begin
    calmxd=tsfull*fitcoeff(2)+fitcoeff(1)
  endif else begin
    calmxd=(tsfull-fitcoeff(1))/fitcoeff(2)  ; rearrange classic regression to get T=f(MXD)
  endelse
  calregts(*,ireg)=calmxd(*)
  filter_cru,/nan,thalf,tsin=tsfull,tshigh=fullhi,tslow=fulllo
  if doclassic eq 0 then begin
    sepmxd=(fullhi*hicoeff(2)+hicoeff(1))+(fulllo*locoeff(2)+locoeff(1))
  endif else begin
    sepmxd=((fullhi-hicoeff(1))/hicoeff(2))+((fulllo-locoeff(1))/locoeff(2))  ; rearrange classic regression to get T=f(MXD)
  endelse
  ;
  ; For overlap period with non-missing data, compute correlation between
  ; the temperature series and the 2-band reconstruction
  ;
  kkk=calkl+(timey(0)-yrmxd(0))
  print,'r(Two-band calibrated series,Temp)='+$
    string(correlate(tstemp(calkl),sepmxd(kkk)),format='(F6.2)')
  ;
  if regname(ireg) eq nhtit then !p.multi=[0,1,2,0,0]
  pause
  plot,timey,tstemp,$
    /xstyle,xrange=[1400,1994],xtitle='Year',$
    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
    /ystyle,yrange=[range0(ireg),range1(ireg)],$
    title=regname(ireg)
  oplot,yrmxd,calmxd,color=10       ;,thick=3
;  oplot,yrmxd,sepmxd,color=11       ;,thick=3
  oplot,!x.crange,[0.,0.],linestyle=1
  filter_cru,/nan,thalf,tsin=tstemp,tslow=ylow
  oplot,timey,ylow,thick=3
  filter_cru,/nan,thalf,tsin=calmxd,tslow=ylow
  oplot,yrmxd,ylow,thick=3,color=10
  ;
  if regname(ireg) eq nhtit then begin
  plot,timey,tstemp,thick=2,$
    /xstyle,xrange=[1850,2000],xtitle='Year',$
    ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
    title=regname(ireg)
  oplot,yrmxd,calmxd,color=10,thick=3
  oplot,yrmxd,sepmxd,color=11,thick=3
  oplot,!x.crange,[0.,0.],linestyle=1
  endif
endif
  ;
endfor
;
; Now output the calibrated series
;
if dosave gt 0 then begin
  tempregts=regtemp
  tempnyr=nyr
  temptimey=timey
  nyr=nfull
  timey=yrmxd
  if doclassic then fnadd='classic' else fnadd=''
  fn='bandtemp_'+fnadd+'calibrated.idlsave'
  print,'Creating: '+fn
  save,filename=fn,$
    nyr,nreg,calregts,regname,timey,tempregts,tempnyr,temptimey
  ;
  if dosave gt 1 then begin
    ;
    missval=-99.9
    ml=where(finite(calregts) eq 0,nmiss)
    if nmiss gt 0 then calregts(ml)=missval
    openw,1,'bandtemp_'+fnadd+'calibrated.dat'
    printf,1,'Calibrated regional reconstructions using tree-ring-density'
    printf,1,'Uses age-banded density records to maintain low frequencies'
    if doclassic then printf,1,'CALIBRATED USING CLASSICAL REGRESSION!!!'
    printf,1
    printf,1,'Calibrated to give estimates of regional-mean April-September'
    printf,1,'temperatures over land areas in some pre-defined regions'
    ;if donh eq 0 then $
    ;  printf,1,'(including region NH which is all land north of 20N).'
    printf,1
    printf,1,'Skill of the reconstruction generally deteriorates back in time,'
    printf,1,'due to fewer chronologies available.'
    printf,1
    printf,1,missval,'  = missing value'
    printf,1,nfull,'  = number of years'
    printf,1
    printf,1,['Year',regname(0:nreg-2)],format='(A4,10A7)'
    for iyr = 0 , nfull-1 do printf,1,yrmxd(iyr),calregts(iyr,0:nreg-2),$
      format='(I4,10F7.2)'
    printf,1
    printf,1,'April-September instrumental temperatures for each region'
    printf,1
    printf,1,['Year',regname(0:nreg-2)],format='(A4,10A7)'
    for iyr = 0 , tempnyr-1 do printf,1,temptimey(iyr),tempregts(iyr,0:nreg-2),$
      format='(I4,10F7.2)'
    close,1
    ;
    ; Also output the NH or ALL series
    ;
    openw,1,'bandtemp_'+fnadd+'calibrated_nh.dat'
    printf,1,'2  # no of columns'
    printf,1,'1  # column with data'
    printf,1,txt(nfull-1)+'  # number of time values'
    printf,1,'2  # no of header lines to skip'
    printf,1
    printf,1,'Year',regname[nreg-1],format='(A4,A7)'
    for iyr = 0 , nfull-2 do printf,1,yrmxd(iyr),calregts(iyr,nreg-1),$
      format='(I4,F7.2)'
    close,1
    ;
    ; Also output the series for MATLAB spectral analysis
    ;
    ;openw,1,'matlab_mxdband_all.dat'
    ;printf,1,calregts(*,nreg-1),format='(F7.2)'
    ;close,1
    ;
  endif
  ;
endif
;
end