; bayesCov.fits produced with .r bayes2fits.idl
; the first 2 columns are Nsample & Chi2 (removed during weight map
; processing)
a=mrdfits('bayesCov.fits',1)
names=tag_names(a)
nx=n_elements(names)
ny=n_elements(a)
d=fltarr(nx-2,ny)
for i=2, nx-1 do for j=0L, ny-1 do d[i-2,j]=a[j].(i)
names=names[2:nx-1]
nx=nx-2 ;; remove Nsample & lnhood columns

cor=correlate(d, /double)
writefits, 'corr.fits', cor

cov=correlate(d, /covariance, /double)
writefits, 'cov.fits', cov

if determ(cov) eq 0 then begin 
    print, "ERROR: null determinant" 
    tmp=cov
    for i=0, nx-1 do tmp[i,i]=0 ; keep only non diag terms
    s = sort(abs(tmp))
    s = reverse(s)
    for i=0, 10, 2 do begin 
        i_x = s[i] / nx
        i_y = s[i] mod nx
        print, format='(A," correlates with ",A," @ cov[",I3,",",I3,"] = ",E11.3)', names[i_x],names[i_y],i_x,i_y,cov[i_x,i_y]
    endfor


     
    stop 
endif
    
weight=la_invert(cov, /double) 
writefits, 'weight.fits', weight 

END