text:
/*
 * BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

/*!
 *  \file reconstruct_con.c
 *  \brief reconstruct continuum.
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_randist.h>

#include "brains.h"

void *best_model_con;   /*!< best model */
void *best_model_std_con;  /*!< standard deviation of the best model */

/*!
 *  this function does postprocess for continuum.
 */
void postprocess_con()
{
  char posterior_sample_file[BRAINS_MAX_STR_LENGTH];
  double *pm, *pmstd;
  int num_ps, i, j;
  void *posterior_sample, *post_model;
  int size_of_modeltype = num_params * sizeof(double);

  best_model_con = malloc(size_of_modeltype);
  best_model_std_con = malloc(size_of_modeltype);

  if(thistask == roottask)
  {
    char fname[200];
    FILE *fp, *fcon;

    /* get file name of posterior sample file */
    //get_posterior_sample_file(dnest_options_file, posterior_sample_file);
    dnest_get_posterior_sample_file(posterior_sample_file);

    /* open file for posterior sample */
    fp = fopen(posterior_sample_file, "r");
    if(fp == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s.\n", posterior_sample_file);
      exit(0);
    }

    /* open file for continuum reconstruction */
    sprintf(fname, "%s/%s", parset.file_dir, "data/con_rec.txt");
    fcon = fopen(fname, "w");
    if(fcon == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
      exit(0);
    }

    /* read number of points in posterior sample */
    if(fscanf(fp, "# %d", &num_ps) < 1)
    {
      fprintf(stderr, "# Error: Cannot read file %s.\n", posterior_sample_file);
      exit(0);
    }
    printf("# Number of points in posterior sample: %d\n", num_ps);

    post_model = malloc(size_of_modeltype);
    posterior_sample = malloc(num_ps * size_of_modeltype);

    which_parameter_update = -1;
    which_particle_update = 0;

    for(i=0; i<num_ps; i++)
    {
      for(j=0; j<num_params; j++)
      {
        if(fscanf(fp, "%lf", (double *)post_model + j) < 1)
        {
          fprintf(stderr, "# Error: Cannot read file %s.\n", posterior_sample_file);
          exit(0);
        }
      }
      fscanf(fp, "\n");

      memcpy(posterior_sample+i*size_of_modeltype, post_model, size_of_modeltype);

      /* not output reconstruction when flag_dim = -1 */
      if(parset.flag_dim != -1)
      {
        calculate_con_from_model_semiseparable(post_model);

        if(gsl_rng_uniform(gsl_r) < 1.0)
        {
          for(j=0; j<parset.n_con_recon; j++)
          {
            fprintf(fcon, "%e %e %f\n", Tcon[j]*(1.0+parset.redshift), Fcon[j]/con_scale, Fcerrs[j]/con_scale);
          }
          fprintf(fcon, "\n");
        }
      }
    }

    fclose(fp);
    fclose(fcon);

    /* calcaulte mean and standard deviation of posterior samples. */
    pm = (double *)best_model_con;
    pmstd = (double *)best_model_std_con;
    for(j=0; j<num_params; j++)
    {
      pm[j] = pmstd[j] = 0.0;
    }
    for(i=0; i<num_ps; i++)
    {
      for(j =0; j<num_params; j++)
        pm[j] += *((double *)posterior_sample + i*num_params + j );
    }

    for(j=0; j<num_params; j++)
      pm[j] /= num_ps;

    for(i=0; i<num_ps; i++)
    {
      for(j=0; j<num_params; j++)
        pmstd[j] += pow( *((double *)posterior_sample + i*num_params + j ) - pm[j], 2.0 );
    }

    for(j=0; j<num_params; j++)
    {
      if(num_ps > 1)
        pmstd[j] = sqrt(pmstd[j]/(num_ps-1.0));
      else
        pmstd[j] = 0.0;
    }

    for(j = 0; j<num_params_var; j++)
      printf("Best params %d %f +- %f\n", j, *((double *)best_model_con + j),
                                             *((double *)best_model_std_con+j) );

    free(post_model);
    free(posterior_sample);
  }
  return;
}

/*!
 * This function runs dnest sampling and postprocessing, reconstructs the continuum
 *  using the best estimates of model parameters.
 */
void reconstruct_con()
{
  int i, argc=0;
  char **argv;

  // configure restart of dnest
  argv = malloc(9*sizeof(char *));
  for(i=0; i<9; i++)
  {
    argv[i] = malloc(BRAINS_MAX_STR_LENGTH*sizeof(char));
  }
  //setup argc and argv
  strcpy(argv[argc++], "dnest");
  strcpy(argv[argc++], "-s");
  strcpy(argv[argc], parset.file_dir);
  strcat(argv[argc++], "/data/restart_dnest.txt");

  if(parset.flag_restart == 1 && parset.flag_dim == 0)
  {
    strcpy(argv[argc++], "-r");
    strcpy(argv[argc], parset.file_dir);
    strcat(argv[argc], "/");
    strcat(argv[argc++], "data/restart_dnest.txt");
  }
  if(parset.flag_postprc == 1 || (parset.flag_dim !=0 && parset.flag_restart == 1))
  {
    strcpy(argv[argc++], "-p");
  }
  if(parset.flag_temp == 1)
  {
    sprintf(argv[argc++], "-t%f", parset.temperature);
  }
  if(parset.flag_sample_info == 1)
  {
    strcpy(argv[argc++], "-c");
  }

  //level-dependent sampling
  {
    strcpy(argv[argc++], "-l");
  }

  reconstruct_con_init();
  logz_con = dnest_con(argc, argv);

  if(parset.flag_para_name != 1)
  {
    if(parset.flag_exam_prior != 1)
    {
      postprocess_con();

      if(thistask == roottask)
      {
        which_parameter_update = -1;
        which_particle_update = 0;

        calculate_con_from_model_semiseparable(best_model_con);

        FILE *fp;
        char fname[200];
        int i;
        sprintf(fname, "%s/%s", parset.file_dir, parset.pcon_out_file);
        fp = fopen(fname, "w");
        if(fp == NULL)
        {
          fprintf(stderr, "# Error: Cannot open file %s\n", fname);
          exit(-1);
        }

        for(i=0; i<parset.n_con_recon; i++)
        {
          fprintf(fp, "%e %e %e\n", Tcon[i]*(1.0+parset.redshift), Fcon[i] / con_scale, Fcerrs[i]/con_scale);
        }
        fclose(fp);

        memcpy(var_param, best_model_con, num_params_var*sizeof(double));
        memcpy(var_param_std, best_model_std_con, num_params_var*sizeof(double));
      }
    }
    else
    {
      for(i=0; i<num_params_var; i++)
      {
        var_param[i] = 0.5*(par_range_model[i][0] + par_range_model[i][1]);
        var_param_std[i] =0.5*(par_range_model[i][1] - par_range_model[i][0])/2.35;
      }
    }

    //use the posterior mean and standard variance as the prior for the 1d and 2d RM.
    //only for variability parameters.
    MPI_Bcast(var_param, num_params_var, MPI_DOUBLE, roottask, MPI_COMM_WORLD);
    MPI_Bcast(var_param_std, num_params_var, MPI_DOUBLE, roottask, MPI_COMM_WORLD);
  }

  reconstruct_con_end();

  //clear up argv
  for(i=0; i<9; i++)
  {
    free(argv[i]);
  }
  free(argv);

  return;
}

/*!
 *  this function calculates continuum light curves from model parameters
 *  and stores it into Fcon.
 */
void calculate_con_from_model(const void *model)
{
  double *ybuf, *y, *yu, *Cq, *yq, *yuq;
  double syserr;

  double *pm = (double *)model;
  double sigma, tau, alpha;
  int i, info;

  syserr = (exp(pm[0]) - 1.0) * con_error_mean;  // systematic error
  tau = exp(pm[2]);
  sigma = exp(pm[1]) * sqrt(tau);
  alpha = 1.0;

  ybuf = workspace;
  y = ybuf + n_con_data*nq;
  Cq = y + n_con_data;
  yq = Cq + nq*nq;
  yu = yq + nq;
  yuq = yu + parset.n_con_recon;

  set_covar_Pmat_data(sigma, tau, alpha, syserr);
  set_covar_Umat(sigma, tau, alpha);

  inverse_pomat(PCmat_data, n_con_data, &info);

  // Cq^-1 = L^TxC^-1xL
  multiply_mat_MN(PCmat_data, Larr_data, ybuf, n_con_data, nq, n_con_data);
  multiply_mat_MN_transposeA(Larr_data, ybuf, Cq, nq, nq, n_con_data);

  // L^TxC^-1xy
  multiply_matvec(PCmat_data, Fcon_data, n_con_data, ybuf);
  multiply_mat_MN_transposeA(Larr_data, ybuf, yq, nq, 1, n_con_data);

  // (hat q) = Cqx(L^TxC^-1xy)
  inverse_pomat(Cq, nq, &info);
  multiply_mat_MN(Cq, yq, ybuf, nq, 1, nq);

  // q = uq + (hat q)
  Chol_decomp_L(Cq, nq, &info);
  multiply_matvec(Cq, &pm[3], nq, yq);
  for(i=0; i<nq; i++)
    yq[i] += ybuf[i];

  // y = yc - Lxq
  multiply_matvec_MN(Larr_data, n_con_data, nq, yq, ybuf);
  for(i=0; i<n_con_data; i++)
  {
    y[i] = Fcon_data[i] - ybuf[i];
  }

  // (hat s) = SxC^-1xy
  multiply_matvec(PCmat_data, y, n_con_data, ybuf);
  multiply_matvec_MN(USmat, parset.n_con_recon, n_con_data, ybuf, Fcon);

  // SxC^-1xS^T
  multiply_mat_MN(USmat, PCmat_data, PEmat1, parset.n_con_recon, n_con_data, n_con_data);
  multiply_mat_MN_transposeB(PEmat1, USmat, PEmat2, parset.n_con_recon, parset.n_con_recon, n_con_data);

  // set DRW covariance at reconstructed points
  set_covar_Pmat(sigma, tau, alpha);

  // S - SxC^-1xS^T; covariance of reconstructed points
  for(i=0; i<parset.n_con_recon*parset.n_con_recon; i++)
    PQmat[i] = PSmat[i] - PEmat2[i];

  // assign errors
  for(i=0; i<parset.n_con_recon; i++)
  {
    Fcerrs[i] = sqrt(PQmat[i*parset.n_con_recon + i]);
  }

  Chol_decomp_L(PQmat, parset.n_con_recon, &info);
  multiply_matvec(PQmat, &pm[num_params_var], parset.n_con_recon, yu);

  // add back long-term trend of continuum
  multiply_matvec_MN(Larr_rec, parset.n_con_recon, nq, yq, yuq);

  for(i=0; i<parset.n_con_recon; i++)
  {
    Fcon[i] += yu[i] + yuq[i];
  }

  return;
}

/*!
 *  this function calculates continuum light curves from model parameters
 *  and stores it into Fcon.
 *
 *  using fast algorithm for semiseparable matrices
 */
void calculate_con_from_model_semiseparable(const void *model)
{
  double *Lbuf, *ybuf, *y, *Cq, *yq, *W, *D, *phi, *u, *v;
  double syserr;

  double *pm = (double *)model;
  double sigma, sigma2, tau, alpha;
  int i, info;

  syserr = (exp(pm[0]) - 1.0) * con_error_mean;  // systematic error
  tau = exp(pm[2]);
  sigma = exp(pm[1]) * sqrt(tau);
  sigma2 = sigma*sigma;
  alpha = 1.0;

  Lbuf = workspace;
  ybuf = Lbuf + n_con_data*nq;
  y = ybuf + n_con_max;
  Cq = y + n_con_max;
  yq = Cq + nq*nq;
  W = yq + nq;
  D = W + n_con_max;
  phi = D + n_con_max;
  u = phi + n_con_max;
  v = u + parset.n_con_recon;

  set_covar_Umat(sigma, tau, alpha);

  compute_semiseparable_drw(Tcon_data, n_con_data, sigma2, 1.0/tau, Fcerrs_data, syserr, W, D, phi);
  // Cq^-1 = L^TxC^-1xL
  multiply_mat_semiseparable_drw(Larr_data, W, D, phi, n_con_data, nq, sigma2, Lbuf);
  multiply_mat_MN_transposeA(Larr_data, Lbuf, Cq, nq, nq, n_con_data);

  // L^TxC^-1xy
  multiply_matvec_semiseparable_drw(Fcon_data, W, D, phi, n_con_data, sigma2, ybuf);
  multiply_mat_MN_transposeA(Larr_data, ybuf, yq, nq, 1, n_con_data);

  // (hat q) = Cqx(L^TxC^-1xy)
  inverse_pomat(Cq, nq, &info);
  multiply_mat_MN(Cq, yq, ybuf, nq, 1, nq);

  // q = uq + (hat q)
  Chol_decomp_L(Cq, nq, &info);
  multiply_matvec(Cq, &pm[3], nq, yq);
  for(i=0; i<nq; i++)
    yq[i] += ybuf[i];

  // y = yc - Lxq
  multiply_matvec_MN(Larr_data, n_con_data, nq, yq, ybuf);
  for(i=0; i<n_con_data; i++)
  {
    y[i] = Fcon_data[i] - ybuf[i];
  }

  // (hat s) = SxC^-1xy
  multiply_matvec_semiseparable_drw(y, W, D, phi, n_con_data, sigma2, ybuf);
  multiply_matvec_MN(USmat, parset.n_con_recon, n_con_data, ybuf, Fcon);

  // SxC^-1xS^T
  multiply_mat_transposeB_semiseparable_drw(USmat, W, D, phi, n_con_data, parset.n_con_recon, sigma2, PEmat1);
  multiply_mat_MN(USmat, PEmat1, PEmat2, parset.n_con_recon, parset.n_con_recon, n_con_data);

  // set DRW covariance at reconstructed points
  set_covar_Pmat(sigma, tau, alpha);

  // S - SxC^-1xS^T; covariance of reconstructed points
  for(i=0; i<parset.n_con_recon*parset.n_con_recon; i++)
    PQmat[i] = PSmat[i] - PEmat2[i];

  // assign errors
  for(i=0; i<parset.n_con_recon; i++)
  {
    Fcerrs[i] = sqrt(PQmat[i*parset.n_con_recon + i]);
  }

  Chol_decomp_L(PQmat, parset.n_con_recon, &info);
  multiply_matvec(PQmat, &pm[num_params_var], parset.n_con_recon, y);

  // add back long-term trend of continuum
  multiply_matvec_MN(Larr_rec, parset.n_con_recon, nq, yq, ybuf);

  for(i=0; i<parset.n_con_recon; i++)
  {
    Fcon[i] += y[i] + ybuf[i];
  }
  return;
}

