WMAP_5yr_tetbeebbeb_pixlike.cc

Go to the documentation of this file.

/****************************************************************
 *   This file is part of the C++ port of the Fortran likelihood
 *   code for the WMAP 5yr release provided by the WMAP team
 *   at http://lambda.gsfc.nasa.gov/ .
 *
 *   The code was ported by Georg Robbers for easier
 *   interfacing with cmbeasy (http://www.cmbeasy.org).
 *   Bugs in this port should be reported to the
 *   cmbeasy authors (bugs@cmbeasy.org).
 ****************************************************************/
#include "WMAP_5yr_tetbeebbeb_pixlike.h"
#include "WMAP_5yr_util.h"
#include "WMAP_5yr_options.h"

#include "read_archive_map.h"
#include "read_fits.h"

#include "gsl/gsl_linalg.h"
#include "gsl/gsl_blas.h"
#include "gsl/gsl_cblas.h"

#include <fstream>
#include <iostream>

// ===========================================================================
namespace  wmap_tetbeebbeb_lowl
{

// This code calculates the likelihood function of Q/U
// polarization for TE/TB/EE/BB/EB signal plus noise.
//
// Use combined 3yr data.
// - Combined Ninv and weighted maps have been generated by
// maps/weight_degrade_coadd_p2.f90
//
// - NinvPlNinv has been generated by
// ninv/generate_ninvplninv_r3_p2.f90
//
// E. Komatsu, October 4, 2005
// Incorporated into likelihood code October 7th RB.
// BB has been added, November 6th EK.
//
// - Implementation of David's likelihood
// Ek, December 16, 2005
//
// << MODIFICATION HISTORY AFTER RELEASE on March 16, 2006 >>
//
// E. Komatsu, June 18, 2006
// -- Use alm^tt from the released 3yr ILC map (wmap_ilc_3yr_v2.fits)
// -- Marginalization over the polarized foreground
//
// E. Komatsu, February 3, 2007
// -- 5-yr version
// -- Use double-precision LAPACK
//
// E. Komatsu, December 25, 2007
// -- TB & EB have been added as optional
// -- Foreground marginalization done to N^{-1} outside of this code
//
// ===========================================================================

int  nsmax,nlmax,mp;
long int np;
#ifdef OPTIMIZE
real8_1d xxx, yyy;
real4_2d ninvplninv2;
real4_3d ninvplninv3;
real8_3d ninvplninv2_double;
real8_3d ninvplninv3_double;
#else
real4_3d ninvplninv1;
real4_3d ninvplninv2;
real8_3d ninvplninv1_double;
real8_3d ninvplninv2_double;
#endif

Array1D<int>  ngood;
Matrix<REAL_8>  Dp0;
Array1D<REAL_8>  m_r3,w_r3,p_r3,f_r3,zzz;
Matrix<COMPLEX>  alm_tt,NinvYe, NinvYb;
Array1D<double>  wl;
//static const REAL sig_temp = 0.01;
static const REAL sig_temp = 0.0; // FG marginalization done outside of this code



//===========================================
void  tetbeebbeb_lowl_like_setup()
{
  using namespace std;
  //===========================================
  using namespace wmap_util;
  string WMAP_data_dir=WMAP_OPTIONS::self()->WMAP_data_dir();

  //string(LEN=2)  rlz;
  static const int ires = 3;
  real_1d  T,N,Mask_R3;
  //string(len=3)  da(10),sband;
  int  ReadStatus;
  string qfile,ufile,maskfile,filename[13], tetbeebbebdir;

  int  nsmax,nlmax,l,ip,jp,m,stat,i,j;
  unsigned int k;
  real_2d NinvQUr3;

#ifdef TIMING
  wmap_timing_start( "tetbeebbeb_lowl_like_setup" );
#endif

  //------------------------------
  // Set initial parameters
  //-----------------------------

  //X   da=(/\"K1\",\"Ka\",\"Q1\",\"Q2\",\"V1\",\"V2\",\"W1\",\"W2\",\"W3\",\"W4\"/)
  //X   da(2)=\"Ka1\";


  nsmax = (int)pow( 2.,ires );
  nlmax = 3*nsmax-1; // using Nyquist sampling rather than 2*nsmax
  np = 12*nsmax*nsmax;

  WMAP_OPTIONS* options = WMAP_OPTIONS::self();

  if ( options->using_gibbs_pol_cleaning() )  {
    tetbeebbebdir = WMAP_data_dir+"lowlP/gibbs/";
    filename[0]=tetbeebbebdir+"gibbs_masked_ee_ninvplninv_qu_r3_5yr.fits";
    filename[1]=tetbeebbebdir+"gibbs_masked_bb_ninvplninv_qu_r3_5yr.fits";
    filename[10]=tetbeebbebdir+"tbeb/gibbs_masked_eb_ninvplninv_qu_r3_5yr.fits";
    filename[11]=tetbeebbebdir+"tbeb/gibbs_masked_be_ninvplninv_qu_r3_5yr.fits";

    filename[2]=tetbeebbebdir+"gibbs_masked_ninv_qu_r3_5yr.fits";
    filename[3]=tetbeebbebdir+"gibbs_wt_r3_5yr.map_q";
    filename[4]=tetbeebbebdir+"gibbs_wt_r3_5yr.map_u";
    filename[6]=tetbeebbebdir+"gibbs_masked_ninvy_e_qu_r3_5yr.fits";
    filename[12]=tetbeebbebdir+"tbeb/gibbs_masked_ninvy_b_qu_r3_5yr.fits";
  } else {
    tetbeebbebdir = WMAP_data_dir+"lowlP/std/";
    filename[0]=tetbeebbebdir+"masked_ee_ninvplninv_qu_r3_corrected_5yr.KaQV.fits";
    filename[1]=tetbeebbebdir+"masked_bb_ninvplninv_qu_r3_corrected_5yr.KaQV.fits";
    filename[10]=tetbeebbebdir+"tbeb/masked_eb_ninvplninv_qu_r3_corrected_5yr.KaQV.fits";
    filename[11]=tetbeebbebdir+"tbeb/masked_be_ninvplninv_qu_r3_corrected_5yr.KaQV.fits";

    filename[2]=tetbeebbebdir+"masked_ninv_qu_r3_corrected_5yr.KaQV.fits";
    filename[3]=tetbeebbebdir+"wt_r3_5yr.KaQV.map_q";
    filename[4]=tetbeebbebdir+"wt_r3_5yr.KaQV.map_u";
    filename[6]=tetbeebbebdir+"masked_ninvy_e_qu_r3_corrected_5yr.KaQV.fits";
    filename[12]=tetbeebbebdir+"tbeb/masked_ninvy_b_qu_r3_corrected_5yr.KaQV.fits";
  }
  filename[5]=WMAP_data_dir+"lowlP/alm_tt_fs_r9_ilc_nopixwin_5yr.dat";
  filename[9]=WMAP_data_dir+"healpix_data/pixel_window_n0008.txt";

  tetbeebbebdir = WMAP_data_dir+"lowlP/std/";

    //------------------------------
    // read in res3 mask
    //------------------------------

  N.resize( Range( 0, np-1 ) );

  maskfile=WMAP_data_dir+"lowlP/mask_r3_p06_jarosik.fits";
  ifstream testStream( maskfile.c_str() );
  if (!testStream.good())
  {
    cout << "tetbeebbeb maskfile not found " << maskfile << endl;
  }
  Mask_R3.resize(np);
  T.resize(np);
  long int dum;
  Read_Archive_Map(maskfile, T, Mask_R3,dum,ReadStatus);
  if (ReadStatus!=0)
  {
    cout << "unable to read in mask" << maskfile << endl;
    exit( -1 );
  }

  ngood.resize( np );
  mp = 0;
  for (ip=0; ip<= np-1; ++ip)
  {
    if (Mask_R3(ip)!=0)
    {
      ngood(mp) = ip;
      mp = mp + 1;
    }
  }

  //------------------------------
  // read in N^{-1}P_l N^{-1} at res3
  //------------------------------

  string errorString;
  try
  {

#ifdef OPTIMIZE
    errorString="ninvplninv2";
    ninvplninv2.resize(Range(1,mp*(2*mp+1)), Range(1,4*(nlmax-1)));

    errorString="ninvplninv3";

   Read_FITS_REAL_3D (filename[0], ninvplninv3, stat);
   if (stat  !=  0)  {
     cout <<  "Error " << stat << " while reading " << filename[0] << endl ;
     exit(-1);
   }



   k = 1;
   for (i = 0; i <= 2*mp-1; ++i ) {
     for (j = i; j <= 2*mp-1; ++j ) {
       ip = ngood(i%mp) + np*(i/mp);
       jp = ngood(j%mp) + np*(j/mp);

       for (int l=1; l <=nlmax-1; ++l) ninvplninv2(k,l) = ninvplninv3(l+1,ip,jp);
       k = k + 1;
     }
   }

  Read_FITS_REAL_3D (filename[1], ninvplninv3, stat);
  if (stat  !=  0)  {
    cout << "Error " << stat << " while reading " << filename[1] << endl;
    exit(-1);
  }

  k = 1;
  for (i = 0; i <= 2*mp-1; ++i ) {
      for (j = i; j <= 2*mp-1; ++j ) {
          ip = ngood((i%mp)) + np*(i/mp);
          jp = ngood((j%mp)) + np*(j/mp);

          for (int l=nlmax; l <=2*(nlmax-1); ++l) ninvplninv2(k,l) = ninvplninv3(l-nlmax+2,ip,jp);
          k = k + 1;
      }
  }

  Read_FITS_REAL_3D (filename[10], ninvplninv3, stat);
  if (stat  !=  0)  {
    cout << "Error " << stat << " while reading " << filename[10] << endl;
    exit(-1);
  }

  k = 1;
  for (i = 0; i <= 2*mp-1; ++i ) {
      for (j = i; j <= 2*mp-1; ++j ) {
          ip = ngood((i%mp)) + np*(i/mp);
          jp = ngood((j%mp)) + np*(j/mp);

          for (int l=2*nlmax-1; l <=3*(nlmax-1); ++l) ninvplninv2(k,l) = ninvplninv3(l-(2*nlmax-1)+2,ip,jp);
          k = k + 1;
      }
  }

  Read_FITS_REAL_3D (filename[11], ninvplninv3, stat);
  if (stat  !=  0)  {
    cout << "Error " << stat << " while reading " << filename[1] << endl;
    exit(-1);
  }

  k = 1;
  for (i = 0; i <= 2*mp-1; ++i ) {
      for (j = i; j <= 2*mp-1; ++j ) {
          ip = ngood((i%mp)) + np*(i/mp);
          jp = ngood((j%mp)) + np*(j/mp);

          for (int l=3*nlmax-2; l <=4*(nlmax-1); ++l) ninvplninv2(k,l) = ninvplninv3(l-(3*nlmax-2)+2,ip,jp);
          k = k + 1;
      }
  }

  ninvplninv3.free();

  xxx.resize(Range(1,4*(nlmax-1)));
  yyy.resize(Range(1,mp*(2*mp+1)));

#else

    errorString = "ninvplninv";

    Read_FITS_REAL_3D (filename[0], ninvplninv1, stat);
    if (stat  !=  0)
    {
      cout <<  "Error " << stat << " while reading " << filename[0] << endl ;
      exit(-1);
    }

    Read_FITS_REAL_3D(filename[1], ninvplninv2, stat);
    if (stat  !=  0)
    {
      cout <<  "Error " << stat << " while reading " << filename[1] << endl;
      exit( -1 );
    }

    Read_FITS_REAL_3D(filename[10], ninvplninv3, stat);
    if (stat  !=  0)
    {
      cout <<  "Error " << stat << " while reading " << filename[10] << endl;
      exit( -1 );
    }

    Read_FITS_REAL_3D(filename[11], ninvplninv4, stat);
    if (stat  !=  0)
    {
      cout <<  "Error " << stat << " while reading " << filename[11] << endl;
      exit( -1 );
    }
#endif

    //------------------------------
    // read in N^{-1} at res3
    //------------------------------
    errorString = "NinvQUr3";

    Read_FITS_REAL_2D (filename[2], NinvQUr3, stat);
    if (stat  !=  0)
    {
      cout <<  "Error " << stat << " while reading " << filename[2] << endl ;
      exit(-1);
    }

    //------------------------------
    // read in maps at res3
    //------------------------------

    errorString = "w_r3, m_r3, p_r3 or zzz";
    w_r3.resize(2*mp);
    m_r3.resize(2*mp);
    p_r3.resize(2*mp);
    zzz.resize(2*mp);
    qfile=filename[3];
    Read_Archive_Map(qfile,T,N,np,ReadStatus);
    if (ReadStatus!=0)
    {
      cout << "unable to read in " << qfile << endl;
      exit( -1 );
    }

    for (ip=0; ip<= mp-1; ++ip)
    {
      w_r3(ip) = T(ngood(ip))*Mask_R3(ngood(ip));
    }

    ufile=filename[4];
    T.resize(np);
    N.resize(np);
    Read_Archive_Map(ufile,T,N,np,ReadStatus);
    if (ReadStatus!=0)
    {
      cout << "unable to read in " << ufile << endl;
      exit( -1 );
    }

    for (ip=0; ip<= mp-1; ++ip)
    {
      w_r3(mp+ip) = T(ngood(ip))*Mask_R3(ngood(ip));
    }

    Dp0.resize(2*mp,2*mp);

    for (ip=0; ip<= mp-1; ++ip)
    {
      for (jp=0; jp<= mp-1; ++jp)
      {
        Dp0(ip,jp) = NinvQUr3(ngood(ip),ngood(jp));
        Dp0(ip,mp+jp) = NinvQUr3(ngood(ip),np+ngood(jp));
        Dp0(mp+ip,jp) = NinvQUr3(np+ngood(ip),ngood(jp));
        Dp0(mp+ip,mp+jp) = NinvQUr3(np+ngood(ip),np+ngood(jp));
      }
    }

  //------------------------------
  // read in FG templates at res3 [NO LONGER USED: 12/25/07, EK]
  //------------------------------

  errorString="f_r3";
  f_r3.resize(Range(0,2*mp-1));

  f_r3 = 0.0;
  if (sig_temp != 0) {
    cout <<  "Using foreground templates!!!" << endl;
    cout <<  "----------------------------------------------------------------------" << endl;
    exit(-1);

    qfile=filename[7];
    Read_Archive_Map(qfile,T,N,np,ReadStatus);
    if (ReadStatus!=0) {
      cout << "unable to read in " << qfile << endl;
      exit(-1);
    }

    for (ip=0; ip<= mp-1; ++ip) {
      f_r3(ip) = T(ngood(ip))*Mask_R3(ngood(ip));
    }

    ufile=filename[8];
    Read_Archive_Map(ufile,T,N,np,ReadStatus);
    if (ReadStatus!=0)  {
      cout << "unable to read in " << ufile << endl;
      exit(-1);
    }

    for (ip=0; ip<= mp-1; ++ip) {
      f_r3(mp+ip) = T(ngood(ip))*Mask_R3(ngood(ip));
    }
  }

  //------------------------------
  // read in alm_tt
  //------------------------------
  errorString = "alm_tt";
  alm_tt.resize( nlmax+1, nlmax+1 );
  ifstream lun( filename[5].c_str() );
  for (l=0; l<= nlmax; ++l)
  {
    for (m=0; m<= l; ++m)
    {
      lun >> alm_tt(l,m);
    }
  }
  lun.close();

  //------------------------------
  // read in Ninv Y
  //------------------------------
  errorString = "NinvYe";
  NinvYe.resize(1536,300);

  Read_FITS_Complex_2D_LM (filename[6], NinvYe, stat, 0 /* IndFmt, was 2*/);
  if (stat  !=  0)
  {
    cout <<  "Error " << stat << " while reading " << filename[6] << endl;
    exit( -1 );
  }
  NinvYe.reindexSelf(0, 1); // instead of the IndFmt=2 above

  errorString = "NinvYb";
  NinvYb.resize(1536,300);

  Read_FITS_Complex_2D_LM (filename[12], NinvYb, stat, 0 /* IndFmt, was 2*/);
  if (stat  !=  0)
  {
    cout <<  "Error " << stat << " while reading " << filename[12] << endl;
    exit( -1 );
  }
  NinvYb.reindexSelf(0, 1); // instead of the IndFmt=2 above

  //
  // read in pixel window
  //
  wl.resize( Range(0, nlmax ) );

  ifstream wlfile(filename[9].c_str());
  for ( int i = 0; i <= nlmax; ++i )
      wlfile >> wl(i);
  }
  catch ( std::bad_alloc )
  {
    cout << "Memory Allocation error for: " << errorString << endl;
    exit( -1 );
  }

#ifdef TIMING
  wmap_timing_end();
#endif

} //  tetbeebbeb_lowl_like_setup


unsigned int tetbeebbeb_pixlike_dof()
{
  return 2*mp;
}

//===========================================
void  tetbeebbeb_lowl_likelihood(int nlmaxin, real8_1d& Clttin, real8_1d& Cltein, real8_1d& Cltbin,
                                 real8_1d& Cleein, real8_1d& Clbbin, real8_1d& Clebin,
                                 REAL_8& chisq_r3, REAL_8& lndet)
{

  //=======================================
  using namespace wmap_util;
  using namespace std;

  const double teeebb_pixlike_lndet_offset= WMAP_OPTIONS::self()->teeebb_pixlike_lndet_offset();
#ifndef OPTIMIZE
  Matrix<REAL_8> NinvSNinv;
  Matrix<REAL_8> CQU;
#endif
  Matrix<REAL_8> Dp;
  Array1D<REAL_8> Clee,Clbb,Cltt,Clte, Cltb, Cleb;
  int  ip,jp,nlmax,info,l,m,i,k;
  REAL Omega_pix;

#ifdef TIMING
  wmap_timing_start( "tetbeebbeb_lowl_likelihood" );
#endif

  nlmax = 23;
  Omega_pix = 3.14159/(3.*8.*8.);
  if (nlmaxin != nlmax)
  {
    cout << "need nlmax" << nlmax << ", in tetbeebbeb likelihood, currently" << nlmaxin;
    exit(0);
  }

  Cltt.resize(Range(2,nlmax));
  Clte.resize(Range(2,nlmax));
  Clee.resize(Range(2,nlmax));
  Clbb.resize(Range(2,nlmax));
  Cltb.resize(Range(2,nlmax));
  Cleb.resize(Range(2,nlmax));
  for (l=2; l<= nlmax; ++l)
  {
    double wl2 = wl(l)*wl(l);
    Cltt(l) = Clttin(l)/double(l*(l+1))*(2.*3.14159)*1.e-6*wl2;
    Clte(l) = Cltein(l)/double(l*(l+1))*(2.*3.14159)*1.e-6*wl2;
    Clee(l) = Cleein(l)/double(l*(l+1))*(2.*3.14159)*1.e-6*wl2;
    Clbb(l) = Clbbin(l)/double(l*(l+1))*(2.*3.14159)*1.e-6*wl2;
    Cltb(l) = Cltbin(l)/double(l*(l+1))*(2.*3.14159)*1.e-6*wl2;
    Cleb(l) = Clebin(l)/double(l*(l+1))*(2.*3.14159)*1.e-6*wl2;
  }

  //------------------------------
  // compute [N^{-1}S N^{-1} + N^{-1}]^{-1}
  //------------------------------


  lndet = 0.;

#ifdef OPTIMIZE
  Dp.resize(Range(0,2*mp-1),Range(0,2*mp-1));


  for (l = 2; l <= nlmax; ++l ) {
    xxx(l-1) = Clee(l)-Clte(l)*Clte(l)/Cltt(l);
    xxx(l+nlmax-2) = Clbb(l)-Cltb(l)*Cltb(l)/Cltt(l);
    xxx(l+2*nlmax-3) = Cleb(l)-Clte(l)*Cltb(l)/Cltt(l);
    xxx(l+3*nlmax-4) = Cleb(l)-Cltb(l)*Clte(l)/Cltt(l);
  }

//X call DGEMV( 'N', mp*(2*mp+1), 2*(nlmax-1), 1.d0, DBLE(ninvplninv2), &
//X              mp*(2*mp+1), xxx, 1, 0.d0, yyy, 1 )
  yyy = 0.;
  for (unsigned int i = 0; i < ninvplninv2.rows(); ++i)
      for (unsigned int j = 0; j < ninvplninv2.cols(); ++j)
          *(yyy.data()+i) += double(ninvplninv2(i+1,j+1))*(*(xxx.data()+j));

  k = 1;
  for (ip = 0; ip <= 2*mp-1; ++ip ) {
      for (jp = ip; jp <= 2*mp-1; ++jp ) {
          Dp(ip,jp) = Dp0(ip,jp) + yyy(k);
              k = k + 1;
      }
  }

#else
  //
  // fill only the upper triangular part of N^{-1}SN^{-1}.
  // add the foreground error term.

  NinvSNinv.resize( 2*mp, 2*mp);
  NinvSNinv = 0.;

  double* eetefact = new double[nlmax+1];
  double* bbtbfact = new double[nlmax+1];
  double* ebtefact = new double[nlmax+1];
  double* ebtbfact = new double[nlmax+1];
  for ( int i = 2; i <= nlmax; ++i ) {
    eetefact[i] = (Clee(i)- Clte(i)*Clte(i)/Cltt(i));
    bbtbfact[i] = (Clbb(i)- Cltb(i)*Cltb(i)/Cltt(i));
    ebtefact[i] = (Cleb(i)- Clte(i)*Cltb(i)/Cltt(i));
    ebtbfact[i] = (Cleb(i)- Cltb(i)*Clte(i)/Cltt(i));
  }
  gsl_vector_view eetevw = gsl_vector_view_array( eetefact+2, nlmax-1 );
  gsl_vector_ *eetev = &eetevw.vector;
  gsl_vector_view bbtbvw = gsl_vector_view_array( bbtbfact+2, nlmax-1 );
  gsl_vector_ *bbtbv = &bbtbvw.vector;
  gsl_vector_view ebtevw = gsl_vector_view_array( ebtefact+2, nlmax-1 );
  gsl_vector_ *ebtev = &ebtevw.vector;
  gsl_vector_view ebtbvw = gsl_vector_view_array( ebtbfact+2, nlmax-1 );
  gsl_vector_ *ebtbv = &ebtbvw.vector;

  //gsl_vector_view bbvw = gsl_vector_view_array( Clbb.data(), Clbb.size() );
  //gsl_vector * bbv = &bbvw.vector;

  gsl_vector_view v1w, v2w, v3w, v4w;
  gsl_vector      *v1, *v2, *v3, *v4;



  for (ip=0; ip<= mp-1; ++ip)
  {
    double ddot1, ddot2, ddot3, ddot4;
    for (jp=ip; jp<= mp-1; ++jp)
    {
      v1w = gsl_vector_view_array_with_stride ( &ninvplninv1_double(2,ngood( ip ), ngood( jp )),
                                                    ninvplninv1_double.stride(firstDim), nlmax-1 );
      v2w = gsl_vector_view_array_with_stride ( &ninvplninv2_double(2,ngood( ip ), ngood( jp ) ),
                                                    ninvplninv2_double.stride(firstDim), nlmax-1 );
      v3w = gsl_vector_view_array_with_stride ( &ninvplninv3_double(2,ngood( ip ), ngood( jp ) ),
                                                    ninvplninv3_double.stride(firstDim), nlmax-1 );
      v4w = gsl_vector_view_array_with_stride ( &ninvplninv4_double(2,ngood( ip ), ngood( jp ) ),
                                                    ninvplninv4_double.stride(firstDim), nlmax-1 );
      v1 = &v1w.vector;
      v2 = &v2w.vector;
      v3 = &v3w.vector;
      v4 = &v4w.vector;
      gsl_blas_ddot( eetevw, v1, &ddot1 );
      gsl_blas_ddot( bbtbvw, v2, &ddot2 );
      gsl_blas_ddot( ebtevw, v3, &ddot3 );
      gsl_blas_ddot( ebtbvw, v4, &ddot4 );
      NinvSNinv(ip,jp) = ddot1+ddot2+ddot3+ddot4;

      v1w = gsl_vector_view_array_with_stride ( &ninvplninv1_double(2,ngood( ip ),( int )np+ngood( jp )),
                                                    ninvplninv1_double.stride(firstDim), nlmax-1 );
      v2w = gsl_vector_view_array_with_stride ( &ninvplninv2_double( 2,ngood( ip ),( int )np+ ngood( jp ) ),
                                                    ninvplninv2_double.stride(firstDim), nlmax-1 );
      v3w = gsl_vector_view_array_with_stride ( &ninvplninv3_double( 2,ngood( ip ),( int )np+ ngood( jp ) ),
                                                    ninvplninv3_double.stride(firstDim), nlmax-1 );
      v4w = gsl_vector_view_array_with_stride ( &ninvplninv4_double( 2,ngood( ip ),( int )np+ ngood( jp ) ),
                                                    ninvplninv4_double.stride(firstDim), nlmax-1 );
      v1 = &v1w.vector;
      v2 = &v2w.vector;
      v3 = &v3w.vector;
      v4 = &v4w.vector;
      gsl_blas_ddot( eetevw, v1, &ddot1 );
      gsl_blas_ddot( bbtbvw, v2, &ddot2 );
      gsl_blas_ddot( ebtevw, v3, &ddot3 );
      gsl_blas_ddot( ebtbvw, v4, &ddot4 );
      NinvSNinv(ip,mp+jp) = ddot1+ddot2+ddot3+ddot4;

      v1w=gsl_vector_view_array_with_stride(&ninvplninv1_double(2,(int)np+ngood(ip),(int)np+ngood(jp)),
                                            ninvplninv1_double.stride(firstDim), nlmax-1 );
      v2w=gsl_vector_view_array_with_stride(&ninvplninv2_double(2,(int)np+ngood(ip),(int)np+ngood(jp)),
                                            ninvplninv2_double.stride(firstDim), nlmax-1 );
      v3w=gsl_vector_view_array_with_stride(&ninvplninv3(2,(int)np+ngood(ip),(int)np+ngood(jp)),
                                            ninvplninv3.stride(firstDim), nlmax-1 );
      v4w=gsl_vector_view_array_with_stride(&ninvplninv4(2,(int)np+ngood(ip),(int)np+ngood(jp)),
                                            ninvplninv4.stride(firstDim), nlmax-1 );
      v1 = &v1w.vector;
      v2 = &v2w.vector;
      v3 = &v3w.vector;
      v4 = &v4w.vector;
      gsl_blas_ddot( eetevw, v1, &ddot1 );
      gsl_blas_ddot( bbtbvw, v2, &ddot2 );
      gsl_blas_ddot( ebtevw, v3, &ddot3 );
      gsl_blas_ddot( ebtbvw, v4, &ddot4 );
      NinvSNinv(mp+ip,mp+jp) = ddot1+ddot2+ddot3+ddot4;
    }

    for (jp=0; jp<= ip-1; ++jp)
    {
      v1w = gsl_vector_view_array_with_stride ( &ninvplninv1_double(2,ngood( ip ),( int )np+ngood( jp )),
                                                 ninvplninv1_double.stride(firstDim), nlmax-1 );
      v2w = gsl_vector_view_array_with_stride ( &ninvplninv2_double(2,ngood( ip ),( int )np+ngood( jp )),
                                                ninvplninv2_double.stride(firstDim), nlmax-1 );
      v3w = gsl_vector_view_array_with_stride ( &ninvplninv3_double( 2,ngood( ip ),( int )np+ ngood( jp ) ),
                                                ninvplninv3_double.stride(firstDim), nlmax-1 );
      v4w = gsl_vector_view_array_with_stride ( &ninvplninv4_double( 2,ngood( ip ),( int )np+ ngood( jp ) ),
                                                ninvplninv4_double.stride(firstDim), nlmax-1 );

      v1 = &v1w.vector;
      v2 = &v2w.vector;
      v3 = &v3w.vector;
      v4 = &v4w.vector;
      gsl_blas_ddot( eetevw, v1, &ddot1 );
      gsl_blas_ddot( bbtbvw, v2, &ddot2 );
      gsl_blas_ddot( ebtevw, v3, &ddot3 );
      gsl_blas_ddot( ebtbvw, v4, &ddot4 );
      NinvSNinv(ip,mp+jp) += ddot1+ddot2+ddot3+ddot4;
    }
  }
  delete[] eetefact;
  delete[] bbtbfact;
  delete[] ebtefact;
  delete[] ebtbfact;

  Dp = Dp0 +NinvSNinv ;

  NinvSNinv.free();
#endif

#ifdef TIMING
  wmap_timing_checkpoint( "finished Dp" );
#endif

  gsl_matrix_view DpView = gsl_matrix_view_array( Dp.data(), 2*mp, 2*mp );
  gsl_matrix* gslDp = &DpView.matrix;



  info = gsl_matrix_transpose( gslDp ); // gsl decomposes into a lower triang. matrix
  info = gsl_linalg_cholesky_decomp( gslDp );

  if (info != 0)
  {
    //wmap_likelihood_error( "tetbeebbeb: bad dpotrf", info );
    //chisq_r3 = 0.0;
    //lndet = 0.0;
    chisq_r3 = 1e10;
    lndet = 0.0;
    cout <<  "dpotrf failed, info = " << info << endl ;
#ifdef TIMING
    wmap_timing_end();
#endif
     ofstream lun;
     lun.open("debug.txt");
     for (l = 2; l <= 23; ++l ) {
       lun << l << " " << Cltt(l) << " " << Clte(l) << " " << Cltb(l) << " "
           << Clee(l) << " " << Clbb(l) << " " << Cleb(l) << endl;
       cout << Cltt(l) << " " << Clte(l) << " " << Cltb(l) << " " << Clee(l)
            << " " << Clbb(l) << " " << Cleb(l) << endl;
     }
     lun.close();
     wmap_likelihood_error( "tetbeebbeb: bad dpotrf", info );
    // exit(-1);
    return;
  }

#ifdef TIMING
  wmap_timing_checkpoint( "finished dpotrf" );
#endif

  for (ip=0; ip< 2*mp; ++ip)
  {
    lndet = lndet + 2.*log( gsl_matrix_get( gslDp, ip, ip ) );
  }


#ifndef OPTIMIZE
  cholesky_invert( gslDp);
#endif

#ifdef TIMING
  wmap_timing_checkpoint( "finished dpotri" );
#endif

  if (info != 0)
  {
    wmap_likelihood_error( "tetbeebbeb: bad dpotri", info );
    chisq_r3 = 0.0;
    lndet = 0.0;
    return;

  }

#ifndef OPTIMIZE
  gsl_matrix *  gslCQU = gslDp ;
  Dp.free();
#endif


  //-----------------------------------
  // calculate the predicted QU at res3
  //-----------------------------------
  p_r3 = 0.;
  for (ip=0; ip<= mp-1; ++ip)
  {
    i = 3;
    for (l=2; l<= 23; ++l)
    {
      i = i+1;
      double tmpe = ( alm_tt(l,0)*NinvYe(ngood(ip),i) ).real();
      double tmpb = ( alm_tt(l,0)*NinvYb(ngood(ip),i) ).real();
      p_r3(ip)    = p_r3(ip)
                   + Clte(l)/Cltt(l)*wl(l)*tmpe
                   + Cltb(l)/Cltt(l)*wl(l)*tmpb;
      tmpe =( alm_tt(l,0).real()*NinvYe(ngood(ip)+np,i) ).real();
      tmpb =( alm_tt(l,0).real()*NinvYb(ngood(ip)+np,i) ).real();
      p_r3(ip+mp) = p_r3(ip+mp)
                    + Clte(l)/Cltt(l)*wl(l)*tmpe
                    + Cltb(l)/Cltt(l)*wl(l)*tmpb;
      for (m=1; m<= l; ++m)
      {
        i = i+1;
        tmpe = conj(alm_tt(l,m)*NinvYe(ngood(ip),i)).real();
        double tmpke =  ( alm_tt(l,m)*NinvYe(ngood(ip),i) ).real();
        tmpb = conj(alm_tt(l,m)*NinvYb(ngood(ip),i)).real();
        double tmpkb =  ( alm_tt(l,m)*NinvYb(ngood(ip),i) ).real();
        p_r3(ip) = p_r3(ip)
                   + Clte(l)/Cltt(l)*wl(l)*( tmpke + tmpe )
                   + Cltb(l)/Cltt(l)*wl(l)*( tmpkb + tmpb );
        tmpe = conj(alm_tt(l,m)*NinvYe(ngood(ip)+np,i)).real();
        tmpke = (alm_tt(l,m)*NinvYe(ngood(ip)+np,i)).real();
        tmpb = conj(alm_tt(l,m)*NinvYb(ngood(ip)+np,i)).real();
        tmpkb = (alm_tt(l,m)*NinvYb(ngood(ip)+np,i)).real();
        p_r3(ip+mp) = p_r3(ip+mp)
                      + Clte(l)/Cltt(l)*wl(l)*(tmpke+tmpe)
                      + Cltb(l)/Cltt(l)*wl(l)*(tmpkb+tmpb);
      }
    }
  }

#ifdef OPTIMIZE
  m_r3 = w_r3-p_r3;
//X   call DPOTRS( "U", 2*mp, 1, Dp, 2*mp, m_r3, 2*mp, info );
  gsl_vector_view m_r3view = gsl_vector_view_array(m_r3.data(), m_r3.size());
  gsl_linalg_cholesky_svx(gslDp, &m_r3view.vector);
  //chisq_r3 = DDOT(2*mp,m_r3,1,w_r3-p_r3,1) 
  zzz = w_r3-p_r3;
//X   chisq_r3 = sum( m_r3*zzz );
  int istart = m_r3.firstIndex();
  int iend = istart + m_r3.size();
  chisq_r3 = 0;
  for (int i=istart; i < iend; ++i) chisq_r3+=m_r3(i)*zzz(i);
#else

  gsl_vector *x = gsl_vector_alloc( 2*mp );
  gsl_vector *y = gsl_vector_alloc( 2*mp );
  for ( i = 0; i < 2*mp; ++i )
  {
    gsl_vector_set(x, i, w_r3(i)-p_r3(i));
    gsl_vector_set(y, i, m_r3(i));
  }

  int status = gsl_blas_dsymv( CblasUpper, 1., gslCQU, x, 0., y );

  status = gsl_blas_ddot(x,y,&chisq_r3);

  CQU.free();
  gsl_vector_free(x);
  gsl_vector_free(y);
#endif

  chisq_r3 = chisq_r3/2.;
  lndet = (lndet - teeebb_pixlike_lndet_offset)/2.0;

#ifdef TIMING
  wmap_timing_end();
#endif

} //  TETBEEBBEB_LOWL_LIKELIHOOD


/* doing the cholesky inversion by hand along the lines of LAPACK */
void cholesky_invert( gsl_matrix* A)
{

  int N = A->size1-1;
  for ( int J = N; J >= 0;  --J )
  {
    gsl_matrix_set(A, J, J, 1. / gsl_matrix_get(A,  J, J ) );
    double AJJ = -gsl_matrix_get( A,  J, J );

    if (  J<N )
    {

      //Compute elements j+1:n of j-th column.
      gsl_matrix_view subA = gsl_matrix_submatrix( A, J+1, J, N-J, N-J );
      gsl_vector_view v = gsl_matrix_column ( &subA.matrix, 0 );
      subA = gsl_matrix_submatrix( A, J+1, J+1, N-J, N-J );

      gsl_blas_dtrmv ( CblasLower, CblasNoTrans, CblasNonUnit, &subA.matrix, &v.vector ) ;
      gsl_blas_dscal ( AJJ, &v.vector );
    }
  }


  //X  Compute the product L' * L.
  N = A->size1-1;

  double* matrixData = A->data;
  unsigned int tda = A->tda;

  for ( int i = 0; i <= N;  ++i )
    for ( int j = 0; j <= i;  ++j )
    {
      double res=0;
      for ( int k = i; k <= N;  ++k )
      {
        res += matrixData[ k*tda+i ] * matrixData[ k*tda+j ];
      }
      matrixData[j*tda+i ]=res;
    }
}


} // end namespace  WMAP_tetbeebbeb_lowl

---