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 dnest_sa1d.c
 *  \brief run dnest sampling for sa and 1d rm analysis.
 */

#ifdef SA

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <stddef.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_interp.h>
#include <mpi.h>

#include "brains.h"

DNestFptrSet *fptrset_sa1d;

/*
 * run denst for SA and 1D RM analysis.
 */
int dnest_sa1d(int argc, char **argv)
{
  int i;
  double logz_sa1d;

  if(parset.flag_sa_par_mutual == 0) /* SA and RM have the same BLR */
  {
    set_blr_model2d();
    num_params_sa_blr_model = 0;
  }
  else /* SA and RM have different BLRs but share the same inc. */
  {
    set_blr_model1d();
    set_sa_blr_model();
  }
  /* RM */
  num_params_blr = num_params_blr_model + 1; /* include line sys err */
  num_params_rm = parset.n_con_recon + num_params_var + num_params_blr;

  /* SA */
  num_params_sa_blr = num_params_sa_blr_model + num_params_sa_extpar;
  num_params_sa = num_params_sa_blr;

  /* total */
  num_params_blr_tot = num_params_blr + num_params_sa_blr;
  num_params = num_params_sa + num_params_rm;
  idx_resp = num_params_blr_tot + num_params_drw + num_params_trend;
  idx_difftrend = idx_resp + num_params_resp;

  par_fix = (int *) malloc(num_params * sizeof(int));
  par_fix_val = (double *) malloc(num_params * sizeof(double));
  par_range_model = malloc( num_params * sizeof(double *));
  par_prior_gaussian = malloc( num_params * sizeof(double *));
  for(i=0; i<num_params; i++)
  {
    par_range_model[i] = malloc(2*sizeof(double));
    par_prior_gaussian[i] = malloc(2*sizeof(double));
  }
  par_prior_model = malloc( num_params * sizeof(int));

  fptrset_sa1d = dnest_malloc_fptrset();
  /* setup functions used for dnest*/
  fptrset_sa1d->from_prior = from_prior_sa1d;
  fptrset_sa1d->print_particle = print_particle_sa1d;
  fptrset_sa1d->restart_action = restart_action_sa;
  fptrset_sa1d->accept_action = accept_action_sa1d;
  fptrset_sa1d->kill_action = kill_action_sa1d;
  fptrset_sa1d->perturb = perturb_sa1d;
  fptrset_sa1d->read_particle = read_particle_sa1d;

  fptrset_sa1d->log_likelihoods_cal_initial = log_likelihoods_cal_initial_sa1d;
  fptrset_sa1d->log_likelihoods_cal_restart = log_likelihoods_cal_restart_sa1d;
  fptrset_sa1d->log_likelihoods_cal = log_likelihoods_cal_sa1d;

  set_par_range_sa1d();
  set_par_fix_blrmodel();
  set_par_fix_sa_blrmodel();

  /* the rest parameters */
  for(i=num_params_blr_model; i<num_params_blr; i++)
  {
    par_fix[i] = 0;
    par_fix_val[i] = -DBL_MAX;
  }
  for(i=num_params_blr + num_params_sa_blr_model; i<num_params; i++)
  {
    par_fix[i] = 0;
    par_fix_val[i] = -DBL_MAX;
  }

  /* fix non-linear response */
  if(parset.flag_nonlinear !=1)
  {
    par_fix[idx_resp+1] = 1;
    par_fix_val[idx_resp+1] = 0.0;
  }

  /* fix systematic error of line */
  if(parset.flag_line_sys_err != 1)
  {
    par_fix[num_params_blr-1] = 1;
    par_fix_val[num_params_blr-1] = log(1.0);
  }

  /* fix inc of SA BLR if SA and RM have different BLR
   * in this case, DA must be fixed.
   */
  if(parset.flag_sa_par_mutual != 0)
  {
    /* get the index for mbh and inclination parameters */
    set_idx_par_mutual();

    /* inc */
    par_fix[num_params_blr + idx_sa_par_mutual[1]] = 1;
    par_fix_val[num_params_blr + idx_sa_par_mutual[1]] = 0.0;

    /* fix DA */
    par_fix[num_params_blr + num_params_sa_blr_model] = 1;
    par_fix_val[num_params_blr + num_params_sa_blr_model] = log(550.0);
  }

  /* fix FA */
  par_fix[num_params_blr + num_params_sa_blr_model+2] = 1;
  par_fix_val[num_params_blr + num_params_sa_blr_model+2] = log(1.0);

  /* fix systematic error of continuum */
  if(parset.flag_con_sys_err != 1)
  {
    par_fix[num_params_blr_tot] = 1;
    par_fix_val[num_params_blr_tot] = log(1.0);
  }

  /* fix continuum variation parameter sigma and tau if flag_fixvar is true */
  if(parset.flag_fixvar == 1)
  {
    par_fix[num_params_blr_tot + 1] = 1;
    par_fix_val[num_params_blr_tot + 1] = var_param[1];
    par_fix[num_params_blr_tot + 2] = 1;
    par_fix_val[num_params_blr_tot + 2] = var_param[2];
  }

  print_par_names_sa1d();

  force_update = parset.flag_force_update;
  if(parset.flag_para_name != 1)
    logz_sa1d = dnest(argc, argv, fptrset_sa1d, num_params, dnest_options_file);

  dnest_free_fptrset(fptrset_sa1d);

  return 0;
}

