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_line2d.c
 *  \brief run dnest for 2d line analysis.
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.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_line2d;

/*!
 * this function does dnest sampling.
 */
int dnest_line2d(int argc, char **argv)
{
  int i;

  set_blr_model2d();

  num_params_blr = num_params_blr_model + num_params_nlr
                 + num_params_res + num_params_linecenter + 1; /* include line sys err */
  num_params_blr_tot = num_params_blr;
  num_params = parset.n_con_recon + num_params_blr_tot + num_params_var;
  idx_resp = num_params_blr_tot + num_params_drw + num_params_trend;
  idx_difftrend = idx_resp + num_params_resp;
  idx_linecenter = num_params_blr_model + num_params_nlr + num_params_res;

  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_line2d = dnest_malloc_fptrset();
  /* setup functions used for dnest*/
  fptrset_line2d->from_prior = from_prior_line2d;
  fptrset_line2d->print_particle = print_particle_line2d;
  fptrset_line2d->restart_action = restart_action_2d;
  fptrset_line2d->accept_action = accept_action_2d;
  fptrset_line2d->kill_action = kill_action_2d;
  fptrset_line2d->perturb = perturb_line2d;
  fptrset_line2d->read_particle = read_particle_line2d;

  if(parset.flag_exam_prior != 1)
  {
    fptrset_line2d->log_likelihoods_cal_initial = log_likelihoods_cal_initial_line2d;
    fptrset_line2d->log_likelihoods_cal_restart = log_likelihoods_cal_restart_line2d;
    fptrset_line2d->log_likelihoods_cal = log_likelihoods_cal_line2d;
  }
  else
  {
    fptrset_line2d->log_likelihoods_cal_initial = log_likelihoods_cal_line2d_exam;
    fptrset_line2d->log_likelihoods_cal_restart = log_likelihoods_cal_line2d_exam;
    fptrset_line2d->log_likelihoods_cal = log_likelihoods_cal_line2d_exam;
  }

  set_par_range_model2d();
  /* setup the fixed parameters */
  set_par_fix_blrmodel();

  /* setup the remaining paramters */
  for(i=num_params_blr_model; i<num_params; i++)
  {
    par_fix[i] = 0;
    par_fix_val[i] = -DBL_MAX;
  }

  /* cope with narrow line */
  if(parset.flag_narrowline == 2)
  {
    /* flux */
    if(parset.flux_narrowline_err == 0.0)
    {
      par_fix[num_params_blr_model] = 1.0;
      par_fix_val[num_params_blr_model] = 0.0;
    }
  }
  if(parset.flag_narrowline >= 2)
  {
    /* width */
    if(parset.width_narrowline_err == 0.0)
    {
      par_fix[num_params_blr_model+1] = 1.0;
      par_fix_val[num_params_blr_model+1] = 0.0;
    }
    /* shift */
    if(parset.shift_narrowline_err == 0.0)
    {
      par_fix[num_params_blr_model+2] = 1.0;
      par_fix_val[num_params_blr_model+2] = 0.0;
    }
  }

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

  /* 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 systematic error of continuum */
  if(parset.flag_con_sys_err != 1)
  {
    par_fix[num_params_blr] = 1;
    par_fix_val[num_params_blr] = log(1.0);
  }

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

  print_par_names_model2d();

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

  dnest_free_fptrset(fptrset_line2d);
  return 0;
}

/*!
 * this function setups parameter ranges.
 *
 * The order of parameters is:                              \n
 *   I.   blr model.............()                          \n
 *   II.  narrow line...........(if flag_narrowline is true)\n
 *   III. spectral broadening...()                          \n
 *   IV.  line center...........(if flag_linecenter is true)\n
 *   V.   systematic error......()                          \n
 *   VI.  variability...........()                          \n
 *   VII. long-term trend.......()                          \n
 *   VIII.response A and Ag.....()                          \n
 *   IX.  different trend.......(if flag_difftend is ture)  \n
 *   X.   continuum light curve.()                          \n
 */
void set_par_range_model2d()
{
  int i, j, i1, i2;

  /* setup parameter range, 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;
  }

  /* cope with narrow line */
  for(i=num_params_blr_model; i<num_params_blr_model + num_params_nlr; i++)
  {
    par_range_model[i][0] = nlr_range_model[i - num_params_blr_model][0];
    par_range_model[i][1] = nlr_range_model[i - num_params_blr_model][1];

    par_prior_model[i] = nlr_prior_model[i - num_params_blr_model];
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0; /* note that for logarithm prior of flux, this value is not used, so does not matter */
  }
  /* cope with spectral broadening */
  for(i=num_params_blr_model + num_params_nlr; i<num_params_blr_model + num_params_nlr + num_params_res; i++)
  {
    par_range_model[i][0] = -10.0;
    par_range_model[i][1] =  10.0;

    par_prior_model[i] = GAUSSIAN;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0;
  }
  /* cope with line center */
  for(i=num_params_blr-num_params_linecenter-1; i< num_params_blr-1; i++)
  {
    par_range_model[i][0] = -10.0;
    par_range_model[i][1] =  10.0;

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

  /* the last is systematic 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;

  /* variability parameters */
   /* first systematic error */
  i = num_params_blr;
  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;

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

      /* 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][0]);
      par_range_model[i][1] = fmin(par_range_model[i][1], var_range_model[i-num_params_blr][1]);

      par_prior_model[i] = GAUSSIAN;
      par_prior_gaussian[i][0] = var_param[i-num_params_blr];
      par_prior_gaussian[i][1] = var_param_std[i-num_params_blr];
    }
    else
    {
      par_range_model[i][0] = var_range_model[i-num_params_blr][0];
      par_range_model[i][1] = var_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;
    }
  }
  /* long-term trend of continuum */
  for(i=num_params_drw + num_params_blr; i< num_params_drw + num_params_trend + num_params_blr; 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 in continuum and line */
  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 ligth curve values */
  for(i=num_params_blr+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_model2d()
{
  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_model2d.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]);
  }

  for(j=0; j<num_params_nlr; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "narrow 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_res; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "line broaden", 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_linecenter; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "line center", 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]);

  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_line2d(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_blr + num_params_var; 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_line2d(const void *model)
{
  double logL;
  logL = prob_initial_line2d(model);
  return logL;
}

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

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

/*!
 * this function prints out parameters.
 */
void print_particle_line2d(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_line2d(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_line2d(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+num_params_var)/num_params))
      which = dnest_rand_int(num_params_blr + num_params_var);
    else
      which = dnest_rand_int(parset.n_con_recon) + num_params_blr + num_params_var;
  }while(par_fix[which] == 1);

  which_parameter_update = which;

  /* level-dependent width */
  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
  {
    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]);
    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]);
  }

  return logH;
}

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

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

  // the parameter previously updated
  param = which_parameter_update;
  // continuum parameter is updated
  if(param >= num_params_blr)
  {
    /*
     *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)
  {
    /* BLR parameter is updated
     * Note a) that the (num_par_blr-1)-th parameter is systematic error of line.
     * when this parameter is updated, Trans2D and Fline are unchanged.
     *      b) Fline is always changed, except for 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 = Trans2D_at_veldata_particles[which_particle_update];
      Trans2D_at_veldata_particles[which_particle_update] = Trans2D_at_veldata_particles_perturb[which_particle_update];
      Trans2D_at_veldata_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;

    }
  }

  return;
}

/*
 * action when particle i is killed in cdnest sampling.
 * particle i_copy's properties is copyed to particle i.
 */
void kill_action_2d(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 * n_vel_data_ext  * sizeof(double));
  memcpy(TransTau_particles[i], TransTau_particles[i_copy], parset.n_tau*sizeof(double));
  memcpy(Trans2D_at_veldata_particles[i], Trans2D_at_veldata_particles[i_copy], parset.n_tau * n_vel_data_ext * sizeof(double));
  return;
}

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS


text:
/*
 * Copyright 2016 Maikel Nadolski
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0

 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef HMM_IO_H_
#define HMM_IO_H_

#include <map>
#include <istream>
#include <Eigen/Dense>
#include <gsl_assert.h>
#include <range/v3/all.hpp>
#include <iostream>

#include "maikel/hmm/hidden_markov_model.h"
#include "maikel/function_profiler.h"

namespace maikel { namespace hmm {

  struct getline_error: public std::runtime_error {
      getline_error(std::string s): std::runtime_error(s) {}
  };

  struct read_ascii_matrix_error: public std::runtime_error {
      read_ascii_matrix_error(std::string s): std::runtime_error(s) {}
  };

  struct read_sequence_error: public std::runtime_error {
      read_sequence_error(std::string s): std::runtime_error(s) {}
  };

  inline std::istream& getline(std::istream& in, std::istringstream& linestream)
  {
    Expects(in);
    std::string line;
    if (!std::getline(in, line))
      throw getline_error("Could not read the line from given stream.");
    linestream.str(line);
    linestream.clear();
    return in;
  }

  template <class size_type>
    std::pair<size_type, size_type> getdims(std::istream& in)
    {
      Expects(in);
      std::pair<size_type, size_type> dim;
      std::istringstream line(ranges::front(ranges::getlines(in)));
      if (!(line >> dim.first >> dim.second))
        throw read_ascii_matrix_error("Could not read dimensions.");
      return dim;
    }

  template <class float_type>
    typename std::enable_if<
        std::is_floating_point<float_type>::value,
    Eigen::Matrix<float_type, Eigen::Dynamic, Eigen::Dynamic>>::type
    read_ascii_matrix(std::istream& in, std::size_t rows, std::size_t cols)
    {
      Expects(in);
      typename Eigen::Matrix<float_type, Eigen::Dynamic, Eigen::Dynamic> matrix(rows, cols);
      std::istringstream line;
      for (std::size_t i = 0; i < rows; ++i) {
        getline(in, line);
        for (std::size_t j = 0; j < cols; ++j)
          if (!(line >> matrix(i,j)))
            throw read_ascii_matrix_error("Could not read entries in line: " + line.str() + ".");
      }
      Ensures(gsl::narrow<std::size_t>(matrix.rows()) == rows &&
              gsl::narrow<std::size_t>(matrix.cols()) == cols);
      return matrix;
    }

  template <class Derived>
    Eigen::DenseBase<Derived>& normalize_rows(Eigen::DenseBase<Derived>& matrix)
    {
      using Index = typename Eigen::DenseBase<Derived>::Index;
      for (Index i = 0; i < matrix.rows(); ++i)
        matrix.row(i) /= matrix.row(i).sum();
      return matrix;
    }

  template <class float_type>
    typename std::enable_if<
        std::is_floating_point<float_type>::value,
    hidden_markov_model<float_type>>::type
    read_hidden_markov_model(std::istream& in)
    {
      Expects(in);
      using matrix     = typename hidden_markov_model<float_type>::matrix;
      using row_vector = typename hidden_markov_model<float_type>::row_vector;
      using index      = typename hidden_markov_model<float_type>::size_type;
      index states;
      index symbols;

      std::tie(states, symbols) = getdims<index>(in);
      matrix A  = read_ascii_matrix<float_type>(in, states, states);
      matrix B  = read_ascii_matrix<float_type>(in, states, symbols);
      row_vector pi = read_ascii_matrix<float_type>(in, 1, states);
      normalize_rows(A);
      normalize_rows(B);
      normalize_rows(pi);

      Ensures(A.rows() == states && A.cols() == states);
      Ensures(B.rows() == states && B.cols() == symbols);
      Ensures(pi.rows() == 1 && pi.cols() == states);
      return hidden_markov_model<float_type>(A, B, pi);
    }

  template <class float_type>
    typename std::enable_if<
        std::is_floating_point<float_type>::value,
    void>::type
    print_model_parameters(std::ostream& out, hidden_markov_model<float_type> const& model)
    {
      out << "N= " << model.states() << "\n";
      out << "M= " << model.symbols() << "\n";
      out << "A:\n" << model.A << "\n";
      out << "B:\n" << model.B << "\n";
      out << "pi:\n" << model.pi << "\n";
      out << std::flush;
    }

  template <class size_type>
    size_type read_sequence_length(std::istream& in)
    {
      std::string line;
      std::getline(in, line);
      std::istringstream linestream(line);
      size_type length;
      linestream >> length;
      return length;
    }

  template <class Integral>
    std::map<std::string, Integral>
    read_symbol_map(std::istream& in)
    {
      std::map<std::string, Integral> symbol_to_index;
      Integral count{};
      std::string buffer;
      std::getline(in, buffer);
      std::istringstream lstream(buffer);
      while (lstream >> buffer)
        if (symbol_to_index.insert(make_pair(buffer, count)).second)
          ++count;
      return symbol_to_index;
    }

  template <class Integral>
    std::vector<Integral>
    read_sequence(std::istream& in)
    {
      std::vector<Integral> sequence;
      std::map<std::string, Integral> symbol_to_index = read_symbol_map<Integral>(in);
      sequence.reserve(read_sequence_length<std::size_t>(in));
      auto symbol_map = [&symbol_to_index] (std::string const& symbol) {
          auto found = symbol_to_index.find(symbol);
          if (found == symbol_to_index.end())
            throw read_sequence_error("Unkown Symbols in Input.");
          return found->second;
      };
      auto sequence_input = ranges::istream_range<std::string>(in);
      ranges::copy(sequence_input | ranges::view::transform(symbol_map), ranges::back_inserter(sequence));
      return sequence;
    }

  template <class Integral, class Symbol>
    std::vector<Integral>
    read_sequence(std::istream& in, std::map<Symbol,Integral>& symbol_to_index)
    {
      MAIKEL_PROFILER;
      std::vector<Integral> sequence;
      sequence.reserve(read_sequence_length<std::size_t>(in));
      auto symbol_map = [&symbol_to_index] (Symbol const& symbol) {
          auto found = symbol_to_index.find(symbol);
          if (found == symbol_to_index.end())
            throw read_sequence_error("Unkown Symbols in Input.");
          return found->second;
      };
      auto sequence_input = ranges::istream_range<Symbol>(in);
      ranges::copy(sequence_input | ranges::view::transform(symbol_map), ranges::back_inserter(sequence));
      return sequence;
    }
} // namespace hmm
} // namespace maikel

#endif /* HMM_IO_H_ */

ext: h

max_stars_repo_name: maikel/hidden-markov-model

max_issues_repo_name: maikel/Hidden-Markov-Model

max_forks_repo_name: maikel/Hidden-Markov-Model


text:
#include <math.h>
#include <gsl/gsl_errno.h>
#include "bicycle.h"
//#include <stdio.h>
//#include <gsl/gsl_matrix.h>
//#include <gsl/gsl_odeiv.h>

void eval_qdots(double qd[], const double q[], const double u[], const Bicycle_t *bike)
{
  qd[0] = (-sin(q[2])*u[0] + cos(q[2])*u[2]) / cos(q[1]);
  qd[1] = cos(q[2])*u[0] + sin(q[2])*u[2];
  qd[2] = sin(q[2])*tan(q[1])*u[0] + u[1] - cos(q[2])*tan(q[1])*u[2];
  qd[3] = u[3] - u[1];
  qd[4] = u[4] - u[2];
  qd[5] = sin(q[4])*u[0] - cos(q[4])*u[1] + u[5];
  qd[6] = -bike->rr*cos(q[0])*(qd[2] + qd[3]) - bike->rrt*(sin(q[0])*qd[1] +
          cos(q[0])*cos(q[1])*(qd[2] + qd[3]));
  qd[7] = -bike->rr*sin(q[0])*(qd[2] + qd[3]) + bike->rrt*(cos(q[0])*qd[1] +
          sin(q[0])*cos(q[1])*(qd[2] + qd[3]));
} // eval_qdots

void Bicycle_to_Bicycle_benchmark(Bicycle_benchmark_t *bike_b, const Bicycle_t *bike)
{
} // Bicycle_to_Bicycle_benchmark

void Bicycle_benchmark_to_Bicycle(Bicycle_t *bike, const Bicycle_benchmark_t *bike_b)
{
} // Bicycle_benchmark_to_Bicycle

void eval_dependent_speeds(double u[], double B[3][6], const int ind[], const int dep[])
{
  double Bi[3][3], Bd[3][3], Bd_adj[3][3], Bd_det, ui[3], ud[3];

  // Independent speeds
  ui[0] = u[ind[0]];
  ui[1] = u[ind[1]];
  ui[2] = u[ind[2]];

  // Dependent speeds
  ud[0] = u[dep[0]];
  ud[1] = u[dep[1]];
  ud[2] = u[dep[2]];

  // Form -Bi
  Bi[0][0] = -B[0][ind[0]];
  Bi[0][1] = -B[0][ind[1]];
  Bi[0][2] = -B[0][ind[2]];

  Bi[1][0] = -B[1][ind[0]];
  Bi[1][1] = -B[1][ind[1]];
  Bi[1][2] = -B[1][ind[2]];

  Bi[2][0] = -B[2][ind[0]];
  Bi[2][1] = -B[2][ind[1]];
  Bi[2][2] = -B[2][ind[2]];

  // Form Bd
  Bd[0][0] = B[0][dep[0]];
  Bd[0][1] = B[0][dep[1]];
  Bd[0][2] = B[0][dep[2]];

  Bd[1][0] = B[1][dep[0]];
  Bd[1][1] = B[1][dep[1]];
  Bd[1][2] = B[1][dep[2]];

  Bd[2][0] = B[2][dep[0]];
  Bd[2][1] = B[2][dep[1]];
  Bd[2][2] = B[2][dep[2]];

  // Forming Bd_adj and Bd_det
  madj3(Bd_adj, Bd);
  mdet3(&Bd_det, Bd);

  // Bd_adj * (-Bi), store the result in Bd instead of creating a new matrix to
  // store the result
  mm3(Bd, Bd_adj, Bi);
  mv3(ud, Bd, ui);

  // Form the dependent speeds
  u[dep[0]] = ud[0] / Bd_det;
  u[dep[1]] = ud[1] / Bd_det;
  u[dep[2]] = ud[2] / Bd_det;
} // eval_dependent_speeds

/*
 * 3x3 Matrix - Matrix multiplication
 */
void mm3(double AB[3][3], double A[][3], double B[][3])
{
  unsigned short i, j, k;
  double sum;
  for (i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) {
      sum = 0.0;
      for (k = 0; k < 3; k++) {
        sum += A[i][k] * B[k][j];
      } // for k
      AB[i][j] = sum;
    } // for j
  } // for i
} // mm3

/*
 * 3x3 Matrix - 3x1 Vector multiplication
 */
void mv3(double Ab[3], double A[3][3], double b[3])
{
  unsigned short i, j;
  double sum;
  for (i = 0; i < 3; i++) {
    sum = 0.0;
    for (j = 0; j < 3; j++) {
      sum += A[i][j]*b[j];
    } // for j
    Ab[i] = sum;
  } // for i
} // mv 3

/*
 * 3x3 Matrix adjugate
 */
void madj3(double adj[3][3], double m[3][3])
{
  adj[0][0] = m[1][1]*m[2][2] - m[1][2]*m[2][1];
  adj[0][1] = m[0][2]*m[2][1] - m[0][1]*m[2][2];
  adj[0][2] = m[0][1]*m[1][2] - m[0][2]*m[1][1];

  adj[1][0] = m[1][2]*m[2][0] - m[1][0]*m[2][2];
  adj[1][1] = m[0][0]*m[2][2] - m[0][2]*m[2][0];
  adj[1][2] = m[0][2]*m[1][0] - m[0][0]*m[1][2];

  adj[2][0] = m[1][0]*m[2][1] - m[1][1]*m[2][0];
  adj[2][1] = m[0][1]*m[2][0] - m[0][0]*m[2][1];
  adj[2][2] = m[0][0]*m[1][1] - m[0][1]*m[1][0];
} // madj3

/*
 * Determinant of a 3x3 matrix
 */
void mdet3(double *det, double m[3][3])
{
  *det = m[0][0]*(m[1][1]*m[2][2] - m[1][2]*m[2][1]) +
         m[0][1]*(m[1][2]*m[2][0] - m[1][0]*m[2][2]) +
         m[0][2]*(m[1][0]*m[2][1] - m[1][1]*m[2][0]);
} // mdet3

void form_constraint_matrix(double B[3][6], const double q[], const Bicycle_t *bike)
{
  double den = sqrt(pow(sin(q[1])*sin(q[4]) - cos(q[1])*sin(q[2])*cos(q[4]), 2)
               + pow(cos(q[1])*cos(q[3]), 2)),
         g31 = (sin(q[1])*sin(q[4]) - cos(q[1])*sin(q[2])*cos(q[5])) / den,
         g33 = cos(q[1])*cos(q[2]) / den;

  B[0][0] = cos(q[4])*sin(q[4])*(g33*bike->rf + cos(q[1])*cos(q[2])*bike->rft);
  B[0][1] = bike->ls +
            pow(sin(q[4]), 2)*(g33*bike->rf + cos(q[1])*cos(q[2])*bike->rft);
  B[0][2] = bike->rrt*sin(q[2]);
  B[0][3] = - cos(q[2])*(bike->rr + cos(q[1])*bike->rrt);
  B[0][4] = - (sin(q[1])*bike->rft + sin(q[4])*(g31*bike->rf + bike->lf));
  B[0][5] = cos(q[4])*(g33*bike->rf + cos(q[1])*cos(q[2])*bike->rft);

  B[1][0] = - bike->ls + cos(q[2])*(bike->rr + cos(q[1])*bike->rrt)
            - pow(cos(q[4]), 2)*(g33*bike->rf + cos(q[1])*cos(q[2])*bike->rft);
  B[1][1] = -cos(q[4])*sin(q[4])*(g33*bike->rf + cos(q[1])*cos(q[2])*bike->rft);
  B[1][2] = bike->lr + sin(q[2])*(bike->rr + cos(q[1])*bike->rrt);
  B[1][3] = 0.0;
  B[1][4] = cos(q[4])*(bike->lf + g31*bike->rf) + cos(q[1])*sin(q[2])*bike->rft;
  B[1][5] = 0.0;

  B[2][0] = 0.0;
  B[2][1] = 0.0;
  B[2][2] = 0.0;
  B[2][3] = 0.0;
  B[2][4] = 0.0;
  B[2][5] = 0.0;
} // form_constraint_matrix

ext: c

max_stars_repo_name: certik/pydy

max_issues_repo_name: certik/pydy

max_forks_repo_name: certik/pydy


text:
/*
 *
 * ML fitting of size (sigma, s) and number of photons (n)
 * of already localized dots using a 2D Gaussian model provided by gaussianInt2.
 *
 * Only one dot at a time is considered, hence the 1 in mlfit1sn
 *
 * Uses 3D images however only a xy-plane at at time is used for the
 * fitting.
 *
 * See also:
 *   mlfit1.c, fitting of xy
 *   df_mlfit1sn.c, a MATLAB interface
 *
 * 2017.04.12. Valgrind ok.
 *
 * Future work:
 * - Return the fitting status
 * - Study the failed cases
 * - More restrictions on the search domain, i.e., cap sigma depending
 *   on the window size

 ==9245== Invalid write of size 8
==9245==    at 0x1094C9: mlfit1sn_dot (mlfit1sn.c:227)
==9245==    by 0x109A9B: mlfit1sn (mlfit1sn.c:278)
==9245==    by 0x109A9B: unit_tests (mlfit1sn.c:346)
==9245==    by 0x57F73F0: (below main) (libc-start.c:291)
==9245==  Address 0x5babf98 is 24 bytes after a block of size 32,000 in arena "client"
==9245==

*/

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

#include <gsl/gsl_multimin.h>
#include <gsl/gsl_statistics.h>

#include "gaussianInt2.h"
#include "mlfit.h"
#include "blit3.h"

//
// Headers
//

// Number of photons in signal
double estimateNphot(double *, size_t);
// Locate multiple dots in a volume
int    mlfit1sn(double * , size_t , size_t , size_t ,
                            double * , size_t , double * );
// Localization for single dot
int    mlfit1sn_dot(double *, size_t, double *,  double *);
// Cost function
double my_f (const gsl_vector *, void *);
// Random number in a range
double rand_range(double, double);
int    unit_tests(void);

// 0: no information, 1: per dot, 2: even more
#ifndef verbose
#define verbose 0
#endif

// When low on bugs,
// replace gsl_vector_get(v,i) by v->data[i*v->stride]
// #define GSL_RANGE_CHECK_OFF

// Globals
#if verbose>1
uint32_t maxInterations = 5;
#else
uint32_t maxInterations = 5000;
#endif

double convCriteria = 1e-6;
// In theory also the window size, denoted Ws in most places

// Optimization constants
typedef struct {
  double * R; // local region
  double * G; // A already allocated temporary space for the Gaussian kernel
  size_t Rw; // R and G has size Rw*Rw
  double sigma; // For constant sigma fitting
  double bg;
  double x;
  double y;
  double midvalue; // value of central pixel in model - bg
} optParams;

double my_f (const gsl_vector *v, void *params)
  // The function to optimize
{
  optParams *p = (optParams *) params;

  size_t Rw = p->Rw;
  double * R = p->R;
  double * GI = p->G;
  double bg = p->bg;
  double x = p->x;
  double y = p->y;

  // Get the other parameters ...
  double sigma = gsl_vector_get(v, 0);
  double Nphot = gsl_vector_get(v, 1);

#if verbose > 1
  printf("Rw: %lu\n", Rw);
  printf("Nphot: %f\n", Nphot);
  printf("sigma: %f\n", sigma);
  printf("bg   : %f\n", bg);
  printf("x: %f, y: %f\n", x, y);
  printf("GI:\n");
  showRegion(GI, Rw);
#endif

  if(sigma<0)
    return(INFINITY);

  if(Nphot<0)
    return(INFINITY);

  // Create Gaussian ...
  double mu[] = {x,y};
  gaussianInt2(GI, mu, &sigma, Rw);
  p->midvalue = Nphot*GI[ (Rw*Rw-1)/2 ];
  for(size_t kk = 0; kk<Rw*Rw; kk++)
    GI[kk] = bg+Nphot*GI[kk];

  /*
   * from LL2PG.m
   *   model = x(3)+x(4)*gaussianInt2([x(1), x(2)], x(5), (size(patch, 1)-1)/2);
   *   mask = disk2d((size(patch,1)-1)/2);
   *   %L = -sum(sum(-(patch-model).^2./model - .5*log(model)));
   *   L = -sum(sum( mask.*( -(patch-model).^2./model - .5*log(model) ) ));
   */

  double E = 0;
  for (size_t kk=0; kk<Rw*Rw; kk++)
  {
    E+= (GI[kk]-R[kk])*(GI[kk]-R[kk])/GI[kk] - .5*log(GI[kk]); // ML
    //E+= (GI[kk]-R[kk])*(GI[kk]-R[kk]); // Quadratic
  }
  //  E = -E;

  #if verbose > 1
  printf("E: %f\n", E);
#endif
  return E;
}

int mlfit1sn_dot(double * V, size_t Vm,
    double * D,  double * F)
  // Localization for a dot roughly centered in V of size Vm x Vm
  // D[0], D[1], D[3] are the global coordinates of the dot
  // F are the fitted coordinates
{
  // Non-optimized parameters
  optParams par;
  par.R = V;
  par.Rw = Vm;
  par.x = D[0]-nearbyint(D[0]);
  par.y = D[1]-nearbyint(D[1]);
  // par.z required for accurate nphot counting, sigma is invariant to
  // this shift
  par.bg = estimateBG(V, Vm);
  par.G = malloc(Vm*Vm*sizeof(double));

#if verbose>0
  printf("localizeDot\n");
  printf("x: %f y: %f\n", par.x, par.y);
  printf("bg: %f\n", par.bg);
#endif

  const gsl_multimin_fminimizer_type *T =
    gsl_multimin_fminimizer_nmsimplex2;
  gsl_multimin_fminimizer *s = NULL;
  gsl_vector *ss, *x;
  gsl_multimin_function minex_func;

  size_t iter = 0;
  int status;
  double size;

  /* Starting point */
  x = gsl_vector_alloc (2);
  gsl_vector_set(x, 0, 1.2); // sigma
  double nphot0 =estimateNphot(V, Vm);
  gsl_vector_set(x, 1, nphot0); // Nphot

  /* Set initial step sizes */
  ss = gsl_vector_alloc(2);
  gsl_vector_set(ss, 0, 0.1); // sigma
  gsl_vector_set(ss, 1, nphot0/100); // Number of photons

  /* Initialize method and iterate */
  minex_func.n = 2;
  minex_func.f = my_f;
  minex_func.params = ∥

  s = gsl_multimin_fminimizer_alloc (T, 2);
  gsl_multimin_fminimizer_set (s, &minex_func, x, ss);

  do
  {
    iter++;
    status = gsl_multimin_fminimizer_iterate(s);

    if (status)
      break;

    size = gsl_multimin_fminimizer_size(s);
    status = gsl_multimin_test_size(size, convCriteria);

    if (status == GSL_SUCCESS)
    {
#if verbose > 0
      printf ("converged to minimum at\n");
      printf ("%5lu sigma:%10.3e NP:%10.3e f() = %7.3f size = %10.3e\n",
          iter,
          gsl_vector_get (s->x, 0),
          gsl_vector_get (s->x, 1),
          s->fval, size);
#endif
    }
  }
  while (status == GSL_CONTINUE && iter < maxInterations);

  F[0] = gsl_vector_get(s->x, 0);
  F[1] = par.midvalue;
  F[2] = s->fval;
  F[3] = status;
  F[4] = par.bg;

  gsl_vector_free(x);
  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free(s);

  free(par.G);

  return status;
}

int mlfit1sn(double * V, size_t Vm, size_t Vn, size_t Vp,
    double * D, size_t Dm, double * F)
  // run the localization routine for a list of dots, D [3xDm]
{

  // Return values have to be checked
  gsl_set_error_handler_off();
  // as a custom error handler can also be used.

  int Ws_pref = 15; // Window size is Ws x Ws
  int Ws = Ws_pref;
  double * W = malloc(Ws_pref*Ws_pref*sizeof(double));

  for(size_t kk =0; kk<Dm; kk++)
  {
    // Copy the neighbourhood around each D into W
    int hasRegion = 0;

    for(Ws = Ws_pref; Ws>7; Ws = Ws-2) {
      if(getRegion(W, Ws, V, Vm, Vn, Vp, D+kk*3) == 0) {
        hasRegion = 1;
        break;
      } else {
#if verbose >0
        printf("Ws: %d not possible\n", Ws);
#endif
      }
    }
#if verbose>0
    printf("Ws: %d\n", Ws);
#endif

    if(hasRegion == 1) {
#if verbose > 1
      printf("W:\n");
      showRegion(W, Ws);
#endif
      // Local fitting in W
#if verbose > 0
      int status = mlfit1sn_dot(W,  Ws, D+kk*3, F+kk*5);
      printf("Status: %d\n", status);
#else
     mlfit1sn_dot(W,  Ws, D+kk*3, F+kk*5);
#endif
    }
    else
    {
      F[kk*5] = -1;
      F[kk*5+1] = -1;
      F[kk*5+2] = -1;
      F[kk*5+3] =  -3;
      F[kk*5+4] = 0;
    }
  }
  free(W);
  return 0;
}

#ifdef standalone
int unit_tests()
{
  double * V; // image
  int Vm = 1024; int Vn = 1024; int Vp = 60;
  double * D; // list of dots
  size_t Dm = 1000; // number of dots
  double * F; // fitted dots

  printf("Image size: %dx%dx%d\n", Vm, Vn, Vp);
  printf("Localizing %lu dots\n", Dm);

  V = malloc(Vm*Vn*Vp*sizeof(double));
  D = malloc(3*Dm*sizeof(double));
  F = malloc(4*Dm*sizeof(double));

  // Initialize the data
  for(int kk=0; kk<Vm*Vn*Vp; kk++)
    V[kk] = rand_range(0,0);

  for(uint32_t kk=0; kk<Dm; kk++)
  {
    size_t pos = kk*3;
    D[pos] =   rand_range(-1, Vm+1);
    D[pos+1] = rand_range(6, Vn-7);
    D[pos+2] = rand_range(6, Vp-7);
    if(D[pos] < 6)
      D[pos] = 6;
    if(D[pos+1] < 6)
      D[pos+1] = 6;
    if(D[pos+2] < 6)
      D[pos+2] = 6;
#if verbose > 0
    printf("D %03d %f %f %f\n", kk, D[pos], D[pos+1], D[pos+2]);
#endif
  }

  for(uint32_t kk=0; kk<Dm; kk++)
  {
    size_t pos = D[kk*3] + D[kk*3+1]*Vm + D[kk*3+2]*Vm*Vn;
    V[pos] = 5;
  }

  D[0] = 100; D[1] = 100; D[2] = 30;
  V[(int) D[0]+(int) D[1]*Vm+(int) D[2]*Vm*Vn] = 7;
  V[(int) D[0]+1+(int) D[1]*Vm+(int) D[2]*Vm*Vn] = 6;

  // Run the optimization
  mlfit1sn(V, Vm, Vn, Vp, D, Dm, F);

  // In next version, also supply clustering information

#if verbose >0
  for(size_t kk = 0; kk<Dm; kk++)
  {
    size_t pos = kk*3;
    size_t dpos = kk*4;
    printf("%6lu (%f, %f, %f) : s: %f, n: %f e: %f status: %d\n", kk,
        D[pos], D[pos+1], D[pos+2],
        F[dpos], F[dpos+1], F[dpos+2], (int) F[dpos+3]);
  }
#endif
  free(F);
  free(D);
  free(V);
  return 0;
}

int main(int argc, char ** argv)
  // For testing, not used any more, see the MATLAB interface in
  // df_mlfit.c
{

  printf("%s\n", argv[0]);

  if(argc == 1) {
    return unit_tests();
  }

  return 1;
}
#endif

max_stars_repo_name: elgw/dotter

max_issues_repo_name: elgw/dotter

max_forks_repo_name: elgw/dotter


text:
#ifdef HAVE_GSL
#include <gsl/gsl_sf_fermi_dirac.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_errno.h>
#else
#include <stddef.h>
#ifdef __cplusplus
extern "C"
{
#endif
   typedef  double gsl_multifit_linear_workspace ;
   typedef  struct { size_t size; double *v;}  gsl_vector ;
   typedef  struct { size_t size_rows; size_t size_columns; double* m;}  gsl_matrix ;
   double gsl_sf_fermi_dirac_mhalf(double x);
   double gsl_sf_fermi_dirac_half(double x);
   double gsl_sf_bessel_Kn(const int n, const double x);

   double gsl_vector_get (const gsl_vector * v, size_t i);
   gsl_matrix *gsl_matrix_alloc (size_t i, size_t j);
   gsl_vector *gsl_vector_alloc (size_t i);
   void gsl_matrix_set (gsl_matrix *m, size_t i, size_t j, double x);
   void gsl_vector_set (gsl_vector *v, size_t i, double x);
   void gsl_vector_free (gsl_vector * v);
   void gsl_matrix_free (gsl_matrix * v);
   gsl_multifit_linear_workspace *gsl_multifit_linear_alloc (size_t m,size_t n);
   void gsl_multifit_linear_free (gsl_multifit_linear_workspace *);
   int gsl_multifit_wlinear (const gsl_matrix * X, const gsl_vector * w, const gsl_vector * y, gsl_vector * c, gsl_matrix * cov, double * chisq, gsl_multifit_linear_workspace * work);
#ifdef __cplusplus
}
#endif
#endif

max_stars_repo_name: gonsie/ddcMD

max_issues_repo_name: gonsie/ddcMD

max_forks_repo_name: gonsie/ddcMD


text:
#include <petsc.h>

EXTERN_C_BEGIN
extern void formInitial(int*,int*,int*,double*,
                        double*,double*);
extern void formFunction(const int*,const int*,const int*,const double*,
                         const double*,const double[],const double[],double[]);
EXTERN_C_END

typedef struct AppCtx {
  PetscInt    nx,ny,nz;
  PetscScalar h[3];
} AppCtx;

#undef __FUNCT__
#define __FUNCT__ "FormInitial"
PetscErrorCode FormInitial(PetscReal t, Vec X, void *ctx)
{
  PetscScalar    *x;
  AppCtx         *app = (AppCtx*) ctx;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = VecGetArray(X,&x);CHKERRQ(ierr);
  /**/
  formInitial(&app->nx,&app->ny,&app->nz,app->h,
              &t,x);
  /**/
  ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "FormFunction"
PetscErrorCode FormFunction(TS ts, PetscReal t, Vec X, Vec Xdot,Vec F, void *ctx)
{
  const PetscScalar *x;
  const PetscScalar    *xdot;
  PetscScalar    *f;
  AppCtx         *app = (AppCtx*) ctx;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
  ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr);
  ierr = VecGetArray(F,&f);CHKERRQ(ierr);
  /**/
  formFunction(&app->nx,&app->ny,&app->nz,app->h,
               &t,x,xdot,f);
  /**/
  ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr);
  ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "RunTest"
PetscErrorCode RunTest(int nx, int ny, int nz, int loops, double *wt)
{
  Vec            x,f;
  TS             ts;
  AppCtx         _app,*app=&_app;
  double         t1,t2;
  PetscErrorCode ierr;
  PetscFunctionBegin;

  app->nx = nx; app->h[0] = 1./(nx-1);
  app->ny = ny; app->h[1] = 1./(ny-1);
  app->nz = nz; app->h[2] = 1./(nz-1);

  ierr = VecCreate(PETSC_COMM_SELF,&x);CHKERRQ(ierr);
  ierr = VecSetSizes(x,nx*ny*nz,nx*ny*nz);CHKERRQ(ierr);
  ierr = VecSetUp(x);CHKERRQ(ierr);
  ierr = VecDuplicate(x,&f);CHKERRQ(ierr);

  ierr = TSCreate(PETSC_COMM_SELF,&ts);CHKERRQ(ierr);
  ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
  ierr = TSSetType(ts,TSTHETA);CHKERRQ(ierr);
  ierr = TSThetaSetTheta(ts,1.0);CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts,0.01);CHKERRQ(ierr);
  ierr = TSSetTime(ts,0.0);CHKERRQ(ierr);
  ierr = TSSetDuration(ts,10,1.0);CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);

  ierr = TSSetSolution(ts,x);CHKERRQ(ierr);
  ierr = TSSetIFunction(ts,f,FormFunction,app);CHKERRQ(ierr);
  ierr = PetscOptionsSetValue(NULL,"-snes_mf","1");CHKERRQ(ierr);
  {
    SNES snes;
    KSP  ksp;
    ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
    ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
    ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
  }
  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
  ierr = TSSetUp(ts);CHKERRQ(ierr);

  *wt = 1e300;
  while (loops-- > 0) {
    ierr = FormInitial(0.0,x,app);CHKERRQ(ierr);
    ierr = PetscTime(&t1);CHKERRQ(ierr);
    ierr = TSSolve(ts,x);CHKERRQ(ierr);
    ierr = PetscTime(&t2);CHKERRQ(ierr);
    *wt = PetscMin(*wt,t2-t1);
  }

  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&f);CHKERRQ(ierr);
  ierr = TSDestroy(&ts);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "GetInt"
PetscErrorCode GetInt(const char* name, PetscInt *v, PetscInt defv)
{
  PetscErrorCode ierr;
  PetscFunctionBegin;
  *v = defv;
  ierr = PetscOptionsGetInt(NULL,NULL,name,v,NULL);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char *argv[])
{
  double         wt;
  PetscInt       n,start,step,stop,samples;
  PetscErrorCode ierr;

  ierr = PetscInitialize(&argc,&argv,NULL,NULL);CHKERRQ(ierr);

  ierr = GetInt("-start",   &start,   12);CHKERRQ(ierr);
  ierr = GetInt("-step",    &step,    4);CHKERRQ(ierr);
  ierr = GetInt("-stop",    &stop,    start);CHKERRQ(ierr);
  ierr = GetInt("-samples", &samples, 1);CHKERRQ(ierr);

  for (n=start; n<=stop; n+=step){
    int nx=n+1, ny=n+1, nz=n+1;
    ierr = RunTest(nx,ny,nz,samples,&wt);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_SELF,
                       "Grid  %3d x %3d x %3d -> %f seconds (%2d samples)\n",
                       nx,ny,nz,wt,samples);CHKERRQ(ierr);
  }

  ierr = PetscFinalize();CHKERRQ(ierr);
  return 0;
}


text:
/**
 *
 * @file core_strasm.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:46 2014
 *
 **/
#include <cblas.h>
#include <math.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_strasm - Computes the sums of the absolute values of elements in a same
 *  row or column in a triangular matrix.
 *  This function is an auxiliary function to triangular matrix norm computations.
 *
 *******************************************************************************
 *
 * @param[in] storev
 *          Specifies whether the sums are made per column or row.
 *          = PlasmaColumnwise: Computes the sum on each column
 *          = PlasmaRowwise:    Computes the sum on each row
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or lower triangular
 *          = PlasmaUpper: Upper triangle of A is referenced;
 *          = PlasmaLower: Lower triangle of A is referenced.
 *
 * @param[in] diag
 *          Specifies whether or not A is unit triangular:
 *          = PlasmaNonUnit: A is non unit;
 *          = PlasmaUnit:    A us unit.
 *
 * @param[in] M
 *          M specifies the number of rows of the matrix A. M >= 0.
 *
 * @param[in] N
 *          N specifies the number of columns of the matrix A. N >= 0.
 *
 * @param[in] A
 *          A is a M-by-N matrix.
 *
 * @param[in] lda
 *          The leading dimension of the array A. lda >= max(1,M).
 *
 * @param[in,out] work
 *          Array of dimension M if storev = PlasmaRowwise; N otherwise.
 *          On exit, contains the sums of the absolute values per column or row
 *          added to the input values.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_strasm = PCORE_strasm
#define CORE_strasm PCORE_strasm
#endif
void CORE_strasm(PLASMA_enum storev, PLASMA_enum uplo, PLASMA_enum diag,
                 int M, int N,
                 const float *A, int lda, float *work)
{
    const float *tmpA;
    int i, j, imax;
    int idiag = (diag == PlasmaUnit) ? 1 : 0;

    /*
     * PlasmaUpper / PlasmaColumnwise
     */
    if  (uplo == PlasmaUpper ) {
        M = min(M, N);

        if (storev == PlasmaColumnwise) {
            for (j = 0; j < N; j++) {
                tmpA = A+(j*lda);
                imax = min(j+1-idiag, M);

                if ( j < M )
                    work[j] += idiag;

                for (i = 0; i < imax; i++) {
                    work[j] += fabsf(*tmpA);
                    tmpA++;
                }
            }
        }
        /*
         * PlasmaUpper / PlasmaRowwise
         */
        else {
            if (diag == PlasmaUnit) {
                for (i = 0; i < M; i++) {
                    work[i] += 1.;
                }
            }
            for (j = 0; j < N; j++) {
                tmpA = A+(j*lda);
                imax = min(j+1-idiag, M);

                for (i = 0; i < imax; i++) {
                    work[i] += fabsf(*tmpA);
                    tmpA++;
                }
            }
        }
    } else {
        N = min(M, N);

        /*
         * PlasmaLower / PlasmaColumnwise
         */
        if (storev == PlasmaColumnwise) {
            for (j = 0; j < N; j++) {
                tmpA = A + j * (lda+1) + idiag;

                work[j] += idiag;
                for (i = j+idiag; i < M; i++) {
                    work[j] += fabsf(*tmpA);
                    tmpA++;
                }
            }
        }
        /*
         * PlasmaLower / PlasmaRowwise
         */
        else {
            if (diag == PlasmaUnit) {
                for (i = 0; i < N; i++) {
                    work[i] += 1.;
                }
            }
            for (j = 0; j < N; j++) {
                tmpA = A + j * (lda+1) + idiag;

                for (i = j+idiag; i < M; i++) {
                    work[i] += fabsf(*tmpA);
                    tmpA++;
                }
            }
        }
    }
}


text:
#ifndef INFMCMC_H
#define INFMCMC_H

#include <complex.h>
#include <fftw3.h>
#include <gsl/gsl_rng.h>

//void sampleRMWH(CHAIN *C);
//void sampleIndependenceSampler(CHAIN *C);
//void updateMean(CHAIN *C);
//void updateVar(CHAIN *C);

struct _CHAIN {
  int nj, nk; // number of Fourier coefficients in x/y direction respectively
  int numKeptSamples;
  int sizeObsVector;
  int currentIter;
  int accepted;
  double _shortTimeAccProbAvg;
  double _bLow, _bHigh;

  double *currentPhysicalState, *avgPhysicalState, *varPhysicalState, *_M2;
  double *proposedPhysicalState;
  double logLHDCurrentState;
  double accProb, avgAccProb;
  double rwmhStepSize, alphaPrior, priorVar, priorStd;

  // -- potentially not used --
  double *currentStateObservations, *proposedStateObservations;
  double *data;
  double currentLSQFunctional;
  double currentStateL2Norm2;
  double obsStdDev;
  // --------------------------

  fftw_complex *currentSpectralState, *avgSpectralState;
  fftw_complex *priorDraw, *proposedSpectralState;

  fftw_plan _c2r;
  fftw_plan _r2c;

  gsl_rng *r;
};

typedef struct _CHAIN CHAIN;
typedef struct _CHAIN INFCHAIN;

void infmcmc_initChain(INFCHAIN *C, const int nj, const int nk);
void infmcmc_freeChain(INFCHAIN *C);
void infmcmc_resetChain(INFCHAIN *C);
void infmcmc_proposeRWMH(INFCHAIN *C);
void infmcmc_updateRWMH(INFCHAIN *C, double logLHDOfProposal);
void infmcmc_seedWithPriorDraw(INFCHAIN *C);
void infmcmc_writeChainInfo(const INFCHAIN *C, FILE *fp);
void infmcmc_writeVFChain(const INFCHAIN *U, const INFCHAIN *V, FILE *fp);
void infmcmc_writeChain(const INFCHAIN *C, FILE *fp);
void infmcmc_printChain(INFCHAIN *C);
void infmcmc_setRWMHStepSize(INFCHAIN *C, double beta);
void infmcmc_adaptRWMHStepSize(INFCHAIN *C, double inc);
void infmcmc_setPriorAlpha(INFCHAIN *C, double alpha);
void infmcmc_setPriorVar(INFCHAIN *C, double var);
double infmcmc_L2Current(INFCHAIN *C);
double infmcmc_L2Proposed(INFCHAIN *C);
double infmcmc_L2Prior(INFCHAIN *C);

void infmcmc_seedWithDivFreePriorDraw(INFCHAIN *C1, INFCHAIN *C2);
void infmcmc_proposeDivFreeRWMH(INFCHAIN *C1, INFCHAIN *C2);
void infmcmc_updateVectorFieldRWMH(INFCHAIN *C1, INFCHAIN *C2, double logLHDOfProposal);

void randomPriorDrawOLD(gsl_rng *r, double PRIOR_ALPHA, fftw_complex *randDrawCoeffs);

#endif

max_stars_repo_name: dmcdougall/mcmclib

max_issues_repo_name: dmcdougall/mcmclib

max_forks_repo_name: dmcdougall/mcmclib


text:
// Copyright (c) 2021 Stig Rune Sellevag
//
// This file is distributed under the MIT License. See the accompanying file
// LICENSE.txt or http://www.opensource.org/licenses/mit-license.php for terms
// and conditions.

#ifndef SCILIB_LINALG_MATRIX_DECOMPOSITION_H
#define SCILIB_LINALG_MATRIX_DECOMPOSITION_H

#ifdef USE_MKL
#include <mkl.h>
#else
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreturn-type-c-linkage"
#endif
#include <lapacke.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#endif

#include <scilib/mdarray.h>
#include <scilib/linalg_impl/lapack_types.h>
#include <experimental/mdspan>
#include <exception>
#include <cassert>
#include <type_traits>

namespace Sci {
namespace Linalg {

namespace stdex = std::experimental;

// LU factorization.
template <class Layout>
inline void lu(Sci::Matrix_view<double, Layout> a,
               Sci::Vector_view<BLAS_INT, Layout> ipiv)
{
    static_assert(a.is_contiguous());
    static_assert(ipiv.is_contiguous());

    const BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT n = static_cast<BLAS_INT>(a.extent(1));

    assert(static_cast<BLAS_INT>(ipiv.size()) >= std::min(m, n));

    auto matrix_layout = LAPACK_ROW_MAJOR;
    BLAS_INT lda = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = LAPACK_COL_MAJOR;
        lda = m;
    }

    BLAS_INT info =
        LAPACKE_dgetrf(matrix_layout, m, n, a.data(), lda, ipiv.data());
    if (info < 0) {
        throw std::runtime_error("dgetrf: illegal input parameter");
    }
    if (info > 0) {
        throw std::runtime_error("dgetrf: U matrix is singular");
    }
}

template <class Layout, class Allocator_a, class Allocator_ipiv>
inline void lu(Sci::Matrix<double, Layout, Allocator_a>& a,
               Sci::Vector<BLAS_INT, Layout, Allocator_ipiv>& ipiv)
{
    lu(a.view(), ipiv.view());
}

// QR factorization.
template <class Layout>
inline void qr(Sci::Matrix_view<double, Layout> a,
               Sci::Matrix_view<double, Layout> q,
               Sci::Matrix_view<double, Layout> r)
{
    assert(q.extent(0) == a.extent(0) && q.extent(1) == a.extent(1));
    assert(r.extent(0) == a.extent(0) && r.extent(1) == a.extent(1));

    const BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT n = static_cast<BLAS_INT>(a.extent(1));

    auto matrix_layout = LAPACK_ROW_MAJOR;
    BLAS_INT lda = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = LAPACK_COL_MAJOR;
        lda = m;
    }
    Sci::copy(a, q);
    Sci::Vector<double, Layout> tau(std::min(m, n));

    // Compute QR factorization:

    BLAS_INT info =
        LAPACKE_dgeqrf(matrix_layout, m, n, q.data(), lda, tau.data());
    if (info != 0) {
        throw std::runtime_error("dgeqrf failed");
    }

    // Compute Q:

    info = LAPACKE_dorgqr(matrix_layout, m, n, n, q.data(), lda, tau.data());
    if (info != 0) {
        throw std::runtime_error("dorgqr failed");
    }

    // Compute R:

    matrix_product(transposed(q), a, r);
    transposed(q);
}

template <class Layout, class Allocator>
inline void qr(Sci::Matrix<double, Layout, Allocator>& a,
               Sci::Matrix<double, Layout, Allocator>& q,
               Sci::Matrix<double, Layout, Allocator>& r)
{
    qr(a.view(), q.view(), r.view());
}

// Singular value decomposition.
template <class Layout>
inline void svd(Sci::Matrix_view<double, Layout> a,
                Sci::Vector_view<double, Layout> s,
                Sci::Matrix_view<double, Layout> u,
                Sci::Matrix_view<double, Layout> vt)
{
    const BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    const BLAS_INT n = static_cast<BLAS_INT>(a.extent(1));
    const BLAS_INT ldu = m;
    const BLAS_INT ldvt = n;

    assert(static_cast<BLAS_INT>(s.extent(0)) == std::min(m, n));
    assert(static_cast<BLAS_INT>(u.extent(0)) == m);
    assert(static_cast<BLAS_INT>(u.extent(1)) == ldu);
    assert(static_cast<BLAS_INT>(vt.extent(0)) == n);
    assert(static_cast<BLAS_INT>(vt.extent(1)) == ldvt);

    auto matrix_layout = LAPACK_ROW_MAJOR;
    BLAS_INT lda = n;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = LAPACK_COL_MAJOR;
        lda = m;
    }

    Sci::Vector<double, Layout> superb(std::min(m, n) - 1);

    BLAS_INT info =
        LAPACKE_dgesvd(matrix_layout, 'A', 'A', m, n, a.data(), lda, s.data(),
                       u.data(), ldu, vt.data(), ldvt, superb.data());
    if (info != 0) {
        throw std::runtime_error("dgesvd failed");
    }
}

template <class Layout, class Allocator>
inline void svd(Sci::Matrix<double, Layout, Allocator>& a,
                Sci::Vector<double, Layout, Allocator>& s,
                Sci::Matrix<double, Layout, Allocator>& u,
                Sci::Matrix<double, Layout, Allocator>& vt)
{
    svd(a.view(), s.view(), u.view(), vt.view());
}

} // namespace Linalg
} // namespace Sci

#endif // SCILIB_LINALG_MATRIX_DECOMPOSITION_H

max_stars_repo_name: stigrs/scilib

max_issues_repo_name: stigrs/scilib

max_forks_repo_name: stigrs/scilib


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

#ifndef SIM_PHYSICS_AVL_AERO_DATABASE_H_
#define SIM_PHYSICS_AVL_AERO_DATABASE_H_

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

#include <stdint.h>

#include <string>

#include "common/c_math/linalg.h"
#include "common/c_math/vec3.h"
#include "common/macros.h"
#include "sim/sim_params.h"
#include "sim/sim_types.h"
#include "system/labels.h"

class AvlAeroDatabase {
  friend class AeroTest;

 public:
  AvlAeroDatabase(const std::string &filename,
                  const AeroSimParams &aero_sim_params);
  ~AvlAeroDatabase();

  // Applies coefficient fudge factors.
  //
  // Args:
  //   coeff_offsets: Offsets to apply.
  //   force_coeff_w_scale_factors: Scale factors [#] to apply to the
  //       force coefficients.  These scale factors should be represented in
  //       the wind axes.
  //   moment_coeff_b_scale_factors: Scale factors [#] to apply to the
  //       moment coefficients.  These scale factors should be represented in
  //       the body axes.
  void ApplyCoeffAdjustments(const AeroCoeffOffsets &coeff_offsets,
                             const AeroCoeffs &force_coeff_w_scale_factors,
                             const AeroCoeffs &moment_coeff_b_scale_factors);

  // Returns the total force-moment coefficient at a specific (alpha,
  // beta, omega_hat, flap deflection) operating point.
  void CalcForceMomentCoeff(double alpha, double beta, const Vec3 &omega_hat,
                            const Vec &flaps, ForceMoment *cfm,
                            double thrust_coeff) const {
    CalcForceMomentCoeff(alpha, beta, omega_hat, flaps, cfm, thrust_coeff,
                         nullptr);
  }

  // See function above.  Optionally returns the breakdown of the
  // force-moment coefficient into control and stability derivatives
  // if aero_coeffs is not null.
  void CalcForceMomentCoeff(double alpha, double beta, const Vec3 &omega_hat,
                            const Vec &flaps, ForceMoment *cfm,
                            double thrust_coeff,
                            AvlAeroCoeffs *aero_coeffs) const;

  // Returns the nominal elevator deflection for a given angle-of-attack
  // and sideslip angle.
  double GetNominalElevatorDeflection(double alpha, double beta) const;

  double reynolds_number() const { return reynolds_number_; }

  int32_t num_flaps() const { return num_flaps_; }
  const Vec3 &omega_hat_0() const { return omega_hat_0_; }

 private:
  // Returns true if the database passes basic sign checks.
  bool IsValid() const;

  // Reduces flap effectiveness at large deflections.
  void AdjustNonlinearFlaps(const Vec &flaps, Vec *flaps_nonlin) const;

  // Returns the nominal force-moment coefficient and the control and
  // stability derivatives.
  void CalcAvlAeroCoeffs(double alpha, double beta,
                         AvlAeroCoeffs *coeffs) const;

  const AeroSimParams &aero_sim_params_;

  double reynolds_number_;
  int32_t num_flaps_;
  gsl_vector *alphads_, *betads_;
  gsl_matrix *eleds_;

  Vec3 omega_hat_0_;
  gsl_matrix *CX_, *CY_, *CZ_, *CL_, *CM_, *CN_;
  gsl_matrix *CXp_, *CYp_, *CZp_, *CLp_, *CMp_, *CNp_;
  gsl_matrix *CXq_, *CYq_, *CZq_, *CLq_, *CMq_, *CNq_;
  gsl_matrix *CXr_, *CYr_, *CZr_, *CLr_, *CMr_, *CNr_;
  gsl_matrix *CXd_[kNumFlaps];
  gsl_matrix *CYd_[kNumFlaps];
  gsl_matrix *CZd_[kNumFlaps];
  gsl_matrix *CLd_[kNumFlaps];
  gsl_matrix *CMd_[kNumFlaps];
  gsl_matrix *CNd_[kNumFlaps];

  DISALLOW_COPY_AND_ASSIGN(AvlAeroDatabase);
};

#endif  // SIM_PHYSICS_AVL_AERO_DATABASE_H_


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

#include "dcp/dcp-types.h"
#include "evp_engine_test.h"
#include "evp_store_single_threaded_test.h"
#include "vbucket_fwd.h"
#include <memcached/engine_error.h>
#include <gsl/gsl>

class Item;
class MockDcpProducer;
class MockActiveStream;
struct dcp_message_producers;

/**
 * Test fixture for unit tests related to DCP.
 */
class DCPTest : public EventuallyPersistentEngineTest {
protected:
    void SetUp() override;

    void TearDown() override;

    // Create a DCP producer; initially with no streams associated.
    void create_dcp_producer(
            int flags = 0,
            IncludeValue includeVal = IncludeValue::Yes,
            IncludeXattrs includeXattrs = IncludeXattrs::Yes,
            std::vector<std::pair<std::string, std::string>> controls = {});

    // Setup a DCP producer and attach a stream and cursor to it.
    void setup_dcp_stream(
            int flags = 0,
            IncludeValue includeVal = IncludeValue::Yes,
            IncludeXattrs includeXattrs = IncludeXattrs::Yes,
            std::vector<std::pair<std::string, std::string>> controls = {});

    void destroy_dcp_stream();

    struct StreamRequestResult {
        ENGINE_ERROR_CODE status;
        uint64_t rollbackSeqno;
    };

    /**
     * Helper function to simplify calling producer->streamRequest() - provides
     * sensible default values for common invocations.
     */
    static StreamRequestResult doStreamRequest(DcpProducer& producer,
                                               uint64_t startSeqno = 0,
                                               uint64_t endSeqno = ~0,
                                               uint64_t snapStart = 0,
                                               uint64_t snapEnd = ~0,
                                               uint64_t vbUUID = 0);

    /**
     * Helper function to simplify the process of preparing DCP items to be
     * fetched from the stream via step().
     * Should be called once all expected items are present in the vbuckets'
     * checkpoint, it will then run the correct background tasks to
     * copy CheckpointManager items to the producers' readyQ.
     */
    static void prepareCheckpointItemsForStep(
            dcp_message_producers& msgProducers,
            MockDcpProducer& producer,
            VBucket& vb);

    /*
     * Creates an item with the key \"key\", containing json data and xattrs.
     * @return a unique_ptr to a newly created item.
     */
    std::unique_ptr<Item> makeItemWithXattrs();

    /*
     * Creates an item with the key \"key\", containing json data and no xattrs.
     * @return a unique_ptr to a newly created item.
     */
    std::unique_ptr<Item> makeItemWithoutXattrs();

    /* Add items onto the vbucket and wait for the checkpoint to be removed */
    void addItemsAndRemoveCheckpoint(int numItems);

    void removeCheckpoint(int numItems);

    void runCheckpointProcessor(dcp_message_producers& producers);

    std::shared_ptr<MockDcpProducer> producer;
    std::shared_ptr<MockActiveStream> stream;
    VBucketPtr vb0;

    /*
     * Fake callback emulating dcp_add_failover_log
     */
    static ENGINE_ERROR_CODE fakeDcpAddFailoverLog(
            vbucket_failover_t* entry,
            size_t nentries,
            gsl::not_null<const void*> cookie) {
        callbackCount++;
        return ENGINE_SUCCESS;
    }

    // callbackCount needs to be static as its used inside of the static
    // function fakeDcpAddFailoverLog.
    static int callbackCount;
};


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

/* Author:  G. Jungman */

#ifndef __GSL_SF_TRIG_H__
#define __GSL_SF_TRIG_H__

#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_types.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/* Sin(x) with GSL semantics. This is actually important
 * because we want to control the error estimate, and trying
 * to guess the error for the standard library implementation
 * every time it is used would be a little goofy.
 */
GSL_EXPORT int gsl_sf_sin_e(double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_sin(const double x);

/* Cos(x) with GSL semantics.
 */
GSL_EXPORT int gsl_sf_cos_e(double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_cos(const double x);

/* Hypot(x,y) with GSL semantics.
 */
GSL_EXPORT int gsl_sf_hypot_e(const double x, const double y, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_hypot(const double x, const double y);

/* Sin(z) for complex z
 *
 * exceptions: GSL_EOVRFLW
 */
GSL_EXPORT int gsl_sf_complex_sin_e(const double zr, const double zi, gsl_sf_result * szr, gsl_sf_result * szi);

/* Cos(z) for complex z
 *
 * exceptions: GSL_EOVRFLW
 */
GSL_EXPORT int gsl_sf_complex_cos_e(const double zr, const double zi, gsl_sf_result * czr, gsl_sf_result * czi);

/* Log(Sin(z)) for complex z
 *
 * exceptions: GSL_EDOM, GSL_ELOSS
 */
GSL_EXPORT int gsl_sf_complex_logsin_e(const double zr, const double zi, gsl_sf_result * lszr, gsl_sf_result * lszi);

/* Sinc(x) = sin(pi x) / (pi x)
 *
 * exceptions: none
 */
GSL_EXPORT int gsl_sf_sinc_e(double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_sinc(const double x);

/* Log(Sinh(x)), x > 0
 *
 * exceptions: GSL_EDOM
 */
GSL_EXPORT int gsl_sf_lnsinh_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_lnsinh(const double x);

/* Log(Cosh(x))
 *
 * exceptions: none
 */
GSL_EXPORT int gsl_sf_lncosh_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_lncosh(const double x);

/* Convert polar to rectlinear coordinates.
 *
 * exceptions: GSL_ELOSS
 */
GSL_EXPORT int gsl_sf_polar_to_rect(const double r, const double theta, gsl_sf_result * x, gsl_sf_result * y);

/* Convert rectilinear to polar coordinates.
 * return argument in range [-pi, pi]
 *
 * exceptions: GSL_EDOM
 */
GSL_EXPORT int gsl_sf_rect_to_polar(const double x, const double y, gsl_sf_result * r, gsl_sf_result * theta);

/* Sin(x) for quantity with an associated error.
 */
GSL_EXPORT int gsl_sf_sin_err_e(const double x, const double dx, gsl_sf_result * result);

/* Cos(x) for quantity with an associated error.
 */
GSL_EXPORT int gsl_sf_cos_err_e(const double x, const double dx, gsl_sf_result * result);

/* Force an angle to lie in the range (-pi,pi].
 *
 * exceptions: GSL_ELOSS
 */
GSL_EXPORT int gsl_sf_angle_restrict_symm_e(double * theta);
GSL_EXPORT double gsl_sf_angle_restrict_symm(const double theta);

/* Force an angle to lie in the range [0, 2pi)
 *
 * exceptions: GSL_ELOSS
 */
GSL_EXPORT int gsl_sf_angle_restrict_pos_e(double * theta);
GSL_EXPORT double gsl_sf_angle_restrict_pos(const double theta);

GSL_EXPORT int gsl_sf_angle_restrict_symm_err_e(const double theta, gsl_sf_result * result);

GSL_EXPORT int gsl_sf_angle_restrict_pos_err_e(const double theta, gsl_sf_result * result);

__END_DECLS

#endif /* __GSL_SF_TRIG_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#ifndef GBNET_MODELORNOR
#define GBNET_MODELORNOR

#include <cmath>
#include <time.h>

#include <gsl/gsl_rstat.h>

#include "ModelBase.h"
#include "GraphORNOR.h"

namespace gbn
{

    class ModelORNOR: public ModelBase
    {
        private:
        protected:
        public:
            ModelORNOR();
            ~ModelORNOR() override;
            ModelORNOR(const network_t, const evidence_dict_t, const prior_active_tf_set_t = prior_active_tf_set_t(),
                       const double [3 * 3] = gbn::SPRIOR, double z_alpha = 25., double z_beta = 25., double z0_alpha = 25., double z0_beta = 25., double t_alpha = 25., double t_beta = 25.,
                       unsigned int = 3, bool = true, bool = true, bool = false,
                       double = 2., double = 2., double = 2., double = 2., double = 8., double = 2.);
    };
}

#endif

max_stars_repo_name: umbibio/gbnet

max_issues_repo_name: umbibio/gbnet

max_forks_repo_name: umbibio/gbnet