/*!
 * This function calculate continuum light curves from variability parameters.
 */
void reconstruct_con_from_varmodel(double sigma_hat, double tau, double alpha, double syserr)
{
  double *Lbuf, *ybuf, *y, *yq, *yuq, *Cq, sigma;
  int i, info;
  double *PEmat3, *PEmat4;

  sigma = sigma_hat * sqrt(tau);

  Lbuf = workspace;
  ybuf = Lbuf + n_con_data * nq;
  y = ybuf + n_con_data;
  Cq = y + n_con_data;
  yq = Cq + nq*nq;
  yuq = yq + nq;

  PEmat3 = malloc(parset.n_con_recon * nq * sizeof(double));
  PEmat4 = malloc(parset.n_con_recon * parset.n_con_recon * sizeof(double));

  set_covar_Pmat_data(sigma, tau, alpha, syserr);
  set_covar_Umat(sigma, tau, alpha);

  inverse_pomat(PCmat_data, n_con_data, &info);

  multiply_mat_MN(PCmat_data, Larr_data, Lbuf, n_con_data, nq, n_con_data);
  multiply_mat_MN_transposeA(Larr_data, Lbuf, Cq, nq, nq, n_con_data);

  multiply_matvec(PCmat_data, Fcon_data, n_con_data, ybuf);
  multiply_mat_MN_transposeA(Larr_data, ybuf, yq, nq, 1, n_con_data);

  inverse_pomat(Cq, nq, &info);
  multiply_mat_MN(Cq, yq, ybuf, nq, 1, nq);

  memcpy(yq, ybuf, nq*sizeof(double));

  multiply_matvec_MN(Larr_data, n_con_data, nq, yq, ybuf);
  for(i=0; i<n_con_data; i++)
  {
    y[i] = Fcon_data[i] - ybuf[i];
  }

  multiply_matvec(PCmat_data, y, n_con_data, ybuf);
  multiply_matvec_MN(USmat, parset.n_con_recon, n_con_data, ybuf, Fcon);

  multiply_matvec_MN(Larr_rec, parset.n_con_recon, nq, yq, yuq);

  for(i=0; i<parset.n_con_recon; i++)
    Fcon[i] += yuq[i];

  // SxC^-1xS^T
  multiply_mat_MN(USmat, PCmat_data, PEmat1, parset.n_con_recon, n_con_data, n_con_data);
  multiply_mat_MN_transposeB(PEmat1, USmat, PEmat2, parset.n_con_recon, parset.n_con_recon, n_con_data);

  multiply_mat_MN(PEmat1, Larr_data, PEmat3, parset.n_con_recon, nq, n_con_data);
  for(i=0; i<parset.n_con_recon*nq; i++)PEmat3[i] -= Larr_rec[i];
  multiply_mat_MN(PEmat3, Cq, PEmat1, parset.n_con_recon, nq, nq);
  multiply_mat_MN_transposeB(PEmat1, PEmat3, PEmat4, parset.n_con_recon, parset.n_con_recon, nq);

  for(i=0; i<parset.n_con_recon; i++)
  {
    Fcerrs[i] = sqrt(sigma*sigma - PEmat2[i*parset.n_con_recon+i] + PEmat4[i*parset.n_con_recon + i]);
  }

  free(PEmat3);
  free(PEmat4);
  return;
}

/*!
 * this function calculates likelihood
 */
double prob_con_variability(const void *model)
{
  double prob = 0.0;
  int i, info;
  double *pm = (double *)model;
  double tau, sigma, alpha, lndet, syserr;
  double *Lbuf, *ybuf, *y, *yq, *Cq;

  which_particle_update = dnest_get_which_particle_update();

  if( which_parameter_update < num_params_var)
  {
    syserr = (exp(pm[0])-1.0)*con_error_mean;
    tau = exp(pm[2]);
    sigma = exp(pm[1]) * sqrt(tau);
    alpha = 1.0;

    Lbuf = workspace;
    ybuf = Lbuf + n_con_data*nq;
    y = ybuf + n_con_data;
    yq = y + n_con_data;
    Cq = yq + nq;

    set_covar_Pmat_data(sigma, tau, alpha, syserr);
    memcpy(IPCmat_data, PCmat_data, n_con_data*n_con_data*sizeof(double));

    lndet = lndet_mat(PCmat_data, n_con_data, &info);

    inverse_pomat(IPCmat_data, n_con_data, &info); /* calculate C^-1 */

    /* calculate L^T*C^-1*L */
    multiply_mat_MN(IPCmat_data, Larr_data, Lbuf, n_con_data, nq, n_con_data);
    multiply_mat_MN_transposeA(Larr_data, Lbuf, Cq, nq, nq, n_con_data);

    /* calculate L^T*C^-1*y */
    multiply_matvec(IPCmat_data, Fcon_data, n_con_data, ybuf);
    multiply_mat_MN_transposeA(Larr_data, ybuf, yq, nq, 1, n_con_data);

    /* calculate (L^T*C^-1*L)^-1 * L^T*C^-1*y */
    inverse_pomat(Cq, nq, &info);
    multiply_mat_MN(Cq, yq, ybuf, nq, 1, nq);

    Chol_decomp_L(Cq, nq, &info);
    multiply_matvec(Cq, &pm[3], nq, yq);
    for(i=0; i<nq; i++)
      yq[i] += ybuf[i];

    multiply_matvec_MN(Larr_data, n_con_data, nq, yq, ybuf);
    for(i=0; i<n_con_data; i++)
    {
      y[i] = Fcon_data[i] - ybuf[i];
    }

    multiply_matvec(IPCmat_data, y, n_con_data, ybuf);
    prob = -0.5 * cblas_ddot(n_con_data, y, 1, ybuf, 1);
    prob += -0.5*lndet;

    prob_con_particles_perturb[which_particle_update] = prob;
  }
  else
  {
    prob = prob_con_particles[which_particle_update];
  }

  return prob;
}

/*!
 * this function calculates likelihood
 *
 * using fast algorithm for semiseparable matrices
 */
double prob_con_variability_semiseparable(const void *model)
{
  double prob = 0.0;
  int i, info;
  double *pm = (double *)model;
  double tau, sigma2, lndet, syserr;
  double *Lbuf, *ybuf, *y, *yq, *Cq, *W, *D, *phi;

  which_particle_update = dnest_get_which_particle_update();

  if( which_parameter_update < num_params_var)
  {

    syserr = (exp(pm[0])-1.0)*con_error_mean;
    tau = exp(pm[2]);
    sigma2 = exp(2.0*pm[1]) * tau;

    Lbuf = workspace;
    ybuf = Lbuf + n_con_data*nq;
    y = ybuf + n_con_data;
    yq = y + n_con_data;
    Cq = yq + nq;
    W = Cq + nq*nq;
    D = W + n_con_data;
    phi = D + n_con_data;

    compute_semiseparable_drw(Tcon_data, n_con_data, sigma2, 1.0/tau, Fcerrs_data, syserr, W, D, phi);
    lndet = 0.0;
    for(i=0; i<n_con_data; i++)
      lndet += log(D[i]);

    /* calculate L^T*C^-1*L */
    multiply_mat_semiseparable_drw(Larr_data, W, D, phi, n_con_data, nq, sigma2, Lbuf);
    multiply_mat_MN_transposeA(Larr_data, Lbuf, Cq, nq, nq, n_con_data);

    /* calculate L^T*C^-1*y */
    multiply_matvec_semiseparable_drw(Fcon_data, W, D, phi, n_con_data, sigma2, ybuf);
    multiply_mat_MN_transposeA(Larr_data, ybuf, yq, nq, 1, n_con_data);

    /* calculate (L^T*C^-1*L)^-1 * L^T*C^-1*y */
    inverse_pomat(Cq, nq, &info);
    multiply_mat_MN(Cq, yq, ybuf, nq, 1, nq);

    Chol_decomp_L(Cq, nq, &info);
    multiply_matvec(Cq, &pm[3], nq, yq);
    for(i=0; i<nq; i++)
      yq[i] += ybuf[i];

    multiply_matvec_MN(Larr_data, n_con_data, nq, yq, ybuf);
    for(i=0; i<n_con_data; i++)
    {
      y[i] = Fcon_data[i] - ybuf[i];
    }

    /* y^T x C^-1 x y*/
    multiply_matvec_semiseparable_drw(y, W, D, phi, n_con_data, sigma2, ybuf);
    prob = -0.5 * cblas_ddot(n_con_data, y, 1, ybuf, 1);
    prob += -0.5*lndet;

    prob_con_particles_perturb[which_particle_update] = prob;
  }
  else
  {
    prob = prob_con_particles[which_particle_update];
  }

  return prob;
}

/*!
 * this function calculates likelihood at initital step.
 */
double prob_con_variability_initial(const void *model)
{
  double prob = 0.0;
  int i, info;
  double *pm = (double *)model;
  double tau, sigma, alpha, lndet, syserr;
  double *Lbuf, *ybuf, *y, *yq, *Cq;

  syserr = (exp(pm[0])-1.0)*con_error_mean;
  tau = exp(pm[2]);
  sigma = exp(pm[1]) * sqrt(tau);
  alpha = 1.0;

  Lbuf = workspace;
  ybuf = Lbuf + n_con_data*nq;
  y = ybuf + n_con_data;
  yq = y + n_con_data;
  Cq = yq + nq;

  set_covar_Pmat_data(sigma, tau, alpha, syserr);
  memcpy(IPCmat_data, PCmat_data, n_con_data*n_con_data*sizeof(double));

  lndet = lndet_mat(PCmat_data, n_con_data, &info);

  inverse_pomat(IPCmat_data, n_con_data, &info);

  multiply_mat_MN(IPCmat_data, Larr_data, Lbuf, n_con_data, nq, n_con_data);
  multiply_mat_MN_transposeA(Larr_data, Lbuf, Cq, nq, nq, n_con_data);

  multiply_matvec(IPCmat_data, Fcon_data, n_con_data, ybuf);
  multiply_mat_MN_transposeA(Larr_data, ybuf, yq, nq, 1, n_con_data);

  inverse_pomat(Cq, nq, &info);
  multiply_mat_MN(Cq, yq, ybuf, nq, 1, nq);

  Chol_decomp_L(Cq, nq, &info);
  multiply_matvec(Cq, &pm[3], nq, yq);
  for(i=0; i<nq; i++)
    yq[i] += ybuf[i];

  multiply_matvec_MN(Larr_data, n_con_data, nq, yq, ybuf);
  for(i=0; i<n_con_data; i++)
  {
    y[i] = Fcon_data[i] - ybuf[i];
  }

  multiply_matvec(IPCmat_data, y, n_con_data, ybuf);
  prob = -0.5 * cblas_ddot(n_con_data, y, 1, ybuf, 1);
  prob += -0.5*lndet;

  prob_con_particles[which_particle_update] = prob;
  return prob;
}

/*!
 * this function calculates likelihood at initital step.
 *
 * using fast algorithm for semiseparable matrices
 */
double prob_con_variability_initial_semiseparable(const void *model)
{
  double prob = 0.0;
  int i, info;
  double *pm = (double *)model;
  double tau, sigma2, lndet, syserr;
  double *Lbuf, *ybuf, *y, *yq, *Cq, *W, *D, *phi;

  syserr = (exp(pm[0])-1.0)*con_error_mean;
  tau = exp(pm[2]);
  sigma2 = exp(2.0*pm[1]) * tau;

  Lbuf = workspace;
  ybuf = Lbuf + n_con_data*nq;
  y = ybuf + n_con_data;
  yq = y + n_con_data;
  Cq = yq + nq;
  W = Cq + nq*nq;
  D = W + n_con_data;
  phi = D + n_con_data;

  compute_semiseparable_drw(Tcon_data, n_con_data, sigma2, 1.0/tau, Fcerrs_data, syserr, W, D, phi);
  lndet = 0.0;
  for(i=0; i<n_con_data; i++)
  lndet += log(D[i]);

  /* calculate L^T*C^-1*L */
  multiply_mat_semiseparable_drw(Larr_data, W, D, phi, n_con_data, nq, sigma2, Lbuf);
  multiply_mat_MN_transposeA(Larr_data, Lbuf, Cq, nq, nq, n_con_data);

  /* calculate L^T*C^-1*y */
  multiply_matvec_semiseparable_drw(Fcon_data, W, D, phi, n_con_data, sigma2, ybuf);
  multiply_mat_MN_transposeA(Larr_data, ybuf, yq, nq, 1, n_con_data);

  inverse_pomat(Cq, nq, &info);
  multiply_mat_MN(Cq, yq, ybuf, nq, 1, nq);

  Chol_decomp_L(Cq, nq, &info);
  multiply_matvec(Cq, &pm[3], nq, yq);
  for(i=0; i<nq; i++)
    yq[i] += ybuf[i];

  multiply_matvec_MN(Larr_data, n_con_data, nq, yq, ybuf);
  for(i=0; i<n_con_data; i++)
  {
    y[i] = Fcon_data[i] - ybuf[i];
  }

  /* y^T x C^-1 x y*/
  multiply_matvec_semiseparable_drw(y, W, D, phi, n_con_data, sigma2, ybuf);
  prob = -0.5 * cblas_ddot(n_con_data, y, 1, ybuf, 1);
  prob += -0.5*lndet;

  prob_con_particles[which_particle_update] = prob;
  return prob;
}

/*!
 * this function sets the covariance matrix at time points for reconstruction.
 * store into PSmat.
 */
void set_covar_Pmat(double sigma, double tau, double alpha)
{
  double t1, t2;
  int i, j;

  for(i=0; i<parset.n_con_recon; i++)
  {
    t1 = Tcon[i];
    for(j=0; j<=i; j++)
    {
      t2 = Tcon[j];
      PSmat[i*parset.n_con_recon+j] = sigma*sigma* exp (- pow (fabs(t1-t2) / tau, alpha));
      PSmat[j*parset.n_con_recon+i] = PSmat[i*parset.n_con_recon+j];
    }
  }
  return;
}

/*!
 * this function sets the covariance matrix at data time points
 */
void set_covar_Pmat_data(double sigma, double tau, double alpha, double syserr)
{
  double t1, t2, err2;
  int i, j;

  for(i=0; i<n_con_data; i++)
  {
    t1 = Tcon_data[i];

    for(j=0; j<i; j++)
    {
      t2 = Tcon_data[j];
      PSmat_data[i*n_con_data+j] = sigma*sigma* exp (- pow (fabs(t1-t2) / tau, alpha));
      PSmat_data[j*n_con_data+i] = PSmat_data[i*n_con_data+j];

      //PNmat_data[i*n_con_data+j] = PNmat_data[j*n_con_data+i] = 0.0;

      PCmat_data[i*n_con_data+j] = PCmat_data[j*n_con_data+i] = PSmat_data[i*n_con_data+j];
    }

    PSmat_data[i*n_con_data+i] = sigma * sigma;
    //PNmat_data[i*n_con_data+i] = Fcerrs_data[i]*Fcerrs_data[i] + syserr*syserr;
    //PCmat_data[i*n_con_data+i] = PSmat_data[i*n_con_data+i] + PNmat_data[i*n_con_data+i];
    err2 = Fcerrs_data[i]*Fcerrs_data[i] + syserr*syserr;
    PCmat_data[i*n_con_data+i] = PSmat_data[i*n_con_data+i] + err2;
  }
  return;
}

/*!
 * this function sets the covariance matrix at time of data points and reconstruction points
 */
void set_covar_Umat(double sigma, double tau, double alpha)
{
  double t1, t2;
  int i, j;

  for(i=0; i<parset.n_con_recon; i++)
  {
    t1 = Tcon[i];
    for(j=0; j<n_con_data; j++)
    {
      t2 = Tcon_data[j];
      USmat[i*n_con_data+j] = sigma*sigma * exp (- pow (fabs(t1-t2) / tau, alpha) );
    }
  }
  return;
}

/*!
 * this function initializes the continuum reconstruction.
 */
void reconstruct_con_init()
{
  int i, j;
  double dT, Tspan;

  Tspan = Tcon_data[n_con_data -1] - Tcon_data[0];

  /* set time array for continuum */
  if(parset.time_back > 0.0)
    Tcon_min = Tcon_data[0] - parset.time_back;
  else
    Tcon_min = Tcon_data[0] - fmax(0.05*Tspan, 10.0);

  Tcon_max = Tcon_data[n_con_data-1] + fmax(0.05*Tspan, 10.0);

  if(thistask == roottask)
    printf("Tcon_min_max: %f %f\n", Tcon_min - Tcon_data[0], Tcon_max - Tcon_data[n_con_data-1]);

  dT = (Tcon_max - Tcon_min)/(parset.n_con_recon -1);

  for(i=0; i<parset.n_con_recon; i++)
  {
    Tcon[i] = Tcon_min + i*dT;
  }

  /* set Larr_rec */
  for(i=0;i<parset.n_con_recon;i++)
  {
    Larr_rec[i*nq + 0]=1.0;
    for(j=1; j<nq; j++)
      Larr_rec[i*nq + j] = pow(Tcon[i], j);
  }

  sprintf(dnest_options_file, "%s/%s", parset.file_dir, "src/OPTIONSCON");
  if(thistask == roottask)
  {
    get_num_particles(dnest_options_file);
  }
  MPI_Bcast(&parset.num_particles, 1, MPI_INT, roottask, MPI_COMM_WORLD);

  Fcon = malloc(parset.n_con_recon * sizeof(double *));
  prob_con_particles = malloc(parset.num_particles * sizeof(double));
  prob_con_particles_perturb = malloc(parset.num_particles * sizeof(double));

  return;
}

/*!
 * this function finalize the continuum reconstruction.
 */
void reconstruct_con_end()
{
  int i;

  for(i=0; i<num_params; i++)
  {
    free(par_range_model[i]);
    free(par_prior_gaussian[i]);
  }
  free(par_range_model);
  free(par_prior_gaussian);
  free(par_prior_model);

  free(Fcon);

  free(prob_con_particles);
  free(prob_con_particles_perturb);

  free(best_model_con);
  free(best_model_std_con);

  return;
}

/*!
 *  create continuum from random number generate based on damped random walk model.
 */

void create_con_from_random(double sigma_hat, double tau, double alpha, double syserr)
{
  int i, info;
  double *Prandvec;
  double sigma = sigma_hat * sqrt(tau);
  double mean;

  Prandvec = malloc(parset.n_con_recon*sizeof(double));

  set_covar_Pmat(sigma, tau, alpha);

  for(i=0; i<parset.n_con_recon; i++)
  {
    Prandvec[i] = gsl_ran_ugaussian(gsl_r);
  }

  Chol_decomp_L(PSmat, parset.n_con_recon, &info);

  multiply_matvec(PSmat, Prandvec, parset.n_con_recon, Fcon);

  mean = 0.0;
  for(i=0; i<parset.n_con_recon; i++)
  {
    Fcon[i] += sigma*3.0; /* shift light curve to have positive fluxes*/
    mean += Fcon[i];
  }
  mean /= parset.n_con_recon;

  for(i=0; i<parset.n_con_recon; i++)
  {
    Fcon[i] /= mean;
  }

  for(i=0; i<parset.n_con_recon; i++)
  {
    Fcon[i] +=  gsl_ran_ugaussian(gsl_r) * con_error_mean;
    Fcerrs[i] = con_error_mean;
  }

  free(Prandvec);
  return;
}

max_stars_repo_name: LiyrAstroph/BRAINS

max_issues_repo_name: LiyrAstroph/BRAINS

max_forks_repo_name: LiyrAstroph/BRAINS


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_core_typeutil_h_
#define SQ_INCLUDE_GUARD_core_typeutil_h_

#include <concepts>
#include <cstddef>
#include <gsl/gsl>
#include <range/v3/range/concepts.hpp>
#include <stdexcept>
#include <string>
#include <system_error>
#include <type_traits>
#include <variant>

// The following macros can't be replaced with constexpr functions or constants
// so tell clang-tidy not to complain

// NOLINTNEXTLINE(cppcoreguidelines-macro-usage)
#define SQ_FWD(x) static_cast<decltype(x) &&>(x)

// NOLINTNEXTLINE(cppcoreguidelines-macro-usage)
#define SQ_ND [[nodiscard]]

// NOLINTNEXTLINE(cppcoreguidelines-macro-usage)
#define SQ_MU [[maybe_unused]]

namespace sq {

/**
 * Concept for types that can be dumped to a std::ostream.
 */
template <typename T> concept Printable = requires(T x, std::ostream &os) {
  os << x;
};

///@{
/**
 * Get whether the type T is a specialization of std::variant.
 */
template <typename T> struct IsVariant : std::false_type {};

template <typename... Types>
struct IsVariant<std::variant<Types...>> : std::true_type {};

template <typename T> inline constexpr bool is_variant_v = IsVariant<T>::value;
///@}

///@{
/**
 * Get whether the type T is one of the alternative types of the variant V.
 */
template <typename T, typename V>
requires is_variant_v<V> struct IsAlternative : std::false_type {};

template <typename T, typename... Types>
struct IsAlternative<T, std::variant<Types...>>
    : std::disjunction<std::is_same<T, Types>...> {};

template <typename T, typename V>
inline constexpr bool is_alternative_v = IsAlternative<T, V>::value;
///@}

/**
 * Concept for types that are alternatives of the variant V.
 */
template <typename T, typename V> concept Alternative = is_alternative_v<T, V>;

///@{
/**
 * Get whether T can be converted to one of the variant V's alternatives.
 */
template <typename T, typename V>
requires is_variant_v<V> struct IsConvertibleToAlternative : std::false_type {};

template <typename T, typename... Types>
struct IsConvertibleToAlternative<T, std::variant<Types...>>
    : std::disjunction<std::is_convertible<T, Types>...> {};

template <typename T, typename V>
inline constexpr bool is_convertible_to_alternative_v =
    IsConvertibleToAlternative<T, V>::value;
///@}

/**
 * Concept for types that are convertible to one of the variant V's
 * alternatives.
 */
template <typename T, typename V>
concept ConvertibleToAlternative = is_convertible_to_alternative_v<T, V>;

///@{
/**
 * A range whose distance can be obtained without invalidating the range.
 *
 * This is different to std::sized_range because it doesn't require:
 * - that size(rng) is available;
 * - that distance(rng) is a constant time operation.
 *
 * Essentially, it just rules out std::input_ranges without a known size.
 */
template <typename T>
concept SlowSizedRange =
    ranges::cpp20::forward_range<T> || ranges::cpp20::sized_range<T>;

/**
 * Get the name of the "base type" of the given expression.
 *
 * The "base type" is the type with no references or cv qualification.
 */
SQ_ND std::string base_type_name(const auto &thing);

/**
 * Create a std::error code from a platform dependent error code.
 */
SQ_ND inline std::error_code make_error_code(int code) {
  return std::make_error_code(static_cast<std::errc>(code));
}

} // namespace sq

namespace ranges {

// Tell the ranges library that gsl::spans with dynamic extent are views. Note
// that the view concept requires a default constructor so sized gsl::spans are
// *not* views...
template <typename T>
inline constexpr bool enable_view<gsl::span<T, gsl::dynamic_extent>> = true;

// ... unless the size is zero.
template <typename T> inline constexpr bool enable_view<gsl::span<T, 0>> = true;

} // namespace ranges

#include "core/typeutil.inl.h"

#endif // SQ_INCLUDE_GUARD_core_typeutil_h_

ext: h

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


text:
#include <string.h>
#include <stdint.h>
#include <mpi.h>

#include <gsl/gsl_rng.h>
#include <fastpm/libfastpm.h>
#include <fastpm/logging.h>
#include <fastpm/string.h>
#include <fastpm/transfer.h>
#include "chealpix.h"
#include <mpsort.h>

#include "pmpfft.h"
#include "pmghosts.h"

static double RNDTABLE[8192];

gsl_rng * random_generator;
void fastpm_utils_init_randtable() {
    random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    gsl_rng_set(random_generator, 37 * (rank + 1) + 1);  /* start-up seed */
    int i;
    for(i = 0; i < 8192; i ++) {
        RNDTABLE[i] = gsl_rng_uniform(random_generator);
    }
    //gsl_rng_free(random_generator);
}

double
fastpm_utils_get_random(uint64_t id)
{
    return gsl_rng_uniform(random_generator);
    uint64_t ind = 0;
    ind = id;

    while(id != 0) {
        ind += id;
        id /= 8192;
    }
    ind %= 8192;
    return RNDTABLE[ind];
}

void
fastpm_utils_dump(PM * pm , const char * filename, FastPMFloat *data)
{
    char fn1[1024];
    char fn2[1024];
    if(pm->NTask > 1) {
        sprintf(fn1, "%s.%03d.geometry", filename, pm->ThisTask);
        sprintf(fn2, "%s.%03d", filename, pm->ThisTask);
    } else {
        sprintf(fn1, "%s.geometry", filename);
        sprintf(fn2, "%s", filename);
    }

    fastpm_path_ensure_dirname(filename);

    FILE * fp;
    fp = fopen(fn2, "w");
    if(pm->allocsize != fwrite(data, sizeof(FastPMFloat), pm->allocsize, fp)) {
        fastpm_raise(-1, "IO failed.\n");
    }
    fclose(fp);
    fp = fopen(fn1, "w");
    fprintf(fp, "# real\n");
    fprintf(fp, "start: %td %td %td\n",
                    pm->IRegion.start[0],
                    pm->IRegion.start[1],
                    pm->IRegion.start[2]);
    fprintf(fp, "size: %td %td %td\n",
                    pm->IRegion.size[0],
                    pm->IRegion.size[1],
                    pm->IRegion.size[2]);
    fprintf(fp, "strides: %td %td %td\n",
                    pm->IRegion.strides[0],
                    pm->IRegion.strides[1],
                    pm->IRegion.strides[2]);

    fprintf(fp, "# complex\n");
    fprintf(fp, "start: %td %td %td\n",
                    pm->ORegion.start[0],
                    pm->ORegion.start[1],
                    pm->ORegion.start[2]);
    fprintf(fp, "size: %td %td %td\n",
                    pm->ORegion.size[0],
                    pm->ORegion.size[1],
                    pm->ORegion.size[2]);
    fprintf(fp, "strides: %td %td %td\n",
                    pm->ORegion.strides[0],
                    pm->ORegion.strides[1],
                    pm->ORegion.strides[2]);

    fclose(fp);
}