void set_par_range_sa1d()
{
  int i, j, i1, i2;

  /* RM BLR parameters first */
  for(i=0; i<num_params_blr_model; i++)
  {
    par_range_model[i][0] = blr_range_model[i][0];
    par_range_model[i][1] = blr_range_model[i][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;
  }

  /* systematic line error */
  i = num_params_blr -1;
  par_range_model[i][0] = sys_err_line_range[0];
  par_range_model[i][1] = sys_err_line_range[1];

  par_prior_model[i] = UNIFORM;
  par_prior_gaussian[i][0] = 0.0;
  par_prior_gaussian[i][1] = 0.0;

  /* SA BLR parameters */
  for(i=num_params_blr; i<num_params_blr + num_params_sa_blr_model; i++)
  {
    par_range_model[i][0] = sa_blr_range_model[i-num_params_blr][0];
    par_range_model[i][1] = sa_blr_range_model[i-num_params_blr][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;
  }
  /* SA extra parameters */
  for(i=num_params_blr + num_params_sa_blr_model; i<num_params_blr_tot; i++)
  {
    par_range_model[i][0] = sa_extpar_range[i-(num_params_blr + num_params_sa_blr_model)][0];
    par_range_model[i][1] = sa_extpar_range[i-(num_params_blr + num_params_sa_blr_model)][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;
  }

  /* variability parameters */
  /* first systematic error */
  i = num_params_blr_tot;
  par_range_model[i][0] = var_range_model[0][0];
  par_range_model[i][1] = var_range_model[0][1];

  par_prior_model[i] = UNIFORM;
  par_prior_gaussian[i][0] = 0.0;
  par_prior_gaussian[i][1] = 0.0;

  /* drw parameters */
  for(i=num_params_blr_tot+1; i<num_params_drw + num_params_blr_tot; i++)
  {
    if(var_param_std[i-num_params_blr_tot] > 0.0)
    {
      par_range_model[i][0] = var_param[i-num_params_blr_tot] - 5.0 * var_param_std[i-num_params_blr_tot];
      par_range_model[i][1] = var_param[i-num_params_blr_tot] + 5.0 * var_param_std[i-num_params_blr_tot];

      /* make sure that the range lies within the initial range */
      par_range_model[i][0] = fmax(par_range_model[i][0], var_range_model[i-num_params_blr_tot][0]);
      par_range_model[i][1] = fmin(par_range_model[i][1], var_range_model[i-num_params_blr_tot][1]);

      par_prior_model[i] = GAUSSIAN;
      par_prior_gaussian[i][0] = var_param[i-num_params_blr_tot];
      par_prior_gaussian[i][1] = var_param_std[i-num_params_blr_tot];
    }
    else
    {
      par_range_model[i][0] = var_range_model[i-num_params_blr_tot][0];
      par_range_model[i][1] = var_range_model[i-num_params_blr_tot][1];

      par_prior_model[i] = UNIFORM;
      par_prior_gaussian[i][0] = 0.0;
      par_prior_gaussian[i][1] = 0.0;
    }
  }
  /* long-term trend */
  for(i=num_params_drw + num_params_blr_tot; i< num_params_drw + num_params_trend + num_params_blr_tot; i++)
  {
    par_range_model[i][0] = var_range_model[3][0];
    par_range_model[i][1] = var_range_model[3][1];

    par_prior_model[i] = GAUSSIAN;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0;
  }

  /* response A and Ag */
  j = 0;
  i1 = idx_resp;
  i2 = idx_resp + num_params_resp;
  for(i=i1; i<i2; i++)
  {
    par_range_model[i][0] = resp_range[j][0];
    par_range_model[i][1] = resp_range[j][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;

    j++;
  }

  /* different trend */
  j = 0;
  i1 = idx_difftrend;
  i2 = idx_difftrend + num_params_difftrend;
  for(i=i1; i< i2; i++)
  {
    par_range_model[i][0] = var_range_model[4 + j][0];
    par_range_model[i][1] = var_range_model[4 + j][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;

    j++;
  }

  /* continuum light curve parameters */
  for(i=num_params_blr_tot+num_params_var; i<num_params; i++)
  {
    par_range_model[i][0] = var_range_model[4+num_params_difftrend][0];
    par_range_model[i][1] = var_range_model[4+num_params_difftrend][1];

    par_prior_model[i] = GAUSSIAN;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0;
  }
  return;
}

/*!
 *  print names and prior ranges for parameters
 *
 */
void print_par_names_sa1d()
{
  if(thistask!= roottask)
    return;

  int i, j;
  FILE *fp;
  char fname[BRAINS_MAX_STR_LENGTH], str_fmt[BRAINS_MAX_STR_LENGTH];

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

  strcpy(str_fmt, "%4d %-15s %10.6f %10.6f %4d %4d %15.6e\n");

  printf("# Print parameter name in %s\n", fname);

  fprintf(fp, "#*************************************************\n");
  fprint_version(fp);
  fprintf(fp, "#*************************************************\n");

  fprintf(fp, "%4s %-15s %10s %10s %4s %4s %15s\n", "#", "Par", "Min", "Max", "Prior", "Fix", "Val");

  i=-1;
  for(j=0; j<num_params_blr_model; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "BLR model", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  i++;
  fprintf(fp, str_fmt, i, "sys_err_line", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  for(j=0; j<num_params_sa_blr_model; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "SA BLR model", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  for(j=0; j<num_params_sa_extpar; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "SA Extra Par", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  i++;
  fprintf(fp, str_fmt, i, "sys_err_con", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  i++;
  fprintf(fp, str_fmt, i, "sigmad", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  i++;
  fprintf(fp, str_fmt, i, "taud", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  for(j=0; j<num_params_trend; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "trend", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  i++;
  fprintf(fp, str_fmt, i, "A", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  i++;
  fprintf(fp, str_fmt, i, "Ag", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  for(j=0; j<num_params_difftrend; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "diff trend", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  for(j=0; j<parset.n_con_recon; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "time series", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  fclose(fp);
}

/*!
 * this function generates a sample from prior.
 */
void from_prior_sa1d(void *model)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<num_params; i++)
  {
    if(par_prior_model[i] == GAUSSIAN)
    {
      pm[i] = dnest_randn()*par_prior_gaussian[i][1] + par_prior_gaussian[i][0];
      dnest_wrap(&pm[i], par_range_model[i][0], par_range_model[i][1]);
    }
    else
    {
      pm[i] = par_range_model[i][0] + dnest_rand() * (par_range_model[i][1] - par_range_model[i][0]);
    }
  }

  /* cope with fixed parameters */
  for(i=0; i<num_params; i++)
  {
    if(par_fix[i] == 1)
      pm[i] = par_fix_val[i];
  }

  which_parameter_update = -1;

  return;
}

/*!
 * this function calculates likelihood at initial step.
 */
double log_likelihoods_cal_initial_sa1d(const void *model)
{
  double logL;
  logL = prob_initial_sa1d(model);
  return logL;
}

/*!
 * this function calculates likelihood at initial step.
 */
double log_likelihoods_cal_restart_sa1d(const void *model)
{
  double logL;
  logL = prob_initial_sa1d(model);
  return logL;
}

/*!
 * this function calculates likelihood.
 */
double log_likelihoods_cal_sa1d(const void *model)
{
  double logL;
  logL = prob_sa1d(model);
  return logL;
}

/*!
 * this function prints out parameters.
 */
void print_particle_sa1d(FILE *fp, const void *model)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<num_params; i++)
  {
    fprintf(fp, "%e ", pm[i] );
  }
  fprintf(fp, "\n");
  return;
}

/*!
 * This function read the particle from the file.
 */
void read_particle_sa1d(FILE *fp, void *model)
{
  int j;
  double *psample = (double *)model;

  for(j=0; j < dnest_num_params; j++)
  {
    if(fscanf(fp, "%lf", psample+j) < 1)
    {
      printf("%f\n", *psample);
      fprintf(stderr, "#Error: Cannot read file %s.\n", options.sample_file);
      exit(0);
    }
  }
  return;
}

/*!
 * this function perturbs parameters.
 */
double perturb_sa1d(void *model)
{
  double *pm = (double *)model;
  double logH = 0.0, limit1, limit2, width, rnd;
  int which, which_level;

  /*
   * fixed parameters need not to update
   * perturb important parameters more frequently
   */
  do
  {
    rnd = dnest_rand();
    if(rnd < fmax(0.2, 1.0*(num_params_blr_tot+num_params_var)/num_params))
      which = dnest_rand_int(num_params_blr_tot + num_params_var);
    else
      which = dnest_rand_int(parset.n_con_recon) + num_params_blr_tot + num_params_var;

  }while(par_fix[which] == 1);

  which_parameter_update = which;

  /* level-dependent width */
  count_saves = dnest_get_count_saves();
  which_level_update = dnest_get_which_level_update();
  which_level = which_level_update > (size_levels-10)?(size_levels-10):which_level_update;

  if( which_level > 0)
  {
    limit1 = limits[(which_level-1) * num_params *2 + which *2];
    limit2 = limits[(which_level-1) * num_params *2 + which *2 + 1];
    width = (limit2 - limit1);
  }
  else
  {
    limit1 = par_range_model[which][0];
    limit2 = par_range_model[which][1];
    width = (par_range_model[which][1] - par_range_model[which][0]);
  }
  width /= (2.35);

  if(par_prior_model[which] == GAUSSIAN)
  {
    logH -= (-0.5*pow((pm[which] - par_prior_gaussian[which][0])/par_prior_gaussian[which][1], 2.0) );
    pm[which] += dnest_randh() * width;
    dnest_wrap(&pm[which], par_range_model[which][0], par_range_model[which][1]);
    //dnest_wrap(&pm[which], limit1, limit2);
    logH += (-0.5*pow((pm[which] - par_prior_gaussian[which][0])/par_prior_gaussian[which][1], 2.0) );
  }
  else
  {
    pm[which] += dnest_randh() * width;
    dnest_wrap(&(pm[which]), par_range_model[which][0], par_range_model[which][1]);
    //dnest_wrap(&pm[which], limit1, limit2);
  }

  return logH;
}

/*!
 * this function calculates likelihood.
 */
double log_likelihoods_cal_sa1d_exam(const void *model)
{
  return 0.0;
}

void accept_action_sa1d()
{
  int param;
  double *ptemp;

  // the parameter previously updated
  param = which_parameter_update;

  /* continuum parameter is updated */
  if( param >= num_params_blr_tot )
  {
    /*
     *note that (response) Fline is also changed as long as Fcon is changed.
     *num_params_blr-the parameter is the systematic error of continuum.
     *the change of this parameter also changes continuum reconstruction.
     */

    ptemp = Fcon_rm_particles[which_particle_update];
    Fcon_rm_particles[which_particle_update] = Fcon_rm_particles_perturb[which_particle_update];
    Fcon_rm_particles_perturb[which_particle_update] = ptemp;

    if(force_update != 1)
    {
      ptemp = Fline_at_data_particles[which_particle_update];
      Fline_at_data_particles[which_particle_update] = Fline_at_data_particles_perturb[which_particle_update];
      Fline_at_data_particles_perturb[which_particle_update] = ptemp;
    }
  }
  else if( param != num_params_blr-1 && force_update != 1)
  {
    /* RM and SA BLR parameter is updated
     * Note a) that the (num_par_blr-1)-th parameter is systematic error of line.
     * when this parameter is updated, Trans1D and Fline, phase_sa, and Fline_sa are unchanged.
     *      b) Fline is always changed, except param = num_params_blr-1.
     */

    {
      ptemp = TransTau_particles[which_particle_update];
      TransTau_particles[which_particle_update] = TransTau_particles_perturb[which_particle_update];
      TransTau_particles_perturb[which_particle_update] = ptemp;

      ptemp = Trans1D_particles[which_particle_update];
      Trans1D_particles[which_particle_update] = Trans1D_particles_perturb[which_particle_update];
      Trans1D_particles_perturb[which_particle_update] = ptemp;

      ptemp = Fline_at_data_particles[which_particle_update];
      Fline_at_data_particles[which_particle_update] = Fline_at_data_particles_perturb[which_particle_update];
      Fline_at_data_particles_perturb[which_particle_update] = ptemp;

      prob_sa_particles[which_particle_update] = prob_sa_particles_perturb[which_particle_update];
    }
  }

  return;
}

/*
 * action when particle i is killed in cdnest sampling.
 * particle i_copy's properties is copyed to particle i.
 */
void kill_action_sa1d(int i, int i_copy)
{
  memcpy(Fcon_rm_particles[i], Fcon_rm_particles[i_copy], parset.n_con_recon * sizeof(double));
  memcpy(Fline_at_data_particles[i], Fline_at_data_particles[i_copy], n_line_data * sizeof(double));
  memcpy(TransTau_particles[i], TransTau_particles[i_copy], parset.n_tau*sizeof(double));
  memcpy(Trans1D_particles[i], Trans1D_particles[i_copy], parset.n_tau * sizeof(double));

  prob_sa_particles[i] = prob_sa_particles[i_copy];

  return;
}

#endif

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS


text:
/**
 *
 * @file qwrapper_dlaset.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:58 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlaset(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int M, int N,
                       double alpha, double beta,
                       double *A, int LDA)
{
    DAG_CORE_LASET;
    QUARK_Insert_Task(quark, CORE_dlaset_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(double),         &alpha, VALUE,
        sizeof(double),         &beta,  VALUE,
        sizeof(double)*LDA*N,    A,      OUTPUT,
        sizeof(int),                        &LDA,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlaset_quark = PCORE_dlaset_quark
#define CORE_dlaset_quark PCORE_dlaset_quark
#endif
void CORE_dlaset_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int M;
    int N;
    double alpha;
    double beta;
    double *A;
    int LDA;

    quark_unpack_args_7(quark, uplo, M, N, alpha, beta, A, LDA);
    LAPACKE_dlaset_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        M, N, alpha, beta, A, LDA);
}

ext: c


text:
/* cdf/gammainv.c
 *
 * Copyright (C) 2003, 2007 Brian Gough
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>

#include <stdio.h>

double
gsl_cdf_gamma_Pinv (const double P, const double a, const double b)
{
  double x;

  if (P == 1.0)
    {
      return GSL_POSINF;
    }
  else if (P == 0.0)
    {
      return 0.0;
    }

  /* Consider, small, large and intermediate cases separately.  The
     boundaries at 0.05 and 0.95 have not been optimised, but seem ok
     for an initial approximation.

     BJG: These approximations aren't really valid, the relevant
     criterion is P*gamma(a+1) < 1. Need to rework these routines and
     use a single bisection style solver for all the inverse
     functions.
  */

  if (P < 0.05)
    {
      double x0 = exp ((gsl_sf_lngamma (a) + log (P)) / a);
      x = x0;
    }
  else if (P > 0.95)
    {
      double x0 = -log1p (-P) + gsl_sf_lngamma (a);
      x = x0;
    }
  else
    {
      double xg = gsl_cdf_ugaussian_Pinv (P);
      double x0 = (xg < -0.5*sqrt (a)) ? a : sqrt (a) * xg + a;
      x = x0;
    }

  /* Use Lagrange's interpolation for E(x)/phi(x0) to work backwards
     to an improved value of x (Abramowitz & Stegun, 3.6.6)

     where E(x)=P-integ(phi(u),u,x0,x) and phi(u) is the pdf.
   */

  {
    double lambda, dP, phi;
    unsigned int n = 0;

  start:
    dP = P - gsl_cdf_gamma_P (x, a, 1.0);
    phi = gsl_ran_gamma_pdf (x, a, 1.0);

    if (dP == 0.0 || n++ > 32)
      goto end;

    lambda = dP / GSL_MAX (2 * fabs (dP / x), phi);

    {
      double step0 = lambda;
      double step1 = -((a - 1) / x - 1) * lambda * lambda / 4.0;

      double step = step0;
      if (fabs (step1) < 0.5 * fabs (step0))
        step += step1;

      if (x + step > 0)
        x += step;
      else
        {
          x /= 2.0;
        }

      if (fabs (step0) > 1e-10 * x || fabs(step0 * phi) > 1e-10 * P)
        goto start;
    }

  end:
    if (fabs(dP) > GSL_SQRT_DBL_EPSILON * P)
      {
        GSL_ERROR_VAL("inverse failed to converge", GSL_EFAILED, GSL_NAN);
      }

    return b * x;
  }
}

double
gsl_cdf_gamma_Qinv (const double Q, const double a, const double b)
{
  double x;

  if (Q == 1.0)
    {
      return 0.0;
    }
  else if (Q == 0.0)
    {
      return GSL_POSINF;
    }

  /* Consider, small, large and intermediate cases separately.  The
     boundaries at 0.05 and 0.95 have not been optimised, but seem ok
     for an initial approximation. */

  if (Q < 0.05)
    {
      double x0 = -log (Q) + gsl_sf_lngamma (a);
      x = x0;
    }
  else if (Q > 0.95)
    {
      double x0 = exp ((gsl_sf_lngamma (a) + log1p (-Q)) / a);
      x = x0;
    }
  else
    {
      double xg = gsl_cdf_ugaussian_Qinv (Q);
      double x0 = (xg < -0.5*sqrt (a)) ? a : sqrt (a) * xg + a;
      x = x0;
    }

  /* Use Lagrange's interpolation for E(x)/phi(x0) to work backwards
     to an improved value of x (Abramowitz & Stegun, 3.6.6)

     where E(x)=P-integ(phi(u),u,x0,x) and phi(u) is the pdf.
   */

  {
    double lambda, dQ, phi;
    unsigned int n = 0;

  start:
    dQ = Q - gsl_cdf_gamma_Q (x, a, 1.0);
    phi = gsl_ran_gamma_pdf (x, a, 1.0);

    if (dQ == 0.0 || n++ > 32)
      goto end;

    lambda = -dQ / GSL_MAX (2 * fabs (dQ / x), phi);

    {
      double step0 = lambda;
      double step1 = -((a - 1) / x - 1) * lambda * lambda / 4.0;

      double step = step0;
      if (fabs (step1) < 0.5 * fabs (step0))
        step += step1;

      if (x + step > 0)
        x += step;
      else
        {
          x /= 2.0;
        }

      if (fabs (step0) > 1e-10 * x)
        goto start;
    }

  }

end:
  return b * x;
}


text:
/************************************************************/
/*********  Interval Censoring data analysis   ***********/
/*********  by Hao Liu on 6/7/07                  **********/
/************************************************************/

/*** accepted data file: first p colmuns are covariates, last 2 columns are observed time intervals, 9999=censored ***/

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sort_double.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_vector.h>

#define MAX_ITER 10000
#define TOL 1.0e-6
#define sim_num 300

/** sim_num is the bootstrap replicates **/

#define machine_p 1.0e-15
#define myPHO 0.5
#define myPI 0.01
#define Sigma 20.0

struct fun_params {
    double *dij;
    double *zij;
    double *pj;
    int nn;
    int np;
    int mk;
};

double fdf_eplike(int, const double *, double *, void *);
void myxij(double *, const double *, const double *, void *);
void deveplike(double *, double *, const double *, double *, double *, void *);

double maxp(double *, double *, double, const double *, double *, void *);
double normf(double *, double *, double, double, const double *, double *, void *);

extern int conmin(int, double *, double *, double *, int *, int *, double, int, double *, int, double, int,
                  double (*)(int, const double *, double *, void *), void *);

int main(int argc, char *argv[])
{
    FILE *fpi, *fpo1, *fpo2;
    if(argc<=4) {
        printf("need: 1) data file name, 2) sample size n, 3) number of covariate p, and 4) the output file1 5) files2...\n");
        return 0;
    }
    fpi=fopen(argv[1], "r");
    fpo1=fopen(argv[4], "w");
    fpo2=fopen(argv[5], "w");

    int nn=0, np=0, i, j;
    char *cs=malloc(1000*sizeof(char));

    sscanf(argv[2], "%i", &nn);
    sscanf(argv[3], "%i", &np);

    if (np>nn) {
        printf("sample size should be larger than the number of covariates.");
        return 0;
    }

    double *li= malloc(nn*sizeof(double)), *lia= malloc(nn*sizeof(double)), *li0=malloc(nn*sizeof(double));
    double *ri= malloc(nn*sizeof(double)), *ria= malloc(nn*sizeof(double)), *ri0=malloc(nn*sizeof(double));
    double *zij= malloc(nn*np*sizeof(double)), *zij0= malloc(nn*np*sizeof(double)), *dij=malloc(nn*nn*sizeof(double));
    double *q=malloc(nn*sizeof(double)), *p=malloc(nn*sizeof(double)), *gam= malloc(nn*nn*sizeof(double));

    double *pj=calloc(nn, sizeof(double)), *mypj=calloc(nn, sizeof(double));

    double alpha0=1.0, theta0=1.0, nvold;
    register double *tau= calloc(nn, sizeof(double));
    double past_time, now_time, mysize, tmp00;

    double curt, wt, u, likval1, likval2, likval;

	int j1, k, mk, ss, iter, status, iter2, cure;
    int *bsi=malloc(nn*sizeof(int)), *bsi2=malloc(nn*sizeof(int));

    double *ww=calloc(2*(np+1)*(np+8), sizeof(double)),
        *xx=calloc(2*(np+1), sizeof(double)), *gg=calloc(2*(np+1), sizeof(double)), fval;
    double acc=1e-11, eps=1e-5;;
    int nflag=-1, ifun=-1, ifunter=-1;

    struct fun_params *mydata=malloc(sizeof(struct fun_params));

    const gsl_rng_type *Ty;
    gsl_rng *rn;

    gsl_rng_env_setup();
    Ty = gsl_rng_default;
    rn = gsl_rng_alloc (Ty);
    gsl_rng_set(rn, 100);

    /** read the data, first line is the hearder and thus removed ***/
    fgets(cs, 1000, fpi);

    for (i = 0; i<nn; i++) {
        for (j=0;j<np;j++) fscanf(fpi, "%lf,", zij0+np*i+j);
        fscanf(fpi, "%lf,%lf", li0+i, ri0+i);
        bsi[i]=i;
    }
    fclose(fpi);
    mydata->nn = nn;
    mydata->np = np;

    /** sim_num is the bootstrap replicates **/

    for (ss=0;ss<sim_num;ss++){

        /** select the bootstrap sample ***/
        /*** random sample with replacement **/

        gsl_ran_sample (rn, bsi2, nn, bsi, nn, sizeof(int));
        for(i=0;i<nn;i++) {
            for (j=0;j<np;j++) zij[np*i+j] = zij0[np*bsi2[i]+j];
            lia[i]=li0[bsi2[i]];
            li[i]=li0[bsi2[i]];
            ria[i]=ri0[bsi2[i]];
            ri[i]=ri0[bsi2[i]];
        }

        mydata->zij = zij;
        gsl_sort(ria,1,nn);
        gsl_sort(lia,1,nn);

        /*** generate p, q data ***/
        p[0]=ria[0];
        i=j=k=0;
        do {
            while (lia[i]<=p[k] && i<nn) i++;
            q[k]=lia[i-1];
            k++;
            while (ria[j]<lia[i] && j<nn) j++;
            if (j<nn) p[k]=ria[j];
            else p[k]=ria[j-1];
        } while (i<nn && j<nn);

        mk=k-1;
        mydata->mk=mk;
        for (i=0; i<nn; i++)
            for (j = 0; j <= mk; j++)
                dij[i*(mk+1)+j]=(li[i]<=q[j] && p[j]<=ri[i]) ? 1.0: 0.0;
        mydata->dij= dij;

        /** Initial values **/

        for (j=0; j<mk; j++) pj[j] = 1.0/mk;
        (mydata->pj)=pj;

        for (j=0; j<=np; j++) xx[j] = 0.0;

        nflag=2;
        fval=9999.0;
        iter2=0;
        do {
            iter2++;
            likval2=fval;

            for (j = 0; j < mk; j++) tau[j]=1.0/(mydata->pj)[j];
            nvold = 0.1;

            /*** calculate expected Xij ***/
            myxij(gam, xx, pj, mydata);

            /*** Primal-dual interior-point method ***/

            nvold = maxp(pj, tau, nvold, xx, gam, mydata);

            /** minimize to get theta **/

            if (nflag!=0) for (j=0; j<=np; j++) xx[j] = 0.1;
            nflag= conmin(np+1, xx, &fval, gg, &ifun, &ifunter, 0.001, 300, ww, 400, 1e-18, 1, fdf_eplike, mydata);

            likval1=fval;

        } while(fabs(likval1-likval2)>0.0001 & iter2 < MAX_ITER);

        now_time = (double) clock ();

        /*** outputs ****/

        for(j=0; j<=np; j++) fprintf(fpo1, "%.10f, ", xx[j]);
        fprintf(fpo1, "%.10f, %.2f \n", fval, now_time/(double) CLOCKS_PER_SEC);
        for(j=0; j<mk; j++) fprintf(fpo2, "%.10f,", (mydata->pj)[j]);
        fprintf(fpo2, "%i\n", mk);
    }

    fclose(fpo1);
    fclose(fpo2);

    free(li);
    free(ri);
    free(zij);

    free(li0);
    free(ri0);
    free(zij0);

    free(lia);
    free(ria);
    free(p);
    free(q);
    free(mypj);

    free(gam);
    free(tau);

    free(xx);
    free(ww);
    free(gg);

    free(mydata);

    return 0;
}

/*** the expected likelihood, 0=alpha, 1=theta and the derivative ***/

double fdf_eplike(int np1, const double *theta, double *df, void * par)
{
    struct fun_params * param = (struct fun_params *) par;
    int i,j,j1, m=(param->mk), nn=(param->nn), tmpi;

    double alp0, sz, likval, msumval, tmpz, tmpz1, tmpa, tmpc, tmpd, tmp0, dfz1, dfz2;
    double *xij, *cij, *pps, pp;

    xij=calloc(m+1, sizeof(double));
    pps=calloc(m+1, sizeof(double));
    cij=calloc(m+1, sizeof(double));

    alp0=theta[0];

    for(j=0;j<m;j++) {
        pps[j]=0.0;
        for(j1=0;j1<j;j1++) pps[j] += ((param->pj)[j1]) ;
    }

    for(j=0;j<np1;j++) df[j]=0.0;

    likval=0.0;
    for(i=0;i<nn;i++) {
        tmpz1=0.0;
        for(j=1;j<np1;j++) tmpz1 += theta[j] * (param->zij)[i*(np1-1)+j-1];
        tmpz=exp(alp0 + tmpz1);

        tmpc=0.0;
        for(j=0;j<m;j++) {
            pp = (param->pj)[j];
            xij[j] = ((param->dij)[i*(m+1)+j]) * exp( -tmpz* pps[j]) *(1.0- exp(-pp*tmpz));
            tmpc += xij[j];
        }
        xij[m]= (param -> dij)[i*(m+1)+m] * exp(-tmpz);
        tmpc += xij[m];

        if (tmpc == 0.0) {
            for(j=0;j<np1;j++) df[j]=-9999999.0;
            return 9999999.0;
        }

        for(j=0;j<=m;j++) {
            cij[j]=0.0;
            for(j1=j+1; j1<m; j1++) cij[j] += xij[j1];
            xij[j] /=tmpc;
        }

        msumval=0.0;
        for(j=0;j<m;j++) {
            cij[j] /=tmpc;
            pp= (param -> pj)[j];
            likval += ( xij[j]* log(1.0- exp(- pp* tmpz )) - cij[j] * pp * tmpz);
            msumval += (( xij[j] *exp(-pp* tmpz) /(1.0- exp(-pp* tmpz)) - cij[j] )*pp);
        }
        df[0] += ((msumval-xij[m])*tmpz);

        for(j=1;j<np1;j++) df[j] += ((msumval-xij[m])*(param->zij)[i*(np1-1)+j-1]*tmpz);
        likval -= xij[m]*tmpz;
    }

    for(j=0;j<np1;j++) df[j] = - df[j]/((double) nn);

    free(cij);
    free(xij);
    free(pps);

    return -likval/((double) nn);
}

/*** Calculate the updated Xij ***/

void myxij(double *lamij, const double *the, const double *mypj, void * parm)
{
    struct fun_params * para = (struct fun_params *) parm;

    double alp0, tmp0, tmpc, pp, tmpz, tmpz1;
    int i,j,j1, tmpi, m=(para->mk), nn=(para->nn), np=(para->np);
    double *gj=calloc(m+1, sizeof(double)), *xij=calloc(m+1, sizeof(double));
    double *dij=(para->dij), *zij=(para->zij);

    alp0=the[0];

    for(j=0;j<m;j++) {
        gj[j]=0.0;
        for(j1=0;j1<j;j1++) gj[j] += mypj[j1];
    }
    for(i=0;i<nn;i++) {
        tmpz1=0.0;
        for(j=0;j<np;j++) tmpz1 += the[j+1] * zij[j+i*np];
        tmpz=exp(alp0+ tmpz1);

        tmpc=0.0;
        for(j=0;j<m;j++) {
            pp = mypj[j];
            xij[j] = dij[i*(m+1)+j] * exp( -tmpz* gj[j]) *(1.0- exp(-pp*tmpz));
            tmpc += xij[j];
        }
        xij[m]= dij[i*(m+1)+m] * exp(-tmpz);
        tmpc += xij[m];

        for(j=0;j<m;j++)  lamij[i*m+j] = xij[j] / tmpc;
    }

    free(xij);
    free(gj);
}

/*** Calculate the first derivative and second derivative wrt p of the expected liklihood ***/
void deveplike(double *df, double *ddf, const double *the, double *gam, double * myp, void *parm)
{
    struct fun_params * para = (struct fun_params *) parm;
    double alp0, tmp1, tmpz1, tmp2, tmpa, pp;
    int i,j,k,j1, m=(para->mk), nn=(para->nn), np=(para->np);
    double *zij=(para->zij);

    alp0=the[0];
    for(j=0;j<m;j++) {

        df[j] = ddf[j] = 0.0;
        pp = myp[j];

        for(i=0;i<nn;i++) {

            tmpz1=0.0;
            for(k=0;k<np;k++) tmpz1 += the[k+1] * zij[k+i*np];
            tmp1=exp(alp0+ tmpz1);
            tmp2=exp(tmp1 * pp);

            tmpa =0.0;
            for(j1=j+1; j1<m; j1++) tmpa += gam[i*m+j1];
            df[j] += tmp1*( gam[i*m+j] /( tmp2- 1.0) - tmpa);
            ddf[j] -= (gam[i*m+j]*tmp1*tmp1*tmp2 / (tmp2 - 1.0) / (tmp2 - 1.0));
        }
    }
}

/*** NormF ***/

double normf(double *p, double *tau, double nv, double eta, const double *the, double *gam, void *parm)
{
    struct fun_params * para = (struct fun_params *) parm;
	double *fj, *ffj, sum1, tmp1, tmp2;
	int j, m=(para->mk);

    fj=calloc(m, sizeof(double));
    ffj=calloc(m, sizeof(double));

    deveplike(fj, ffj, the, gam, p, parm);

    sum1=0.0;
    for(j=0;j<m;j++) {
        tmp1 = (fj[j] + tau[j] - nv);
        tmp2 = p[j] * tau[j] - 1.0/eta;
        sum1 += (tmp1*tmp1 + tmp2*tmp2);
    }

    free(fj);
    free(ffj);
    return sqrt(sum1);
}

/**** Maximize for p ****/

double maxp(double *pj, double *tau, double nvold, const double *the, double *gam, void *parm)
{
    struct fun_params * para = (struct fun_params *) parm;

	int j, iter, m=(para->mk), nn=(para->nn), np=(para->np);
    double eta, nv, curnv, steps, sum1, sum2, tmp1, tmp2;
    double *fj, *ffj, *delp, *curp, *taunew, *curtau;

    fj = calloc(m, sizeof(double));
    ffj = calloc(m, sizeof(double));
    delp = calloc(m, sizeof(double));
    curp= calloc(m, sizeof(double));
    taunew = calloc(m, sizeof(double));
    curtau = calloc(m, sizeof(double));

    iter=0;
    do {
        iter++;

        // calculate the df=fj and diag(ddf)=ffj

        for(j=0;j<m;j++) eta += pj[j]*tau[j];
        eta = m*Sigma/eta;
        deveplike(fj, ffj, the, gam, pj, parm);

        sum1 = sum2 = 0.0;
        for(j=0;j<m;j++) {
            sum1 += ((1.0/eta + pj[j] * fj[j]) / (tau[j] - pj[j] * ffj[j]));
            sum2 += (pj[j] / (tau[j] - pj[j] * ffj[j]));
        }
        nv=sum1/sum2;

        // backtracking: find the largest positive step = steps

        steps=1.0;
        for(j=0;j<m;j++) {
            delp[j] =  (pj[j] * (fj[j]-nv) + 1.0/eta) / (tau[j] - pj[j] * ffj[j]);
            taunew[j] =   nv - ffj[j] * delp[j] - fj[j];
            tmp1 = taunew[j] - tau[j];
            if (tmp1 <0) steps = fmin(steps, -tau[j]/tmp1);
        }
        steps *= 0.99;

        // backtracking: to keep the minimal p > 0;
        tmp2 = 1.0;
        for(j=0;j<m;j++) tmp2=fmin(tmp2,  pj[j] + steps * delp[j]);
        while(tmp2 <= 0) {
            steps *= myPHO;
            tmp2 = 1.0;
            for(j=0;j<m;j++) tmp2=fmin(tmp2, pj[j] + steps * delp[j]);
        }

        // backtracking:

        for(j=0;j<m;j++) {
            curp[j] = pj[j] + steps * delp[j];
            curtau[j] = tau[j] + steps * (taunew[j]-tau[j]);
        }
        curnv=nvold + steps * (nv-nvold);
        tmp2 = normf(curp, curtau, curnv, eta, the, gam, parm);
        tmp1 = normf(pj, tau, nvold, eta, the, gam, parm);

        while(tmp2 > (1.0-myPI*steps)*tmp1) {
            steps *= myPHO;
            for(j=0;j<m;j++) {
                curp[j] = pj[j] + steps * delp[j];
                curtau[j] = tau[j] + steps * (taunew[j]-tau[j]);
            }
            curnv=nvold + steps * (nv-nvold);
            tmp2=normf(curp, curtau, curnv, eta, the, gam, parm);
        }

        deveplike(fj, ffj, the, gam, curp, parm);

        sum1=sum2=0.0;
        for(j=0;j<m;j++) {
            tmp1 = (fj[j] + curtau[j] - curnv);
            sum1 += (tmp1*tmp1);
            sum2 += curp[j] * curtau[j];
            pj[j] = curp[j];
            tau[j] = curtau[j];
        }

        nvold=curnv;
//         printf("iter=%i: sum1=%.20f sum2=%f \n", iter, sqrt(sum1), sum2/m);

    } while ( (sqrt(sum1)> TOL)& iter < MAX_ITER*10 );

    free(fj);
    free(ffj);
    free(delp);
    free(curtau);
    free(taunew);
    free(curp);

    return nvold;
}


text:
/** @file */

#ifndef __CCL_F1D_H_INCLUDED__
#define __CCL_F1D_H_INCLUDED__

#include <gsl/gsl_spline.h>

CCL_BEGIN_DECLS

/*
 * Spline wrapper
 * Used to take care of evaluations outside the supported range
 */
typedef struct {
  gsl_spline *spline;
  double x0,xf; //Interpolation limits
  double y0,yf; //Constant values to use beyond interpolation limit
} ccl_f1d_t;

ccl_f1d_t *ccl_f1d_t_new(int n,double *x,double *y,double y0,double yf);

double ccl_f1d_t_eval(ccl_f1d_t *spl,double x);

void ccl_f1d_t_free(ccl_f1d_t *spl);

CCL_END_DECLS

#endif

ext: h

max_stars_repo_name: benediktdiemer/CCL

max_issues_repo_name: benediktdiemer/CCL

max_forks_repo_name: benediktdiemer/CCL


text:
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include "bstrlib/bstrlib.h"
#include "DosValues.h"

int write_dos_vals(char *outPath, int num_dos, double *Es, double *dos_vals);

int main(int argc, char *argv[]) {
    if (argc < 15) {
        printf("Usage: DosValues.out 'out_path' 'n' 'num_dos' 'Tae' 'Tce' 'Tbe' 'Tao' 'Tco' 'Tbo' 'EpsilonR' 'EpsilonM' 'M' 'W' 'Mu'\n");
        printf("Example: DosValues.out 'dos_vals' '8' '500' '0.2' '1.0' '0.4' '0.08' '0.4' '0.16' '0.05' '0.05' '1.0' '1.0' '-2.0'\n");
        return 2;
    }
    // Parse arguments.
    char *out_path = argv[1];
    int n = atoi(argv[2]);
    int num_dos = atoi(argv[3]);
    double Tae = atof(argv[4]);
    double Tce = atof(argv[5]);
    double Tbe = atof(argv[6]);
    double Tao = atof(argv[7]);
    double Tco = atof(argv[8]);
    double Tbo = atof(argv[9]);
    double EpsilonR = atof(argv[10]);
    double EpsilonM = atof(argv[11]);
    double M = atof(argv[12]);
    double W = atof(argv[13]);
    double Mu = atof(argv[14]);

    // Set up data.
    double *Es = malloc(num_dos * sizeof(double));
    Environment *env = malloc(sizeof(Environment));
    env->Tae = Tae;
    env->Tce = Tce;
    env->Tbe = Tbe;
    env->Tao = Tao;
    env->Tco = Tco;
    env->Tbo = Tbo;
    env->EpsilonM = EpsilonM;
    env->EpsilonR = EpsilonR;
    env->M = M;
    env->W = W;
    env->Mu = Mu;

    // Calculate DOS.
    double *dos_vals = DosValues(env, n, Es, num_dos);

    // Write out DOS.
    int write_err = write_dos_vals(out_path, num_dos, Es, dos_vals);
    if (write_err != 0) {
        printf("Error writing DOS values\n");
    }

    // Cleanup.
    free(dos_vals);
    free(Es);
    return 0;
}

int write_dos_vals(char *outPath, int num_dos, double *Es, double *dos_vals) {
    FILE *fp = fopen(outPath, "w");
    if (fp == NULL) {
        return 1;
    }

    // Write header.
    fprintf(fp, "E\tDOS\n");
    // Write DOS data.
    int i;
    for (i = 0; i < num_dos; i++) {
        fprintf(fp, "%.10f\t%.10f\n", Es[i], dos_vals[i]);
    }

    fclose(fp);
    return 0;
}


text:
#ifndef LAPACKWRAPPER_H
#define LAPACKWRAPPER_H

#include <vector>

extern "C" {
    // void dbdsqr_(char *uplo, int *n, int *ncvt, int *nru, int *ncc, double *D, double *E, double *VT, int *ldvt, double *U, int *ldu, double *C, int *ldc, double *work, int *info);
    // void dgetrf_(int *m, int *n, double *A, int *lda, int *ipiv, int *info);
    // void dgetrs_(char *trans, int *n, int *nrhs, double *A, int *lda, int *ipiv, double *B, int *ldb, int *info);
    // void dgetri_(int *n, double *A, int *lda, int *ipiv, double *work, int *lwork, int *info);
    // void dgehrd_(int *n, int*ilo, int *ihi, double *A, int *lda, double *tau, double *work, int *lwork, int *info);
    // void dhseqr_(char *job, char *compz, int* n, int* ilo, int *ihi, double *H, int *ldh, double *WR, double *WI, double *Z, int *ldz, double *work, int *lwork, int *info);

    // void dorghr_(int *n, int *ilo, int *ihi, double *A, int *lda, double *tau, double *work, int *lwork, int *info);

    // void dtrsyl_(char *trana, char *tranb, int *isgn, int *m, int *n, double *A, int *lda, double *B, int *ldb, double *C, int *ldc, double *scale, int *info);
    #include <lapack.h>
    #include <cblas.h>
}

//---------------------------------------------------------------------------
// LAPACK FUNCTIONS
//---------------------------------------------------------------------------

// Compute the bidiagonal SVD
int cdbdsqr(char uplo, int n, int ncvt, int nru, int ncc,
            double *D, double *E, double *VT, int ldvt,
            double *U, int ldu, double *C, int ldc);

// Compute LU factorization of A
int cdgetrf(int m, int n, double *A, int lda, int *ipiv);
// Solve linear system from LU factorization from cdgetrf
int cdgetrs(char trans, int n, int nrhs, double *A, int lda, int *ipiv, double *B, int ldb);
// Invert matrix given the LU factorization from dgetrf_
int cdgetri(int n, double *A, int lda, int *ipiv);

//-----------------------------------------------------------------
// Computing the Schur Decompotision to apply preconditioner
// --- 1. Reduce to upper Hessenberg using cdgehrd
// --- 2. Generate Orthogonal Matrix using cdorghr
// --- 3. Compute Schur decomposition of Hessenberg Matrix using cdorghr
//----------------------------------------------------------------

// Reduce matrix to upper Hessenberg form
int cdgehrd(int n, int ilo, int ihi, double *A, int lda, double *tau);
// Schur decomposition of a Hessenberg matrix
int cdhseqr(char job, char compz, int n, int ilo, int ihi, double *H, int ldh, double *WR, double *WI, double *Z, int ldz);
// Generate the orthogonal matrix Q implictly created in dgehrd
int cdorghr(int n, int ilo, int ihi, double *A, int lda, double *tau);

// Solve Sylvester System
int cdtrsyl(char trana, char tranb, int isgn, int m, int n, double *A, int lda, double *B, int ldb, double *C, int ldc, double scale);

//-----------------------------------------------------------------------------------
// C++ Class interfaces for LAPACK FUNCTIONS
//-----------------------------------------------------------------------------------

void Print(const std::string str, const std::vector<int> &size, const double *A);
void Transpose(std::vector<int> &size, double *A, double *AT);
// LU decomposition
class LU{
private:
    // Data is stored column-major fashion
    double *data;
    int *ipiv;
public:
    std::vector<int> size; // (number of rows, number of cols)
    LU(const int m_, const int n_, const double *data_);
    ~LU();
    void Factor();
    void SetData(const double *data_);
    void Solve(const int mb, const int nb, double *B, bool trans);
    void GetInverse(const int m, const int n, double *invA);
    void PrintLU();
};

// Schur Decomposition
class Schur{
private:
    double *Qdata, *Tdata;
public:
    std::vector<int> size;
    Schur(const int m_, const int n_, const double *data);
    ~Schur();
    void Factor();
    double* Q() const;
    double* T() const;
    void PrintQ();
    void PrintT();
};

#endif

max_stars_repo_name: mjrodriguez/qwiklapack

max_issues_repo_name: mjrodriguez/qwiklapack

max_forks_repo_name: mjrodriguez/qwiklapack


text:
/* -*- linux-c -*- */
/* cluster.c

   Copyright (C) 2002-2004 John M. Fregeau

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <getopt.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_roots.h>
#include "fewbody.h"
#include "cluster.h"

/* print the usage */
void print_usage(FILE *stream)
{
	fprintf(stream, "USAGE:\n");
	fprintf(stream, "  cluster [options...]\n");
	fprintf(stream, "\n");
	fprintf(stream, "OPTIONS:\n");
	fprintf(stream, "  -n --n <n>                    : set number of stars in cluster [%d]\n", FB_N);
	fprintf(stream, "  -m --m <m/MSUN>               : set mass of each star [%.6g]\n", FB_M/FB_CONST_MSUN);
	fprintf(stream, "  -r --r <r/RSUN>               : set radius of each star [%.6g]\n", FB_R/FB_CONST_RSUN);
	fprintf(stream, "  -S --sigma <sigma/(km/s)>     : set central velocity dispersion [%.6g]\n", FB_SIGMA/1.0e5);
	fprintf(stream, "  -T --rmax <rmax/PARSEC>       : set truncation radius [%.6g]\n", FB_RMAX/FB_CONST_PARSEC);
	fprintf(stream, "  -t --tstop <tstop/t_dyn>      : set stopping time [%.6g]\n", FB_TSTOP);
	fprintf(stream, "  -P --tphysstop <tphysstop/yr> : set physical stopping time [%.6g]\n", FB_TPHYSSTOP/FB_CONST_YR);
	fprintf(stream, "  -D --dt <dt/t_dyn>            : set approximate output dt [%.6g]\n", FB_DT);
	fprintf(stream, "  -c --tcpustop <tcpustop/sec>  : set cpu stopping time [%.6g]\n", FB_TCPUSTOP);
	fprintf(stream, "  -A --absacc <absacc>          : set integrator's absolute accuracy [%.6g]\n", FB_ABSACC);
	fprintf(stream, "  -R --relacc <relacc>          : set integrator's relative accuracy [%.6g]\n", FB_RELACC);
	fprintf(stream, "  -N --ncount <ncount>          : set number of integration steps between calls\n");
	fprintf(stream, "                                  to fb_classify() [%d]\n", FB_NCOUNT);
	fprintf(stream, "  -z --tidaltol <tidaltol>      : set tidal tolerance [%.6g]\n", FB_TIDALTOL);
	fprintf(stream, "  -x --fexp <f_exp>             : set expansion factor of merger product [%.6g]\n", FB_FEXP);
	fprintf(stream, "  -k --ks                       : turn K-S regularization on or off [%d]\n", FB_KS);
	fprintf(stream, "  -s --seed                     : set random seed [%ld]\n", FB_SEED);
	fprintf(stream, "  -d --debug                    : turn on debugging\n");
	fprintf(stream, "  -V --version                  : print version info\n");
	fprintf(stream, "  -h --help                     : display this help text\n");
}

/* calculate the units used (N-body units) */
void calc_units(fb_hier_t hier, fb_units_t *units)
{
	int i, j, k;
	double petot, ketot, r[3];

	/* the unit of mass is defined so that M_tot=1 */
	units->m = 0.0;
	for (i=0; i<hier.nstar; i++) {
		units->m += hier.hier[hier.hi[1]+i].m;
	}

	/* calculate total potential and kinetic energy */
	petot = 0.0;
	for (i=0; i<hier.nstar-1; i++) {
		for (j=i+1; j<hier.nstar; j++) {
			for (k=0; k<3; k++) {
				r[k] = hier.hier[hier.hi[1]+j].x[k] - hier.hier[hier.hi[1]+i].x[k];
			}
			petot += - FB_CONST_G * hier.hier[hier.hi[1]+i].m * hier.hier[hier.hi[1]+j].m / fb_mod(r);
		}
	}

	ketot = 0.0;
	for (i=0; i<hier.nstar; i++) {
		ketot += 0.5 * hier.hier[hier.hi[1]+i].m * fb_dot(hier.hier[hier.hi[1]+i].v, hier.hier[hier.hi[1]+i].v);
	}

	/* the unit of energy is defined so that E=-1/4 */
	units->E = -4.0 * (petot + ketot);

	/* with G=1, all other units are derived */
	units->t = FB_CONST_G * pow(units->m, 2.5) * pow(units->E, -1.5);
	units->l = FB_CONST_G * fb_sqr(units->m) / units->E;
	units->v = units->l / units->t;
}

/* f_0-f(v/v_esc) for the Plummer model */
double fv(double v, void *params)
{
	double f;

	f = ((double *)params)[0];
	return(f - 512.0/(7.0*FB_CONST_PI) * (sqrt(1.0-fb_sqr(v))*(-7.0*v/256.0 + 121.0*fb_cub(v)/384.0 - 263.0*fb_sqr(v)*fb_cub(v)/480.0 + 31.0*v*fb_sqr(fb_cub(v))/80.0 - fb_cub(fb_cub(v))/10.0) + 7.0*asin(v)/256.0));
}

/* get speed from Plummer distribution */
double vf(double f)
{
	int status, iter;
	double v, params[1];
	gsl_function F;
	const gsl_root_fsolver_type *T;
	gsl_root_fsolver *s;

	/* set up the root solver */
	F.function = &fv;
	F.params = &params;

	/* set the parameters */
	params[0] = f;

	T = gsl_root_fsolver_brent;
	s = gsl_root_fsolver_alloc(T);
	gsl_root_fsolver_set(s, &F, 0.0, 1.0);

	/* get v/v_esc by root-finding */
	iter = 0;
	do {
		iter++;
		gsl_root_fsolver_iterate(s);
		status = gsl_root_test_interval(gsl_root_fsolver_x_lower(s), gsl_root_fsolver_x_upper(s), \
						FB_ROOTSOLVER_ABS_ACC, FB_ROOTSOLVER_REL_ACC);
	} while (status == GSL_CONTINUE && iter < FB_ROOTSOLVER_MAX_ITER);

	if (iter >= FB_ROOTSOLVER_MAX_ITER) {
		fprintf(stderr, "Root finder failed to converge.\n");
		exit(1);
	}

	/* we've got the root */
	v = gsl_root_fsolver_root(s);

	/* free memory associated with root solver */
	gsl_root_fsolver_free(s);

	return(v);
}

/* the main attraction */
int main(int argc, char *argv[])
{
	int i, j, n;
	unsigned long int seed;
	double tphysstop, m, r, Ei, Li[3], Lint[3], t, M, a, sigma, rmax, Xmax, radius, v, theta, phi, xcm[3], vcm[3];
	fb_hier_t hier;
	fb_input_t input;
	fb_ret_t retval;
	fb_units_t units;
	char string1[FB_MAX_STRING_LENGTH], string2[FB_MAX_STRING_LENGTH];
	gsl_rng *rng;
	const gsl_rng_type *rng_type=gsl_rng_mt19937;
	const char *short_opts = "n:m:r:S:T:t:P:D:c:A:R:N:z:x:k:s:dVh";
	const struct option long_opts[] = {
		{"n", required_argument, NULL, 'n'},
		{"m", required_argument, NULL, 'm'},
		{"r", required_argument, NULL, 'r'},
		{"sigma", required_argument, NULL, 'S'},
		{"rmax", required_argument, NULL, 'T'},
		{"tstop", required_argument, NULL, 't'},
		{"tphysstop", required_argument, NULL, 'P'},
		{"dt", required_argument, NULL, 'D'},
		{"tcpustop", required_argument, NULL, 'c'},
		{"absacc", required_argument, NULL, 'A'},
		{"relacc", required_argument, NULL, 'R'},
		{"ncount", required_argument, NULL, 'N'},
		{"tidaltol", required_argument, NULL, 'z'},
		{"fexp", required_argument, NULL, 'x'},
		{"ks", required_argument, NULL, 'k'},
		{"seed", required_argument, NULL, 's'},
		{"debug", no_argument, NULL, 'd'},
		{"version", no_argument, NULL, 'V'},
		{"help", no_argument, NULL, 'h'},
		{NULL, 0, NULL, 0}
	};

	/* set parameters to default values */
	n = FB_N;
	m = FB_M;
	r = FB_R;
	sigma = FB_SIGMA;
	rmax = FB_RMAX;
	input.ks = FB_KS;
	input.tstop = FB_TSTOP;
	tphysstop = FB_TPHYSSTOP;
	input.Dflag = 0;
	input.dt = FB_DT;
	input.tcpustop = FB_TCPUSTOP;
	input.absacc = FB_ABSACC;
	input.relacc = FB_RELACC;
	input.ncount = FB_NCOUNT;
	input.tidaltol = FB_TIDALTOL;
	input.fexp = FB_FEXP;
	seed = FB_SEED;
	fb_debug = FB_DEBUG;

	while ((i = getopt_long(argc, argv, short_opts, long_opts, NULL)) != -1) {
		switch (i) {
		case 'n':
			n = atoi(optarg);
			break;
		case 'm':
			m = atof(optarg) * FB_CONST_MSUN;
			break;
		case 'r':
			r = atof(optarg) * FB_CONST_RSUN;
			break;
		case 'S':
			sigma = atof(optarg) * 1.0e5;
			break;
		case 'T':
			rmax = atof(optarg) * FB_CONST_PARSEC;
			break;
		case 't':
			input.tstop = atof(optarg);
			break;
		case 'P':
			tphysstop = atof(optarg) * FB_CONST_YR;
			break;
		case 'D':
			input.Dflag = 1;
			input.dt = atof(optarg);
			break;
		case 'c':
			input.tcpustop = atof(optarg);
			break;
		case 'A':
			input.absacc = atof(optarg);
			break;
		case 'R':
			input.relacc = atof(optarg);
			break;
		case 'N':
			input.ncount = atoi(optarg);
			break;
		case 'z':
			input.tidaltol = atof(optarg);
			break;
		case 'x':
			input.fexp = atof(optarg);
			break;
		case 'k':
			input.ks = atoi(optarg);
			break;
		case 's':
			seed = atol(optarg);
			break;
		case 'd':
			fb_debug = 1;
			break;
		case 'V':
			fb_print_version(stdout);
			return(0);
		case 'h':
			fb_print_version(stdout);
			fprintf(stdout, "\n");
			print_usage(stdout);
			return(0);
		default:
			break;
		}
	}

	/* check to make sure there was nothing crazy on the command line */
	if (optind < argc) {
		print_usage(stdout);
		return(1);
	}

	/* set up parameters for Plummer model */
	M = ((double) n) * m;
	a = FB_CONST_G * M / (6.0*fb_sqr(sigma));
	Xmax = pow(1.0+fb_sqr(a/rmax), -1.5);

	/* initialize a few things for integrator */
	t = 0.0;
	hier.nstarinit = n;
	hier.nstar = n;
	fb_malloc_hier(&hier);
	fb_init_hier(&hier);

	/* put stuff in log entry */
	snprintf(input.firstlogentry, FB_MAX_LOGENTRY_LENGTH, "  command line:");
	for (i=0; i<argc; i++) {
		snprintf(&(input.firstlogentry[strlen(input.firstlogentry)]),
			 FB_MAX_LOGENTRY_LENGTH-strlen(input.firstlogentry), " %s", argv[i]);
	}
	snprintf(&(input.firstlogentry[strlen(input.firstlogentry)]),
		 FB_MAX_LOGENTRY_LENGTH-strlen(input.firstlogentry), "\n");

	/* print out values of paramaters */
	fprintf(stderr, "PARAMETERS:\n");
	fprintf(stderr, "  ks=%d  seed=%ld\n", input.ks, seed);
	fprintf(stderr, "  M=%.6g MSUN  sigma=%.6g km/s  a=%.6g PARSEC  rmax=%.6g PARSEC  rmax/a=%.6g\n",
		M/FB_CONST_MSUN, sigma/1.0e5, a/FB_CONST_PARSEC, rmax/FB_CONST_PARSEC, rmax/a);
	fprintf(stderr, "  n=%d  m=%.6g MSUN  r=%.6g RSUN  tstop=%.6g  tphysstop=%.6g yr  tcpustop=%.6g\n", \
		n, m/FB_CONST_MSUN, r/FB_CONST_RSUN, input.tstop, tphysstop/FB_CONST_YR, input.tcpustop);
	fprintf(stderr, "  tidaltol=%.6g  abs_acc=%.6g  rel_acc=%.6g  ncount=%d  fexp=%.6g\n\n", \
		input.tidaltol, input.absacc, input.relacc, input.ncount, input.fexp);

	/* initialize GSL rng */
	gsl_rng_env_setup();
	rng = gsl_rng_alloc(rng_type);
	gsl_rng_set(rng, seed);

	/* prepare to calculate the center of mass position and velocity */
	for (i=0; i<3; i++) {
		xcm[i] = 0.0;
		vcm[i] = 0.0;
	}

	/* give the objects some properties */
	for (j=0; j<hier.nstar; j++) {
		hier.hier[hier.hi[1]+j].ncoll = 1;
		hier.hier[hier.hi[1]+j].id[0] = j;
		snprintf(hier.hier[hier.hi[1]+j].idstring, FB_MAX_STRING_LENGTH, "%d", j);
		hier.hier[hier.hi[1]+j].n = 1;
		hier.hier[hier.hi[1]+j].obj[0] = NULL;
		hier.hier[hier.hi[1]+j].obj[1] = NULL;
		hier.hier[hier.hi[1]+j].Eint = 0.0;
		hier.hier[hier.hi[1]+j].Lint[0] = 0.0;
		hier.hier[hier.hi[1]+j].Lint[1] = 0.0;
		hier.hier[hier.hi[1]+j].Lint[2] = 0.0;
		hier.hier[hier.hi[1]+j].m = m;
		hier.hier[hier.hi[1]+j].R = r;

		/* draw radius from Plummer distribution */
		radius = a / sqrt(pow(Xmax*gsl_rng_uniform(rng), -2.0/3.0)-1.0);
		theta = acos(2.0 * gsl_rng_uniform(rng) - 1.0);
		phi = 2.0 * FB_CONST_PI * gsl_rng_uniform(rng);
		hier.hier[hier.hi[1]+j].x[0] = radius * sin(theta) * cos(phi);
		hier.hier[hier.hi[1]+j].x[1] = radius * sin(theta) * sin(phi);
		hier.hier[hier.hi[1]+j].x[2] = radius * cos(theta);

		/* draw velocity from Plummer distribution */
		v = vf(gsl_rng_uniform(rng)) * sqrt(2.0*FB_CONST_G*M/(a*sqrt(1.0+fb_sqr(radius/a))));
		theta = acos(2.0 * gsl_rng_uniform(rng) - 1.0);
		phi = 2.0 * FB_CONST_PI * gsl_rng_uniform(rng);
		hier.hier[hier.hi[1]+j].v[0] = v * sin(theta) * cos(phi);
		hier.hier[hier.hi[1]+j].v[1] = v * sin(theta) * sin(phi);
		hier.hier[hier.hi[1]+j].v[2] = v * cos(theta);

		/* calculate center of mass position and velocity */
		for (i=0; i<3; i++) {
			xcm[i] += hier.hier[hier.hi[1]+j].x[i] / ((double) n);
			vcm[i] += hier.hier[hier.hi[1]+j].v[i] / ((double) n);
		}
	}

	/* transform to center of mass frame */
	for (j=0; j<hier.nstar; j++) {
		for (i=0; i<3; i++) {
			hier.hier[hier.hi[1]+j].x[i] -= xcm[i];
			hier.hier[hier.hi[1]+j].v[i] -= vcm[i];
		}
	}

	/* get the units and normalize */
	calc_units(hier, &units);
	fb_normalize(&hier, units);

	/* reset input.tstop if necessary, based on the value of tphysstop */
	if (tphysstop/units.t < input.tstop) {
		input.tstop = tphysstop/units.t;
		fb_dprintf("decreasing tstop in accord with tphysstop: tstop=%.6g\n", input.tstop);
	}

	fb_dprintf("virial ratio: q=%.6g\n", -fb_ketot(&(hier.hier[hier.hi[1]]), hier.nstar)/fb_petot(&(hier.hier[hier.hi[1]]), hier.nstar));

	fprintf(stderr, "UNITS:\n");
	fprintf(stderr, "  v=%.6g km/s  l=%.6g AU  t=t_dyn=%.6g yr\n", \
		units.v/1.0e5, units.l/FB_CONST_AU, units.t/FB_CONST_YR);
	fprintf(stderr, "  M=%.6g M_sun  E=%.6g erg\n\n", units.m/FB_CONST_MSUN, units.E);

	/* calculate the initial energy and angular momentum, for bookkeeping */
	Ei = fb_petot(&(hier.hier[hier.hi[1]]), hier.nstar) + fb_ketot(&(hier.hier[hier.hi[1]]), hier.nstar) +
		fb_einttot(&(hier.hier[hier.hi[1]]), hier.nstar);
	fb_angmom(&(hier.hier[hier.hi[1]]), hier.nstar, Li);
	fb_angmomint(&(hier.hier[hier.hi[1]]), hier.nstar, Lint);
	for (j=0; j<3; j++) {
		Li[j] += Lint[j];
	}

	/* integrate along */
	fb_dprintf("calling fewbody()...\n");

	/* call fewbody! */
	retval = fewbody(input, &hier, &t);

	/* print information to screen */
	fprintf(stderr, "OUTCOME:\n");
	if (retval.retval == 1) {
		fprintf(stderr, "  encounter complete:  t=%.6g (%.6g yr)  %s  (%s)\n\n",
			t, t * units.t/FB_CONST_YR,
			fb_sprint_hier(hier, string1),
			fb_sprint_hier_hr(hier, string2));
	} else {
		fprintf(stderr, "  encounter NOT complete:  t=%.6g (%.6g yr)  %s  (%s)\n\n",
			t, t * units.t/FB_CONST_YR,
			fb_sprint_hier(hier, string1),
			fb_sprint_hier_hr(hier, string2));
	}

	fb_dprintf("there were %ld integration steps\n", retval.count);
	fb_dprintf("fb_classify() was called %ld times\n", retval.iclassify);

	fprintf(stderr, "FINAL:\n");
	fprintf(stderr, "  t_final=%.6g (%.6g yr)  t_cpu=%.6g s\n", \
		t, t*units.t/FB_CONST_YR, retval.tcpu);

	fprintf(stderr, "  L0=%.6g  DeltaL/L0=%.6g  DeltaL=%.6g\n", fb_mod(Li), retval.DeltaLfrac, retval.DeltaL);
	fprintf(stderr, "  E0=%.6g  DeltaE/E0=%.6g  DeltaE=%.6g\n", Ei, retval.DeltaEfrac, retval.DeltaE);
	fprintf(stderr, "  Rmin=%.6g (%.6g RSUN)  Rmin_i=%d  Rmin_j=%d\n", \
		retval.Rmin, retval.Rmin*units.l/FB_CONST_RSUN, retval.Rmin_i, retval.Rmin_j);
	fprintf(stderr, "  Nosc=%d (%s)\n", retval.Nosc, (retval.Nosc>=1?"resonance":"non-resonance"));

	/* free GSL stuff */
	gsl_rng_free(rng);

	/* free our own stuff */
	fb_free_hier(hier);

	/* done! */
	return(0);
}

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos


text:
/*
 *  mixmax.c
 *  MIXMAX, A Pseudo-Random Number Generator
 *
 *  Copyright Konstantin Savvidy.
 *
 *  Free to use, academic or commercial. Do not redistribute without permission.
 *
 *	G.K.Savvidy and N.G.Ter-Arutyunian,
 *  On the Monte Carlo simulation of physical systems,
 *	J.Comput.Phys. 97, 566 (1991);
 *  Preprint EPI-865-16-86, Yerevan, Jan. 1986
 *
 *  K.Savvidy
 *  The MIXMAX random number generator
 *  Comp. Phys. Commun. 196 (2015), pp 161–165
 *  http://dx.doi.org/10.1016/j.cpc.2015.06.003
 *
 */

#include <stdio.h>
#include <stdint.h>

#ifndef MIXMAX_H_
#define MIXMAX_H_

#define USE_INLINE_ASM

#ifdef __cplusplus
extern "C" {
#endif

#ifndef _N
#define N 240
/* The currently recommended generator is the three-parameter MIXMAX with
N=240, s=487013230256099140, m=2^51+1

   Other newly recommended N are N=8, 17 and 240,
   as well as the ordinary MIXMAX with N=256 and s=487013230256099064

   Since the algorithm is linear in N, the cost per number is almost independent of N.
 */
#else
#define N _N
#endif

#ifndef __LP64__
typedef uint64_t myuint;
//#warning but no problem,  'myuint' is 'uint64_t'
#else
typedef unsigned long long int myuint;
//#warning but no problem,  'myuint' is 'unsigned long long int'
#endif

struct rng_state_st
{
    myuint V[N];
    myuint sumtot;
    int counter;
    FILE* fh;
};

typedef struct rng_state_st rng_state_t; // C struct alias

int  rng_get_N(void); // get the N programmatically, useful for checking the value for which the library was compiled

rng_state_t  *rng_alloc();                 /* allocate the state */
int           rng_free(rng_state_t* X);    /* free memory occupied by the state */
rng_state_t  *rng_copy(myuint *Y);         /* init from vector, takes the vector Y,
                                               returns pointer to the newly allocated and initialized state */
void read_state(rng_state_t* X, const char filename[] );
void print_state(rng_state_t* X);
    int iterate(rng_state_t* X);
    myuint iterate_raw_vec(myuint* Y, myuint sumtotOld);

//   FUNCTIONS FOR SEEDING

typedef uint32_t myID_t;

void seed_uniquestream(rng_state_t* X, myID_t clusterID, myID_t machineID, myID_t runID, myID_t  streamID );
/*
 best choice: will make a state vector from which you can get at least 10^100 numbers
 guaranteed mathematically to be non-colliding with any other stream prepared from another set of 32bit IDs,
 so long as it is different by at least one bit in at least one of the four IDs
			-- useful if you are running a parallel simulation with many clusters, many CPUs each
 */

void seed_spbox(rng_state_t* X, myuint seed);    // non-linear method, makes certified unique vectors,  probability for streams to collide is < 1/10^4600

void seed_vielbein(rng_state_t* X, unsigned int i); // seeds with the i-th unit vector, i = 0..N-1,  for testing only

//   FUNCTIONS FOR GETTING RANDOM NUMBERS

#ifdef __MIXMAX_C
	myuint get_next(rng_state_t* X);         // returns 64-bit int, which is between 1 and 2^61-1 inclusive
	double get_next_float(rng_state_t* X);   // returns double precision floating point number in (0,1]
#endif  //__MIXMAX_C

void fill_array(rng_state_t* X, unsigned int n, double *array); // fastest method: set n to a multiple of N (e.g. n=256)

void iterate_and_fill_array(rng_state_t* X, double *array); // fills the array with N numbers

myuint precalc(rng_state_t* X);
/* needed if the state has been changed by something other than  iterate, but no worries, seeding functions call this for you when necessary */
myuint apply_bigskip(myuint* Vout, myuint* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t  streamID );
// applies a skip of some number of steps calculated from the four IDs
void branch_inplace( rng_state_t* Xin, myID_t* ID ); // almost the same as apply_bigskip, but in-place and from a vector of IDs

#define BITS  61

/* magic with Mersenne Numbers */

#define M61   2305843009213693951ULL

    myuint modadd(myuint foo, myuint bar);
    myuint modmulM61(myuint s, myuint a);
    myuint fmodmulM61(myuint cum, myuint s, myuint a);

#define MERSBASE M61 //xSUFF(M61)
#define MOD_PAYNE(k) ((((k)) & MERSBASE) + (((k)) >> BITS) )  // slightly faster than my old way, ok for addition
#define MOD_REM(k) ((k) % MERSBASE )  // latest Intel CPU is supposed to do this in one CPU cycle, but on my machines it seems to be 20% slower than the best tricks
#define MOD_MERSENNE(k) MOD_PAYNE(k)

//#define INV_MERSBASE (0x1p-61)
#define INV_MERSBASE (0.4336808689942017736029811203479766845703E-18)
//const double INV_MERSBASE=(0.4336808689942017736029811203479766845703E-18); // gives "duplicate symbol" error

// the charpoly is irreducible for the combinations of N and SPECIAL

#if (N==256)
#define SPECIALMUL 0
#define SPECIAL 487013230256099064 // s=487013230256099064, m=1 -- good old MIXMAX
#define MOD_MULSPEC(k) fmodmulM61( 0, SPECIAL , (k) )

#elif (N==8)
#define SPECIALMUL 53 // m=2^53+1
#define SPECIAL 0

#elif (N==17)
#define SPECIALMUL 36 // m=2^36+1, other valid possibilities are m=2^13+1, m=2^19+1, m=2^24+1
#define SPECIAL 0

#elif (N==40)
#define SPECIALMUL 42 // m=2^42+1
#define SPECIAL 0

#elif (N==60)
#define SPECIALMUL 52 // m=2^52+1
#define SPECIAL 0

#elif (N==96)
#define SPECIALMUL 55 // m=2^55+1
#define SPECIAL 0

#elif (N==120)
#define SPECIALMUL 51   // m=2^51+1 and a SPECIAL=+1 (!!!)
#define SPECIAL 1
#define MOD_MULSPEC(k) (k)

#elif (N==240)
#define SPECIALMUL 51   // m=2^51+1 and a SPECIAL=487013230256099140
#define SPECIAL 487013230256099140ULL
#define MOD_MULSPEC(k) fmodmulM61( 0, SPECIAL , (k) )

#elif (N==44851)
#define SPECIALMUL 0
#define SPECIAL -3
#define MOD_MULSPEC(k) MOD_MERSENNE(3*(MERSBASE-(k)))

#else
#warning Not a verified N, you are on your own!
#define SPECIALMUL 58
#define SPECIAL 0

#endif // list of interesting N for modulus M61 ends here

        inline 	myuint get_next_inlined(rng_state_t* X) { // forward output
            int i;
            i=X->counter;

            if (i<=(N-1) ){
                X->counter++;
                return X->V[i];
            }else{
                X->sumtot = iterate_raw_vec(X->V, X->sumtot);
                X->counter=2;
                return X->V[1];
            }
        }

//    inline 	myuint get_next_inlined(rng_state_t* X) { // backward output, slightly faster
//        int i;
//        i=X->counter;
//
//        if ( i!=0 ){
//            X->counter--;
//            return X->V[i];
//        }else{
//            X->sumtot = iterate_raw_vec(X->V, X->sumtot);
//            X->counter=N-2;
//            return X->V[N-1];
//        }
//    }

    inline double get_next_float_inlined(rng_state_t* X){
        /* cast to signed int trick suggested by Andrzej Görlich     */
        int64_t Z=(int64_t)get_next_inlined(X);
        double F;
#if defined(__GNUC__) && (__GNUC__ < 5) && (!defined(__ICC)) && defined(__x86_64__) && defined(__SSE2_MATH__) && defined(USE_INLINE_ASM)
#warning Using the inline assembler
        /* using SSE inline assemly to zero the xmm register, just before int64 -> double conversion,
         not really necessary in GCC-5 or better, but huge penalty on earlier compilers
         */
        __asm__  __volatile__("pxor %0, %0; "
                              :"=x"(F)
                              );
#endif
        F=Z;
        return F*INV_MERSBASE;
    }

#ifndef __MIXMAX_C // c++ can put code into header files, why cant we? (with the inline declaration, should be safe from duplicate-symbol error)

#define get_next(X) get_next_inlined(X)
#define get_next_float(X) get_next_float_inlined(X) // get_next_float_packbits(X) //

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wstrict-aliasing"
    inline double convert1double(uint64_t u){
        const double one = 1;
        const uint64_t onemask = *(uint64_t*)&one;
        uint64_t tmp = (u>>9) | onemask; // must not have bits higher than 52 !
        double d = *(double*)&tmp;
        return d-1.0;
    }

    inline double get_next_float_packbits(rng_state_t* X){
        uint64_t Z=get_next_inlined(X); // must not have bits higher than 52 !
        return convert1double(Z);
    }
#pragma GCC diagnostic pop

#endif // __MIXMAX_C

// ERROR CODES - exit() is called with these
#define ARRAY_INDEX_OUT_OF_BOUNDS   0xFF01
#define SEED_WAS_ZERO               0xFF02
#define ERROR_READING_STATE_FILE    0xFF03
#define ERROR_READING_STATE_COUNTER       0xFF04
#define ERROR_READING_STATE_CHECKSUM      0xFF05

#ifdef __cplusplus
}
#endif

//#define HOOKUP_GSL 1

#ifdef HOOKUP_GSL // if you need to use mixmax through GSL, pass -DHOOKUP_GSL=1 to the compiler
#ifndef __MIXMAX_C

#include <gsl/gsl_rng.h>
unsigned long gsl_get_next(void *vstate);
double gsl_get_next_float(void *vstate);
void seed_for_gsl(void *vstate, unsigned long seed);

static const gsl_rng_type mixmax_type =
{"MIXMAX",                        /* name */
    MERSBASE,                         /* RAND_MAX */
    1,                                /* RAND_MIN */
    sizeof (rng_state_t),
    &seed_for_gsl,
    &gsl_get_next,
    &gsl_get_next_float
};

unsigned long gsl_get_next(void *vstate) {
    rng_state_t* X= (rng_state_t*)vstate;
    return (unsigned long)get_next(X);
}

double gsl_get_next_float(void *vstate) {
    rng_state_t* X= (rng_state_t*)vstate;
    return ( (double)get_next(X)) * INV_MERSBASE;
}

void seed_for_gsl(void *vstate, unsigned long seed){
    rng_state_t* X= (rng_state_t*)vstate;
    seed_spbox(X,(myuint)seed);
}

const gsl_rng_type *gsl_rng_mixmax = &mixmax_type;

#endif // HOOKUP_GSL
#endif // not inside __MIXMAX_C
#endif // closing MIXMAX_H_

max_stars_repo_name: ruicamposcolabpt/MontyCarlo

max_issues_repo_name: ruicamposcolabpt/MontyCarlo

max_forks_repo_name: ruicamposcolabpt/MontyCarlo


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

/* solve_quadratic.c - finds the real roots of a x^2 + b x + c = 0 */

#include <config.h>
#include <math.h>

#include <gsl/gsl_poly.h>

int
gsl_poly_solve_quadratic (double a, double b, double c,
                          double *x0, double *x1)
{
  double disc = b * b - 4 * a * c;

  if (a == 0) /* Handle linear case */
    {
      if (b == 0)
        {
          return 0;
        }
      else
        {
          *x0 = -c / b;
          return 1;
        };
    }

  if (disc > 0)
    {
      if (b == 0)
        {
          double r = fabs (0.5 * sqrt (disc) / a);
          *x0 = -r;
          *x1 =  r;
        }
      else
        {
          double sgnb = (b > 0 ? 1 : -1);
          double temp = -0.5 * (b + sgnb * sqrt (disc));
          double r1 = temp / a ;
          double r2 = c / temp ;

          if (r1 < r2)
            {
              *x0 = r1 ;
              *x1 = r2 ;
            }
          else
            {
              *x0 = r2 ;
              *x1 = r1 ;
            }
        }
      return 2;
    }
  else if (disc == 0)
    {
      *x0 = -0.5 * b / a ;
      *x1 = -0.5 * b / a ;
      return 2 ;
    }
  else
    {
      return 0;
    }
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>

#include "su2_matrix.c"

#ifndef gauge_config
#include "config.h"
#define gauge_config 1

#define sq(x) ((x)*(x))
#define cb(x) ((x)*(x)*(x))
#endif

int acc, tot; // acceptance rate

#define FORWARD 1 //true
#define BACKWARD 0 //false

// Access pattern methods --------------------------------------------------

// index of link site=(t,x,y,z) in direction dir
int index(const int* site, int dir){
  return dir+site[0]*4+site[1]*N*4+site[2]*N*N*4+site[3]*N*N*N*4;
}

// shift site array by one to direction dir
void shift(int* site, int dir, int forward){
  if(forward == 1){
    site[dir] = (site[dir]+1)%N;
  }else{
    site[dir] = (site[dir]-1+N)%N;
  }
}

// end of access pattern methods -------------------------------------------

// Plaquette methods -------------------------------------------------------

// plaquette operator at site (t,x,y,z) in plane (mu,nu)
// normally mu>nu
double plaquette(const su2_matrix* links, int mu, int nu, int* site){
  su2_matrix A,B,C,D;
  A = links[index(site,mu)];
  shift(site,mu,FORWARD);
  B = links[index(site,nu)];
  shift(site,mu,BACKWARD);
  shift(site,nu,FORWARD);
  C = su2_inv(links[index(site,mu)]);
  shift(site,nu,BACKWARD); // restore site position
  D = su2_inv(links[index(site,nu)]);

  A = su2_mul(A,B);
  A = su2_mul(A,C);
  A = su2_mul(A,D);

  return su2_trace(A)/2.0;
}

double plaquette_mean(const su2_matrix* links, int mu, int nu){
  double mean = 0;
  int site[4];
  for(int z=0;z<N;z++){
    for(int y=0;y<N;y++){
      for(int x=0;x<N;x++){
        for(int t=0;t<N;t++){
            site[0]=t; site[1]=x; site[2]=y; site[3]=z;
            mean += plaquette(links,mu,nu,site);
        }
      }
    }
  }
  return mean / (N*N*N*N);
}

// end of plaquette methods ------------------------------------------------

// Metropolis update --------------------------------------------------------

void do_hb_update(su2_matrix* links, const gsl_rng* r, const su2_matrix* rands, int* site, int dir);
double action(su2_matrix* links, int* site, int dir);

// the big loop
void hb_update(su2_matrix* links, const gsl_rng* r, const su2_matrix* rands){
  int site[4];
  for(int z=0;z<N;z++){
    for(int y=0;y<N;y++){
      for(int x=0;x<N;x++){
        for(int t=0;t<N;t++){
          for(int dir=0;dir<4;dir++){
            site[0]=t; site[1]=x; site[2]=y; site[3]=z;
            do_hb_update(links,r,rands,site,dir);
          }
        }
      }
    }
  }
}

// update
void do_hb_update(su2_matrix* links, const gsl_rng* r, const su2_matrix* rands, int* site, int dir){
  if(verbose) {printf("Updating {%d, %d, %d, %d} direction %d",site[0],site[1],site[2],site[3],dir);getchar();}
  double dS,action_old;
  su2_matrix link_old;
  for(int times=0; times<Nk; times++){
    action_old = action(links,site,dir);
    link_old = links[index(site,dir)];
    links[index(site,dir)]=su2_mul(links[index(site,dir)],rands[gsl_rng_uniform_int(r, 50)]);

    dS = action(links,site,dir)-action_old;

    if(verbose) {printf("Old action = %f, ",action_old); su2_print(link_old); getchar();}
    if(verbose) {printf("New action = %f, ",action_old+dS); su2_print(links[index(site,dir)]); getchar();}
    acc++;tot++;
    if(dS>0){
      if(exp(-dS)<gsl_rng_uniform(r)){
        acc--;
        links[index(site,dir)]=link_old; // revert update
        if(verbose) {printf("Reverted"); getchar();}
      }
    }
  }

}

// calculating the action
// for each direction we need to find plaquette at 2 sites in each of the 3 planes
// mu > nu
double action(su2_matrix* links, int* site, int dir){
  int mu,nu;
  double result=0;
  for(int k=0; k<4; k++){
    if(k==dir) continue; // loop over planes
    if(k>dir){mu=k; nu=dir;} // mu=larger number
    if(k<dir){mu=dir; nu=k;}

    result += plaquette(links, mu,nu,site);
    shift(site,k,BACKWARD); // shift to the perpendicular direction backwards
    result += plaquette(links, mu,nu,site);
    shift(site,k,FORWARD);  // shift it back
  }

  return -beta*result;
}

// end of heat bath update -------------------------------------------------

// Wilson loop -------------------------------------------------------------

double wloop(int r, int t, int mu,const su2_matrix* links, int* site){
  su2_matrix A[(r+t)*2];

  // bottom side
  for(int i=0;i<t;i++){
    A[i]=links[index(site,0)];
    shift(site,0,FORWARD);
  }

  // right side
  for(int i=0;i<r;i++){
    A[t+i]=links[index(site,mu)];
    shift(site,mu,FORWARD);
  }

  // top side
  for(int i=0;i<t;i++){
    shift(site,0,BACKWARD);
    A[t+r+i]=su2_inv(links[index(site,0)]);
  }

  // left side
  for(int i=0;i<r;i++){
    shift(site,mu,BACKWARD);
    A[t+r+t+i]=su2_inv(links[index(site,mu)]);
  }

  // multiply
  for(int i=1;i<2*(r+t);i++){
    A[0]=su2_mul(A[0],A[i]);
  }

  return su2_trace(A[0])/2.0;
}

double wloop_mean(int wR, int wT, const su2_matrix* links){
  double mean = 0;
  int site[4];
  for(int z=0;z<N;z++){
    for(int y=0;y<N;y++){
      for(int x=0;x<N;x++){
        for(int t=0;t<N;t++){
            site[0]=t; site[1]=x; site[2]=y; site[3]=z;
            for(int mu=1;mu<4;mu++){
              mean += wloop(wR,wT,mu,links,site);
            }
        }
      }
    }
  }
  return mean/ (3*N*N*N*N); // 3 directions & N^3 sites
}

// end of Wilson loop ------------------------------------------------------

max_stars_repo_name: uchuutamashi/lqcd

max_issues_repo_name: uchuutamashi/lqcd

max_forks_repo_name: uchuutamashi/lqcd


text:
/**
 *
 * @file qwrapper_dpotrf.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:56 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlauum(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int n, int nb,
                       double *A, int lda)
{
    DAG_CORE_LAUUM;
    QUARK_Insert_Task(quark, CORE_dlauum_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &n,     VALUE,
        sizeof(double)*nb*nb,    A,             INOUT,
        sizeof(int),                        &lda,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlauum_quark = PCORE_dlauum_quark
#define CORE_dlauum_quark PCORE_dlauum_quark
#endif
void CORE_dlauum_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int N;
    double *A;
    int LDA;

    quark_unpack_args_4(quark, uplo, N, A, LDA);
    LAPACKE_dlauum_work(LAPACK_COL_MAJOR, lapack_const(uplo), N, A, LDA);
}


text:
#include <gsl/gsl_complex.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_statistics_float.h>

#include <stdio.h>

int vector_float_convolve(int cs, const float* c, int as, const float* a, int rs, float* r)
{
  int h = cs / 2;
  int li, ri;
  int i,j;

  for (i = 0; i < cs; i++) {
    li = i - h;
    ri = i + h;
    r[i] = 0;
    for (j = (li >= 0 ? li : 0); j < (ri < as ? ri : (as - 1)); j++) {
      r[i] += a[j]*c[j+h+1];
    }
  }
  return 0;
}

int vector_double_convolve(int cs, const double* c, int as, const double* a, int rs, double* r)
{
  int h = cs / 2;
  int li, ri;
  int i,j;

  for (i = 0; i < cs; i++) {
    li = i - h;
    ri = i + h;
    r[i] = 0;
    for (j = (li >= 0 ? li : 0); j < (ri < as ? ri : (as - 1)); j++) {
      r[i] += a[j]*c[j+h+1];
    }
  }
  return 0;
}

int vector_complex_convolve(int cs, const gsl_complex* c, int as, const gsl_complex* a, int rs, gsl_complex* r)
{
  int h = cs / 2;
  int li, ri;
  int i,j;

  for (i = 0; i < cs; i++) {
    li = i - h;
    ri = i + h;
    r[i].dat[0] = 0;
    r[i].dat[1] = 0;
    for (j = (li >= 0 ? li : 0); j < (ri < as ? ri : (as - 1)); j++) {
      r[i].dat[0] += a[j].dat[0]*c[j+h+1].dat[0]-a[j].dat[1]*c[j+h+1].dat[1];
      r[i].dat[1] += a[j].dat[0]*c[j+h+1].dat[1]+a[j].dat[1]*c[j+h+1].dat[0];
    }
  }
  return 0;
}

int filter_double(int ls, const double* l, int ks, const double* k, int vs, const double* v, int rs, double* r)
{
  if (ls > vs || ks > vs) return 2000; // BAD_SIZE

  int i,j;

  double L = l[0];
  double K = k[0];

  int N = ls - 1;
  int M = ks - 1;

  for (i = 0; i < vs; i++) {
    r[i] = 0;
    for (j = 0; j < N; j++) {
      if (i - j > 0) r[i] -= (l[j+1])*v[i-j];
    }
    for (j = 0; j < M; j++) {
      if (i - j > 0) r[i] += (k[j+1])*r[i-j];
    }
  }
  return 0;
}

int filter_float(int ls, const float* l, int ks, const float* k, int vs, const float* v, int rs, float* r)
{
  if (ls > vs || ks > vs) return 2000; // BAD_SIZE

  int i,j;

  float L = l[0];
  float K = k[0];

  int N = ls - 1;
  int M = ks - 1;

  for (i = 0; i < vs; i++) {
    r[i] = 0;
    for (j = 0; j < N; j++) {
      if (i - j > 0) r[i] -= (l[j+1])*v[i-j];
    }
    for (j = 0; j < M; j++) {
      if (i - j > 0) r[i] += (k[j+1])*r[i-j];
    }
  }
  return 0;
}

int hilbert(int rs, gsl_complex* r)
{
  int s = rs;

  gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s);
  gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s);

  // forward fourier transform
  gsl_fft_complex_forward ((double*)r, 1, s, wavetable, workspace);
  // zero negative coefficients and double positive

  int i;
  int m = s/2;
  for (i = 1; i < s; i++) {
    if (i <= m) {
      r[i].dat[0] *= 2;
      r[i].dat[1] *= 2;
    }
    else if (s % 2 == 0 && i == m+1) {
    }
    else {
      r[i].dat[0] = 0;
      r[i].dat[1] = 0;
    }
  }

  // inverse fourier transform
  gsl_fft_complex_inverse ((double*)r, 1, s, wavetable, workspace);

  gsl_fft_complex_wavetable_free (wavetable);
  gsl_fft_complex_workspace_free (workspace);

  return 0;
}

int pwelch(int w, int vs, const gsl_complex* v, int rs, double* r)
{
  if (w > vs) return 2000; // BAD_SIZE

  int i,j;

  int fs = w;

  int num_windows = vs / fs; // ignore end

  double s[fs];
  for (i = 0; i < fs; i++) s[i] = 0;

  gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (fs);
  gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (fs);

  gsl_complex* f = malloc(sizeof(gsl_complex)*fs);
  gsl_vector_view F = gsl_vector_view_array((double*)f, 2*fs);

  gsl_vector_view X;

  for (i = 0; i < num_windows; i++) {
    X = gsl_vector_view_array((double*)(&v[i*fs]), 2*fs); // v is gsl_complex*
    gsl_blas_dcopy(&X.vector,&F.vector);
    gsl_fft_complex_forward ((double*)f, 1, fs, wavetable, workspace);
    for (j = 0; j < fs; j++) s[j] += f[j].dat[0]*f[j].dat[0] + f[j].dat[1]*f[j].dat[1];
  }
  for (j = 0; j < rs; j++) {
    if (j == 0) r[j] = s[j];
    else if (j == (rs-1)) r[j] = s[j];
    else r[j] = s[j] + s[fs-j+1];

    r[j] /= num_windows;
    r[j] = sqrt(r[j]);
  }
  gsl_fft_complex_wavetable_free (wavetable);
  gsl_fft_complex_workspace_free (workspace);

  free(f);

  return 0;
}

int hamming_double(int rs, double* r)
{
  int i;

  for (i = 0; i < rs; i++) r[i] = 0.54 - 0.46 * cos(2*M_PI*i/rs);

  return 0;
}

int hamming_float(int rs, float* r)
{
  int i;

  for (i = 0; i < rs; i++) r[i] = 0.54 - 0.46 * cos(2*M_PI*i/rs);

  return 0;
}

int real_poly_complex_eval(int cs, const double* c, int zs, const gsl_complex* z, int rs, gsl_complex* r)
{
  int i;

  for (i = 0; i < zs; i++)
    r[i] = gsl_poly_complex_eval(c,cs,z[i]);

  return 0;
}

int complex_power_double(int cs, const gsl_complex* c, int rs, double* r)
{
  if (rs != cs) return 2000; // BAD_SIZE

  int i;

  for (i = 0; i < cs; i++)
    r[i] = c[i].dat[0]*c[i].dat[0] + c[i].dat[1]*c[i].dat[1];

  return 0;
}

int complex_power_float(int cs, const gsl_complex* c, int rs, float* r)
{
  if (rs != cs) return 2000; // BAD_SIZE

  int i;

  for (i = 0; i < cs; i++)
    r[i] = c[i].dat[0]*c[i].dat[0] + c[i].dat[1]*c[i].dat[1];

  return 0;
}

int downsample_double(int n, int xs, const double* x, int rs, double* r)
{
  if (rs != xs/n) return 2000; // BAD_SIZE

  int i;

  for (i = 0; i < rs; i++)
    r[i] = x[i*n];

  return 0;
}

int downsample_float(int n, int xs, const float* x, int rs, float* r)
{
  if (rs != xs/n) return 2000; // BAD_SIZE

  int i;

  for (i = 0; i < rs; i++)
    r[i] = x[i*n];

  return 0;
}

int vector_diff_double(int xs, const double* x, int rs, double* r)
{
  if (rs != xs - 1) return 2000; // BAD_SIZE

  int i;

  for (i = 0; i < rs; i++)
    r[i] = x[i+1] - x[i];

  return 0;
}

int vector_diff_float(int xs, const float* x, int rs, float* r)
{
  if (rs != xs - 1) return 2000; // BAD_SIZE

  int i;

  for (i = 0; i < rs; i++)
    r[i] = x[i+1] - x[i];

  return 0;
}

int unwrap_double(int xs, const double* x, int rs, double* r)
{
  if (rs != xs) return 2000; // BAD_SIZE

  int i;

  r[0] = x[0];

  double c = 0;

  int tmp;

  for (i = 1; i < rs; i++) {
    tmp = x[i-1] - x[i];
    if (tmp > M_PI) {
      r[i] = 2*M_PI;
    }
    else if (tmp < (-M_PI)) {
      r[i] = -2*M_PI;
    }
    else {
      r[i] = 0;
    }
    c += r[i];
    r[i] = c + x[i];
  }

  return 0;
}

int unwrap_float(int xs, const float* x, int rs, float* r)
{
  if (rs != xs) return 2000; // BAD_SIZE

  int i;

  r[0] = x[0];

  double c = 0;

  int tmp;

  for (i = 1; i < rs; i++) {
    tmp = x[i-1] - x[i];
    if (tmp > M_PI) {
      r[i] = 2*M_PI;
    }
    else if (tmp < (-M_PI)) {
      r[i] = -2*M_PI;
    }
    else {
      r[i] = 0;
    }
    c += r[i];
    r[i] = c + x[i];
  }

  return 0;
}

int cross_covariance_double(int max_lag,
			    double* sx,
			    double* sy,
			    int xs, const double* x,
			    int ys, const double* y,
			    int rs, double* r)
{
  if (xs != ys)        return 2000; // BAD_SIZE
  if (rs != 2*max_lag) return 2000; // BAD_SIZE

  double mx = gsl_stats_mean(x,1,xs);
  double my = gsl_stats_mean(y,1,ys);

  *sx = gsl_stats_sd(x,1,xs);
  *sy = gsl_stats_sd(y,1,ys);

  int delay;
  int i, j;

  double sxy;

  for (delay=-max_lag; delay < max_lag; delay++) {
    sxy = 0;
    for (i=0;i<xs;i++) {
      j = i+delay;
      if (j < 0 || j >= xs) sxy += (x[i]-mx)*(-my);
      else sxy += (x[i]-mx)*(y[j]-my);
      /* or should it be:
      if (j < 0 || j >= xs) continue;
      else sxy += (x[i]-mx)*(y[j]-my);
      */
    }
    r[delay+max_lag] = sxy/xs;
  }

  return 0;
}

int cross_covariance_float(int max_lag,
			    float* sx,
			    float* sy,
			    int xs, const float* x,
			    int ys, const float* y,
			    int rs, double* r)
{
  if (xs != ys)        return 2000; // BAD_SIZE
  if (rs != 2*max_lag) return 2000; // BAD_SIZE

  float mx = gsl_stats_float_mean(x,1,xs);
  float my = gsl_stats_float_mean(y,1,ys);

  *sx = gsl_stats_float_sd(x,1,xs);
  *sy = gsl_stats_float_sd(y,1,ys);

  int delay;
  int i, j;

  float sxy;

  for (delay=-max_lag; delay < max_lag; delay++) {
    sxy = 0;
    for (i=0;i<xs;i++) {
      j = i+delay;
      if (j < 0 || j >= xs) sxy += (x[i]-mx)*(-my);
      else sxy += (x[i]-mx)*(y[j]-my);
      /* or should it be:
      if (j < 0 || j >= xs) continue;
      else sxy += (x[i]-mx)*(y[j]-my);
      */
    }
    r[delay+max_lag] = sxy/xs;
  }

  return 0;
}

int cum_sum_double(int xs, const double* x, int rs, double* r)
{
  int i;

  r[0] = x[0];

  for (i=1;i<xs;i++) {
    r[i] = r[i-1]+x[i];
  }

  return 0;
}

int cum_sum_float(int xs, const float* x, int rs, float* r)
{
  int i;

  r[0] = x[0];

  for (i=1;i<xs;i++) {
    r[i] = r[i-1]+x[i];
  }

  return 0;
}

max_stars_repo_name: amcphail/hsignal

max_issues_repo_name: amcphail/hsignal

max_forks_repo_name: amcphail/hsignal


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C code headers for the implementation of the Fourier domain response for LISA-like detectors
 *
 */

#ifndef _LISAFDRESPONSE_H
#define _LISAFDRESPONSE_H

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "waveform.h"
#include "EOBNRv2HMROMstruct.h"
#include "LISAgeometry.h"

#if defined(__cplusplus)
extern "C" {
#elif 0
} /* so that editors will match preceding brace */
#endif

/**************************************************/
/**************** Prototypes **********************/

/* Core function processing a signal (in the form of a list of modes) through the Fourier-domain LISA response, for given values of the inclination, position in the sky and polarization angle - here simplified version for just the y_21 observable */
int LISASimFDResponse21(
  LISAconstellation *variant,                                 /* Provides specifics on the variant of LISA */
  struct tagListmodesCAmpPhaseFrequencySeries **list,  /* Input/Output: list of modes in Frequency-domain amplitude and phase form as produced by the ROM, and output after FD response processing */
  const double inclination,                                   /* Inclination of the source */
  const double lambda,                                        /* First angle for the position in the sky */
  const double beta,                                          /* Second angle for the position in the sky */
  const double psi);                                          /* Polarization angle */

//WARNING: tRef is ignored for now in the response - i.e. set to 0
/* Core function processing a signal (in the form of a list of modes) through the Fourier-domain LISA response, for given values of the inclination, position in the sky and polarization angle */
int LISASimFDResponsey12(
  LISAconstellation *variant,                                 /* Provides specifics on the variant of LISA */
  struct tagListmodesCAmpPhaseFrequencySeries **list,      /* Input: list of modes in Frequency-domain amplitude and phase form as produced by the ROM */
  struct tagListmodesCAmpPhaseFrequencySeries **listy12,   /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the y12 observable */
  const double torb,                                       /* Reference orbital time - tf as read from the hlm gives t-tinj, this arg allows to pass tinj to the response */
  const double lambda,                                        /* First angle for the position in the sky */
  const double beta,                                          /* Second angle for the position in the sky */
  const double inclination,                                   /* Inclination of the source */
  const double psi,                                           /* Polarization angle */
  const int tagfrozenLISA,                                    /* Tag to treat LISA as frozen at its torb configuration  */
  const ResponseApproxtag responseapprox);                    /* Tag to select possible low-f approximation level in FD response */

//WARNING: tRef is ignored for now in the response - i.e. set to 0
/* Core function processing a signal (in the form of a list of modes) through the Fourier-domain LISA response, for given values of the inclination, position in the sky and polarization angle */
int LISASimFDResponseTDI3Chan(
  int tagtRefatLISA,                                          /* 0 to measure Tref from SSB arrival, 1 at LISA guiding center */
  LISAconstellation *variant,                                 /* Provides specifics on the variant of LISA */
  struct tagListmodesCAmpPhaseFrequencySeries **list,      /* Input: list of modes in Frequency-domain amplitude and phase form as produced by the ROM */
  struct tagListmodesCAmpPhaseFrequencySeries **listTDI1,  /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the TDI channel 1 */
  struct tagListmodesCAmpPhaseFrequencySeries **listTDI2,  /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the TDI channel 2 */
  struct tagListmodesCAmpPhaseFrequencySeries **listTDI3,  /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the TDI channel 3 */
  const double torb,                                       /* Reference orbital time - tf as read from the hlm gives t-tinj, this arg allows to pass tinj to the response */
  const double lambda,                                        /* First angle for the position in the sky */
  const double beta,                                          /* Second angle for the position in the sky */
  const double inclination,                                   /* Inclination of the source */
  const double psi,                                           /* Polarization angle */
  const double m1,                                            /* m1 in solar masses - used for resampling */
  const double m2,                                            /* m2 in solar masses - used for resampling */
  const double maxf,                                          /* Maximal frequency to consider - used to ignore hard-to-resolve response at f>1Hz - NOTE: for now, no recomputation of the boundary, so when not resampling can lose a bit of support between the last frequency point covered and maxf */
  const TDItag tditag,                                        /* Selector for the set of TDI observables */
  const int tagfrozenLISA,                                    /* Tag to treat LISA as frozen at its torb configuration  */
  const ResponseApproxtag responseapprox);                    /* Tag to select possible low-f approximation level in FD response */

// int LISASimFDResponseTDI1Chan(
//   struct tagListmodesCAmpPhaseFrequencySeries **list,      /* Input: list of modes in Frequency-domain amplitude and phase form as produced by the ROM */
//   struct tagListmodesCAmpPhaseFrequencySeries **listTDI,   /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the TDI channel 1 */
//   const double tRef,                                          /* Time at coalescence */
//   const double lambda,                                        /* First angle for the position in the sky */
//   const double beta,                                          /* Second angle for the position in the sky */
//   const double inclination,                                   /* Inclination of the source */
//   const double psi,                                           /* Polarization angle */
//   const double maxf,                                          /* Maximal frequency to consider - used to ignore hard-to-resolve response at f>1Hz - NOTE: for now, no recomputation of the boundary, so when not resampling can lose a bit of support between the last frequency point covered and maxf */
//   const TDItag tditag);                                       /* Selector for the set of TDI observables */

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _LISAFDRESPONSE_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare