text:
/**
 * The driver function and some auxillaries
 * to compute the spectrum from an image and
 * information on the location and distortion of the spectra.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <limits.h>

#include <gsl/gsl_roots.h>
#include <gsl/gsl_vector.h>

#include "aXe_grism.h"
#include "aXe_utils.h"
#include "spce_sect.h"
#include "spce_is_in.h"
#include "spce_pathlength.h"

#define DEBUG_ME 0x10

/**
  fill in the path length field (xi in ap_pixel) from the abscissas (xs in
  ap_pixel) (helper function for spc_extract).

  @param func the spectrum trace to use for the transformation from xs to xi
  @param table a pointer to the table of aperture pixels, terminated with
  x=-1
*/
static int
transform_to_pathlen (const trace_func * const func, ap_pixel * const table)
{
  int i, len;
  ap_pixel *cur_p = table;
  gsl_vector *section_points;

  while (cur_p->p_x != -1)
    cur_p++;
  len = cur_p - table;
  if (len==0) return -1; // exit now if the table is empty

  section_points = gsl_vector_alloc (len);
  for (i = 0; i < len; i++)
    {
      gsl_vector_set (section_points, i, table[i].xs);
    }

  if (abscissa_to_pathlength (func, section_points))
    {
      gsl_vector_free (section_points);
      return -1;
    }

  for (i = 0; i < len; i++)
    {
      table[i].xi = gsl_vector_get (section_points, i);
    }

  gsl_vector_free (section_points);

  return 0;
}

/**
   creates an ap_pixel.

   @param x x coordinate of pixel relative to beam's reference point
   @param y y coordinate of pixel relative to beam's reference point
   @param px absolute x coordinate of the pixel
   @param py absolute y coordinate of the pixel
   @param ob the observation to get the pixels from
   @param sf pointer to sectionfun structure for this beam
   @param cur_ap a pointer to the next free aperture pixel
   @return a pointer to the next free aperture pixel after the new pixels
    have been added.

*/
static ap_pixel *
handle_one_pixel (const double x, const double y, const int px, const int py,
		  const observation * const obs, sectionfun * const sf,
		  const trace_func * const tracefun, ap_pixel * cur_ap)
{
  double res;
  double sect_y;
  double tmp;
  double phi_trace;

  if (find_section_point (sf, x, y, &res))
    {
      return cur_ap;	/* FIXME: Issue warning here */
    }
  /* set the extraction weight to 1. */
  cur_ap->weight = 1.;
  cur_ap->xs = res;
  phi_trace = atan (tracefun->deriv (cur_ap->xs, tracefun->data));
  cur_ap->dxs = phi_trace;
  cur_ap->ys = tracefun->func (cur_ap->xs, tracefun->data);
  sect_y = tracefun->func (res, tracefun->data);
  tmp = tracefun->func (x, tracefun->data);

  cur_ap->contam = -1.0;
  cur_ap->model = 0.0;
  cur_ap->p_x = px;
  cur_ap->p_y = py;
  cur_ap->x = x;
  cur_ap->y = y;
  cur_ap->dist =
    sqrt ((sect_y - y) * (sect_y - y) +
	  (cur_ap->xs - x) * (cur_ap->xs - x));

  if ( y < tmp ) {
    cur_ap->dist = cur_ap->dist * (-1.); /* If pixel is bwlow the trace, dist is neg. */
  }
  cur_ap->count = gsl_matrix_get (obs->grism, px, py);
  cur_ap->error = gsl_matrix_get (obs->pixerrs, px, py);
  if (obs->dq != NULL) cur_ap->dq = (long) gsl_matrix_get (obs->dq, px, py);
  else cur_ap->dq = 0;

  cur_ap++;
  return cur_ap;
}

/**
   Does some sanity checks on make_spc_table's input.

   @param ob the object to check
   @return 0 if everything is ok, -1 otherwise
*/
static int
sanitycheck (object * const ob)
{
  int i;

  for (i = 0; i < ob->nbeams; i++)
    {
      while (ob->beams[i].orient < 0)
	{
	  ob->beams[i].orient += M_PI;
	}
      while (ob->beams[i].orient > M_PI)
	{
	  ob->beams[i].orient -= M_PI;
	}
    }
  return 0;
}

#define MIN(x,y) (((x)<(y))?(x):(y))
#define MAX(x,y) (((x)>(y))?(x):(y))
/**
   computes a table of aperture pixels, i.e. of tuples containing the
   source coordinates, the distance to the spectrum trace, the path
   length along the trace, and the intensity.  The end of the table
   is marked with an ap_pixel->x==-1.  This routine can do subsampling,
   whereby each pixel is divided into n_sub*n_sub smaller pixels.
   If n_sub==1, a special (faster) handling is enabled.

   @param ob the object struct to examine
   @param beamorder the order of the spectrum to examine
   @param flags warning flags like in the warning field of the
     spectrum structure.  This has to be passed in initialized.
*/
ap_pixel *
make_spc_table (object * const ob, const int beamorder,
		int *const flags)
{
  int bb_x, bb_y, bb_w, bb_h;
  int x, y;
  beam *curbeam = ob->beams + beamorder;
  double dx, dy;
  ap_pixel *table, *cur_ap;
  is_in_descriptor iid;
  sectionfun sf;
  trace_func *tracefun = curbeam->spec_trace;

  quad_to_bbox (curbeam->corners, curbeam->bbox, curbeam->bbox + 1);
  bb_x = curbeam->bbox[0].x;
  bb_y = curbeam->bbox[0].y;
  bb_w = curbeam->bbox[1].x - curbeam->bbox[0].x + 1;
  bb_h = curbeam->bbox[1].y - curbeam->bbox[0].y + 1;

  if (sanitycheck (ob))
    {
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "Input data failed sanity check");
      return NULL;
    }

  if (fill_is_in_descriptor (&iid, curbeam->corners))
    return NULL;

  if (fill_in_sectionfun (&sf, curbeam->orient, curbeam))
    return NULL;

  if (!(table =
	malloc ((bb_w * bb_h + 1) * sizeof (ap_pixel))))
    return NULL;

  /* We have a little coordinate confusion here.  There are three
     systems:
     (a) the absolute one, rooted in (0,0) of the image.
     (b) the relative one, rooted in the reference point and used
     when evaluating the spectrum trace.  This one is used
     for finding the section points, etc.
     (c) one relative to the bounding box; x and y below live in this
     system and have to be converted into whatever system is
     required.  dx, dy are used for the conversion to (b)
  */
  dx = bb_x - curbeam->refpoint.x;
  dy = bb_y - curbeam->refpoint.y;
  cur_ap = table;
  for (y = 0; y <= bb_h; y++)
    {
      for (x = 0; x < bb_w; x++)
	{
	  if ((bb_x + x < 0) || (bb_y + y < 0)
	      || (bb_x + x >= (int)ob->grism_obs->grism->size1)
	      || (bb_y + y >= (int)ob->grism_obs->grism->size2))
	    {
	      *flags |= SPC_W_BORDER;
	      continue;
	    }

	  // check whether the trace description has
	  // an order higher than linear
	  if (curbeam->spec_trace->type > 1)
	    {
	      // new criteria based on the true trace distance
	      // which means the true distance from the section point
	      if (!tracedist_criteria(x + dx, y + dy, &sf, tracefun, curbeam->width+2.0))
		{
		  continue;
		}
	    }
	  else
	    {
	      // old criteria based on the box model
	      if (!is_in (bb_x + x, bb_y + y, &iid))
		{
		  continue;
		}
	    }

	  if (isnan
	      (gsl_matrix_get
	       (ob->grism_obs->grism, bb_x + x, bb_y + y)))
	    {
	      continue;
	    }
	  //	  if (ob->ID == 11 && x + bb_x == 59)
	    //	    fprintf(stdout, "xx: %i, yy: %i: %i\n", x + bb_x, y + bb_y, is_in (x + bb_x, y + bb_y, &iid));
	  cur_ap =
	    handle_one_pixel (x + dx, y + dy, x + bb_x,
			      y + bb_y, ob->grism_obs, &sf,
			      tracefun, cur_ap);

	}
    }

  cur_ap->p_x = -1;
  cur_ap->p_y = -1;
  cur_ap->count = -1;

  free_sectionfun (&sf);

  if (transform_to_pathlen (curbeam->spec_trace, table))
    {
      //   free (table);
      //   table = NULL;
      //   return NULL;
    }
  return table;
}

/**
   Similar to make_spc_table, it just works only
   on one spot given in the parameters.
*/
ap_pixel *
make_gps_table (object * const ob, const int beamorder,
		int *const flags, int xval, int yval)
{
  int bb_x, bb_y, bb_w, bb_h;
  int x, y;
  beam *curbeam = ob->beams + beamorder;
  double dx, dy;
  ap_pixel *table, *cur_ap;
  is_in_descriptor iid;
  sectionfun sf;
  trace_func *tracefun = curbeam->spec_trace;

  quad_to_bbox (curbeam->corners, curbeam->bbox, curbeam->bbox + 1);
  bb_x = curbeam->bbox[0].x;
  bb_y = curbeam->bbox[0].y;
  bb_w = curbeam->bbox[1].x - curbeam->bbox[0].x + 1;
  bb_h = curbeam->bbox[1].y - curbeam->bbox[0].y + 1;

  if (sanitycheck (ob))
    {
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "Input data failed sanity check");
      return NULL;
    }

  if (fill_is_in_descriptor (&iid, curbeam->corners))
    return NULL;

  if (fill_in_sectionfun (&sf, curbeam->orient, curbeam))
    return NULL;

  if (!
      (table =
       malloc (2 * sizeof (ap_pixel))))
    return NULL;

  /* We have a little coordinate confusion here.  There are three
     systems:
     (a) the absolute one, rooted in (0,0) of the image.
     (b) the relative one, rooted in the reference point and used
     when evaluating the spectrum trace.  This one is used
     for finding the section points, etc.
     (c) one relative to the bounding box; x and y below live in this
     system and have to be converted into whatever system is
     required.  dx, dy are used for the conversion to (b)
  */
  dx = bb_x - curbeam->refpoint.x;
  dy = bb_y - curbeam->refpoint.y;
  x  = xval - bb_x;
  y  = yval - bb_y;

  cur_ap = table;

  cur_ap =
    handle_one_pixel (x + dx, y + dy, x + bb_x,
		      y + bb_y, ob->grism_obs, &sf,
		      tracefun, cur_ap);

  cur_ap->p_x = -1;
  cur_ap->p_y = -1;
  cur_ap->count = -1;

  free_sectionfun (&sf);

  transform_to_pathlen (curbeam->spec_trace, table);

  return table;
}

void
print_ap_pixel_table (const ap_pixel * ap_p)
{
  printf ("# x y pathlen distance lambda count error\n");
  while (ap_p->p_x != -1)
    {
      printf ("%d %d %f %f %f %f %f %f %f %ld\n", ap_p->p_x, ap_p->p_y,
	      ap_p->x, ap_p->y, ap_p->xi, ap_p->dist, ap_p->lambda,
	      ap_p->count, ap_p->error, ap_p->dq);
      ap_p++;
    }
}

/*
int
tracedist_criteria(const double x, const double y, sectionfun * const sf,
		   const trace_func *tracefun, const double width)
{
  double x_sect, y_sect;
  double dist, max_dist;

  int ireturn=0;

  find_section_point (sf, x, y, &x_sect);
  y_sect = tracefun->func (x_sect, tracefun->data);

  dist =  sqrt((y_sect-y)*(y_sect-y)+(x_sect-x)*(x_sect-x));

  if (dist < width)
    ireturn=1;

  return ireturn;
}
*/

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#ifndef INITIAL_POINT_H
#define INITIAL_POINT_H

#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>

int pdip_initial_point_feasible_x( const gsl_matrix * A, const gsl_vector *b, gsl_vector *x );
int pdip_initial_point_feasible_s( const gsl_matrix * C, const gsl_vector *d, const gsl_vector *x, gsl_vector *s );
int pdip_initial_point_y( const gsl_matrix *Q, const gsl_vector *q, const gsl_matrix *A, const gsl_vector *x, gsl_vector *y );
int pdip_initial_point_z( gsl_vector *z );
int pdip_initial_point_strict_feasible( gsl_vector *x, gsl_vector *s );

#endif

ext: h

max_stars_repo_name: wmotte/toolkid

max_issues_repo_name: wmotte/toolkid

max_forks_repo_name: wmotte/toolkid


text:
#include "common.h"
#include "dMatrix.h"
#include "options.h"
#include "linefile.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_statistics_double.h>

struct rnaBinder
/* Name and correlation of rna binding proteins. */
{
    struct rnaBinder *next; /* Next in list. */
    char *psName;           /* Name of probe set. Something like
		 	     * E@NM_XXXXXX_a_at or GXXXXXXX_a_at.*/
    char *geneName;         /* Name of gene probed by probe
			     * set. Hopefully something readable. */
    char *pfamAcc;          /* Name of pfam accession. */
    char *pfamName;         /* Name of pfam domain. */
    double corr;           /* Correlation between rnaBinder and cluster of interest. */
};

struct clusterMember
/* Member of the cluster of interest. */
{
    struct clusterMember *next; /* Next in list. */
    char *geneId;               /* Gene id in cluster. */
    char *psName;               /* Name of the probe set. */
    char *desc;                 /* Description of the probe set. */
};

static struct optionSpec optionSpecs[] =
/* Our acceptable options to be called with. */
{
    {"help", OPTION_BOOLEAN},
    {"sjIndexFile", OPTION_STRING},
    {"psFile", OPTION_STRING},
    {"clusterFile", OPTION_STRING},
    {"rnaBindingFile", OPTION_STRING},
    {"outputFile", OPTION_STRING},
    {"antiCorrelation", OPTION_BOOLEAN},
    {NULL, 0}
};

static char *optionDescripts[] =
/* Description of our options for usage summary. */
{
    "Display this messge.",
    "File with splice junction indexes.",
    "File with probe set intensities.",
    "File with probe sets in cluster of interest.",
    "File with rna binding gene sets.",
    "File to output correlation results to.",
    "Sort by anti-correlation rather than correlation."
};

void usage()
/* Print usage and quit. */
{
int i=0;
warn("rnaBindingClusterCorr - Program to calculate the correlation\n"
     "between rna binding protein expression data and splice junctions\n"
     "co-regulated as determined by clustering splice junction indexes.\n"
     "options are:");
for(i=0; i<ArraySize(optionSpecs) -1; i++)
    fprintf(stderr, "  -%s -- %s\n", optionSpecs[i].name, optionDescripts[i]);
errAbort("\nusage:\n"
	 "   rnaBindingClusterCorr -sjIndexFile=input/sjIndexTest.tab -psFile=input/psIntenTest.tab \\ \n"
         "      -clusterFile=input/clusterTest.tab -rnaBindingFile=input/rnaBinderTest.tab \\ \n"
         "      -outputFile=output/test.out.tab\n");
}

double correlation(double *X, double *Y, int count)
/* Compute the correlation between X and Y
   correlation(X,Y) = cov(X,Y)/ squareRt(var(X)var(Y))
   page 332 Sheldon Ross "A First Course in Probability" 1998
*/
{
double varX = gsl_stats_variance(X, 1, count);
double varY = gsl_stats_variance(Y, 1, count);
double covXY = gsl_stats_covariance(X, 1, Y, 1, count);

double correlation = covXY / sqrt(varX *varY);
return correlation;
}

double calcDistanceFromCluster(struct rnaBinder *rb, struct clusterMember *cmList,
			       struct dMatrix *sjIndex, struct dMatrix *psInten)
/* Calculate the distance from the rnaBinder intensity measurement to
   the sjIndexes of the cluster members. If no intensity present use
   0 as it will fall in the middle of [-1,1]. */
{
double sum = 0;
int count = 0;
int sjIx = 0, gsIx = 0;
struct clusterMember *cm = NULL;
double corr = 0;
if(sjIndex->colCount != psInten->colCount)
    errAbort("Splice Junction and Intensity files must have same number of columns.");

/* Get the index of the gene set in the intensity file. */
gsIx = hashIntValDefault(psInten->nameIndex, rb->psName, -1);
if(gsIx == -1)
    {
/*     warn("Probe Set %s not found in intensitiy file."); */
    return 0;
    }
for(cm = cmList; cm != NULL; cm = cm->next)
    {
    /* For each member get the index in the splice junction file. */
    sjIx = hashIntValDefault(sjIndex->nameIndex, cm->psName, -1);
    if(sjIx == -1)
	errAbort("Probe Set %s not found in SJ index file.");
    corr = correlation(psInten->matrix[gsIx], sjIndex->matrix[sjIx], sjIndex->colCount);
    sum += corr;
    count++;
    }
if(count == 0)
    errAbort("No junctions in cluster.");
sum = sum / (double) count;
return sum;
}

int rnaBinderCmp(const void *va, const void *vb)
/* Compare to sort based correlation. */
{
const struct rnaBinder *a = *((struct rnaBinder **)va);
const struct rnaBinder *b = *((struct rnaBinder **)vb);
return a->corr < b->corr;
}

struct rnaBinder *loadRnaBinders()
/* Load the probe sets that encode genes thought to
   bind rnas. Expected order is probeSet, geneName, pfamAcc, pfamName */
{
struct rnaBinder *rbList = NULL, *rb = NULL;
char *words[4];
struct lineFile *lf = NULL;
char *inputFile = optionVal("rnaBindingFile", NULL);

assert(inputFile);
lf = lineFileOpen(inputFile, TRUE);
while(lineFileChopCharNext(lf, '\t', words, ArraySize(words)))
    {
    AllocVar(rb);
    rb->psName = cloneString(words[0]);
    rb->geneName = cloneString(words[1]);
    rb->pfamAcc = cloneString(words[2]);
    rb->pfamName = cloneString(words[3]);
    slAddHead(&rbList, rb);
    }
lineFileClose(&lf);
slReverse(&rbList);
return rbList;
}

struct clusterMember *loadClusterMembers()
/* Load the probe sets that are in our cluster of interest. */
{
struct clusterMember *cmList = NULL, *cm = NULL;
char *words[3];
struct lineFile *lf = NULL;
char *inputFile = optionVal("clusterFile", NULL);
int wordCount = 0;
assert(inputFile);
lf = lineFileOpen(inputFile, TRUE);
while((wordCount = lineFileChopCharNext(lf, '\t', words, ArraySize(words))) != 0)
    {
    AllocVar(cm);
    if(wordCount == 3)
	{
	cm->geneId = cloneString(words[0]);
	cm->psName = cloneString(words[1]);
	cm->desc = cloneString(words[2]);
	}
    else if(wordCount == 2)
	{
	cm->psName = cloneString(words[0]);
	cm->desc = cloneString(words[1]);
	}
    else
	errAbort("Got %d words at line %d", wordCount, lf->lineIx);
    slAddHead(&cmList, cm);
    }
lineFileClose(&lf);
slReverse(&cmList);
return cmList;
}

void reportRnaBinders(struct rnaBinder *rbList)
/* Report the rnaBinder distances. */
{
struct rnaBinder *rb = NULL;
char *fileName = optionVal("outputFile", NULL);
FILE *out = NULL;

if(fileName == NULL)
    errAbort("Must specify an outputFile");
out = mustOpen(fileName, "w");
for(rb = rbList; rb != NULL; rb = rb->next)
    {
    fprintf(out, "%s\t%s\t%s\t%s\t%.4f\n",
	    rb->psName, rb->geneName, rb->pfamAcc, rb->pfamName, rb->corr);
    }
carefulClose(&out);
}

void rnaBindingClusterCorr()
/* Top level function to calculate correlations. */
{
char *sjIndexFile = optionVal("sjIndexFile", NULL);
char *psFile = optionVal("psFile", NULL);
struct dMatrix *sjIndex = NULL, *psInten = NULL;
struct rnaBinder *rbList = NULL, *rb = NULL;
struct clusterMember *cmList = NULL, *cm = NULL;

warn("Loading Files.");
sjIndex = dMatrixLoad(sjIndexFile);
psInten = dMatrixLoad(psFile);
rbList = loadRnaBinders();
cmList = loadClusterMembers();

warn("Calculating Distance.");
/* Calculate distance. */
for(rb = rbList; rb != NULL; rb = rb->next)
    rb->corr = calcDistanceFromCluster(rb, cmList, sjIndex, psInten);

/* Sort by distance (correlati2n) */
slSort(&rbList, rnaBinderCmp);

/* If we want anti-correlation reverse. */
if(optionExists("antiCorrelation"))
    slReverse(&rbList);

warn("Writing out results.");
/* Do some outputting. */
reportRnaBinders(rbList);

/* Cleanup. */
dMatrixFree(&sjIndex);
dMatrixFree(&psInten);
warn("Done.");
}

int main(int argc, char *argv[])
{
if(argc == 1)
    usage();
optionInit(&argc, argv, optionSpecs);
if(optionExists("help"))
    usage();
rnaBindingClusterCorr();
return 0;
}

ext: c

max_stars_repo_name: andypohl/kent

max_issues_repo_name: andypohl/kent

max_forks_repo_name: andypohl/kent


text:
#ifndef _linalg_h_included_
#define _linalg_h_included_

#include <boost/config.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/operation.hpp>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

#include <jsc/util/log.hpp>

extern "C" {
	extern void dgetrf_(int *,int *,double *,int *,double *,int *);
	extern void dgetri_(int *,double *,int *,double *,double *,int *,int *);
	extern void dpotrf_(char *,int *,double *,int *,int *);
	extern void dpotri_(char *,int *,double *,int *,int *);
};

using namespace boost;
using namespace boost::numeric;
using namespace jsc::util;

class linalg {
	public:
		static void invert_matrix(const ublas::matrix<double> & A,
				ublas::matrix<double> & inv) {
			invert_matrix_gsl_lu(A, inv);
			// invert_matrix_lapack_lu(A, inv);
			// invert_matrix_lapack_cholesky(A, inv);
			unsigned long K = A.size1() + 1;
			ublas::matrix<double> rst(K-1, K-1);
			ublas::matrix<double> ident = ublas::identity_matrix<double>(K-1);
			axpy_prod(A, inv, rst);
			for (unsigned long p = 0; p < K - 1; ++p) {
				for (unsigned long q = 0; q < K - 1; ++q) {
					if (abs(rst(p,q) - ident(p,q)) > 1E-5) {
						L_(debug) << "A*(invA)[" << p << "," << q
							<< "] = " << rst(p,q);
					}
				}
			}
		}
	private:
		static void invert_matrix_gsl_lu(const ublas::matrix<double> & A,
				ublas::matrix<double> & inv) {
			unsigned long K = A.size1() + 1;

			gsl_matrix* m = gsl_matrix_alloc(K - 1, K - 1);
			for (unsigned long p = 0; p < K - 1; ++p) {
				for (unsigned long q = 0; q < K - 1; ++q) {
					gsl_matrix_set(m, p, q, A(p,q));
				}
			}

			gsl_matrix* inverse = gsl_matrix_alloc(K - 1, K - 1);
			gsl_permutation*perm = gsl_permutation_alloc(K-1);
			int s=0;
			gsl_linalg_LU_decomp(m, perm, &s);
			gsl_linalg_LU_invert(m, perm, inverse);
			for (unsigned long p = 0; p < K - 1; ++p) {
				for (unsigned long q = 0; q < K - 1; ++q) {
					inv(p,q) = gsl_matrix_get(inverse, p, q);
				}
			}
			gsl_permutation_free(perm);
			gsl_matrix_free(inverse);
			gsl_matrix_free(m);
		}

		static void invert_matrix_lapack_lu(const ublas::matrix<double> & A,
				ublas::matrix<double> & inv) {
			unsigned long K = A.size1() + 1;

			double * m = (double *)malloc((K - 1)*(K - 1)*sizeof(double));
			for (unsigned long p = 0; p < K - 1; ++p) {
				for (unsigned long q = 0; q < K - 1; ++q) {
					m[p + q*(K-1)] = A(p,q);
				}
			}
			double * ipiv = (double *)malloc((K - 1)*sizeof(double));
			int N = K - 1;
			int lda = N;
			int info;
			dgetrf_(&N, &N, &m[0], &lda, &ipiv[0], &info);
			assert(info == 0);

			double * work = (double *)malloc((K - 1)*sizeof(double));
			int lwork = N;
			dgetri_(&N, &m[0], &lda, &ipiv[0], &work[0], &lwork, &info);
			assert(info == 0);

			for (unsigned long p = 0; p < K - 1; ++p) {
				for (unsigned long q = 0; q < K - 1; ++q) {
					inv(p,q) = m[p + q*(K-1)];
				}
			}

			free(m);
			free(ipiv);
			free(work);
		}

		static void invert_matrix_lapack_cholesky(const ublas::matrix<double> & A,
				ublas::matrix<double> & inv) {
			unsigned long K = A.size1() + 1;

			double * m = (double *)malloc((K - 1)*(K - 1)*sizeof(double));
			for (unsigned long p = 0; p < K - 1; ++p) {
				for (unsigned long q = 0; q < K - 1; ++q) {
					m[p + q*(K-1)] = A(p,q);
				}
			}
			int N = K - 1;
			int lda = N;
			int info;
			char UPLO = 'U';
			dpotrf_(&UPLO, &N, &m[0], &lda, &info);
			assert(info == 0);
			dpotrf_(&UPLO, &N, &m[0], &lda, &info);
			assert(info == 0);

			for (unsigned long p = 0; p < K - 1; ++p) {
				for (unsigned long q = 0; q < K - 1; ++q) {
					inv(p,q) = m[p + q*(K-1)];
				}
			}

			free(m);
		}
};

#endif

max_stars_repo_name: gersteinlab/LESSeq

max_issues_repo_name: gersteinlab/LESSeq

max_forks_repo_name: gersteinlab/LESSeq


text:
// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//

#pragma once
#define _CRT_SECURE_NO_WARNINGS
#define _SCL_SECURE_NO_WARNINGS
#define _CRT_NON_CONFORMING_SWPRINTFS
#define HPX_COMPONENT_EXPORTS
#include "targetver.h"

#include <gsl.h>
#include <stdio.h>
#include <tchar.h>
#include <iostream>
#include <string>

// TODO: reference additional headers your program requires here
#include <hpx/hpx_init.hpp>
#include <hpx/include/actions.hpp>
#include <hpx/include/lcos.hpp>
#include <hpx/include/components.hpp>
#include <hpx/include/serialization.hpp>

#include <boost/any.hpp>
#include <hpx/include/iostreams.hpp>

#include "functions.hpp"
#include "components.hpp"
#include "stubs.h"
#include "clients.h"

// SmallServer references
#include "../../Libs/SmallServer/SmallServer.h"


text:
#include <cblas.h>
#include <stdio.h>
#include <stdlib.h>

int
main()
{
	unsigned int const len = 1000;
	double x[len];
	double y[len];
	for (int i = 0; i < len; i++)
	{
		x[i] = rand() / (double)RAND_MAX;
		y[i] = rand() / (double)RAND_MAX;
	}
	double ret;
	int inc = 1;
	for (int i = 0; i < len*len; i++)
		ret = cblas_ddot(len, x, inc, y, inc);
	return (int)ret;
}

max_stars_repo_name: charles-cooper/idylblas

max_issues_repo_name: charles-cooper/idylblas

max_forks_repo_name: charles-cooper/idylblas


text:
/* ode-initval/gear2.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Gear 2 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "odeiv_util.h"
#include <gsl/gsl_odeiv.h>

/* gear2 state object */
typedef struct
{
  int primed;                   /* flag indicating that yim1 is ready */
  double t_primed;              /* system was primed for this value of t */
  double last_h;                /* last step size */
  gsl_odeiv_step *primer;       /* stepper to use for priming */
  double *yim1;                 /* y_{i-1}    */
  double *k;                    /* work space */
  double *y0;                   /* work space */
  double *y0_orig;
  double *y_onestep;
  int stutter;
}
gear2_state_t;

static void *
gear2_alloc (size_t dim)
{
  gear2_state_t *state = (gear2_state_t *) malloc (sizeof (gear2_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for gear2_state", GSL_ENOMEM);
    }

  state->yim1 = (double *) malloc (dim * sizeof (double));

  if (state->yim1 == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for yim1", GSL_ENOMEM);
    }

  state->k = (double *) malloc (dim * sizeof (double));

  if (state->k == 0)
    {
      free (state->yim1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k", GSL_ENOMEM);
    }

  state->y0 = (double *) malloc (dim * sizeof (double));

  if (state->y0 == 0)
    {
      free (state->k);
      free (state->yim1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
    }

  state->y0_orig = (double *) malloc (dim * sizeof (double));

  if (state->y0_orig == 0)
    {
      free (state->y0);
      free (state->k);
      free (state->yim1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0_orig", GSL_ENOMEM);
    }

  state->y_onestep = (double *) malloc (dim * sizeof (double));

  if (state->y_onestep == 0)
    {
      free (state->y0_orig);
      free (state->y0);
      free (state->k);
      free (state->yim1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0_orig", GSL_ENOMEM);
    }

  state->primed = 0;
  state->primer = gsl_odeiv_step_alloc (gsl_odeiv_step_rk4imp, dim);

  if (state->primer == 0)
    {
      free (state->y_onestep);
      free (state->y0_orig);
      free (state->y0);
      free (state->k);
      free (state->yim1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for primer", GSL_ENOMEM);
    }

  state->last_h = 0.0;

  return state;
}

static int
gear2_step (double *y, gear2_state_t * state,
            const double h, const double t,
            const size_t dim, const gsl_odeiv_system * sys)
{
  /* Makes a Gear2 advance with step size h.
     y0 is the initial values of variables y.
     The implicit matrix equations to solve are:
     k = y0 + h * f(t + h, k)
     y = y0 + h * f(t + h, k)
   */

  const int iter_steps = 3;
  int nu;
  size_t i;
  double *y0 = state->y0;
  double *yim1 = state->yim1;
  double *k = state->k;

  /* Iterative solution of k = y0 + h * f(t + h, k)
     Note: This method does not check for convergence of the
     iterative solution!
   */

  for (nu = 0; nu < iter_steps; nu++)
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, k);

      if (s != GSL_SUCCESS)
        {
          return s;
        }

      for (i = 0; i < dim; i++)
        {
          y[i] = ((4.0 * y0[i] - yim1[i]) + 2.0 * h * k[i]) / 3.0;
        }
    }

  return GSL_SUCCESS;
}

static int
gear2_apply (void *vstate,
             size_t dim,
             double t,
             double h,
             double y[],
             double yerr[],
             const double dydt_in[],
             double dydt_out[], const gsl_odeiv_system * sys)
{
  gear2_state_t *state = (gear2_state_t *) vstate;

  state->stutter = 0;

  if (state->primed == 0 || t == state->t_primed || h != state->last_h)
    {
      /* Execute a single-step method to prime the process.  Note that
       * we do this if the step size changes, so frequent step size
       * changes will cause the method to stutter.
       *
       * Note that we reuse this method if the time has not changed,
       * which can occur when the adaptive driver is attempting to find
       * an appropriate step-size on its first iteration */

      int status;
      DBL_MEMCPY (state->yim1, y, dim);

      status =
        gsl_odeiv_step_apply (state->primer, t, h, y, yerr, dydt_in, dydt_out,
                              sys);

      /* Make note of step size and indicate readiness for a Gear step. */

      state->primed = 1;
      state->t_primed = t;
      state->last_h = h;
      state->stutter = 1;

      return status;
    }
  else
    {
      /* We have a previous y value in the buffer, and the step
       * sizes match, so we go ahead with the Gear step.
       */

      double *const k = state->k;
      double *const y0 = state->y0;
      double *const y0_orig = state->y0_orig;
      double *const yim1 = state->yim1;
      double *y_onestep = state->y_onestep;

      int s;
      size_t i;

      DBL_MEMCPY (y0, y, dim);

      /* iterative solution */

      if (dydt_out != NULL)
        {
          DBL_MEMCPY (k, dydt_out, dim);
        }

      /* First traverse h with one step (save to y_onestep) */

      DBL_MEMCPY (y_onestep, y, dim);

      s = gear2_step (y_onestep, state, h, t, dim, sys);

      if (s != GSL_SUCCESS)
        {
          return s;
        }

      /* Then with two steps with half step length (save to y) */

      s = gear2_step (y, state, h / 2.0, t, dim, sys);

      if (s != GSL_SUCCESS)
        {
          /* Restore original y vector */
          DBL_MEMCPY (y, y0_orig, dim);
          return s;
        }

      DBL_MEMCPY (y0, y, dim);

      s = gear2_step (y, state, h / 2.0, t + h / 2.0, dim, sys);

      if (s != GSL_SUCCESS)
        {
          /* Restore original y vector */
          DBL_MEMCPY (y, y0_orig, dim);
          return s;
        }

      /* Cleanup update */

      if (dydt_out != NULL)
        {
          s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);

          if (s != GSL_SUCCESS)
            {
              /* Restore original y vector */
              DBL_MEMCPY (y, y0_orig, dim);
              return s;
            }
        }

      /* Estimate error and update the state buffer. */

      for (i = 0; i < dim; i++)
        {
          yerr[i] = 4.0 * (y[i] - y_onestep[i]);
          yim1[i] = y0[i];
        }

      /* Make note of step size. */
      state->last_h = h;

      return 0;
    }
}

static int
gear2_reset (void *vstate, size_t dim)
{
  gear2_state_t *state = (gear2_state_t *) vstate;

  DBL_ZERO_MEMSET (state->yim1, dim);
  DBL_ZERO_MEMSET (state->k, dim);
  DBL_ZERO_MEMSET (state->y0, dim);

  state->primed = 0;
  state->last_h = 0.0;
  return GSL_SUCCESS;
}

static unsigned int
gear2_order (void *vstate)
{
  gear2_state_t *state = (gear2_state_t *) vstate;
  state = 0;                    /* prevent warnings about unused parameters */
  return 3;
}

static void
gear2_free (void *vstate)
{
  gear2_state_t *state = (gear2_state_t *) vstate;

  free (state->yim1);
  free (state->k);
  free (state->y0);
  free (state->y0_orig);
  free (state->y_onestep);
  gsl_odeiv_step_free (state->primer);

  free (state);
}

static const gsl_odeiv_step_type gear2_type = { "gear2",        /* name */
  1,                            /* can use dydt_in */
  0,                            /* gives exact dydt_out */
  &gear2_alloc,
  &gear2_apply,
  &gear2_reset,
  &gear2_order,
  &gear2_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_gear2 = &gear2_type;

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*
**
** G.Lohmann, MPI-KYB, May 2018
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "viaio/Vlib.h"
#include "viaio/VImage.h"
#include "viaio/mu.h"
#include "viaio/option.h"

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>

extern void VCheckImage(VImage src);
extern VImage VFilterGauss3d (VImage src,VImage dest,double *sigma);

void Gauss4d(VAttrList list,double *sigma)
{
  /* get pointers to image data */
  double u=0;
  int j,b,r,c;
  int nrows=0,ncols=0,nt=0;
  int nslices = VAttrListNumImages(list);
  VImage *src = VAttrListGetImages(list,nslices);
  VImageDimensions(src,nslices,&nt,&nrows,&ncols);
  fprintf(stderr," image dimensions:  %d x %d x %d,  nt: %d\n",nslices,nrows,ncols,nt);

  VImage tmp = VCreateImage(nslices,nrows,ncols,VFloatRepn);
  VFillImage(tmp,VAllBands,0);
  VImage dest = VCreateImage(nslices,nrows,ncols,VFloatRepn);
  VFillImage(dest,VAllBands,0);

  for (j=0; j<nt; j++) {
    if (j%5==0) fprintf(stderr," %5d  of  %d\r",j,nt);
    for (b=0; b<nslices; b++) {
      for (r=0; r<nrows; r++) {
	for (c=0; c<ncols; c++) {
	  u = VGetPixel(src[b],j,r,c);
	  VSetPixel(tmp,b,r,c,u);
	}
      }
    }
    dest = VFilterGauss3d (tmp,dest,sigma);

    for (b=0; b<nslices; b++) {
      for (r=0; r<nrows; r++) {
	for (c=0; c<ncols; c++) {
	  u = VGetPixel(dest,b,r,c);
	  VSetPixel(src[b],j,r,c,u);
	}
      }
    }
  }
  fprintf(stderr,"\n");
  VDestroyImage(tmp);
  VDestroyImage(dest);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include "utilities.h"
#include "compute_signal.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_histogram.h>
#include <fftw3.h>
#include <stdlib.h>
#include <stdio.h>

#include "parameters.h"

int main(void){
  gsl_rng * r = init_rng();

  double * signal, * spectrum;
  allocate_arrays(&signal, &spectrum, LENGTH);

  gsl_histogram * hist = init_histogram(HIST_MIN, HIST_MAX, BIN_N);

  load_wisdom("wisdom");

  fftw_plan plan = make_plan(signal, LENGTH);

  for(size_t i = 0; i < COUNT; i++){
    printf("---- pass %zu ----\n", i);

    compute_signal(signal, LENGTH, DELTA, r, hist);
    if(i==0) output_signal("signal.dat", signal, LENGTH, DELTA, N_POINTS);

    remove_mean(signal, LENGTH);
    make_fft(plan);
    compute_spectrum(signal, spectrum, LENGTH);

    puts("Done.");
  }

  final_spectrum(spectrum, LENGTH, COUNT, DELTA, SMOOTH_NEIGHBORS);

  output_spectrum("spectrum.dat", spectrum, LENGTH, DELTA, N_POINTS);
  output_distribution("distribution.dat", hist);

  cleanup(signal, spectrum, plan, r, hist);
  return EXIT_SUCCESS;
}


text:
#include <math.h>
#include <stdlib.h>
#if !defined(__APPLE__)
#include <malloc.h>
#endif
#include <stdio.h>
#include <assert.h>
#include <time.h>
#include <string.h>

#include <fftw3.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_deriv.h>

#include "../cosmolike_core/theory/basics.c"
#include "../cosmolike_core/theory/structs.c"
#include "../cosmolike_core/theory/parameters.c"
#include "../cosmolike_core/emu17/P_cb/emu.c"
#include "../cosmolike_core/theory/recompute.c"
#include "../cosmolike_core/theory/cosmo3D.c"
#include "../cosmolike_core/theory/redshift.c"
#include "../cosmolike_core/theory/halo.c"
#include "../cosmolike_core/theory/HOD.c"
#include "../cosmolike_core/theory/cosmo2D_fourier.c"
#include "../cosmolike_core/theory/IA.c"
#include "../cosmolike_core/theory/cluster.c"
#include "../cosmolike_core/theory/BAO.c"
#include "../cosmolike_core/theory/external_prior.c"
#include "../cosmolike_core/theory/covariances_3D.c"
#include "../cosmolike_core/theory/covariances_fourier.c"
#include "../cosmolike_core/theory/covariances_cluster.c"
#include "init_SRD.c"

void run_cov_N_N (char *OUTFILE, char *PATH, int nzc1, int nzc2,int start);
void run_cov_cgl_N (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster,int N1, int nzc2, int start);
void run_cov_cgl_cgl (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int start);
void run_cov_cgl_cgl_all (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster);
void run_cov_shear_N (char *OUTFILE, char *PATH, double *ell, double *dell, int N1, int nzc2, int start);
void run_cov_shear_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start);
void run_cov_ggl_N (char *OUTFILE, char *PATH, double *ell, double *dell, int N1, int nzc2, int start);
void run_cov_ggl_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start);
void run_cov_cl_N (char *OUTFILE, char *PATH, double *ell, double *dell, int N1, int nzc2, int start);
void run_cov_cl_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start);

void run_cov_ggl_shear(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_clustering_shear(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_clustering_ggl(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_clustering(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_ggl(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);
void run_cov_shear_shear(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start);

void run_cov_N_N (char *OUTFILE, char *PATH, int nzc1, int nzc2,int start)
{
  int nN1, nN2,i,j;
  double cov;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  for (nN1 = 0; nN1 < Cluster.N200_Nbin; nN1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
      i += Cluster.N200_Nbin*nzc1+nN1;
      j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
      j += Cluster.N200_Nbin*nzc2+nN2;

      cov =cov_N_N(nzc1,nN1, nzc2, nN2);
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j,0.0,0.0, nzc1, nN1, nzc2, nN2,cov,0.0);
    }
  }
  fclose(F1);
}

void run_cov_cgl_N (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster,int N1, int nzc2, int start)
{
  int nN1, nN2, nl1, nzc1, nzs1,i,j;
  double cov;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  nzc1 = ZC(N1);
  nzs1 = ZSC(N1);
  for (nN1 = 0; nN1 < Cluster.N200_Nbin; nN1 ++){
    for( nl1 = 0; nl1 < Cluster.lbin; nl1 ++){
     for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
       i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
       i += (N1*Cluster.N200_Nbin+nN1)*Cluster.lbin +nl1;
       j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
       j += Cluster.N200_Nbin*nzc2+nN2;

       cov =cov_cgl_N(ell_Cluster[nl1],nzc1,nN1, nzs1, nzc2, nN2);
       fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j, ell_Cluster[nl1], 0., nzc1, nzs1, nzc2, nN2,cov,0.);
     }
   }
 }
 fclose(F1);
}

void run_cov_cgl_cgl (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int start)
{
  int nN1, nN2, nl1, nzc1, nzs1, nl2, nzc2, nzs2,i,j;
  double c_g, c_ng;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  nzc1 = ZC(N1);
  nzs1 = ZSC(N1);
  nzc2 = ZC(N2);
  nzs2 = ZSC(N2);
  for (nN1 = 0; nN1 < Cluster.N200_Nbin; nN1 ++){
    for( nl1 = 0; nl1 < Cluster.lbin; nl1 ++){
      for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
        for( nl2 = 0; nl2 < Cluster.lbin; nl2 ++){
          i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
          i += (N1*Cluster.N200_Nbin+nN1)*Cluster.lbin +nl1;
          j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
          j += (N2*Cluster.N200_Nbin+nN2)*Cluster.lbin +nl2;

          c_g = 0;
          c_ng = cov_NG_cgl_cgl(ell_Cluster[nl1],ell_Cluster[nl2],nzc1,nN1, nzs1, nzc2, nN2,nzs2);
          if (nl2 == nl1){c_g =cov_G_cgl_cgl(ell_Cluster[nl1],dell_Cluster[nl1],nzc1,nN1, nzs1, nzc2, nN2,nzs2);}
          fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j,ell_Cluster[nl1],ell_Cluster[nl2], nzc1, nzs1, nzc2, nzs2,c_g, c_ng);
        }
      }
    }
  }
  fclose(F1);
}

void run_cov_cgl_cgl_all (char *OUTFILE, char *PATH, double *ell_Cluster, double *dell_Cluster)
{
  int nN1, nN2, nl1, nzc1, nzs1, nl2, nzc2, nzs2,i,j, N1,N2;
  double c_g, c_ng;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s",PATH,OUTFILE);
  F1 =fopen(filename,"w");
  for (N1 = 0; N1 < tomo.cgl_Npowerspectra; N1 ++){
    for (N2 = 0; N2 < tomo.cgl_Npowerspectra; N2 ++){
      nzc1 = ZC(N1);
      nzs1 = ZSC(N1);
      nzc2 = ZC(N2);
      nzs2 = ZSC(N2);
      for (nN1 = 0; nN1 < Cluster.N200_Nbin; nN1 ++){
        for( nl1 = 0; nl1 < Cluster.lbin; nl1 ++){
          for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
            for( nl2 = 0; nl2 < Cluster.lbin; nl2 ++){
              i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
              i += (N1*Cluster.N200_Nbin+nN1)*Cluster.lbin +nl1;
              j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
              j += (N2*Cluster.N200_Nbin+nN2)*Cluster.lbin +nl2;
              c_g = 0;
              c_ng = cov_NG_cgl_cgl(ell_Cluster[nl1],ell_Cluster[nl2],nzc1,nN1, nzs1, nzc2, nN2,nzs2);
              if (nl2 == nl1){c_g =cov_G_cgl_cgl(ell_Cluster[nl1],dell_Cluster[nl1],nzc1,nN1, nzs1, nzc2, nN2,nzs2);}
              fprintf(F1,"%d %d %e %e %d %d %d  %d %d %d  %e %e\n",i,j,ell_Cluster[nl1],ell_Cluster[nl2], nzc1, nN1,nzs1, nzc2, nN2,nzs2,c_g, c_ng);
            }
          }
        }
      }
    }
  }
  fclose(F1);
}

void run_cov_shear_N (char *OUTFILE, char *PATH, double *ell, double *dell, int N1, int nzc2, int start)
{
  int nz1,nz2, nN2, nl1, nzc1, i,j;
  double cov;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  nz1 = Z1(N1);
  nz2 = Z2(N1);
  for( nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      cov = 0.;
      i = like.Ncl*N1+nl1;
      j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
      j += Cluster.N200_Nbin*nzc2+nN2;

      if (ell[nl1] < like.lmax_shear){cov =cov_shear_N(ell[nl1],nz1,nz2, nzc2, nN2);}
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j, ell[nl1], 0., nz1, nz2, nzc2, nN2,cov,0.);
    }
  }
  fclose(F1);
}

void run_cov_shear_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start)
{
  int nN1, nN2, nzs1, nzs2,nl2, nzc2, nzs3,i,j;
  double c_g, c_ng;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  nzs1 = Z1(N1);
  nzs2 = Z2(N1);
  nzc2 = ZC(N2);
  nzs3 = ZSC(N2);
  for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      for( nl2 = 0; nl2 < Cluster.lbin; nl2 ++){
        i = like.Ncl*N1+nl1;
        j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
        j += (N2*Cluster.N200_Nbin+nN2)*Cluster.lbin +nl2;
        c_g = 0.;
        c_ng = 0.;
        if (ell[nl1] < like.lmax_shear){
          c_ng = cov_NG_shear_cgl(ell[nl1],ell_Cluster[nl2],nzs1, nzs2, nzc2, nN2,nzs3);
          if (fabs(ell[nl1]/ell_Cluster[nl2] -1.) < 0.001){
            c_g =cov_G_shear_cgl(ell[nl1],dell_Cluster[nl2],nzs1,nzs2, nzc2, nN2,nzs3);
          }
        }
        fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",i,j,ell[nl1],ell_Cluster[nl2], nzs1, nzs2, nzc2, nzs3,c_g, c_ng);
      }
    }
  }
  fclose(F1);
}

void run_cov_ggl_N (char *OUTFILE, char *PATH, double *ell, double *dell, int N1, int nzc2, int start)
{
  int zl,zs, nN2, nl1, nzc1, i,j;
  double cov,weight;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  zl = ZL(N1);
  zs = ZS(N1);
  for( nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      i = like.Ncl*(tomo.shear_Npowerspectra+N1)+nl1;
      j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
      j += Cluster.N200_Nbin*nzc2+nN2;
      cov = 0.;
      weight = test_kmax(ell[nl1],zl);
      if (weight){
        cov =cov_ggl_N(ell[nl1],zl,zs, nzc2, nN2);
      }
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j, ell[nl1], 0., zl, zs, nzc2, nN2,cov,0.);
    }
  }
  fclose(F1);
}

void run_cov_ggl_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start)
{
  int nN2, zl, zs, nzs1, nl2, nzc2, nzs3,i,j;
  double c_g, c_ng,weight;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  zl = ZL(N1);
  zs = ZS(N1);
  nzc2 = ZC(N2);
  nzs3 = ZSC(N2);
  for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      for( nl2 = 0; nl2 < Cluster.lbin; nl2 ++){
        i = like.Ncl*(tomo.shear_Npowerspectra+N1)+nl1;
        j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
        j += (N2*Cluster.N200_Nbin+nN2)*Cluster.lbin +nl2;

        c_g = 0; c_ng = 0.;
        weight = test_kmax(ell[nl1],zl);
        if (weight){
          c_ng = cov_NG_ggl_cgl(ell[nl1],ell_Cluster[nl2],zl,zs, nzc2, nN2,nzs3);
          if (fabs(ell[nl1]/ell_Cluster[nl2] -1.) < 0.1){c_g =cov_G_ggl_cgl(ell[nl1],dell_Cluster[nl2],zl,zs, nzc2, nN2,nzs3);}
        }
        fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j,ell[nl1],ell_Cluster[nl2], zl, zs, nzc2, nzs3,c_g, c_ng);
      }
    }
  }
  fclose(F1);
}

void run_cov_cl_N (char *OUTFILE, char *PATH, double *ell, double *dell,int N1, int nzc2, int start)
{
  int zl,zs, nN2, nl1, nzc1, i,j;
  double cov,weight;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  for( nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+N1)+nl1;
      j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra);
      j += Cluster.N200_Nbin*nzc2+nN2;
      cov = 0.;
      weight = test_kmax(ell[nl1],N1);
      if (weight){
        cov =cov_cl_N(ell[nl1],N1,N1,nzc2,nN2);
      }
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j, ell[nl1], 0., N1, N1, nzc2, nN2,cov,0.);
    }
  }
  fclose(F1);
}

void run_cov_cl_cgl (char *OUTFILE, char *PATH, double *ell, double *dell, double *ell_Cluster, double *dell_Cluster,int N1, int N2, int nl1, int start)
{
  int nN2,nzc2, nzs3,i,j,nl2;
  double c_g, c_ng,weight;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  nzc2 = ZC(N2);
  nzs3 = ZSC(N2);
  for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nN2 = 0; nN2 < Cluster.N200_Nbin; nN2 ++){
      for( nl2 = 0; nl2 < Cluster.lbin; nl2 ++){
        i = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+N1)+nl1;
        j = like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+tomo.clustering_Npowerspectra)+Cluster.N200_Nbin*tomo.cluster_Nbin;
        j += (N2*Cluster.N200_Nbin+nN2)*Cluster.lbin +nl2;

        c_g = 0; c_ng = 0.;
        weight = test_kmax(ell[nl1],N1);
        if (weight){
          c_ng = cov_NG_cl_cgl(ell[nl1],ell_Cluster[nl2],N1,N1, nzc2, nN2,nzs3);
          if (fabs(ell[nl1]/ell_Cluster[nl2] -1.) < 0.1){
            c_g =cov_G_cl_cgl(ell[nl1],dell_Cluster[nl2],N1,N1, nzc2, nN2,nzs3);
          }
        }
        fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n",i,j,ell[nl1],ell_Cluster[nl2], N1,N1, nzc2, nzs3,c_g, c_ng);
        //printf("%d %d %e %e %d %d %d %d  %e %e\n",i,j,ell[nl1],ell_Cluster[nl2], N1,N1, nzc2, nzs3,c_g, c_ng);
      }
    }
  }
  fclose(F1);
}

void run_cov_ggl_shear(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int zl,zs,z3,z4,nl1,nl2,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  zl = ZL(n1); zs = ZS(n1);
  printf("\nN_ggl = %d (%d, %d)\n", n1,zl,zs);
  z3 = Z1(n2); z4 = Z2(n2);
  printf("N_shear = %d (%d, %d)\n", n2,z3,z4);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 <like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],zl);
      if (weight && ell[nl2] < like.lmax_shear){
        if (test_zoverlap(zl,z3)*test_zoverlap(zl,z4)){
          c_ng = cov_NG_gl_shear_tomo(ell[nl1],ell[nl2],zl,zs,z3,z4);
        }
        if (nl1 == nl2){
          c_g =  cov_G_gl_shear_tomo(ell[nl1],dell[nl1],zl,zs,z3,z4);
        }
      }
      fprintf(F1, "%d %d %e %e %d %d %d %d  %e %e\n", like.Ncl*(tomo.shear_Npowerspectra+n1)+nl1,like.Ncl*(n2)+nl2, ell[nl1],ell[nl2],zl,zs,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_clustering_shear(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int z1,z2,z3,z4,nl1,nl2,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  z1 = n1; z2 = n1;
  printf("\nN_cl = %d \n", n1);
  z3 = Z1(n2); z4 = Z2(n2);
  printf("N_shear = %d (%d, %d)\n", n2,z3,z4);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 <like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],z1);
      if (weight && ell[nl2] < like.lmax_shear){
        if (test_zoverlap(z1,z3)*test_zoverlap(z1,z4)){
          c_ng = cov_NG_cl_shear_tomo(ell[nl1],ell[nl2],z1,z2,z3,z4);
        }
        if (nl1 == nl2){
          c_g =  cov_G_cl_shear_tomo(ell[nl1],dell[nl1],z1,z2,z3,z4);
        }
      }
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n", like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1,like.Ncl*(n2)+nl2, ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_clustering_ggl(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int z1,z2,zl,zs,nl1,nl2,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  z1 = n1; z2 = n1;
  printf("\nN_cl_1 = %d \n", n1);
  zl = ZL(n2); zs = ZS(n2);
  printf("N_tomo_2 = %d (%d, %d)\n", n2,zl,zs);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      if (z1 == zl){
        weight = test_kmax(ell[nl1],z1)*test_kmax(ell[nl2],zl);
        if (weight){
          c_ng = cov_NG_cl_gl_tomo(ell[nl1],ell[nl2],z1,z2,zl,zs);
          if (nl1 == nl2){
            c_g =  cov_G_cl_gl_tomo(ell[nl1],dell[nl1],z1,z2,zl,zs);
          }
        }
      }
      fprintf(F1, "%d %d %e %e %d %d %d %d  %e %e\n", like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1,like.Ncl*(tomo.shear_Npowerspectra+n2)+nl2, ell[nl1],ell[nl2],z1,z2,zl,zs,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_clustering(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int z1,z2,z3,z4,nl1,nl2,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  z1 = n1; z2 = n1;
  printf("\nN_cl_1 = %d \n", n1);
  z3 = n2; z4 = n2;
  printf("N_cl_2 = %d\n", n2);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      if (z1 == z3){
        weight = test_kmax(ell[nl1],z1)*test_kmax(ell[nl2],z3);
        if (weight) {
          c_ng = cov_NG_cl_cl_tomo(ell[nl1],ell[nl2],z1,z2,z3,z4);
        }
        if (nl1 == nl2){
          c_g =  cov_G_cl_cl_tomo(ell[nl1],dell[nl1],z1,z2,z3,z4);
        }
      }
      fprintf(F1, "%d %d %e %e %d %d %d %d %e %e\n", like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra+n1)+nl1,like.Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra + n2)+nl2, ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_ggl(char *OUTFILE, char *PATH, double *ell, double *dell, int n1, int n2,int start)
{
  int zl1,zl2,zs1,zs2,nl1,nl2, weight;
  double c_ng, c_g;
  double fsky = survey.area/41253.0;
  FILE *F1;
  char filename[300];
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");

  zl1 = ZL(n1); zs1 = ZS(n1);
  printf("\nN_tomo_1 = %d (%d, %d)\n", n1,zl1,zs1);
  zl2 = ZL(n2); zs2 = ZS(n2);
  printf("N_tomo_2 = %d (%d, %d)\n", n2,zl2,zs2);
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      weight = test_kmax(ell[nl1],zl1)*test_kmax(ell[nl2],zl2);
      if (weight && zl1 == zl2) {
        c_ng = cov_NG_gl_gl_tomo(ell[nl1],ell[nl2],zl1,zs1,zl2,zs2);
      }
      if (nl1 == nl2){
        c_g =  cov_G_gl_gl_tomo(ell[nl1],dell[nl1],zl1,zs1,zl2,zs2);
      }
      if (weight ==0 && n2 != n1){
        c_g = 0;
      }
      fprintf(F1,"%d %d %e %e %d %d %d %d  %e %e\n", like.Ncl*(tomo.shear_Npowerspectra+n1)+nl1,like.Ncl*(tomo.shear_Npowerspectra+n2)+nl2, ell[nl1],ell[nl2],zl1,zs1,zl2,zs2,c_g,c_ng);
    }
  }
  fclose(F1);
}

void run_cov_shear_shear(char *OUTFILE, char *PATH, double *ell, double *dell,int n1, int n2,int start)
{
  int z1,z2,z3,z4,nl1,nl2,weight;
  double c_ng, c_g;
  FILE *F1;
  char filename[300];
  z1 = Z1(n1); z2 = Z2(n1);
  printf("N_shear = %d\n", n1);
  z3 = Z1(n2); z4 = Z2(n2);
  printf("N_shear = %d (%d, %d)\n",n2,z3,z4);
  sprintf(filename,"%s%s_%d",PATH,OUTFILE,start);
  F1 =fopen(filename,"w");
  for (nl1 = 0; nl1 < like.Ncl; nl1 ++){
    for (nl2 = 0; nl2 < like.Ncl; nl2 ++){
      c_ng = 0.; c_g = 0.;
      if (ell[nl1] < like.lmax_shear && ell[nl2] < like.lmax_shear){
        c_ng = cov_NG_shear_shear_tomo(ell[nl1],ell[nl2],z1,z2,z3,z4);
      }
      if (nl1 == nl2){
        c_g =  cov_G_shear_shear_tomo(ell[nl1],dell[nl1],z1,z2,z3,z4);
        if (ell[nl1] > like.lmax_shear && n1!=n2){c_g = 0.;}
      }
      fprintf(F1,"%d %d %e %e %d %d %d %d %e %e\n",like.Ncl*n1+nl1,like.Ncl*(n2)+nl2,ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
      //printf("%d %d %e %e %d %d %d %d %e %e\n", like.Ncl*n1+nl1,like.Ncl*(n2)+nl2, ell[nl1],ell[nl2],z1,z2,z3,z4,c_g,c_ng);
    }
  }
  fclose(F1);
}

int main(int argc, char** argv)
{

  int i,l,m,n,o,s,p,nl1,t,k;
  char OUTFILE[400],filename[400],arg1[400],arg2[400];

  int N_scenarios=12;
  double area_table[12]={7500.0,13000.0,16000.0,10000.0,15000.0,20000.0,10000.0,15000.0,20000.0,10000.0,15000.0,20000.0};
  double nsource_table[12]={9.8,12.1,15.1,15.1,18.9,23.5,20.3,23.5,26.9,26.9,30.8,35.0};
  double nlens_table[12]={15.0,20.0,25.0,25.0,32.0,41.0,35.0,41.0,48.0,48.0,57.0,67.0};

  char survey_designation[12][200]={"LSST_Y1","LSST_Y1","LSST_Y1","LSST_Y3","LSST_Y3","LSST_Y3","LSST_Y6","LSST_Y6","LSST_Y6","LSST_Y10","LSST_Y10","LSST_Y10"};
  char source_zfile[12][400]={"WL_zdistri_model0_z0=1.940000e-01_alpha=8.830000e-01","WL_zdistri_model1_z0=1.900000e-01_alpha=8.620000e-01","WL_zdistri_model2_z0=1.860000e-01_alpha=8.410000e-01","WL_zdistri_model3_z0=1.860000e-01_alpha=8.410000e-01","WL_zdistri_model4_z0=1.830000e-01_alpha=8.210000e-01","WL_zdistri_model5_z0=1.790000e-01_alpha=8.000000e-01","WL_zdistri_model6_z0=1.810000e-01_alpha=8.140000e-01","WL_zdistri_model7_z0=1.790000e-01_alpha=8.000000e-01","WL_zdistri_model8_z0=1.760000e-01_alpha=7.860000e-01","WL_zdistri_model9_z0=1.760000e-01_alpha=7.860000e-01","WL_zdistri_model10_z0=1.740000e-01_alpha=7.720000e-01","WL_zdistri_model11_z0=1.710000e-01_alpha=7.590000e-01"};
  char lens_zfile[12][400]={"LSS_zdistri_model0_z0=2.590000e-01_alpha=9.520000e-01","LSS_zdistri_model1_z0=2.610000e-01_alpha=9.370000e-01","LSS_zdistri_model2_z0=2.640000e-01_alpha=9.250000e-01","LSS_zdistri_model3_z0=2.640000e-01_alpha=9.250000e-01","LSS_zdistri_model4_z0=2.680000e-01_alpha=9.150000e-01","LSS_zdistri_model5_z0=2.740000e-01_alpha=9.070000e-01","LSS_zdistri_model6_z0=2.700000e-01_alpha=9.120000e-01","LSS_zdistri_model7_z0=2.740000e-01_alpha=9.070000e-01","LSS_zdistri_model8_z0=2.780000e-01_alpha=9.030000e-01","LSS_zdistri_model9_z0=2.780000e-01_alpha=9.030000e-01","LSS_zdistri_model10_z0=2.830000e-01_alpha=9.000000e-01","LSS_zdistri_model11_z0=2.880000e-01_alpha=8.980000e-01"};

  int Ntomo_lens[12]={5,5,5,7,7,7,9,9,9,10,10,10};
    //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    // use this remapping only if files fail !!!!!!!!!
    //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    // int fail[5]={1134, 259, 497, 623, 718};
    // hit=fail[hit-1];
    //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  int hit=atoi(argv[1]);
  Ntable.N_a=20;
  k=1;

  for(t=0;t<1;t++){

    //RUN MODE setup
    init_cosmo();
    init_binning_fourier(20,20.0,15000.0,3000.0,21.0,5,Ntomo_lens[t]);
    init_survey(survey_designation[t]);
    sprintf(arg1,"zdistris/%s",source_zfile[t]);
    sprintf(arg2,"zdistris/%s",lens_zfile[t]);
    init_galaxies(arg1,arg2,"none","none","source");
    init_clusters();
    init_IA("none", "GAMA");
    init_probes("3x2pt_clusterN_clusterWL");

    //set l-bins for shear, ggl, clustering, clusterWL
    double logdl=(log(like.lmax)-log(like.lmin))/like.Ncl;
    double *ell, *dell, *ell_Cluster, *dell_Cluster;
    ell=create_double_vector(0,like.Ncl-1);
    dell=create_double_vector(0,like.Ncl-1);
    ell_Cluster=create_double_vector(0,Cluster.lbin-1);
    dell_Cluster=create_double_vector(0,Cluster.lbin-1);
    int j=0;
    for(i=0;i<like.Ncl;i++){
      ell[i]=exp(log(like.lmin)+(i+0.5)*logdl);
      dell[i]=exp(log(like.lmin)+(i+1)*logdl)-exp(log(like.lmin)+(i*logdl));
      if(ell[i]<like.lmax_shear) printf("%le\n",ell[i]);
      if(ell[i]>like.lmax_shear){
        ell_Cluster[j]=ell[i];
        dell_Cluster[j]=dell[i];
        printf("%le %le\n",ell[i],ell_Cluster[j]);
        j++;
      }
    }

    printf("----------------------------------\n");
    survey.area=area_table[t];
    survey.n_gal=nsource_table[t];
    survey.n_lens=nlens_table[t];

    sprintf(survey.name,"%s_area%le_ng%le_nl%le",survey_designation[t],survey.area,survey.n_gal,survey.n_lens);
    printf("area: %le n_source: %le n_lens: %le\n",survey.area,survey.n_gal,survey.n_lens);

    sprintf(covparams.outdir,"/home/u17/timeifler/covparallel/");

    printf("----------------------------------\n");
    sprintf(OUTFILE,"%s_ssss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.shear_Npowerspectra; l++){
      for (m=l;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_shear_shear(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        k=k+1;
        //printf("%d\n",k);
      }
    }

    sprintf(OUTFILE,"%s_lsls_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=l;m<tomo.ggl_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_ggl(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
       //printf("%d\n",k);
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_llll_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){ //auto bins only for now!
      for (m=l;m<tomo.clustering_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_clustering(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        k=k+1;
        //printf("%d %d %d\n",l,m,k);
      }
    }
    sprintf(OUTFILE,"%s_llss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_clustering_shear(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        k=k+1;
        //printf("%d\n",k);
      }
    }
    sprintf(OUTFILE,"%s_llls_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.ggl_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_clustering_ggl(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        k=k+1;
        //printf("%d\n",k);
      }
    }
    sprintf(OUTFILE,"%s_lsss_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_ggl_shear(OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        k=k+1;
        //printf("%d\n",k);
      }
    }

    //******************************
    //******cluster covariance******
    //******************************

    sprintf(OUTFILE,"%s_nn_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.cluster_Nbin; l++){
      for (m=0;m<tomo.cluster_Nbin; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_N_N (OUTFILE,covparams.outdir,l,m,k);
          }
        }
        k=k+1;
        //printf("%d\n",k);
      }
    }
    sprintf(OUTFILE,"%s_cscs_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.cgl_Npowerspectra; l++){
      for (m=0;m<tomo.cgl_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_cgl_cgl (OUTFILE,covparams.outdir,ell_Cluster,dell_Cluster,l,m,k);
          }
        }
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_csn_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.cgl_Npowerspectra; l++){
      for (m=0;m<tomo.cluster_Nbin; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
            run_cov_cgl_N (OUTFILE,covparams.outdir,ell_Cluster,dell_Cluster,l,m,k);
          }
        }
        k=k+1;
      }
    }
    //shear X cluster
    sprintf(OUTFILE,"%s_ssn_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.shear_Npowerspectra; l++){
      for (m=0;m<tomo.cluster_Nbin; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
           run_cov_shear_N (OUTFILE,covparams.outdir,ell,dell,l,m,k);
         }
       }
       k=k+1;
     }
   }
   sprintf(OUTFILE,"%s_sscs_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
   for (l=0;l<tomo.shear_Npowerspectra; l++){
      for (m=0;m<tomo.cgl_Npowerspectra; m++){
        //for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
          if(k==hit){
            sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
            if (fopen(filename, "r") != NULL){exit(1);}
            else {
              run_cov_shear_cgl (OUTFILE,covparams.outdir,ell,dell,ell_Cluster,dell_Cluster,l,m,nl1,k);
            }
          }
          k=k+1;
        //}
      }
    }
    // ggl X cluster
    sprintf(OUTFILE,"%s_lsn_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=0;m<tomo.cluster_Nbin; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
           run_cov_ggl_N (OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_lscs_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=0;m<tomo.cgl_Npowerspectra; m++){
        //for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
          if(k==hit){
            sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
            if (fopen(filename, "r") != NULL){exit(1);}
            else {
              run_cov_ggl_cgl (OUTFILE,covparams.outdir,ell,dell,ell_Cluster,dell_Cluster,l,m,nl1,k);
            }
          }
          k=k+1;
        //}
      }
    }
    // clustering X cluster
    sprintf(OUTFILE,"%s_lln_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.cluster_Nbin; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {
           run_cov_cl_N (OUTFILE,covparams.outdir,ell,dell,l,m,k);
          }
        }
        //printf("%d\n",k);
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_llcs_cov_Ncl%d_Ntomo%d",survey.name,like.Ncl,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.cgl_Npowerspectra; m++){
        //for(nl1 = 0; nl1 < like.Ncl; nl1 ++){
          if(k==hit){
            sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
            if (fopen(filename, "r") != NULL){exit(1);}
            else {
              run_cov_cl_cgl (OUTFILE,covparams.outdir,ell,dell,ell_Cluster,dell_Cluster,l,m,nl1,k);
            }
          }
          k=k+1;
        //}
      }
    }
  }
  printf("number of cov blocks for parallelization: %d\n",k-1);
  printf("-----------------\n");
  printf("PROGRAM EXECUTED\n");
  printf("-----------------\n");
  return 0;
}