void
fastpm_utils_load(PM * pm , const char * filename, FastPMFloat *data)
{
    char fn1[1024];
    char fn2[1024];
    if(pm->NTask > 1) {
        sprintf(fn1, "%s.%03d.geometry", filename, pm->ThisTask);
        sprintf(fn2, "%s.%03d", filename, pm->ThisTask);
    } else {
        sprintf(fn1, "%s.geometry", filename);
        sprintf(fn2, "%s", filename);
    }
    FILE * fp;
    fp = fopen(fn2, "r");
    if(pm->allocsize != fread(data, sizeof(FastPMFloat), pm->allocsize, fp)) {
        fastpm_raise(-1, "File was bad\n");
    };
    fclose(fp);
}

static double
tk_eh(double k, struct fastpm_powerspec_eh_params * params)		/* from Martin White */
{
    double q, theta, ommh2, a, s, gamma, L0, C0;
    double tmp;
    double omegam, ombh2, hubble;

    /* other input parameters */
    hubble = params->hubble_param;
    omegam = params->omegam;
    ombh2 = params->omegab * hubble * hubble;

    theta = 2.728 / 2.7;
    ommh2 = omegam * hubble * hubble;
    s = 44.5 * log(9.83 / ommh2) / sqrt(1. + 10. * exp(0.75 * log(ombh2))) * hubble;
    a = 1. - 0.328 * log(431. * ommh2) * ombh2 / ommh2
        + 0.380 * log(22.3 * ommh2) * (ombh2 / ommh2) * (ombh2 / ommh2);
    gamma = a + (1. - a) / (1. + exp(4 * log(0.43 * k * s)));
    gamma *= omegam * hubble;
    q = k * theta * theta / gamma;
    L0 = log(2. * exp(1.) + 1.8 * q);
    C0 = 14.2 + 731. / (1. + 62.5 * q);
    tmp = L0 / (L0 + C0 * q * q);
    return (tmp);
}

double
fastpm_utils_powerspec_eh(double k, struct fastpm_powerspec_eh_params * param)	/* Eisenstein & Hu */
{
    return param->Norm * k * pow(tk_eh(k, param), 2);
}

double
fastpm_utils_powerspec_white(double k, double * amplitude) 	/* White Noise*/
{
    return *amplitude;
}

static void
_sort_pix(const void * ptr, void * radix, void * arg)
{
    memcpy(radix, ptr, 8);
}

/* generate evenly balanced healpix positions that are inside the lightcone. */
void
fastpm_utils_healpix_ra_dec (
                int nside,
                double **ra,
                double **dec,
                uint64_t **pix,
                size_t * n,
                FastPMLightCone * lc,
                MPI_Comm comm
            )
{
    const double rad_to_degree = 180./M_PI;

    size_t npix = nside2npix (nside);
    int ThisTask, NTask;

    MPI_Comm_rank(comm, &ThisTask);
    MPI_Comm_size(comm, &NTask);

    //fastpm_info("healpix npix %ld \n",*npix);
    size_t i = 0;

    size_t pix_start = ThisTask * npix / NTask;
    size_t pix_end = (ThisTask + 1) * npix / NTask;

    uint64_t local_npix = 0;

    uint64_t * pixels = NULL;

    /* two iterations; estimate and fill */
    while(1) {
        size_t j = 0;

        for (i = pix_start; i < pix_end; i++)
        {
            double vec[3];

            pix2vec_ring(nside, i, vec);

            if(!fastpm_lc_inside(lc, vec)) continue;

            if(pixels != NULL) {
                pixels[j] = i;
            }

            j ++;
        }
        if(pixels == NULL) {
            local_npix = j;
            pixels = malloc(sizeof(uint64_t) * local_npix);

        } else {
            break;
        }
    }

    /* count total */
    uint64_t valid_npix = local_npix;
    MPI_Allreduce(MPI_IN_PLACE, &valid_npix, 1, MPI_LONG, MPI_SUM, comm);

    /* redistribute / balance */

    size_t localsize =
        (ThisTask + 1) * valid_npix / NTask
    -
        (ThisTask    ) * valid_npix / NTask;

    uint64_t * recv_buffer = malloc(sizeof(uint64_t) * localsize);

    mpsort_mpi_newarray(pixels, local_npix, recv_buffer, localsize, sizeof(uint64_t), _sort_pix, 8, NULL, comm);

    free(pixels);

    *ra = malloc(sizeof(double) * localsize);
    *dec = malloc(sizeof(double) * localsize);

    for(i = 0; i < localsize; i ++) {
        double phi, theta;

        pix2ang_ring(nside, recv_buffer[i], &theta, &phi);

        phi *= rad_to_degree;
        theta*= rad_to_degree;
        (*ra)[i] = phi;
        (*dec)[i]= 90 - theta;
        recv_buffer[i] += (((uint64_t)nside) << 48);
    }
    *n = localsize;
    *pix = recv_buffer;
}

void
fastpm_gldot(double glmatrix[4][4], double xi[4], double xo[4])
{
    int i, j;
    for(i = 0; i < 4; i ++) {
        xo[i] = 0;
        for(j = 0; j < 4; j ++) {
            xo[i] += glmatrix[i][j] * xi[j];
        }
    }
}

void
fastpm_gldotf(double glmatrix[4][4], float vi[4], float vo[4])
{
    int i, j;
    for(i = 0; i < 4; i ++) {
        vo[i] = 0;
        for(j = 0; j < 4; j ++) {
            vo[i] += glmatrix[i][j] * vi[j];
        }
    }
}

/* returns MPI_LOR reduction */
int
MPIU_Any(MPI_Comm comm, int value)
{
    MPI_Allreduce(MPI_IN_PLACE, &value, 1, MPI_INT, MPI_LOR, comm);
    return value;
}

/* returns MPI_LAND reduction */
int
MPIU_All(MPI_Comm comm, int value)
{
    MPI_Allreduce(MPI_IN_PLACE, &value, 1, MPI_INT, MPI_LAND, comm);
    return value;
}

/*
 * Compute the standard statistics of a value over communicator.
 * fmt is a string, describing the list of stats to return:
 *  < : minimum  (double)
 *  , : rank of minimum (int)
 *  > : maximum  (double)
 *  . : rank of maximum (int)
 *  - : mean  (double)
 *  s : std (biased standard derivation)
 *  v : variance (biased variance)
 *  S : unbiased std
 *  V : unbiased variance
 *
 * e.g.
 *
 *    double min, max;
 *    MPIU_stats(comm, r, "<>", &min, &max);
 *
 * returns number of items converted;
 * if a fmt char is unknown, the input variable is unmodified.
 *
 * */
int
MPIU_stats(MPI_Comm comm,
    const double value,
    const char * fmt, ...)
{

    va_list va;
    va_start(va, fmt);

    int NTask;
    MPI_Comm_size(comm, &NTask);
    int ThisTask;
    MPI_Comm_rank(comm, &ThisTask);

    double sum = value;
    struct {
        double value;
        int rank;
    } min = {value, ThisTask},
      max = {value, ThisTask};

    double sumsq = value * value;

    MPI_Allreduce(MPI_IN_PLACE, &sum, 1, MPI_DOUBLE, MPI_SUM, comm);
    MPI_Allreduce(MPI_IN_PLACE, &sumsq, 1, MPI_DOUBLE, MPI_SUM, comm);
    MPI_Allreduce(MPI_IN_PLACE, &min, 1, MPI_DOUBLE_INT, MPI_MINLOC, comm);
    MPI_Allreduce(MPI_IN_PLACE, &max, 1, MPI_DOUBLE_INT, MPI_MAXLOC, comm);

    double avg = (sum) / NTask;
    double avg_sq = (sumsq) / NTask;
    double var = avg_sq - avg * avg;
    double std = sqrt(avg_sq - avg * avg);

    int i;
    int c = 0;
    for(i = 0; i < strlen(fmt); i ++) {
        void * r = va_arg(va, void *);
        int * ir = (int * ) r;
        double * dr = (double * ) r;

        c++;
        switch(fmt[i]) {
            case '-':
                *dr = avg;
                break;
            case '<':
                *dr = min.value;
                break;
            case '>':
                *dr = max.value;
                break;
            case ',':
                *ir = min.rank;
                break;
            case '.':
                *ir = max.rank;
                break;
            case 's':
                *dr = std;
                break;
            case 'S':
                *dr = std * sqrt(1.0 * NTask / (NTask - 1.0));
                break;
            case 'v':
                *dr = var;
                break;
            case 'V':
                *dr = var * 1.0 * NTask / (NTask - 1.0);
                break;
            default:
                c--;
        }
    }
    va_end(va);

    return c;
}

/*
 *
 * create a copy of string src and broadcast it to all ranks.
 *
 * Every rank receives a copy of src as the return value, including the root.
 *
 * if free is not NULL, the function is called as free(src)
 * */
char *
MPIU_Bcast_string(MPI_Comm comm, char * src, int root, void(*free)(void * ptr))
{
    int srcNULL = src == NULL;
    MPI_Bcast(&srcNULL, 1, MPI_INT, root, comm);

    if(srcNULL) return NULL;

    int rank;
    int N;
    MPI_Comm_rank(comm, &rank);

    if(rank == root) {
        N = strlen(src) + 1;
    }
    MPI_Bcast(&N, 1, MPI_INT, root, comm);

    char * ret = malloc(N);
    if(rank == root) {
        strcpy(ret, src);
    }

    if(free && rank == root) free(src);
    MPI_Bcast(ret, N, MPI_BYTE, root, comm);

    return ret;
}

ext: c

max_stars_repo_name: sbird/FastPMRunner

max_issues_repo_name: sbird/FastPMRunner

max_forks_repo_name: sbird/FastPMRunner


text:
/* filter/test_impulse.c
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * 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 3 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.
 */

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_filter.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

static int
vector_sum(const gsl_vector_int * v)
{
  size_t i;
  int sum = 0;

  for (i = 0; i < v->size; ++i)
    {
      int vi = gsl_vector_int_get(v, i);
      sum += vi;
    }

  return sum;
}

static void
test_impulse_proc(const double tol, const size_t n, const size_t K, const double nsigma,
                  const gsl_filter_end_t etype, const gsl_filter_scale_t stype,
                  const double outlier_percentage, gsl_rng * r)
{
  gsl_vector * x = gsl_vector_alloc(n);
  gsl_vector * y = gsl_vector_alloc(n);
  gsl_vector * z = gsl_vector_alloc(n);
  gsl_vector * y_med = gsl_vector_alloc(n);
  gsl_vector * xmedian = gsl_vector_alloc(n);
  gsl_vector * xsigma = gsl_vector_alloc(n);
  size_t noutlier;
  gsl_vector_int * ioutlier = gsl_vector_int_alloc(n);
  gsl_vector_int * ioutlier_exact = gsl_vector_int_alloc(n);
  size_t i;
  gsl_filter_impulse_workspace *impulse_p = gsl_filter_impulse_alloc(K);
  gsl_filter_median_workspace *median_p = gsl_filter_median_alloc(K);
  size_t noutlier_exact = 0;
  char buf[1024];

  gsl_vector_int_set_zero(ioutlier_exact);

  for (i = 0; i < n; ++i)
    {
      double xi = gsl_ran_gaussian(r, 1.0);
      double vi = gsl_rng_uniform(r);

      if (vi <= outlier_percentage)
        {
          xi += 15.0 * GSL_SIGN(xi);
          ++noutlier_exact;
          gsl_vector_int_set(ioutlier_exact, i, 1);
        }

      gsl_vector_set(x, i, xi);
    }

  /* first test that median filter is equal to impulse filter with nsigma = 0 */

  gsl_filter_median(etype, x, y_med, median_p);
  gsl_filter_impulse(etype, stype, 0.0, x, y, xmedian, xsigma, &noutlier, ioutlier, impulse_p);

  sprintf(buf, "impulse nsigma=0 smf comparison, etype=%u stype=%u", etype, stype);
  compare_vectors(tol, y, y_med, buf);

  /* second test: filter y = impulse(x) with given nsigma */

  gsl_filter_impulse(etype, stype, nsigma, x, y, xmedian, xsigma, &noutlier, ioutlier, impulse_p);

  /* test correct number of outliers detected */
  gsl_test(noutlier != noutlier_exact, "impulse [n=%zu,K=%zu,nsigma=%g,outlier_percentage=%g] noutlier=%zu exact=%zu",
           n, K, nsigma, outlier_percentage, noutlier, noutlier_exact);

#if 0
  {
    for (i = 0; i < n; ++i)
      {
        printf("%.12e %.12e %d %.12e %.12e\n",
               gsl_vector_get(x, i),
               gsl_vector_get(y, i),
               gsl_vector_int_get(ioutlier, i),
               gsl_vector_get(xmedian, i) + nsigma * gsl_vector_get(xsigma, i),
               gsl_vector_get(xmedian, i) - nsigma * gsl_vector_get(xsigma, i));
      }
  }
#endif

  /* test outliers found in correct locations */
  for (i = 0; i < n; ++i)
    {
      int val = gsl_vector_int_get(ioutlier, i);
      int val_exact = gsl_vector_int_get(ioutlier_exact, i);

      gsl_test(val != val_exact, "test_impulse: outlier mismatch [i=%zu,K=%zu,nsigma=%g,outlier_percentage=%g] ioutlier=%d ioutlier_exact=%d",
               i, K, nsigma, outlier_percentage, val, val_exact);
    }

  /* test noutlier = sum(ioutlier) */
  {
    size_t iout_sum = vector_sum(ioutlier);
    gsl_test(noutlier != iout_sum, "impulse [K=%zu,nsigma=%g,outlier_percentage=%g] noutlier=%zu sum(ioutlier)=%zu",
             K, nsigma, outlier_percentage, noutlier, iout_sum);
  }

  /* third test: test in-place filter z = impulse(z) */

  gsl_vector_memcpy(z, x);
  gsl_filter_impulse(etype, stype, nsigma, z, z, xmedian, xsigma, &noutlier, ioutlier, impulse_p);

  sprintf(buf, "impulse in-place nsigma=%g,n=%zu,K=%zu,etype=%u stype=%u", nsigma, n, K, etype, stype);
  compare_vectors(GSL_DBL_EPSILON, z, y, buf);

  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(z);
  gsl_vector_free(y_med);
  gsl_vector_free(xmedian);
  gsl_vector_free(xsigma);
  gsl_vector_int_free(ioutlier);
  gsl_vector_int_free(ioutlier_exact);
  gsl_filter_impulse_free(impulse_p);
  gsl_filter_median_free(median_p);
}

static void
test_impulse(gsl_rng * r)
{
  const double tol = 1.0e-10;

  test_impulse_proc(tol, 1000, 21, 6.0, GSL_FILTER_END_TRUNCATE, GSL_FILTER_SCALE_QN, 0.05, r);
  test_impulse_proc(tol, 1000, 21, 6.0, GSL_FILTER_END_TRUNCATE, GSL_FILTER_SCALE_SN, 0.05, r);
  test_impulse_proc(tol, 1000, 21, 6.0, GSL_FILTER_END_TRUNCATE, GSL_FILTER_SCALE_MAD, 0.05, r);
  test_impulse_proc(tol, 1000, 21, 6.0, GSL_FILTER_END_TRUNCATE, GSL_FILTER_SCALE_IQR, 0.05, r);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/**
 */
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp.h>
#include "aXe_utils.h"
#include "aXe_errors.h"
#include "fringe_conf.h"
#include "fringe_model.h"

#define MAX(x,y) (((x)>(y))?(x):(y))
#define MIN(x,y) (((x)<(y))?(x):(y))

/**
 * Function: compute_fringe_amplitude
 * The function computes the fringe image for a CCD setup stored
 * in a fringe configuration structure. This fringe configuration
 * structure completely describes the problem, and this function
 * executes the loops for every pixel over the wavelength
 * range spanned by the filter.
 *
 * Parameters
 * @param fconf - the fringe configuration structure
 *
 * Returns:
 * @return fringe_image - the image with the computed fringe amplitudes
 */
gsl_matrix  *
compute_fringe_amplitude(fringe_conf *fconf)
{
  gsl_matrix *fringe_image;

  gsl_vector **filter_vectors;

  int index=0;
  int ii=0;
  int jj=0;

  double lambda_mean;
  double pixel_ampl;
  //double phase_number;

  optical_property *optprops;

  // allocate the fringe image
  fringe_image = alloc_fringe_image(fconf->opt_layers);

  // allocate memory for the optical property structure
  optprops = alloc_optprops_list(fconf);

  // find the exact wavelength range,
  // and define the wavelength values and
  // normalized fiter throughput there
  filter_vectors = evaluate_wavelength_steps(fconf);
  //filter_vectors = get_PET_calibration_data();

  // compute the mean wavelength
  lambda_mean = gsl_vector_get(filter_vectors[0],filter_vectors[0]->size-1)/2.0
		 + gsl_vector_get(filter_vectors[0],0)/2.0;

  // initialize some values in the optical property list
  init_optprops_list(fconf, lambda_mean, optprops);

    for (ii=0; ii < (int)fringe_image->size1; ii++)
      //    for (ii=0; ii < 2; ii++)
    {
      fprintf(stderr, "Computing row No.: %i\n", ii);
            for (jj=0; jj < (int)fringe_image->size2; jj++)
	      //for (jj=0; jj < 2; jj++)
	{
	  // fill the optical thickness of the layers
	  // into the structure
	  fill_optprops_thickness(fconf->opt_layers, ii, jj, optprops);

	  pixel_ampl = 0.0;
	  for (index=0; index < (int)filter_vectors[0]->size; index++)
	    {
	      // fill all information in the optical
	      // property list
	      fill_optprops_all(fconf->opt_layers,
				gsl_vector_get(filter_vectors[0],index),
				optprops);

	      // compute and add the contribution at a wavelength
	      pixel_ampl += gsl_vector_get(filter_vectors[1],index)*
		fringe_contrib_single(optprops, fconf);
	    }

	  // finally set the pixel value
	  // in the output image
	  gsl_matrix_set(fringe_image, ii, jj,
			 fconf->fringe_amp * pixel_ampl + 1.0);
	}
    }

  // release the memory in the vectors
  gsl_vector_free(filter_vectors[0]);
  gsl_vector_free(filter_vectors[1]);
  free(filter_vectors);

  // free the optical property structure
  free_optprops_list(optprops);

  // return the fringe image
  return fringe_image;
}

/**
 * Function: alloc_fringe_image
 * The function browses through the the structure for the CCD layers
 * and extracts all information on image sizes stored there in
 * one or several thickness images. This information is checked
 * for consistency.
 * Finally, a matrix for a fringe image is allocated, initialized
 * and returned.
 *
 * Parameters:
 * @param opt_layers - the optical layers in the CCD
 *
 * Returns:
 * @return fringe_image - the allocated matrix for the fringe image
 */
gsl_matrix  *
alloc_fringe_image(const ccd_layers *opt_layers)
{
  int n1         = 0;
  int n2         = 0;
  int n1_new     = 0;
  int n2_new     = 0;
  int index      = 0;
  int is_defined = 0;

  gsl_matrix *fringe_image=NULL;
  // go through all CCD layers
  for (index=0; index < opt_layers->num_layers; index++)
    {
      // check whether the thickness information
      // is represented by an image
      if (opt_layers->opt_layer[index]->thickness2D != NULL)
	{

	  // mark that an image was found
	  is_defined=1;

	  // get the dimension of the image
	  n1_new = opt_layers->opt_layer[index]->thickness2D->size1;
	  n2_new = opt_layers->opt_layer[index]->thickness2D->size2;

	  // check the frst axix value against
	  // previous values, if possible
	  if (n1 && n1_new != n1)
	    // give an error if the new value is different
	    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			 "Thickness image of layer %i "
			 "has %i pixels in first axis."
			 "This differs from the previous value %i\n",
			 index, n1_new, n1);
	  else
	    // store the new value
	    n1 = n1_new;

	  // check the frst axix value against
	  // previous values, if possible
	  if (n2 && n2_new != n2)
	    // give an error if the new value is different
	    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
			 "Thickness image of layer %i "
			 "has %i pixels in second axis."
			 "This differs from the previous value %i\n",
			 index, n2_new, n2);
	  else
	    // store the new value
	    n2 = n2_new;
	}
    }

#ifdef DEBUGFCONF
  fprintf(stdout, "Allocating an image with size: (%i, %i)\n", n1, n2);
#endif

  // check whether at least one image was
  // found as thickness information
  if (is_defined)
    {
      // allocate the image
      fringe_image =  gsl_matrix_alloc(n1, n2);

      // set all vallues to zero
      gsl_matrix_set_all(fringe_image, 0.0);
    }
  else
    {
      // no image found; report the error
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "The size of the fringe image is unknown\n"
		   "since none of the layers has an image "
		   "to specify its thickness!\n");
    }

  // return the image
  return fringe_image;
}

/**
 * Function: get_calibration_data()
 * This function is a data storage for filter throughput values
 * as computed by the orioginal program 'final.f' of J.Walsh.
 * The cubic spline routine used there turned out to give
 * (slightly) different results than the gsl-routine used
 * here. To be able making a in-depth comparison the throughuts
 * were computed for the filter 'ccd_ref.fits' and stored here.
 *
 * Parameters:
 *
 * Returns:
 * @return double_vector - interpolated standard filter throughputs
 */
gsl_vector **
get_calibration_data()
{
  gsl_vector *lambda_values;
  gsl_vector *through_values;

  gsl_vector **double_vector;

  // allocate space for the return vector
  double_vector  = (gsl_vector **)malloc(2*sizeof (gsl_vector *));

  // allocate the space for the vectors
  lambda_values  = gsl_vector_alloc(43);
  through_values = gsl_vector_alloc(43);

  // what follows is the wavelength steps as well as the normalized
  // filter throughputs which were computed in the original
  // program 'final.f' of J.Walsh.
  // The gsl-cubic spline routine yields slightly different values.
  // The values fixed here are good for a detailed cross-check
  // of the two programs.
  gsl_vector_set(through_values, 0, 0.0000000000000000);gsl_vector_set(lambda_values, 0,   0.9179);
  gsl_vector_set(through_values, 1, 0.0095066184313808);gsl_vector_set(lambda_values, 1,   0.9180);
  gsl_vector_set(through_values, 2, 0.0176006475945001);gsl_vector_set(lambda_values, 2,   0.9181);
  gsl_vector_set(through_values, 3, 0.0240259215044616);gsl_vector_set(lambda_values, 3,   0.9182);
  gsl_vector_set(through_values, 4, 0.0289185490230508);gsl_vector_set(lambda_values, 4,   0.9183);
  gsl_vector_set(through_values, 5, 0.0324146730103703);gsl_vector_set(lambda_values, 5,   0.9184);
  gsl_vector_set(through_values, 6, 0.0346504181940869);gsl_vector_set(lambda_values, 6,   0.9185);
  gsl_vector_set(through_values, 7, 0.0357619115684218);gsl_vector_set(lambda_values, 7,   0.9186);
  gsl_vector_set(through_values, 8, 0.0358852710613785);gsl_vector_set(lambda_values, 8,   0.9187);
  gsl_vector_set(through_values, 9, 0.0351566282002872);gsl_vector_set(lambda_values, 9,   0.9188);
  gsl_vector_set(through_values,10, 0.0337121099793692);gsl_vector_set(lambda_values,10,   0.9189);
  gsl_vector_set(through_values,11, 0.0316878433928458);gsl_vector_set(lambda_values,11,   0.9190);
  gsl_vector_set(through_values,12, 0.0292199509018293);gsl_vector_set(lambda_values,12,   0.9191);
  gsl_vector_set(through_values,13, 0.0264445572339866);gsl_vector_set(lambda_values,13,   0.9192);
  gsl_vector_set(through_values,14, 0.0234977916500934);gsl_vector_set(lambda_values,14,   0.9193);
  gsl_vector_set(through_values,15, 0.0205157811443709);gsl_vector_set(lambda_values,15,   0.9194);
  gsl_vector_set(through_values,16, 0.0176346504444861);gsl_vector_set(lambda_values,16,   0.9195);
  gsl_vector_set(through_values,17, 0.0149905288112146);gsl_vector_set(lambda_values,17,   0.9196);
  gsl_vector_set(through_values,18, 0.0127195375723916);gsl_vector_set(lambda_values,18,   0.9197);
  gsl_vector_set(through_values,19, 0.0109578003224067);gsl_vector_set(lambda_values,19,   0.9198);
  gsl_vector_set(through_values,20, 0.0098414542549762);gsl_vector_set(lambda_values,20,   0.9199);
  gsl_vector_set(through_values,21, 0.0095066184313808);gsl_vector_set(lambda_values,21,   0.9200);
  gsl_vector_set(through_values,22, 0.0100435357173845);gsl_vector_set(lambda_values,22,   0.9201);
  gsl_vector_set(through_values,23, 0.0113589226644173);gsl_vector_set(lambda_values,23,   0.9202);
  gsl_vector_set(through_values,24, 0.0133136264280558);gsl_vector_set(lambda_values,24,   0.9203);
  gsl_vector_set(through_values,25, 0.0157684646986687);gsl_vector_set(lambda_values,25,   0.9204);
  gsl_vector_set(through_values,26, 0.0185842800987238);gsl_vector_set(lambda_values,26,   0.9205);
  gsl_vector_set(through_values,27, 0.0216219095843027);gsl_vector_set(lambda_values,27,   0.9206);
  gsl_vector_set(through_values,28, 0.0247421844451007);gsl_vector_set(lambda_values,28,   0.9207);
  gsl_vector_set(through_values,29, 0.0278059450370313);gsl_vector_set(lambda_values,29,   0.9208);
  gsl_vector_set(through_values,30, 0.0306740203832353);gsl_vector_set(lambda_values,30,   0.9209);
  gsl_vector_set(through_values,31, 0.0332072463065172);gsl_vector_set(lambda_values,31,   0.9210);
  gsl_vector_set(through_values,32, 0.0352664608962358);gsl_vector_set(lambda_values,32,   0.9211);
  gsl_vector_set(through_values,33, 0.0367124999751954);gsl_vector_set(lambda_values,33,   0.9212);
  gsl_vector_set(through_values,34, 0.0374061970996460);gsl_vector_set(lambda_values,34,   0.9213);
  gsl_vector_set(through_values,35, 0.0372083835592830);gsl_vector_set(lambda_values,35,   0.9214);
  gsl_vector_set(through_values,36, 0.0359798997100196);gsl_vector_set(lambda_values,36,   0.9215);
  gsl_vector_set(through_values,37, 0.0335815745749970);gsl_vector_set(lambda_values,37,   0.9216);
  gsl_vector_set(through_values,38, 0.0298742553097917);gsl_vector_set(lambda_values,38,   0.9217);
  gsl_vector_set(through_values,39, 0.0247187709375447);gsl_vector_set(lambda_values,39,   0.9218);
  gsl_vector_set(through_values,40, 0.0179759414151792);gsl_vector_set(lambda_values,40,   0.9219);
  gsl_vector_set(through_values,41, 0.0095066184313808);gsl_vector_set(lambda_values,41,   0.9220);
  gsl_vector_set(through_values,42, 0.0000000000000000);gsl_vector_set(lambda_values,42,   0.9221);

  // build up the output array
  double_vector[0] = lambda_values;
  double_vector[1] = through_values;

  // return the two vectors
  return double_vector;
}

/**
 * Function: get_PET_calibration_data()
 */
gsl_vector **
get_PET_calibration_data()
{
  gsl_vector *lambda_values;
  gsl_vector *through_values;

  gsl_vector **double_vector;

  // allocate space for the return vector
  double_vector  = (gsl_vector **)malloc(2*sizeof (gsl_vector *));

  // allocate the space for the vectors
  lambda_values  = gsl_vector_alloc(21);
  through_values = gsl_vector_alloc(21);

  gsl_vector_set(lambda_values, 0,1.069460);  gsl_vector_set(through_values, 0, 0.000000);
  gsl_vector_set(lambda_values, 1,1.069653);  gsl_vector_set(through_values, 1, 0.012630);
  gsl_vector_set(lambda_values, 2,1.069846);  gsl_vector_set(through_values, 2, 0.025261);
  gsl_vector_set(lambda_values, 3,1.070039);  gsl_vector_set(through_values, 3, 0.037891);
  gsl_vector_set(lambda_values, 4,1.070232);  gsl_vector_set(through_values, 4,   0.050522);
  gsl_vector_set(lambda_values, 5,1.070424);  gsl_vector_set(through_values, 5,   0.063152);
  gsl_vector_set(lambda_values, 6,1.070617);  gsl_vector_set(through_values, 6,   0.069010);
  gsl_vector_set(lambda_values, 7,1.070810);  gsl_vector_set(through_values, 7,   0.069010);
  gsl_vector_set(lambda_values, 8,1.071003);  gsl_vector_set(through_values, 8,   0.069010);
  gsl_vector_set(lambda_values, 9,1.071196);  gsl_vector_set(through_values, 9,   0.069010);
  gsl_vector_set(lambda_values, 10,1.071389); gsl_vector_set(through_values, 10,   0.069010);
  gsl_vector_set(lambda_values, 11,1.071582); gsl_vector_set(through_values, 11,   0.069010);
  gsl_vector_set(lambda_values, 12,1.071775); gsl_vector_set(through_values, 12,   0.069010);
  gsl_vector_set(lambda_values, 13,1.071967); gsl_vector_set(through_values, 13,   0.069010);
  gsl_vector_set(lambda_values, 14,1.072160); gsl_vector_set(through_values, 14,   0.069010);
  gsl_vector_set(lambda_values, 15,1.072353); gsl_vector_set(through_values, 15,   0.063152);
  gsl_vector_set(lambda_values, 16,1.072546); gsl_vector_set(through_values, 16,   0.050522);
  gsl_vector_set(lambda_values, 17,1.072739); gsl_vector_set(through_values, 17,   0.037891);
  gsl_vector_set(lambda_values, 18,1.072932); gsl_vector_set(through_values, 18,   0.025261);
  gsl_vector_set(lambda_values, 19,1.073125); gsl_vector_set(through_values, 19,  0.012630);
  gsl_vector_set(lambda_values, 20,1.073318); gsl_vector_set(through_values, 20,   0.000000);

  //Wavelength: 10713.888672, fringe factor: 1.097945, (x,y): (18, 745)

  // build up the output array
  double_vector[0] = lambda_values;
  double_vector[1] = through_values;

  // return the two vectors
  return double_vector;
}

/**
 * Function: evaluate_wavelength_steps
 * The function defines and computes the wavelength steps and the
 * filter throughputs which are used to compute the fringe amplitude
 * for all pixels. The basis for the wavelength data are the
 * filter data and the general are where fringing is significant
 * and computed, which is defined in the fringe
 * configuration structure.
 * First all throughput values below a certain value are discarded
 * at the short and long wavelength edge. Then the final wavelength
 * range is determined using the fringing range. Then all wavelength
 * values and the filter throughputs are determined and returned
 * as two gsl-vectors.
 *
 * Parameters:
 * @param fconf - the fringe configuration structure
 *
 * Returns:
 * @return double_vector - a set of two gsl-vectors with wavlength
                           and throughput
 */
gsl_vector **
evaluate_wavelength_steps(fringe_conf *fconf)
{
  int index=0;
  int lower=0;
  int upper=0;

  int nsteps;

  double lambda_min;
  double lambda_max;
  double lambda_act;
  double through_act=0.0;
  double through_tot=0.0;

  gsl_vector *lambda_values;
  gsl_vector *through_values;

  gsl_vector **double_vector;

  // allocate space for the return vector
  double_vector  = (gsl_vector **)malloc(2*sizeof (gsl_vector *));

  // search through the filter data
  // to narrow the usable range from below
  lower=0;
  while (lower < fconf->filter_through->nvals
	 && fconf->filter_through->yvals[lower] < FILTER_THRESHOLD)
    lower++;
  lower = MAX(0,lower-1);

  // search through the filter data
  // to narrow the usable range from above
  upper=fconf->filter_through->nvals - 1;
  while (upper > -1
	 && fconf->filter_through->yvals[upper] < FILTER_THRESHOLD)
    upper--;
  upper = MIN(fconf->filter_through->nvals - 1, upper+1);

#ifdef DEBUGFCONF
  fprintf(stderr, "New lower range index: %i, new upper: %i, old number: %i\n",
	  lower, upper, fconf->filter_through->nvals);
#endif

  // replace the old filter interpolator
  // with a new one, if necessary
  if (lower != 0 || upper != (fconf->filter_through->nvals - 1))
    fconf->filter_through = redefine_filter_throughput(lower, upper, fconf);

  // compute the lower wavelength range by comparing
  // the filter throughput with the range important
  // for fringing; the latter is given in AA!!
  lambda_min = MAX(gsl_vector_get(fconf->fringe_range, 0),
		   fconf->filter_through->xmin);

  // compute the upper wavelength range by comparing
  // the filter throughput with the range important
  // for fringing; the latter is given in AA!!
  lambda_max = MIN(gsl_vector_get(fconf->fringe_range, 1),
		   fconf->filter_through->xmax);

  // find the number of steps to cover the wavelength range
  nsteps = (int)ceil((lambda_max-lambda_min)/(fconf->fringe_step)) + 1;

  // allocate the space for the vectors
  lambda_values  = gsl_vector_alloc(nsteps);
  through_values = gsl_vector_alloc(nsteps);

  // fill the vectors with the wavelength steps and the
  // filter throughput at those wavelengths.
  lambda_act = lambda_min;
  for (index=0; index < nsteps-1; index+=1)
    {
      gsl_vector_set(lambda_values, index, lambda_act);
      through_act = eval_interp(fconf->filter_through, lambda_act);
      through_tot += through_act;
      gsl_vector_set(through_values, index, through_act);
      lambda_act += fconf->fringe_step;
    }

  // add also the values at lambda_amax
  gsl_vector_set(lambda_values, nsteps-1, lambda_max);
  through_act = eval_interp(fconf->filter_through, lambda_max);
  through_tot += through_act;
  gsl_vector_set(through_values, nsteps-1, through_act);

  for (index=0; index < (int)through_values->size; index++)
    {
      // normalize the filter throughput values
      gsl_vector_set(through_values, index,
		     gsl_vector_get(through_values, index)/through_tot);

      // convert the wavelength from AA to micron
      gsl_vector_set(lambda_values, index,
		     gsl_vector_get(lambda_values, index)*1.0e-04);
    }
#ifdef DEBUGFCONF
  through_tot = 0.0;
  for (index=0; index < through_values->size; index++)
    {
      fprintf(stderr, "Wavelength: %f, Throughput: %f\n",
	      gsl_vector_get(lambda_values, index),
	      gsl_vector_get(through_values, index));
      through_tot += gsl_vector_get(through_values, index);
    }
  fprintf(stderr, "Total throughput: %f\n", through_tot);
#endif

  // build up the output array
  double_vector[0] = lambda_values;
  double_vector[1] = through_values;

  return double_vector;
}

/**
 * Function: redefine_filter_throughput
 * The function re-defines an existing interpolator using only the
 * existing data values in an index range given as parameters.
 *
 * Parameters:
 * @param lower - index with the lowest data point
 * @param upper - index with the highest data point
 *
 * Returns:
 * @return new_interp - the new interpolator
 */
interpolator *
redefine_filter_throughput(const int lower, const int upper,
			   fringe_conf *fconf)
{
  int new_nvals;
  int index;

  double *new_x;
  double *new_y;

  interpolator *new_interp;

   // allocate space for the return structure;
  // complain if this fails
  new_interp = (interpolator *)malloc (sizeof (interpolator));
  if (new_interp == NULL)
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "Could not allocate memory for interpolator");

  // compute the new length of the data arrays
  new_nvals = upper-lower+1;

  // allocate space for the new data arrays
  new_x = (double *)malloc(new_nvals*sizeof(double));
  new_y = (double *)malloc(new_nvals*sizeof(double));

  // transfer the data
  for (index=lower; index < upper+1; index++)
    {
      new_x[index-lower] = fconf->filter_through->xvals[index];
      new_y[index-lower] = fconf->filter_through->yvals[index];
    }

  // free the old structure
  free_interp(fconf->filter_through);

  // build up the interpolator structure
  new_interp->xmin    = new_x[0];
  new_interp->xmax    = new_x[new_nvals-1];
  new_interp->nvals   = new_nvals;
  new_interp->xvals   = new_x;
  new_interp->yvals   = new_y;
  new_interp->acc     = gsl_interp_accel_alloc();
  new_interp->interp  = gsl_interp_alloc(FILTER_INTERP_TYPE,
					 new_interp->nvals);

  // initialize the cubic spline
  gsl_interp_init(new_interp->interp, new_interp->xvals, new_interp->yvals,
		  new_interp->nvals);

  // return the new structure
  return new_interp;
}

/**
 * Function: eval_linear_interp
 * The function computes and returns the interpolated value
 * at a given position for a linear interpolator.
 *
 * Parameters:
 * @param lin_int - the linear interpolator
 * @param xval    - the position to evaluate the linear interpolator
 *
 * Returns:
 * @return (value) - the interpolated data value
 */
double
eval_linear_interp(linear_interp *lin_int, const double xval)
{

  double factor=0.0;

  // check whether the x-value is within
  // the range spanned by the data;
  // complain if the x-value is outside
  if (xval < lin_int->xmin || xval > lin_int->xmax)
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		 "independent linear interpolation value %f "
		 "is outside interval (%f, %f)\n", xval,
		 lin_int->xmin, lin_int->xmax);

  // check whether you have to search upwards or downwards
  if (xval >= gsl_vector_get(lin_int->xvals, lin_int->act_index))
    {

      // in case that you search upwards, go up
      // the independent values until you find the right interval
      lin_int->act_index += 1;
      while(xval > gsl_vector_get(lin_int->xvals, lin_int->act_index))
	lin_int->act_index++;
    }
  else
    {

      // in case that you search downwards, go down
      // the independent values  until you find the right interval
      //      while(wavelength < resp->spec[nact-1].lambda_mean)
      while(xval < gsl_vector_get(lin_int->xvals, lin_int->act_index-1))
	lin_int->act_index--;
    }

  // interpolate within the interval to calculate the
  // sensitivity
  factor = (xval - gsl_vector_get(lin_int->xvals, lin_int->act_index-1))/
    (gsl_vector_get(lin_int->xvals, lin_int->act_index) - gsl_vector_get(lin_int->xvals, lin_int->act_index-1));

  // compute and terutn the interpolated value
    return (gsl_vector_get(lin_int->yvals, lin_int->act_index-1)
    + factor * (gsl_vector_get(lin_int->yvals, lin_int->act_index)-gsl_vector_get(lin_int->yvals, lin_int->act_index-1)));
}

/**
 * Function: get_layer_thickness
 * The function determines the thickness of an individual layer
 * at a given pixel position.
 *
 * Parameters:
 * @param opt_layer - the optical layer
 * @param ii        - the first pixel index
 * @param jj        - the second pixel index
 *
 * Returns:
 * @return (value) - the thickness of the layer at the requested position
 */
double
get_layer_thickness(const ccd_layer *opt_layer, const int ii, const int jj)
{
  // check whether the thickness is 2D
  if (opt_layer->thickness2D)
    // return the 2D value
    return gsl_matrix_get(opt_layer->thickness2D, ii, jj);
  else
    // return the constant value
    return opt_layer->thickness;
}

/**
 * Function: get_complex_refindex
 * The function computes and returns the complex refraction index of
 * a given CCD layer at a given wavelength.
 *
 * Parameters:
 * @param opt_layer - the optical layer
 * @param lambda    - the wavelength
 *
 * Returns:
 * @return (value) - the complex refraction index
 */
gsl_complex
get_complex_refindex(const ccd_layer *opt_layer, const double lambda)
{
  double re;
  double im;

  // compute the real part
  re = eval_interp(opt_layer->re_refraction, lambda);

  // comute the imaginary part
  im = eval_interp(opt_layer->im_refraction, lambda);

  // return the complex number
  return gsl_complex_rect(re, im);
}

/**
 * Function: compute_reflection
 * The function computes and returns the reflection index (in %)
 * for two complex refraction indices  at a given wavelength.
 *
 * Parameters:
 * @param refract_l1 - refraction index of one layer
 * @param refract_l2 - refraction index of the second layer
 *
 * Returns:
 * @return (value) - the reflection index
 */
double
compute_reflection(const gsl_complex refract_l1,
		   const gsl_complex refract_l2)
{
  gsl_complex conj_l1;
  gsl_complex conj_l2;

  gsl_complex refl_compl;

  // get the comlex conjugated values of the input
  conj_l1 = gsl_complex_conjugate(refract_l1);
  conj_l2 = gsl_complex_conjugate(refract_l2);

  // compute the complex reflection
  refl_compl = gsl_complex_div(gsl_complex_sub(conj_l1, conj_l2),
			       gsl_complex_add(conj_l1, conj_l2));

  // return the magnitude of the complex reflection
  return gsl_complex_abs2(refl_compl);
}

/**
 * Function: compute_transmission
 * The function computes and returns the transmission index (in %)
 * for two complex refraction indices at a given wavelength.
 *
 * Parameters:
 * @param refract_l1 - refraction index of one layer
 * @param refract_l2 - refraction index of the second layer
 *
 * Returns:
 * @return (value) - the trnasmission index
 */
double
compute_transmission(const gsl_complex refract_l1,
		     const gsl_complex refract_l2)
{
  gsl_complex conj_l1;
  gsl_complex conj_l2;

  gsl_complex trans_compl;

  // get the comlex conjugated values of the input
  conj_l1 = gsl_complex_conjugate(refract_l1);
  conj_l2 = gsl_complex_conjugate(refract_l2);

  // compute the complex transmission
  trans_compl =
    gsl_complex_div(gsl_complex_sqrt(gsl_complex_mul(conj_l1, conj_l2)),
		    gsl_complex_add(conj_l1, conj_l2));

  // finish the computation of the complex transmission
  trans_compl = gsl_complex_mul_real (trans_compl, 2.0);

  // return the magnitude of the complex transmission
  return gsl_complex_abs2(trans_compl);
}

/**
 * Function: compute_attenuation
 * The function computes the attenuation/damping of a plane wave
 * with a given inverse wavelelength (phase number) traversing
 * TWO TIMES a layer of a given thickness and complex refraction
 * index.
 *
 * Parameters:
 * @param refract      - the complex refraction index
 * @param thickness    - the thickness of the layer
 * @param phase_number - the inverse wavelength
 *
 * Returns:
 * @return (value) - the attenuation
 */
double
compute_attenuation(const gsl_complex refract, const double thickness,
		    const double phase_number)
{
  // check if there is something to compute
  if (GSL_IMAG(refract))
    // compute and return the value
    return exp(-2.0*GSL_IMAG(refract)*thickness*phase_number);
  else
    // return the value for exp(0.0)
    return 1.0;
}

/**
 * Function: compute_pshift
 * The function computes the phase shift of a plane wave
 * with a given inverse wavelelength (phase number) traversing
 * TWO TIMES a layer of a given thickness and complex refraction
 * index.
 *
 * Parameters:
 * @param refract      - the complex refraction index
 * @param thickness    - the thickness of the layer
 * @param phase_number - the inverse wavelength
 *
 * Returns:
 * @return (value) - the complex phase shift
 */
gsl_complex
compute_pshift(const gsl_complex refract, const double thickness,
	       const double phase_number)
{
  // all in one line:
  return gsl_complex_exp(gsl_complex_rect(0.0,2.0*GSL_REAL(refract)*thickness*phase_number));
}

/**
 * Function: alloc_optprops_list
 * The function allocates space for a list of optical property
 * elements. the size of the list is derived from a fringe
 * configuration structure given in the input.
 *
 * Parameters:
 * @param fconf - the fringe configuration file
 *
 * Returns:
 * @return optprops - the list of optical properties
 */
optical_property *
alloc_optprops_list(const fringe_conf *fconf)
{
  optical_property *optprops;

  // allocate large enough memory
  optprops =
    (optical_property *)malloc(fconf->opt_layers->num_layers*sizeof(ap_pixel));

  // return the pointer
  return optprops;
}

/**
 * Function: print_optprops_list
 * The function prints the content of an optical property
 * list onto the screen.
 *
 * Parameters:
 * @param optprops - the optical property list
 *
 * Returns:
 * @return -
 */
void
print_optprops_list(const optical_property *optprops, const int num_entries)
{
  int index;

  for (index = 0; index < num_entries; index++)
    {
      fprintf(stdout, "#%i: %f %f %f %f %f %f \n",index,
	      optprops[index].trans_upper, optprops[index].reflect_upper,
	      GSL_REAL(optprops[index].double_phshift),
	      GSL_IMAG(optprops[index].double_phshift),
	      optprops[index].trans_lower, optprops[index].reflect_lower);
    }
}

/**
 * Function: free_optprops_list
 * The function releases the space allocated in an
 * optical property list.
 *
 * Parameters:
 * @param optprops - the optical property list
 *
 * Returns:
 * @return -
 */
void
free_optprops_list(optical_property *optprops)
{
  if (optprops != NULL)
    {
      // free the structure
      free(optprops);

      // set it to NULL
      optprops = NULL;
    }
}

/*
 * Function: init_optprops_list
 * The function initializes an optical property list. This means
 * it computes and stores values for the list elements. Here
 * only elements are fixed, which will not change during the
 * whole fringe computation and remain constant for all pixels
 * at all wavelengths.
 *
 * Parameters:
 * @params fconf       - the fringe configuration structure
 * @params lambda_mean - the mean wavelength of analyzed waveband
 * @params optprops    - the list of optical property elements
 *
 * Returns:
 * @return -
 */
void
init_optprops_list(const fringe_conf *fconf, const double lambda_mean,
		   optical_property *optprops)
{
  int num_layers;
  int index;

  double n1;
  double n2;

  // determine the number of layers
  num_layers = fconf->opt_layers->num_layers;

  if (num_layers > 0)
    {
      // set the transmission and reflection
      // vacuum -- first CCD layer
      // that's a bit heuristic, since
      // the true value using the true
      // refration index for vacuum would
      // give a different result
      optprops[0].trans_upper   = 1.0;
      optprops[0].reflect_upper = 0.0;

      // also that one is a bit heuristic
      // but assuming n(vacuum)=1.0 its also true
      optprops[0].sign_upper = 1.0;

      // set the sign of the lowest layer -- substrate
      // this is purely heuristic and works
      // only for the HRC/WFC fringing models!!
      optprops[num_layers-1].sign_lower = -1.0;
    }

  for (index=0; index < num_layers-1; index++)
    {
      // get the real part of the refraction index
      // of the upper layer
      n1 = GSL_REAL(get_complex_refindex(fconf->opt_layers->opt_layer[index],
					 lambda_mean));

      // get the real part of the refraction index
      // of the lower layer
      n2 = GSL_REAL(get_complex_refindex(fconf->opt_layers->opt_layer[index+1],
					 lambda_mean));

      // distribute the signs
      // for the upper and lower layer
      if (n1 < n2)
	{
	  optprops[index].sign_lower   =  1.0;
	  optprops[index+1].sign_upper = -1.0;
	}
      else
	{
	  optprops[index].sign_lower   = -1.0;
	  optprops[index+1].sign_upper =  1.0;
	}
    }
}

/*
 * Function: fill_optprops_thickness
 * The function fills the thickness of the optival layers
 * at a given pixel into a list of optical properties.
 *
 * Parameters:
 * @params opt_layers - the structure for the optical layers
 * @params ii         - the mean wavelength of analyzed waveband
 * @params jj         - the list of optical property elements
 * @params optprops   - the list of optical property elements
 *
 * Returns:
 * @return -
 */
void
fill_optprops_thickness(const ccd_layers *opt_layers, const int ii,
			const int jj, optical_property *optprops)
{
  int index;

  // go over all layers
  for (index=0; index < opt_layers->num_layers; index++)
    // get and store the thickness
    optprops[index].thickness =
      get_layer_thickness(opt_layers->opt_layer[index], ii, jj);
}

/*
 * Function: fill_optprops_all
 * This function computes and fills in all pixel AND wavelength
 * dependent values into a list of optical properties. The completely
 * filled list can then be used to determine the fringe contribution
 * at that wavelength according to a certain model which defines
 * the various light-paths.
 *
 * Parameters:
 * @params opt_layers - the structure for the optical layers
 * @params lambda     - the actual wavelength
 * @params optprops   - the list of optical property elements
 *
 * Returns:
 * @return -
 */
void
fill_optprops_all(const ccd_layers *opt_layers, const double lambda,
		  optical_property *optprops)
{
  int index;

  double phase_number;

  gsl_complex compl_upp;
  gsl_complex compl_low;

  int num_layers;

  // store the number of layers
  num_layers = opt_layers->num_layers;

  // compute the phase number
  phase_number = 2.0*M_PI / lambda;

  // get the complex refractive index for the first layer
  compl_upp =
    get_complex_refindex(opt_layers->opt_layer[0],lambda);

  for (index=0; index < num_layers-1; index++)
    {

      // get the refractive index for the next layer
      compl_low =
      	get_complex_refindex(opt_layers->opt_layer[index+1],lambda);

      // compute the reflection and transmission
      // to the lower layers
      optprops[index].reflect_lower =
	compute_reflection(compl_upp, compl_low);
      optprops[index].trans_lower   = 1.0 - optprops[index].reflect_lower;

      // compute the reflection and transmission
      // to the upper layers
      optprops[index+1].reflect_upper = optprops[index].reflect_lower;
      optprops[index+1].trans_upper   = 1.0 - optprops[index+1].reflect_upper;

      // compute the attenuation
      optprops[index].double_attenuation =
      	compute_attenuation(compl_upp,optprops[index].thickness,phase_number);

      // compute the phase shift
      optprops[index].double_phshift =
	compute_pshift(compl_upp, optprops[index].thickness, phase_number);

      // prepare the next iteration:
      // the now upper layer will then be the lower
      compl_upp = compl_low;
    }
  // compute the attenuation for the last layer
  optprops[num_layers-1].double_attenuation =
    compute_attenuation(compl_upp,optprops[num_layers-1].thickness,phase_number);

  // compute the phase shift for the last layer
  optprops[num_layers-1].double_phshift =
    compute_pshift(compl_upp, optprops[num_layers-1].thickness, phase_number);

  // compute the transmission of the last layer
  // towards the substrate
  optprops[num_layers-1].trans_lower =
    eval_interp(opt_layers->substrate, lambda);

  // compute the reflection of the last layer
  // towards the substrate
  optprops[num_layers-1].reflect_lower =
    1.0 - optprops[num_layers-1].trans_lower;

}

/*
 * Function: fill_contrib_single
 * The function computes the contribution to the fringe amplitude
 * on the basis of an appropriately filled list of optical properties.
 * In the properties the ingredients such as phase shifts and
 * attenuation and transmissions and reflections are defined.
 * This function just combines the ingredients using a model
 * which includes all single reflected beam combinations
 * as contributors to the fringing.
 *
 * Parameters:
 * @params optprops   - the list of optical property elements
 * @params fconf      - the fringe configuration structure
 *
 * Returns:
 * @return (value)    - the fringe contribution
 */
double
fringe_contrib_single(const optical_property *optprops,
		      const fringe_conf *fconf)
{
  int index=2;

  gsl_complex amp_tmp;
  gsl_complex amp_tot;

  // compute the initial phase at the top of the
  // first layer
  amp_tot =
    gsl_complex_mul_real(gsl_complex_exp(gsl_complex_rect(0.0,fconf->fringe_phase*M_PI/180.0)), 1.0);

  // transfer the intitial value
  // to the running variable
  amp_tmp = amp_tot;

  // pile up the amplitude,
  // also putting in the transmission
  amp_tmp =
    gsl_complex_mul(amp_tmp,
		    gsl_complex_mul_real(optprops[0].double_phshift,
					 1.0*1.0*optprops[0].double_attenuation));

  // compute the intermediate result
  // for the first layer
  amp_tot =
    gsl_complex_add(amp_tot,
		    gsl_complex_mul_real(amp_tmp,
					 optprops[0].sign_lower*optprops[0].reflect_lower));

  // go over layer 2 to the end
  for (index=1; index < fconf->opt_layers->num_layers; index++)
    {
      // pile up the amplitude
      // to that layer
      amp_tmp =
	gsl_complex_mul(amp_tmp,
			gsl_complex_mul_real(optprops[index].double_phshift,
					   optprops[index-1].trans_lower*optprops[index].trans_upper*optprops[index].double_attenuation));

      // compute the intermediate result
      // for the that layer
      amp_tot =
	gsl_complex_add(amp_tot,
			gsl_complex_mul_real(amp_tmp,
					     optprops[index].sign_lower*optprops[index].reflect_lower));
      //      if (index==3)
      //	fprintf(stdout,"ccc:%f, %f\n", GSL_REAL(amp_tmp),GSL_IMAG(amp_tmp) );

    }

  // return the final value
  return (GSL_REAL(amp_tot) - 1.0);
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#pragma once

#include <bitset>
#include <type_traits>

#include <thrust/execution_policy.h>
#include <thrust/gather.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/scatter.h>

#include <cuda/define_specifiers.hpp>

#include <gsl-lite/gsl-lite.hpp>

#include <thrustshift/bit.h>
#include <thrustshift/gather.h>
#include <thrustshift/managed-vector.h>
#include <thrustshift/not-a-vector.h>
#include <thrustshift/scatter.h>

namespace thrustshift {

namespace permutation {

/*! \brief Return \f$ P_b P_a^T \f$
 *
 *   \param result `result[i]` is the index of element \f$ a_i \f$
 *     after permutation \f$ P_b \f$ is applied to the natural order.
 */
template <class PermutationA,
          class PermutationB,
          class ResultPermutation,
          class MemoryResource>
void multiply(PermutationA&& a,
              PermutationB&& b,
              ResultPermutation&& result,
              MemoryResource& memory_resource) {

	gsl_Expects(a.size() == b.size());
	gsl_Expects(result.size() == a.size());

	using ResultIndex =
	    typename std::remove_reference<ResultPermutation>::type::value_type;

	auto [__, tmp] =
	    make_not_a_vector_and_span<ResultIndex>(result.size(), memory_resource);

	auto cit_begin = thrust::make_counting_iterator(ResultIndex(0));
	auto cit_end = thrust::make_counting_iterator(
	    gsl_lite::narrow<ResultIndex>(result.size()));

	// Rather use `.data()` in case the iterator class is not defined
	// in device code
	thrust::scatter(thrust::device, cit_begin, cit_end, b.data(), tmp.data());
	thrust::gather(thrust::device,
	               a.data(),
	               a.data() + a.size(),
	               tmp.data(),
	               result.data());
}

/*! \brief Multiply successive permutations
 *
 *  Assume you have \f$N\f$ permutations \f$P_{N-1},...,P_0\f$
 *  and \f$N\f$ linear maps \f$M_{N-1},...,M_0\f$ where each
 *  linear map \f$M_i\f$ operates on the new order described by
 *  \f$P_i\f$ relatively to the natural order. If every
 *  linear map is applied to a vector \f$x\f$ in natural order
 *  it yields to:
 *
 *  \f[
 *      P_{N-1}^T M_{N-1} P_{N-1} P_{N-2}^T M^{N-2} ... P_0^T M_0 P_0 x
 *  \f]
 *
 *  The \f$P_i\f$ permutation matrices are given in vectorial form \f$p_i\f$
 *  such that \f$(P_i)_mn = 1\f$ if \f$(p_i)_m = n\f$.
 *
 *  \param input_permutations Range of size N with ranges of length L each
 *  \param merged_permutations Range of size N - 1 with ranges of length L each.
 *    merged_permutations[i] = \f$P_{i+1} P_{i}^T\f$
 */
template <class InputPermutations,
          class MergedPermutations,
          class MemoryResource>
void multiply_successive(InputPermutations&& input_permutations,
                         MergedPermutations&& multiplied_permutations,
                         MemoryResource& memory_resource) {

	const auto N = input_permutations.size();

	if (N == 0) {
		return;
	}

	gsl_Expects(N > 0);
	gsl_Expects(N == multiplied_permutations.size() + 1);

	for (size_t i = 0; i < input_permutations.size() - 1; ++i) {
		multiply(input_permutations[i],
		         input_permutations[i + 1],
		         multiplied_permutations[i],
		         memory_resource);
	}
}

/*! \brief Bitoptimized successive pairwise permutation
 *
 *  This class saves a permutation of `I = {0,...,N-1}`, which
 *  is created by successive pairwise swaps of the sequence `I`:
 *
 *  ```cpp
 *  for (int j = 1; j < N; ++j) {
 *      if (do_swap[j - 1]) {
 *          swap(I[j - 1], I[j]);
 *      }
 *  }
 *  ```
 *
 *  The array `do_swap[]` saves if two values are swapped and is saved in
 *  a bit pattern within this class.
 *  This can also be seen as a path in a binary tree with `N` levels because
 *  the internal representation is a bit pattern, which sets the bit if the
 *  swap is not done. `N` is limited by the amount of Bits in `BitPatternT`.
 *  If `operator[]` is called once, no function to change the `do_swap`
 *  array should be called by the user. The state of the class is then a
 *  read-only mode. Moreover the array `I` can only be read successively
 *  by the `operator[]`. Thus reading the pattern starts with the last element
 *  `i = N - 1` and in the next call `i = i - 1` must be used as a function argument.
 *  This is an highly unsafe class design, which is chosen to avoid the
 *  usage of additional registers. E.g. the safety of the class can be
 *  greatly improved if one register is added, which is avoided here
 *  due to performance reasons.
 *
 *  \note When `std::bitset` is supported for device code. This should be
 *        used to save the bit pattern
 *  \note When compiled for device code CUDA currently supports built-in
 *        function to count leading zeros `clz` for `int` and `long long int`, which
 *        then must be used for `BitPatternT`.
 */
template <class BitPatternT>
class bit_successive_permutation_t {

   public:
	CUDA_FHD bit_successive_permutation_t(int N)
	    : bit_pattern_(BitPatternT(1) << (N - 1)) {
		gsl_Expects(N <= sizeof(BitPatternT) * 8);
#ifndef NDEBUG
		N_ = N;
#endif
	}

	CUDA_FHD void do_swap(int j) {
		gsl_Expects(j < sizeof(BitPatternT) * 8);
#ifndef NDEBUG
		gsl_Expects(!read_only_mode_);
		gsl_Expects(j < N_ - 1);
#endif
		bit_pattern_ &= ~(BitPatternT(1) << j);
	}

	CUDA_FHD void do_not_swap(int j) {
		gsl_Expects(j < sizeof(BitPatternT) * 8);
#ifndef NDEBUG
		gsl_Expects(!read_only_mode_);
		gsl_Expects(j < N_ - 1);
#endif
		bit_pattern_ |= (BitPatternT(1) << j);
	}

	CUDA_FHD void set(int j, bool do_swap_) {
		if (do_swap_) {
			do_swap(j);
		}
		else {
			do_not_swap(j);
		}
	}

	CUDA_FHD int operator[](int i) {
		gsl_Expects(i < sizeof(BitPatternT) * 8);
#ifndef NDEBUG
		gsl_Expects(i < N_);
		if (!read_only_mode_) {
			read_only_mode_ = true;
			last_i_ = N_;
		}
		gsl_Expects(i == last_i_ - 1);
		last_i_ = i;
#endif
		gsl_Expects(i >= 0);
		const bool flag = (1ll << i) & bit_pattern_;
		bit_pattern_ &= ~(1ll << i);

		if (flag) {
			const int lz = count_leading_zeros(bit_pattern_);
			return sizeof(BitPatternT) * 8 - lz;
		}
		else {
			return i + 1;
		}
	}

   private:
#ifndef NDEBUG
	bool read_only_mode_ = false;
	int N_;
	int last_i_;
#endif
	BitPatternT bit_pattern_;
};

} // namespace permutation

} // namespace thrustshift

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
/*
 * halo.h
 *
 *  Created on: Oct 23, 2012
 *      Author: cgiocoli
 *
 *      This class implement the non-linear power spectrum using the Halo Model
 */

#ifndef POWERCDMHM_H_
#define POWERCDMHM_H_
#include <math.h>
#include <cosmo.h>
#include <halo.h>
#include <utilities.h>

#ifdef ENABLE_GSL

#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_rng.h>

/**
 * \brief Class for calculating the non-linear power spectrum using ****** Halo Model.
 *  To use this class you must enable GSL.
 */
class POWERCDMHM{
  public:
  POWERCDMHM(COSMOLOGY *co, double redshift, double minHalomass, double sigma=0.25,int cmRelation=0,double slope=-0.1);

  double nonlinpowerCDMHM(double k);
  double nonlinpowerCDMHM1Halo(double k);
  double nonlinpowerCDMHM2Halo(double k);
  double nonlinKAPPApowerCDMHM(double l, double zs);
  double nonlinKAPPApowerCDMHM1Halo(double l, double zs);
  double nonlinKAPPApowerCDMHM2Halo(double l, double zs);
  double linKAPPApowerCDMHM(double l, double zs);
  double nonlinfitKAPPApowerCDMHM(double l, double zs);
  virtual ~POWERCDMHM ();

  protected:
  gsl_function intPk1,intPk2;
  double Pklin,Pk1,Pk2;
  double Pk20;
  double weight (double z1, double z2);
  double weight (double z0);

  struct weightfnc{
    std::vector<double> ai;
    std::vector<double> wi;
  };
  weightfnc wgf; // weight function
  void Initweight(double zs);
  private:
  float *xf,*wf;
};
#endif

#endif /* POWERCDMHM_H_ */

max_stars_repo_name: glenco/CosmoLib

max_issues_repo_name: glenco/CosmoLib

max_forks_repo_name: glenco/CosmoLib


text:
#if !defined(PETIGAFTN_H)
#define PETIGAFTN_H

#include <petsc.h>

EXTERN_C_BEGIN
extern void IGA_Quadrature_1D(PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
                              PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_Quadrature_2D(PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
                              PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
                              PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_Quadrature_3D(PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
                              PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
                              PetscInt,const PetscReal[],const PetscReal[],const PetscReal*,
                              PetscReal[],PetscReal[],PetscReal[]);
EXTERN_C_END

EXTERN_C_BEGIN
extern void IGA_BasisFuns_1D(PetscInt,
                             PetscInt,PetscInt,const PetscReal[],
                             PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_BasisFuns_2D(PetscInt,
                             PetscInt,PetscInt,const PetscReal[],
                             PetscInt,PetscInt,const PetscReal[],
                             PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_BasisFuns_3D(PetscInt,
                             PetscInt,PetscInt,const PetscReal[],
                             PetscInt,PetscInt,const PetscReal[],
                             PetscInt,PetscInt,const PetscReal[],
                             PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
EXTERN_C_END

EXTERN_C_BEGIN
extern void IGA_Rationalize_1D(PetscInt,PetscInt,PetscInt,const PetscReal[],
                               PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_Rationalize_2D(PetscInt,PetscInt,PetscInt,const PetscReal[],
                               PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_Rationalize_3D(PetscInt,PetscInt,PetscInt,const PetscReal[],
                               PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
EXTERN_C_END

EXTERN_C_BEGIN
extern void IGA_GeometryMap_1D(PetscInt,PetscInt,PetscInt,const PetscReal[],
                               const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                               PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_GeometryMap_2D(PetscInt,PetscInt,PetscInt,const PetscReal[],
                               const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                               PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_GeometryMap_3D(PetscInt,PetscInt,PetscInt,const PetscReal[],
                               const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                               PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_GeometryMap   (PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,const PetscReal[],
                               const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                               PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
EXTERN_C_END

EXTERN_C_BEGIN
extern void IGA_InverseMap_1D(PetscInt,PetscInt,
                              const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                              PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_InverseMap_2D(PetscInt,PetscInt,
                              const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                              PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_InverseMap_3D(PetscInt,PetscInt,
                              const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                              PetscReal[],PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
EXTERN_C_END

EXTERN_C_BEGIN
extern void IGA_ShapeFuns_1D(PetscInt,PetscInt,PetscInt,
                             const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                             const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                             PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_ShapeFuns_2D(PetscInt,PetscInt,PetscInt,
                             const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                             const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                             PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
extern void IGA_ShapeFuns_3D(PetscInt,PetscInt,PetscInt,
                             const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                             const PetscReal[],const PetscReal[],const PetscReal[],const PetscReal[],
                             PetscReal[],PetscReal[],PetscReal[],PetscReal[]);
EXTERN_C_END

#endif/*PETIGAFTN_H*/

max_stars_repo_name: dalcinl/PetIGA

max_issues_repo_name: dalcinl/PetIGA

max_forks_repo_name: dalcinl/PetIGA


text:
/**
* Copyright 2016 BitTorrent Inc.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
*    http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#pragma once

#include <scraps/config.h>

#if SCRAPS_APPLE
    #include <scraps/apple/HMAC.h>
    namespace scraps { using HMACSHA256 = apple::HMAC<apple::HMACAlgorithmType::SHA256>; }
#else
    #include <scraps/sodium/HMACSHA256.h>
    namespace scraps { using HMACSHA256 = sodium::HMACSHA256; }
#endif

#include <gsl.h>

namespace scraps {

template <typename BaseByteType>
struct HMACSHA256ByteTag {};

template <typename BaseByteType>
using HMACSHA256Byte = StrongByte<HMACSHA256ByteTag<BaseByteType>>;

template <typename KeyByteType, std::ptrdiff_t KeySize, typename ByteT, std::ptrdiff_t DataSize>
std::array<HMACSHA256Byte<std::remove_const_t<ByteT>>, HMACSHA256::kResultSize>
GetHMACSHA256(gsl::span<KeyByteType, KeySize> key, gsl::span<ByteT, DataSize> data) {
    std::array<HMACSHA256Byte<std::remove_const_t<ByteT>>, HMACSHA256::kResultSize> ret;

    HMACSHA256 hmac(key.data(), key.size());
    hmac.update(data.data(), data.size());
    hmac.finish(ret.data());

    return ret;
}

} // namespace scraps

max_stars_repo_name: carlbrown/scraps

max_issues_repo_name: carlbrown/scraps

max_forks_repo_name: carlbrown/scraps


text:
/**
 * author: Jochen K"upper
 * created: Jan 2002
 * file: pygsl/src/statistics/longmodule.c
 * $Id: intmodule.c,v 1.6 2004/03/24 08:40:45 schnizer Exp $
 *
 *
 *"
 */

#include <Python.h>
#include <pygsl/error_helpers.h>
#include <pygsl/block_helpers.h>
#include <gsl/gsl_statistics.h>

/* include real functions for different data-types */

#define STATMOD_APPEND_PY_TYPE(X) X ## Int
#define STATMOD_APPEND_PYC_TYPE(X) X ## INT
#define STATMOD_FUNC_EXT(X, Y) X ## _int ## Y
#define STATMOD_PY_AS_C PyInt_AsLong
#define STATMOD_C_TYPE int
#define PyGSL_STATISTICS_IMPORT_API
#include "functions.c"

/* initialization */

PyGSL_STATISTICS_INIT(int, "int")

/*
 * Local Variables:
 * mode: c
 * c-file-style: "Stroustrup"
 * End:
 */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod