;
; Creates regional timeseries (not composite regions) and adjusts them
; for changing sample size using a time-dependent rbar.
; It also computes EPS(t) for each regional timeseries, to use for later
; weighting.
;
plot,[0,1]
multi_plot,nrow=5,layout='large'
if !d.name eq 'X' then begin
  window,ysize=900
  initx
endif
;
ncid=ncdf_open('tree_dens_nh.nc')
ncdf_diminq,ncid,'time',dummy,ntime
ncdf_diminq,ncid,'station',dummy,nstat
ncdf_varget,ncid,'year',x
ncdf_varget,ncid,'density',density
ncdf_attget,ncid,'density','missing_value',valmiss
ncdf_varget,ncid,'fraction',weight
ncdf_close,ncid
;
misslist=where(density eq valmiss,nmiss)
density(misslist)=!values.f_nan
;
; Get regional tree lists
;
restore,filename='reglists.idlsave'
;
for i = 3 , nreg-3 do begin      ; ignore 0, 1 & 2 (ALL, NORTH & SOUTH)
                                 ; & ignore 8 & 9 (ARCTIC & EXTRA-ARCTIC)
  ;
  ; Get grand correlation matrix for this region
  ;
  restore,filename='rbartd_dens_'+regname(i)+'.idlsave'
  ;
  ; Now extract the chronologies used for this region
  ;
  dens=density(*,treelist(0:ntree(i)-1,i))
  wt=weight(*,treelist(0:ntree(i)-1,i))
  ;
  ; Make a time-dependent rbar next
  ;
  nt=n_elements(x)
  print,'Time-dependent rbar for ',nt,'  years'
  rbar_td=fltarr(nt)*!values.f_nan
  for it = 0 , nt-1 do begin
    print,it,format='($,I4)'
    mask = where( finite(dens(it,*)) , nsamp)
    if nsamp gt 0 then begin
      if nsamp eq 1 then begin
        rbar_td(it)=1.               ; if n=1 then n'=1 regardless of rbar
      endif else begin
        ; get reduced grandr
        redr=grandr(mask,*,0)
        redr=redr(*,mask)
        ; Total it and average it, but only do i<>j.
        totr=0.
        totn=0.
        for j = 1 , nsamp-1 do begin
          for k = 0 , j-1 do begin
            totr=totr+redr(j,k)        ; should have no missing values I hope
            totn=totn+1.
            if redr(j,k) eq 1. then message,'Ooops!'
          endfor
        endfor
        rbar_td(it)=totr/totn
      endelse
    endif
  endfor
  print
  ;
  ; Also compute rbar from the full matrix
  ;  
  ; Total it and average it, but only do i<>j.
  totr=0.
  totn=0.
  nsamp=n_elements(grandr(*,0,0))
  for j = 1 , nsamp-1 do begin
    for k = 0 , j-1 do begin
      rval=grandr(j,k,0)
      if finite(rval) then begin
        totr=totr+rval        ; should have no missing values I hope
        totn=totn+1.
        if rval eq 1. then message,'Ooops!'
      endif
    endfor
  endfor
  sigrbar=totr/totn
  print,sigrbar
  ;
  n=total(dens*0.+1.,2,/nan)       ; compute timeseries of number of chronols
  ncore=total(dens*0.+wt,2,/nan)   ; compute timeseries of summed core fraction
  totwdens=total(dens*wt,2,/nan)
  totwvals=total(float(finite(dens))*wt,2)
  wdens=totwdens/float(totwvals)               ; this is the weighted anomaly
  ;
  ; Compute the EPS timeseries
  ;
  denseps=sigrbar / (sigrbar + (1-sigrbar)/n)
  ;
  ; Now need to adjust for varying sample size
  ;
  xadj=mkeffective(wdens,n,rbar=rbar_td,neff=neff)
  wdens=xadj
  ;
  ; Now normalise w.r.t. 1881-1960
  ;
  mknormal,wdens,x,refperiod=[1881,1960],refmean=refmean,refsd=refsd
  ;
  ; Now plot them
  ;
  pause
  filter_cru,20,tsin=wdens,tslow=tslow,/nan 
  cpl_barts,x,wdens,title='Weighted mean normalised max density anomaly'+$ 
    ' for region: '+regname(i),xtitle='Year',/xstyle,$ 
    zeroline=tslow,yrange=[-8,4]
  oplot,x,tslow,thick=2 
  oplot,!x.crange,[0.,0.]
  oplot,!x.crange,[-2.,-2.],linestyle=1
  ;
  plot,x,n,ytitle='Sample size',/xstyle
  ;
  plot,x,rbar_td,/xstyle,ytitle='rbar'
  ;
  plot,x,neff,/xstyle,ytitle='Effective sample size'
  ;
  plot,x,denseps,/xstyle,ytitle='Expressed population signal'
  ;
  ; Now save the weighted mean timeseries, for further analysis
  ;
  densadj=wdens
  save,filename='densadj_'+regname(i)+'.idlsave',x,densadj,n,neff,ncore,$
    denseps
  ;
endfor
;
end