text:
#include <math.h>
#include <stdlib.h>
#if !defined(__APPLE__)
#include <malloc.h>
#endif
#include <stdio.h>
#include <assert.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sf_expint.h>

#include "../cosmolike_light/theory/baryons.h"
#include "../cosmolike_light/theory/basics.c"
#include "../cosmolike_light/theory/structs.c"
#include "../cosmolike_light/theory/parameters.c"
#include "../cosmolike_light/emu17/P_cb/emu.c"
#include "../cosmolike_light/theory/recompute.c"
#include "../cosmolike_light/theory/cosmo3D.c"
#include "../cosmolike_light/theory/redshift_spline.c"
#include "../cosmolike_light/theory/halo.c"
#include "../cosmolike_light/theory/HOD.c"
#include "../cosmolike_light/theory/cosmo2D_fourier.c"
#include "../cosmolike_light/theory/IA.c"
#include "../cosmolike_light/theory/BAO.c"
#include "../cosmolike_light/theory/external_prior.c"
#include "../cosmolike_light/theory/covariances_3D.c"
#include "../cosmolike_light/theory/covariances_fourier.c"
#include "../cosmolike_light/theory/covariances_real.c"
#include "../cosmolike_light/theory/run_covariances_real.c"
#include "../cosmolike_light/theory/init.c"

int main(int argc, char** argv)
{
  int hit=atoi(argv[1]);
  FILE *F1,*F2;
  int i,l,m,n,o,s,p,output;
  double ktmp;
  char OUTFILE[400],filename[400];

  printf("-----------------\n");
  printf("This is not production mode covariance code\n");
  printf("Please only use with care and at your own risk\n");
  printf("-----------------\n");

  Ntable.N_a=20;
  set_cov_parameters_to_("cov_Y1/cov_y1_mcal_revision.ini",1);
  //here: setting values internally

  // set this to zero to quickly run Gaussian-only covariances for testing
  if (covparams.ng==1){
    NG = 1;
  }
  else {
    NG = 0;
  }

  // set this to one to output details about inputs for diagnostics
  output = 0;
  FILE *F;
  printf("running multi_covariance_real with NG = %d\n",NG);

  set_cosmological_parameters_to_("cov_Y1/cov_y1_mcal_revision.ini",1);
  //here: setting values internally
  // cosmology.Omega_m   = 0.286;
  // cosmology.Omega_v   = 1.0-cosmology.Omega_m;
  // cosmology.sigma_8   = 0.82;
  // cosmology.n_spec    = 0.96;

  // cosmology.w0=-1.;
  // cosmology.wa=0.;
  // cosmology.omb=0.04868;
  // cosmology.h0=0.673;
  // cosmology.coverH0= 2997.92458;

  // cosmology.rho_crit = 7.4775e+21;
  // cosmology.f_NL = 0.0;
  // pdeltaparams.runmode="Halofit"

  set_survey_parameters_to_("cov_Y1/cov_y1_mcal_revision.ini",1);
  //here: setting values internally
  //survey.area=1000.0;
  //survey.m_lim=
  //survey.name=
  //survey.Kcorrect_File=
  // redshift.shear_REDSHIFT_FILE=source.nz
  // redshift.clustering_REDSHIFT_FILE=lens.nz
  // survey.sourcephotoz="multihisto";
  // survey.lensphotoz="multihisto";
  // survey.galsample=
  // tomo.shear_Nbin=5;
  // redshift.clustering_histogram_zbins =
  // tomo.clustering_Nbin=5;
  // tomo.clustering_Npowerspectra=tomo.clustering_Nbin;
  // survey.sigma_e=0.24;
  // survey.ggl_overlap_cut=
  // tomo.n_source[0]=1.07535691109;
  // tomo.n_source[1]=1.1606931903;
  // tomo.n_source[2]=1.2275933389;
  // tomo.n_source[3]=1.23453951309;
  // tomo.n_source[4]=1.3387876584;

  // tomo.n_lens[0]=0.02143277;
  // tomo.n_lens[1]=0.05426833;
  // tomo.n_lens[2]=0.08070083;
  // tomo.n_lens[3]=0.05130111;
  // tomo.n_lens[4]=0.01515694;
  // for (i=0;i<tomo.clustering_Nbin:i++) survey.n_lens+=tomo.n_lens[i];
  // for (i=0;i<tomo.shear_Nbin:i++) survey.n_gal+=tomo.n_source[i];

  init_source_sample_();
  init_lens_sample_();

  //init_clusters();
  //init_IA("none", "GAMA");
  //printf("test values: %d, %d, %s",redshift.clustering_photoz,tomo.clustering_Nbin,redshift.clustering_REDSHIFT_FILE);
  // printf("end of setup in main\n");

  double theta_min= covparams.tmin;
  double theta_max= covparams.tmax;
  int Ntheta=covparams.ntheta;

  double logdt=(log(theta_max)-log(theta_min))/Ntheta;
  double *theta,*thetamin,*thetamax, *dtheta;
  theta=create_double_vector(0,Ntheta-1);
  thetamin=create_double_vector(0,Ntheta);
  thetamax=create_double_vector(0,Ntheta-1);
  dtheta=create_double_vector(0,Ntheta-1);
  for(i=0; i<Ntheta ; i++){
    thetamin[i]=exp(log(theta_min)+(i+0.0)*logdt);
    thetamax[i]=exp(log(theta_min)+(i+1.0)*logdt);
    theta[i] = 2./3.*(pow(thetamax[i],3.)-pow(thetamin[i],3.))/(pow(thetamax[i],2.)-pow(thetamin[i],2.));
    dtheta[i]=thetamax[i]-thetamin[i];
  }
  thetamin[Ntheta] = thetamax[Ntheta-1];

  int k=1;
  if (strcmp(covparams.ss,"true")==0)
  {
    sprintf(OUTFILE,"%s_ssss_++_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);

    for (l=0;l<tomo.shear_Npowerspectra; l++){
      for (m=l;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {run_cov_shear_shear_real_bin(OUTFILE,covparams.outdir,thetamin,Ntheta,l,m,1,1,k);}
        }
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_ssss_--_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.shear_Npowerspectra; l++){
      for (m=l;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {run_cov_shear_shear_real_bin(OUTFILE,covparams.outdir,thetamin,Ntheta,l,m,0,0,k);}
        }
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_ssss_+-_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.shear_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {run_cov_shear_shear_real_bin(OUTFILE,covparams.outdir,thetamin,Ntheta,l,m,1,0,k);}
        }
        k=k+1;
      }
    }
  }
  if (strcmp(covparams.ll,"true")==0)
  {
    sprintf(OUTFILE,"%s_llll_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=l;m<tomo.clustering_Npowerspectra; m++){
        if(k==hit){
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {run_cov_clustering_real_bin(OUTFILE,covparams.outdir,thetamin,Ntheta,l,m,k);}
        }
        k=k+1;
      }
    }
  }
  if (strcmp(covparams.ls,"true")==0)
  {
    sprintf(OUTFILE,"%s_lsls_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=l;m<tomo.ggl_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {run_cov_ggl_real_bin(OUTFILE,covparams.outdir,thetamin,Ntheta,l,m,k);}
        }
        k=k+1;
      }
    }
  }
  if (strcmp(covparams.ls,"true")==0 && strcmp(covparams.ss,"true")==0)
  {
    sprintf(OUTFILE,"%s_lsss_+_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {run_cov_ggl_shear_real_bin(OUTFILE,covparams.outdir,thetamin,Ntheta,l,m,1,k);}
        }
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_lsss_-_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.ggl_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {run_cov_ggl_shear_real_bin(OUTFILE,covparams.outdir,thetamin,Ntheta,l,m,0,k);}
        }
        k=k+1;
      }
    }
  }
  if (strcmp(covparams.ll,"true")==0 && strcmp(covparams.ss,"true")==0)
  {
    sprintf(OUTFILE,"%s_llss_+_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {run_cov_clustering_shear_real_bin(OUTFILE,covparams.outdir,thetamin,Ntheta,l,m,1,k);}
        }
        k=k+1;
      }
    }
    sprintf(OUTFILE,"%s_llss_-_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.shear_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {run_cov_clustering_shear_real_bin(OUTFILE,covparams.outdir,thetamin,Ntheta,l,m,0,k);}
        }
        k=k+1;
      }
    }
  }
  if (strcmp(covparams.ll,"true")==0 && strcmp(covparams.ls,"true")==0)
  {
    sprintf(OUTFILE,"%s_llls_cov_Ntheta%d_Ntomo%d",covparams.filename,Ntheta,tomo.shear_Nbin);
    for (l=0;l<tomo.clustering_Npowerspectra; l++){
      for (m=0;m<tomo.ggl_Npowerspectra; m++){
        if(k==hit) {
          sprintf(filename,"%s%s_%d",covparams.outdir,OUTFILE,k);
          if (fopen(filename, "r") != NULL){exit(1);}
          else {run_cov_clustering_ggl_real_bin(OUTFILE,covparams.outdir,thetamin,Ntheta,l,m,k);}
        }
        k=k+1;
      }
    }
  }

  if (hit==0)
  {
    sprintf(OUTFILE,"%s%s",covparams.outdir,covparams.filename);
    write_gglensing_zbins(OUTFILE);

    sprintf(OUTFILE,"%s%s.blocks",covparams.outdir,covparams.filename);
    F1 = fopen(OUTFILE,"w");
    fprintf(F1,"%d\n",k-1);
    fclose(F1);
  }

  printf("number of cov blocks for parallelization: %d\n",k-1);
  printf("-----------------\n");
  printf("PROGRAM EXECUTED\n");
  printf("-----------------\n");
  return 0;
}


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */
/*
 *
 *
 *       Filename:  stationary.c
 *
 *    Description:  Find the stationary state for the given generator
 *
 *        Version:  1.0
 *        Created:  05/05/2014 15:27:48
 *       Revision:  none
 *        License:  BSD
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 *
 */

#include "funcs.h"
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>

/*
 *      FUNCTION
 *         Name:  stationary
 *  Description:  Given the dissipator in Bloch form, reduce to a 3x3 problem and store
 *  			the stationary state in the 3x1 vector *X
 *
 * 			M X = 0
 *
 * 	        	|  0    0    0    0  |  | 1  |    0
 * 		        | M10  M11  M12  M13 |  | X1 |    0
 * 		        | M20  M21  M22  M23 |  | X2 | =  0
 * 			| M30  M31  M32  M33 |  | X3 |    0
 *
 *
 * 			A x = b
 *
 * 			| M11  M12  M13 |  | X1 |   | -M10 |
 * 			| M21  M22  M23 |  | X2 | = | -M20 |
 * 			| M31  M32  M33 |  | X3 |   | -M30 |
 */
int stationary ( const gsl_matrix* M, gsl_vector* stat_state )
{
	/* Store space for the stationary state */
	gsl_vector* req = gsl_vector_calloc ( 4 ) ;
	gsl_vector_set ( req, 0, 1 ) ;

	/* Copy the dissipator matrix in a temporary local matrix m
	 * (because the algorithm destroys it...) */
	gsl_matrix* m = gsl_matrix_calloc ( 4, 4 ) ;
	gsl_matrix_memcpy ( m, M ) ;

	/* Create a view of the spatial part of vector req */
	gsl_vector_view x = gsl_vector_subvector ( req, 1, 3 ) ;

	/* Create a submatrix view of the spatial part of m and a vector view
	 * of the spatial part of the 0-th column, which goes into -b in the system
	 * A x = b */
	gsl_matrix_view A = gsl_matrix_submatrix ( m, 1, 1, 3, 3 ) ;
	gsl_vector_view b = gsl_matrix_subcolumn ( m, 0, 1, 3 ) ;
	int status1 = gsl_vector_scale ( &b.vector, -1.0 ) ;

	/* Solve the system A x = b using Householder transformations.
	 * Changing the view x of req => also req is changed, in the spatial part */
	int status2 = gsl_linalg_HH_solve ( &A.matrix, &b.vector, &x.vector ) ;

	/* Set the returning value for the state stat_state */
	*stat_state = *req ;

	/* Free memory */
	gsl_matrix_free(m) ;

	return status1 + status2 ;
}		/* -----  end of function stationary  ----- */

/*
 *      FUNCTION
 *         Name:  J
 *  Description:  Ohmic spectral density (divided by alpha)
 *
 */
double J ( double w, double oc )
{
	double W = w*exp(-w/oc) ;
	return (W);
}		/* -----  end of function J  ----- */

/*
 *      FUNCTION
 *         Name:  polarization
 *  Description:  Polarization P through the CP formula
 *
 */
double polarization ( void* params, double oc )
{
	struct f_params* pars = (struct f_params*) params ;
	double Omega, omega_1, b ;
	Omega   = pars->Omega ;
	omega_1 = pars->omega_1 ;
	b       = pars->beta ;

	double P, num, den, add1, add2, Jplus, Jminus ;
	Jplus  = J(omega_1 + Omega, oc) ;
	Jminus = J(omega_1 - Omega, oc) ;
	add1 = POW_2(omega_1 - Omega)*Jplus ;
	add2 = POW_2(omega_1 + Omega)*Jminus ;
	num = add1 + add2 ;
	den = add1/(tanh((omega_1+Omega)*b/2)) +
		add2/(tanh((omega_1-Omega)*b/2)) ;
	P = num/den ;

	return (P);
}		/* -----  end of function polarization  ----- */


text:
/* interpolation/test.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stddef.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>

int
test_bsearch(void)
{
  double x_array[5] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
  size_t index_result;
  int status = 0;
  int s;

  /* check an interior point */
  index_result = gsl_interp_bsearch(x_array, 1.5, 0, 4);
  s = (index_result != 1);
  status += s;
  gsl_test (s, "simple bsearch");

  /* check that we get the last interval if x == last value */
  index_result = gsl_interp_bsearch(x_array, 4.0, 0, 4);
  s = (index_result != 3);
  status += s;
  gsl_test (s, "upper endpoint bsearch");

  /* check that we get the first interval if x == first value */
  index_result = gsl_interp_bsearch(x_array, 0.0, 0, 4);
  s = (index_result != 0);
  status += s;
  gsl_test (s, "lower endpoint bsearch");

  /* check that we get correct interior boundary behaviour */
  index_result = gsl_interp_bsearch(x_array, 2.0, 0, 4);
  s = (index_result != 2);
  status += s;
  gsl_test (s, "degenerate bsearch");

  /* check out of bounds above */
  index_result = gsl_interp_bsearch(x_array, 10.0, 0, 4);
  s = (index_result != 3);
  status += s;
  gsl_test (s, "out of bounds bsearch +");

  /* check out of bounds below */
  index_result = gsl_interp_bsearch(x_array, -10.0, 0, 4);
  s = (index_result != 0);
  status += s;
  gsl_test (s, "out of bounds bsearch -");

  return status;
}

typedef double TEST_FUNC (double);
typedef struct _xy_table xy_table;

struct _xy_table
  {
    double * x;
    double * y;
    size_t n;
  };

xy_table make_xy_table (double x[], double y[], size_t n);

xy_table
make_xy_table (double x[], double y[], size_t n)
{
  xy_table t;
  t.x = x;
  t.y = y;
  t.n = n;
  return t;
}

static int
test_interp (
  const xy_table * data_table,
  const gsl_interp_type * T,
  xy_table * test_table,
  xy_table * test_d_table,
  xy_table * test_i_table
  )
{
  int status = 0;
  size_t i;

  gsl_interp_accel *a = gsl_interp_accel_alloc ();
  gsl_interp *interp = gsl_interp_alloc (T, data_table->n);

  gsl_interp_init (interp, data_table->x, data_table->y, data_table->n);

  for (i = 0; i < test_table->n; i++)
    {
      double x = test_table->x[i];
      double y;
      double deriv;
      double integ;
      double diff_y, diff_deriv, diff_integ;
      gsl_interp_eval_e (interp, data_table->x, data_table->y, x, a, &y);
      gsl_interp_eval_deriv_e (interp, data_table->x, data_table->y, x, a, &deriv);
      gsl_interp_eval_integ_e (interp, data_table->x, data_table->y, 0.0, x, a, &integ);
      diff_y = y - test_table->y[i];
      diff_deriv = deriv - test_d_table->y[i];
      diff_integ = integ - test_i_table->y[i];
      if (fabs (diff_y) > 1.e-10 || fabs(diff_deriv) > 1.0e-10 || fabs(diff_integ) > 1.0e-10) {
	status++;
      }
    }

  gsl_interp_accel_free (a);
  gsl_interp_free (interp);

  return status;
}

static int
test_linear (void)
{
  int s;

  double data_x[4] = { 0.0, 1.0, 2.0, 3.0 };
  double data_y[4] = { 0.0, 1.0, 2.0, 3.0 };
  double test_x[6] = { 0.0, 0.5, 1.0, 1.5, 2.5, 3.0 };
  double test_y[6] = { 0.0, 0.5, 1.0, 1.5, 2.5, 3.0 };
  double test_dy[6] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
  double test_iy[6] = { 0.0, 0.125, 0.5, 9.0/8.0, 25.0/8.0, 9.0/2.0 };

  xy_table data_table = make_xy_table(data_x, data_y, 4);
  xy_table test_table = make_xy_table(test_x, test_y, 6);
  xy_table test_d_table = make_xy_table(test_x, test_dy, 6);
  xy_table test_i_table = make_xy_table(test_x, test_iy, 6);

  s = test_interp (&data_table, gsl_interp_linear, &test_table, &test_d_table, &test_i_table);
  gsl_test (s, "linear interpolation");
  return s;
}

static int
test_cspline (void)
{
  int s;

  double data_x[3] = { 0.0, 1.0, 2.0 };
  double data_y[3] = { 0.0, 1.0, 2.0 };
  double test_x[4] = { 0.0, 0.5, 1.0, 2.0 };
  double test_y[4] = { 0.0, 0.5, 1.0, 2.0 };
  double test_dy[4] = { 1.0, 1.0, 1.0, 1.0 };
  double test_iy[4] = { 0.0, 0.125, 0.5, 2.0 };

  xy_table data_table = make_xy_table(data_x, data_y, 3);
  xy_table test_table = make_xy_table(test_x, test_y, 4);
  xy_table test_d_table = make_xy_table(test_x, test_dy, 4);
  xy_table test_i_table = make_xy_table(test_x, test_iy, 4);

  s = test_interp (&data_table, gsl_interp_cspline, &test_table, &test_d_table, &test_i_table);
  gsl_test (s, "cspline interpolation");
  return s;
}

static int
test_akima (void)
{
  int s;

  double data_x[5] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
  double data_y[5] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
  double test_x[4] = { 0.0, 0.5, 1.0, 2.0 };
  double test_y[4] = { 0.0, 0.5, 1.0, 2.0 };
  double test_dy[4] = { 1.0, 1.0, 1.0, 1.0 };
  double test_iy[4] = { 0.0, 0.125, 0.5, 2.0 };

  xy_table data_table = make_xy_table(data_x, data_y, 5);
  xy_table test_table = make_xy_table(test_x, test_y, 4);
  xy_table test_d_table = make_xy_table(test_x, test_dy, 4);
  xy_table test_i_table = make_xy_table(test_x, test_iy, 4);

  s = test_interp (&data_table, gsl_interp_akima, &test_table, &test_d_table, &test_i_table);
  gsl_test (s, "akima interpolation");
  return s;
}

int
main (int argc, char **argv)
{
  int status = 0;

  argc = 0;    /* prevent warnings about unused parameters */
  argv = 0;

  status += test_bsearch();
  status += test_linear();
  status += test_cspline();
  status += test_akima();

  exit (gsl_test_summary());
}