text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for the geometric coefficients entering the response for LIGO-VIRGO detectors.
 *
 *
 */

#ifndef _LLVGEOMETRY_H
#define _LLVGEOMETRY_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"

/********************************************************/
/**************** Type definitions **********************/

/* Enumerator for the detector selector and network selector */
typedef enum {LHO, LLO, VIRGO} Detectortag;
typedef enum {L, H, V, LH, LV, HV, LHV} Networktag;

/*************************************************************/
/**************** Geometrical constants **********************/

/**
 * \name LIGO Hanford Observatory 4km Interferometric Detector constants
 * The following constants describe the location and geometry of the
 * LIGO Hanford Observatory 4km Interferometric Detector.
 */
#define LAL_LHO_4K_DETECTOR_NAME               	"LHO_4k"	/**< LHO_4k detector name string */
#define LAL_LHO_4K_DETECTOR_PREFIX             	"H1"	/**< LHO_4k detector prefix string */
#define LAL_LHO_4K_DETECTOR_LONGITUDE_RAD      	-2.08405676917	/**< LHO_4k vertex longitude (rad) */
#define LAL_LHO_4K_DETECTOR_LATITUDE_RAD       	0.81079526383	/**< LHO_4k vertex latitude (rad) */
#define LAL_LHO_4K_DETECTOR_ELEVATION_SI       	142.554	/**< LHO_4k vertex elevation (m) */
#define LAL_LHO_4K_DETECTOR_ARM_X_AZIMUTH_RAD  	5.65487724844	/**< LHO_4k x arm azimuth (rad) */
#define LAL_LHO_4K_DETECTOR_ARM_Y_AZIMUTH_RAD  	4.08408092164	/**< LHO_4k y arm azimuth (rad) */
#define LAL_LHO_4K_DETECTOR_ARM_X_ALTITUDE_RAD 	-0.00061950000	/**< LHO_4k x arm altitude (rad) */
#define LAL_LHO_4K_DETECTOR_ARM_Y_ALTITUDE_RAD 	0.00001250000	/**< LHO_4k y arm altitude (rad) */
#define LAL_LHO_4K_DETECTOR_ARM_X_MIDPOINT_SI  	1997.54200000000	/**< LHO_4k x arm midpoint (m) */
#define LAL_LHO_4K_DETECTOR_ARM_Y_MIDPOINT_SI  	1997.52200000000	/**< LHO_4k y arm midpoint (m) */
#define LAL_LHO_4K_VERTEX_LOCATION_X_SI        	-2.16141492636e+06	/**< LHO_4k x-component of vertex location in Earth-centered frame (m) */
#define LAL_LHO_4K_VERTEX_LOCATION_Y_SI        	-3.83469517889e+06	/**< LHO_4k y-component of vertex location in Earth-centered frame (m) */
#define LAL_LHO_4K_VERTEX_LOCATION_Z_SI        	4.60035022664e+06	/**< LHO_4k z-component of vertex location in Earth-centered frame (m) */
#define LAL_LHO_4K_ARM_X_DIRECTION_X           	-0.22389266154	/**< LHO_4k x-component of unit vector pointing along x arm in Earth-centered frame */
#define LAL_LHO_4K_ARM_X_DIRECTION_Y           	0.79983062746	/**< LHO_4k y-component of unit vector pointing along x arm in Earth-centered frame */
#define LAL_LHO_4K_ARM_X_DIRECTION_Z           	0.55690487831	/**< LHO_4k z-component of unit vector pointing along x arm in Earth-centered frame */
#define LAL_LHO_4K_ARM_Y_DIRECTION_X           	-0.91397818574	/**< LHO_4k x-component of unit vector pointing along y arm in Earth-centered frame */
#define LAL_LHO_4K_ARM_Y_DIRECTION_Y           	0.02609403989	/**< LHO_4k y-component of unit vector pointing along y arm in Earth-centered frame */
#define LAL_LHO_4K_ARM_Y_DIRECTION_Z           	-0.40492342125	/**< LHO_4k z-component of unit vector pointing along y arm in Earth-centered frame */

/**
 * \name LIGO Livingston Observatory 4km Interferometric Detector constants
 * The following constants describe the location and geometry of the
 * LIGO Livingston Observatory 4km Interferometric Detector.
 */
#define LAL_LLO_4K_DETECTOR_NAME               	"LLO_4k"	/**< LLO_4k detector name string */
#define LAL_LLO_4K_DETECTOR_PREFIX             	"L1"	/**< LLO_4k detector prefix string */
#define LAL_LLO_4K_DETECTOR_LONGITUDE_RAD      	-1.58430937078	/**< LLO_4k vertex longitude (rad) */
#define LAL_LLO_4K_DETECTOR_LATITUDE_RAD       	0.53342313506	/**< LLO_4k vertex latitude (rad) */
#define LAL_LLO_4K_DETECTOR_ELEVATION_SI       	-6.574	/**< LLO_4k vertex elevation (m) */
#define LAL_LLO_4K_DETECTOR_ARM_X_AZIMUTH_RAD  	4.40317772346	/**< LLO_4k x arm azimuth (rad) */
#define LAL_LLO_4K_DETECTOR_ARM_Y_AZIMUTH_RAD  	2.83238139666	/**< LLO_4k y arm azimuth (rad) */
#define LAL_LLO_4K_DETECTOR_ARM_X_ALTITUDE_RAD 	-0.00031210000	/**< LLO_4k x arm altitude (rad) */
#define LAL_LLO_4K_DETECTOR_ARM_Y_ALTITUDE_RAD 	-0.00061070000	/**< LLO_4k y arm altitude (rad) */
#define LAL_LLO_4K_DETECTOR_ARM_X_MIDPOINT_SI  	1997.57500000000	/**< LLO_4k x arm midpoint (m) */
#define LAL_LLO_4K_DETECTOR_ARM_Y_MIDPOINT_SI  	1997.57500000000	/**< LLO_4k y arm midpoint (m) */
#define LAL_LLO_4K_VERTEX_LOCATION_X_SI        	-7.42760447238e+04	/**< LLO_4k x-component of vertex location in Earth-centered frame (m) */
#define LAL_LLO_4K_VERTEX_LOCATION_Y_SI        	-5.49628371971e+06	/**< LLO_4k y-component of vertex location in Earth-centered frame (m) */
#define LAL_LLO_4K_VERTEX_LOCATION_Z_SI        	3.22425701744e+06	/**< LLO_4k z-component of vertex location in Earth-centered frame (m) */
#define LAL_LLO_4K_ARM_X_DIRECTION_X           	-0.95457412153	/**< LLO_4k x-component of unit vector pointing along x arm in Earth-centered frame */
#define LAL_LLO_4K_ARM_X_DIRECTION_Y           	-0.14158077340	/**< LLO_4k y-component of unit vector pointing along x arm in Earth-centered frame */
#define LAL_LLO_4K_ARM_X_DIRECTION_Z           	-0.26218911324	/**< LLO_4k z-component of unit vector pointing along x arm in Earth-centered frame */
#define LAL_LLO_4K_ARM_Y_DIRECTION_X           	0.29774156894	/**< LLO_4k x-component of unit vector pointing along y arm in Earth-centered frame */
#define LAL_LLO_4K_ARM_Y_DIRECTION_Y           	-0.48791033647	/**< LLO_4k y-component of unit vector pointing along y arm in Earth-centered frame */
#define LAL_LLO_4K_ARM_Y_DIRECTION_Z           	-0.82054461286	/**< LLO_4k z-component of unit vector pointing along y arm in Earth-centered frame */

/**
 * \name VIRGO 3km Interferometric Detector constants
 * The following constants describe the location and geometry of the
 * VIRGO 3km Interferometric Detector.
 */
#define LAL_VIRGO_DETECTOR_NAME               	"VIRGO"	/**< VIRGO detector name string */
#define LAL_VIRGO_DETECTOR_PREFIX             	"V1"	/**< VIRGO detector prefix string */
#define LAL_VIRGO_DETECTOR_LONGITUDE_RAD      	0.18333805213	/**< VIRGO vertex longitude (rad) */
#define LAL_VIRGO_DETECTOR_LATITUDE_RAD       	0.76151183984	/**< VIRGO vertex latitude (rad) */
#define LAL_VIRGO_DETECTOR_ELEVATION_SI       	51.884	/**< VIRGO vertex elevation (m) */
#define LAL_VIRGO_DETECTOR_ARM_X_AZIMUTH_RAD  	0.33916285222	/**< VIRGO x arm azimuth (rad) */
#define LAL_VIRGO_DETECTOR_ARM_Y_AZIMUTH_RAD  	5.05155183261	/**< VIRGO y arm azimuth (rad) */
#define LAL_VIRGO_DETECTOR_ARM_X_ALTITUDE_RAD 	0.00000000000	/**< VIRGO x arm altitude (rad) */
#define LAL_VIRGO_DETECTOR_ARM_Y_ALTITUDE_RAD 	0.00000000000	/**< VIRGO y arm altitude (rad) */
#define LAL_VIRGO_DETECTOR_ARM_X_MIDPOINT_SI  	1500.00000000000	/**< VIRGO x arm midpoint (m) */
#define LAL_VIRGO_DETECTOR_ARM_Y_MIDPOINT_SI  	1500.00000000000	/**< VIRGO y arm midpoint (m) */
#define LAL_VIRGO_VERTEX_LOCATION_X_SI        	4.54637409900e+06	/**< VIRGO x-component of vertex location in Earth-centered frame (m) */
#define LAL_VIRGO_VERTEX_LOCATION_Y_SI        	8.42989697626e+05	/**< VIRGO y-component of vertex location in Earth-centered frame (m) */
#define LAL_VIRGO_VERTEX_LOCATION_Z_SI        	4.37857696241e+06	/**< VIRGO z-component of vertex location in Earth-centered frame (m) */
#define LAL_VIRGO_ARM_X_DIRECTION_X           	-0.70045821479	/**< VIRGO x-component of unit vector pointing along x arm in Earth-centered frame */
#define LAL_VIRGO_ARM_X_DIRECTION_Y           	0.20848948619	/**< VIRGO y-component of unit vector pointing along x arm in Earth-centered frame */
#define LAL_VIRGO_ARM_X_DIRECTION_Z           	0.68256166277	/**< VIRGO z-component of unit vector pointing along x arm in Earth-centered frame */
#define LAL_VIRGO_ARM_Y_DIRECTION_X           	-0.05379255368	/**< VIRGO x-component of unit vector pointing along y arm in Earth-centered frame */
#define LAL_VIRGO_ARM_Y_DIRECTION_Y           	-0.96908180549	/**< VIRGO y-component of unit vector pointing along y arm in Earth-centered frame */
#define LAL_VIRGO_ARM_Y_DIRECTION_Z           	0.24080451708	/**< VIRGO z-component of unit vector pointing along y arm in Earth-centered frame */

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

#endif /* _LLVGEOMETRY_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
/* igraph library headers */

#include <igraph/igraph.h>

/* gsl headers */
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_histogram2d.h>

/*headers for hashing */
#include "uthash.h"		/* hash table macros */
#include <stddef.h>		/* offset of */

/* headers for argument parsing */
#include <stdlib.h>
#include <error.h>
#include <argp.h>

/* errno.h has declaration of errno variable */
#include <errno.h>

/** more stuff for argument parsing **/

const char *argp_program_version = "dynet 0.2";
const char *argp_program_bug_address = "<ebo@mail.utexas.edu>";

static char doc[] = "DyNet -- simulator for disease spreading\
 on dynamic networks\
\vThe simulator allows the network to relax from an initial\
 degree distribuion toward a target degree distribution. The\
 condition that degree correlations are close to zero is used\
 to determine the types of edges that are added and deleted.\
 No further documentation is available yet.";

/* A description of the arguments we accept. */
static char args_doc[] = " ";

/* Keys for options without short-options */
#define OPT_ABORT 1

/* macros and functions for network generation*/
#define DYNET_EXP (int) 1
#define DYNET_REG (int) 2
#define DYNET_PARETO (int) 3
#define DYNET_ER (int) 4

/* The options we understand. */
static struct argp_option options[] = {
  {0, 0, 0, 0, "General Options"},
  {"verbose", 'v', 0, 0, "Produce verbose output"},
  {"quiet", 'q', 0, 0, "Don't produce any outpu"},
  {"silent", 's', 0, OPTION_ALIAS},
  {"output", 'o', "FILE", 0,
   "Output to FILE instead of standard output"},
  {0, 0, 0, 0, "Network Model Options"},
  {"network_rate", 'n', "RATE", 0,
   "Network degree distribution approaches equilibirium "
   "at rate RATE (default 0.4)"},
  {"network_tension", 'k', "RATE", 0,
   "Network degree correlations are removed in proportion to "
   "thier size at rate RATE (default 50)"},
  {"network_size", 'N', "COUNT", 0,
   "the network is made to have about COUNT (default 10000) nodes "},
  {"network_type", 'y', "INTEGER_CODE", 0,
   "INTEGER_CODE (default 2) determines the limiting degree distribution:"
    " 1==exponential, 2==regular, 3==Pareto, 4==Poisson"},
  {"par1", 'm', "RATIONAL", 0,
   "If network_type is in {1, 2, 4}, RATIONAL (default 4) is "
    " the mean of the degree distribution. If network_type is 3,"
   " RATIONAL is A in definition of Pareto pdf. in GSL "
   "documentation"},
  {"par2", 'b', "RATIONAL", 0,
   "If network_type is PARETO, RATIONAL is B in the definition of"
   " the Pareto pdf in the GSL documentation"},
  {0, 0, 0, 0, "Disease Model Options"},
  {"trans_rate", 't', "RATE", 0,
   "Disease moves at rate RATE (default 2) across an edge"},
  {"recov_rate", 'r', "RATE", 0,
   "Infected nodes reacover at rate RATE (default 1)"},
  {"interval", 'i', "LENGTH", 0,
   "State variables are printed to output at every LENGTH\
 (default 0.05) time units"},
  {"epsilon", 'e', "FRACTION", 0, "Initial fraction FRACTION "
   "(default 0.1) of edges from suscetible nodes pointing "
   "to infected nodes"},
  {"finish_time", 'f', "LENGTH", 0, "Simulation time will start "
   "at zero and run LENGTH (default 5) time units"},
  {0, 0, 0, 0, "Evolutionary Model Options"},
  {"sample_fraction", 'p', "FRACTION", 0,
    "FRACTION (default 0.001) nodes in virus gene tree are selected"
      "for inclusion in newick output"},
  {0}
};

/* maximum rate allowed in arguments */
#define MAX_RATE 10000.

/* Used by `main' to communicate with `parse_opt'. */
struct arguments
{
  int silent, verbose;		/* `-s', `-v' */
  char *output_file;		/* FILE arg to `--output' */
  double trans_rate;		/* RATE arg to `--trans_rate' */
  double recov_rate;		/* RATE arg to `--recov_rate' */
  double interval;		/* LENGTH arg to `--interval' */
  double epsilon;		/* FRACTION arg to `--epsilon' */
  double finish_time;		/* LENGTH arg to `--finish_time' */
  double network_rate;		/* RATE arg to `--network_rate' */
  double network_tension;	/* RATE arg to `--network_tension' */
  int network_size;		/* COUNT arg to `--network_size' */
  int network_type;             /* INTEGER_CODE arg to `--network_type' */
  double par1;                  /* RATIONAL arg to `--par1' */
  double par2;                  /* RATIONAL arg to `--par2' */
  double sample_fraction;	/* FRACTION arg to `--sample_fraction' */

};

/* Parse a single option. */
static error_t
parse_opt (int key, char *arg, struct argp_state *state)
{
  /* Get the INPUT argument from `argp_parse', which we
     know is a pointer to our arguments structure. */
  struct arguments *arguments = (struct arguments *) state->input;

  switch (key)
    {
    case 'q':
    case 's':
      arguments->silent = 1;
      break;
    case 'v':
      arguments->verbose = 1;
      break;
    case 'o':
      arguments->output_file = arg;
      break;
    case 't':
      arguments->trans_rate = strtod (arg, NULL);
      if (arguments->trans_rate < 0 || arguments->trans_rate > MAX_RATE)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "trans_rate=%g, should be in [0, %g]",
			 arguments->trans_rate, MAX_RATE);

	}
      break;
    case 'r':
      arguments->recov_rate = strtod (arg, NULL);
      if (arguments->recov_rate < 0 || arguments->recov_rate > MAX_RATE)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "recov_rate=%g, should be in [0, %g]",
			 arguments->recov_rate, MAX_RATE);

	}
      break;
    case 'i':
      arguments->interval = strtod (arg, NULL);
      if (arguments->interval <= 0 || arguments->interval > 1e2)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "interval=%g, should be in [0, %g]",
			 arguments->recov_rate, 1e2);

	}
      break;
    case 'e':
      arguments->epsilon = strtod (arg, NULL);
      if (arguments->epsilon < 0 || arguments->epsilon > 1)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "epsilon=%g, should be in [0, 1]",
			 arguments->epsilon);

	}
      break;
    case 'f':
      arguments->finish_time = strtod (arg, NULL);
      if (arguments->finish_time < 0 || arguments->finish_time > 100)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "finish_time=%g, should be in [0, 100]",
			 arguments->finish_time);

	}
      break;
    case 'n':
      arguments->network_rate = strtod (arg, NULL);
      if (arguments->network_rate < 0 || arguments->network_rate > MAX_RATE)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "network_rate=%g, should be in [0, %g]",
			 arguments->network_rate, MAX_RATE);

	}
      break;
    case 'k':
      arguments->network_tension = strtod (arg, NULL);
      if (arguments->network_tension < 0
	  || arguments->network_tension > MAX_RATE)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "network_tension=%g, should be in [0, %g]",
			 arguments->network_tension, MAX_RATE);

	}
      break;
    case 'N':
      arguments->network_size = atoi (arg);
      if (arguments->network_size < 1 || arguments->network_size > 1e6)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "network_size=%d, should be in [1, %d]",
			 arguments->network_size, (int) 1e6);

	}
      break;
    case 'p':
      arguments->sample_fraction = strtod (arg, NULL);
      if (arguments->sample_fraction < 0 || arguments->sample_fraction > 1)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "sample_fraction =%g, should be in [0, 1]",
			 arguments->sample_fraction);

	}
      break;
    case 'y':
      arguments->network_type =  atoi (arg);
      if (arguments->network_type != DYNET_REG
          && arguments->network_type != DYNET_ER
          && arguments->network_type != DYNET_PARETO
          && arguments->network_type != DYNET_EXP)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "network type =%d, should be in {%d, %d, %d, %d}",
			 arguments->network_type, DYNET_EXP, DYNET_REG,
                         DYNET_PARETO, DYNET_ER);

	}
      break;
    case 'm':
      arguments->par1 = strtod (arg, NULL);
      if (arguments->par1< 0 || arguments->par1 > 100)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "par1 =%g, should be in [0, 100]",
			 arguments->par1);

	}
      break;
    case 'b':
      arguments->par2 = strtod (arg, NULL);
      if (arguments->par2< 1 || arguments->par2 > 100)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "par2 =%g, should be in [1, 100]",
			 arguments->par2);

	}
      break;
    case ARGP_KEY_ARG:
      /* Too many arguments. */
      argp_usage (state);
      break;
    default:
      return ARGP_ERR_UNKNOWN;
    }
  return 0;
}

/* Our argp parser. */
static struct argp argp = { options, parse_opt, args_doc, doc };

/*structure for events in simulation */

#define MUTATE (char) 'm'
#define INFECT (char) 'i'
#define RECOVER (char) 'r'
#define EDGE_ADD (char) 'a'
#define EDGE_DEL (char) 'd'

struct event
{
  /* key is an aggregate of event code, ego_id, and alter_id */
  /* TODO use separate structures for different event_codes */
  char event_code;
  int ego_id;
  /* The type of the last part of the key must be consistent
   * with the keylen calculation below */
  int alter_id;

  int i_deg;
  int j_deg;
  int N_ij;			/* number of edges between degree i and j nodes */
  double rate;
  int strain_id;
  int phylo_id;

  /* UT_hash_handle member makes structure a hashtable */
  UT_hash_handle hh;
};

/** genealogy struct **/
struct node
{
  int abv;
  int ndes;
  float time;
  int flag;
};

/*hash table */
struct event *event_table = NULL;

/*keylen holds the size of the hash key */
unsigned keylen;

/* sum of rate of all events in hash table */
double rate_sum;

/* array to track state of the each node */
char *node_states;

/* arrays to track states of nodes of each degree class */
int *S_k;
int *I_k;

/* maximum number of stochastic sim steps */
#define STEPMAX 50000

/*maximum degree correlations allowed */
#define MAXREDGE 0.05

/* SIR model parameters */
double trans_rate;
double recov_rate;

/* OUTSIDE is the ID of any host outside of population that transmit
 * the initial infections */
#define OUTSIDE -99

/* rounding doubles to nearest integer */
/* source http://www.cs.tut.fi/~jkorpela/round.html */
#define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))

/** function prototypes **/

int infect (int infector_id, int infectee_id, igraph_t *g,
            igraph_vector_t *degree);
int recover (int recoverer_id, igraph_t *g, igraph_vector_t *degree);
void delete_all_events ();
void delete_event (struct event *ev);
void delete_event_by_key (char event_code, int ego_id, int alter_id);
void print_rates (struct event *event_table);
int check_event_by_key (char event_code, int ego_id, int alter_id);
int check_infector_events (int infector_id, igraph_t *g);
double get_rate_sum ();
double get_infection_rate_sum ();
int neighbor_state_count (int susceptible_id, int *num_S_nb, int *num_I_nb,
                          igraph_t *g);
int get_num_infected (int max_degree);
int get_approx_case_prod_variance (double trans, double gen_time, int num_infected,
                               double *icpm, double *icpv, igraph_t *g);
double get_approx_case_prod_mean (double gen_time, int num_infected);
FILE *open_output (const char *filename);
int close_file (FILE *f);
void prtree3 (struct node *ptree, int num_nodes,
              FILE** f_tree, FILE** f_coal_times, FILE** f_samp_times);
int sequence_to_graph(int *sequence, int num_nodes, igraph_t *g);
int get_degree_seq (int *resulting_degree_seq, int num_nodes,
                    double  par1, double par2, int net_type);

     /** functions **/

void
delete_all_events ()
{
  struct event *current_event;

  while (event_table)
    {
      current_event = event_table;
      HASH_DEL (event_table, current_event);
      free (current_event);
    }
}

void
print_rates (struct event *event_table)
{
  struct event *ev;
  unsigned z = 0;

  for (ev = event_table; ev != NULL; ev = (struct event *) ev->hh.next)
    {
      printf ("rate is %g\n", ev->rate);
    }
}

void
delete_event (struct event *ev)
{
  HASH_DEL (event_table, ev);
  free (ev);
}

int
check_event_by_key (char event_code, int ego_id, int alter_id)
{

  struct event *ev1, ev2;

  memset (&ev2, 0, sizeof (struct event));
  ev2.event_code = event_code;
  ev2.ego_id = ego_id;
  ev2.alter_id = alter_id;
  HASH_FIND (hh, event_table, &ev2.event_code, keylen, ev1);
  if (!ev1)
    {
      error_at_line (0, errno, __FILE__, __LINE__,
		     "failed to find hash with key %c_%d_%d",
		     event_code, ego_id, alter_id);
    }
  else if (event_code == INFECT && ev1->rate != trans_rate)
    {
      error_at_line (0, 0, __FILE__, __LINE__,
		     "infection event with rate %g, should be %g",
		     ev1->rate, trans_rate);
    }
  else if (event_code == RECOVER && ev1->rate != recov_rate)
    {
      error_at_line (0, 0, __FILE__, __LINE__,
		     "recovery event with rate %g, should be %g",
		     ev1->rate, recov_rate);
    }
  return 0;
}

void
delete_event_by_key (char event_code, int ego_id, int alter_id)
{

  struct event *ev1, ev2;

  memset (&ev2, 0, sizeof (struct event));
  ev2.event_code = event_code;
  ev2.ego_id = ego_id;
  ev2.alter_id = alter_id;
  HASH_FIND (hh, event_table, &ev2.event_code, keylen, ev1);
  if (!ev1)
    {
      fprintf (stderr, "error: %s: %d: failed to find hash key \n",
	       __FILE__, __LINE__);
      exit (1);
    }
  HASH_DEL (event_table, ev1);
  free (ev1);
}

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

  /* index for loops */
  size_t i;

  /* value returned by functions called */
  int ret = 0;

  /* contact network */
  igraph_t g;

  /* number of hosts of degree i*/
  int N_i;

  /* vector with degree of each host  */
  igraph_vector_t degree;

  /* maximum degree of nodes in contact network */
  int max_degree;

  /* initial degree sequence of nodes in contact network */
  int *initial_degree_seq;

  /* element I of degree_dist contains the number of nodes
   * with degree I */
  gsl_histogram *degree_dist;

  /* counter of number of errors */
  errno = 0;

  /* epsilon is the initial fraction of infected nodes */
  double epsilon;

  /* I_initial holds the number of initially infected hosts in
   * based on epsilon */
  int I_initial;

  /* array holding IDS of initially infected hosts */
  int *I_initial_IDs;
  int *all_nodes_IDs;

  /* time in simulation */
  double time = 0;

  /* count of steps in simulation */
  int step_count = 0;

  double rand;
  double cum_density, p_I, p_S;

  struct arguments arguments;

/* Default values for variable set by arguments */
  arguments.silent = 0;
  arguments.verbose = 0;
  arguments.output_file = (char *) "-";
  arguments.trans_rate = 2;
  arguments.recov_rate = 1;
  arguments.interval = 0.05;
  arguments.epsilon = 0.1;
  arguments.finish_time = 5;
  arguments.network_rate = 0.4;
  arguments.network_tension = 50;
  arguments.network_size = 10000;
  arguments.network_type = DYNET_REG;
  arguments.par1 = 4;
  arguments.par2 = 1;
  arguments.sample_fraction = 0.001;

  /*TODO add option to use high tension to eliminate degree correlations
   * in initial network before beginning simulation */

  struct event *ev, *ev1, ev2;
  unsigned z;

  argp_parse (&argp, argc, argv, 0, 0, &arguments);

  fprintf (stdout, "trans_rate = %g\nrecov_rate = %g\ninterval = %g\n"
	   "epsilon = %g\nfinish_time = %g\nnetwork_rate = %g\n"
	   "network_tension = %g\nnetwork_size = %d\nnetwork_type = %d\n"
           "par1 = %g\npar2 = %g\n"
           "sample_fraction= %g\nOUTPUT_FILE = %s\nVERBOSE = %s\nSILENT = %s\n",
	   arguments.trans_rate,
	   arguments.recov_rate,
	   arguments.interval,
	   arguments.epsilon,
	   arguments.finish_time,
	   arguments.network_rate,
	   arguments.network_tension,
	   arguments.network_size,
           arguments.network_type,
           arguments.par1,
           arguments.par2,
           arguments.sample_fraction,
	   arguments.output_file,
	   arguments.verbose ? "yes" : "no", arguments.silent ? "yes" : "no");
#ifdef    DYNET_DEBUG
  fprintf (stdout, "dynet debugging level = %d\n", DYNET_DEBUG);
#endif

  /* node_hoste[I] contains the id of the genealogy node in host  with id I having the greatest time */
  int *node_host;

  struct node *ptree;

  ptree =
    (struct node *) calloc ((unsigned) (2 * arguments.network_size),
			    sizeof (struct node));
  if (!ptree)
    {
      fprintf (stderr, "Error: %s: %d: Malloc of array failed\n",
	       __FILE__, __LINE__);
      return (1);
    }

  node_host =
    (int *) calloc ((unsigned) (arguments.network_size), sizeof (int));

  if (!node_host)
    {
      fprintf (stderr, "Error: %s: %d: Malloc of array failed\n",
	       __FILE__, __LINE__);
      return (1);
    }

  int gen_ind = 0;

  trans_rate = arguments.trans_rate;
  recov_rate = arguments.recov_rate;
  epsilon = arguments.epsilon;

  keylen =
    offsetof (struct event, alter_id) + sizeof (int) - offsetof (struct event,
								 event_code);

  initial_degree_seq = (int *) malloc ((arguments.network_size) * sizeof (int));

  if (!initial_degree_seq)
    {
      fprintf (stderr, "Error: %s: %d: Malloc failed\n",
	       __FILE__, __LINE__);
      return (1);
    }

  if (arguments.network_type != DYNET_ER)
    {
  get_degree_seq (initial_degree_seq, arguments.network_size,
                 arguments.par1, arguments.par2, arguments.network_type);
  sequence_to_graph (initial_degree_seq, arguments.network_size, &g);
    }
  else
    {
      igraph_erdos_renyi_game (&g, IGRAPH_ERDOS_RENYI_GNM,
                               arguments.network_size,
                               round (arguments.network_size * arguments.par1 /2),
                               0, 0);
    }

  /* get the degree distribution */
  igraph_vector_init (&degree, 0);
  igraph_degree (&g, &degree, igraph_vss_all(), IGRAPH_ALL, IGRAPH_LOOPS);
  max_degree = (int) igraph_vector_max (&degree);
  degree_dist = gsl_histogram_alloc (max_degree + 1);
  gsl_histogram_set_ranges_uniform (degree_dist, -0.5, max_degree + 0.5);
  for (i = 0; i < igraph_vector_size (&degree); i++)
    {
      gsl_histogram_increment (degree_dist, VECTOR (degree)[i]);
    }

  // open filestreams and print headers
  FILE *fp, *of;

  if (arguments.output_file == "-")
    {
      of = stdout;
    }
  else
    {
      of = open_output (arguments.output_file);
    }
  fprintf (of, "        time   ");
  for (i = 0; i <= max_degree; i++)
    {
      N_i = round(gsl_histogram_get (degree_dist, i));
      if (N_i)
        {
          fprintf (of, "N_%05u S_%05u I_%05u   ", i + 1, i + 1, i + 1);
        }
    }
fprintf (of, "   incidence  num_infected      gen_time          icpv          icpm      vm_ratio\n");

/*  fprintf (of, "r_edge p_I p_S\n"); */

  const gsl_rng_type *T;
  gsl_rng *rng;

  /* create a generator by the environment variable
   * GSL_RNG_TYPE */

  gsl_rng_env_setup ();

  T = gsl_rng_default;
  rng = gsl_rng_alloc (T);

  double write_point = 0;

  /* put intial infection and recovery event into table */

  node_states = (char *) malloc ((arguments.network_size + 1) * sizeof (char));
  if (!node_states)
    {
      fprintf (stderr, "Error: %s: %d: Malloc of array failed\n",
	       __FILE__, __LINE__);
      return (1);
    }
  memset (node_states, 's', arguments.network_size * sizeof (char));
  node_states[arguments.network_size * sizeof (char)] = '\0';

  S_k = (int *) malloc ((max_degree + 1) * sizeof (int));
  if (!S_k)
    {
      fprintf (stderr, "Error: %s: %d: Malloc of array failed\n",
	       __FILE__, __LINE__);
      return (1);
    }
  for (i = 0; i <= max_degree; i++)
    {
      S_k[i] = gsl_histogram_get (degree_dist, i);
    }

  I_k = (int *) calloc ((max_degree + 1), sizeof (int));
  if (!I_k)
    {
      fprintf (stderr, "Error: %s: %d: Malloc of array failed\n",
	       __FILE__, __LINE__);
      return (1);
    }

  /* TODO select host at random from network */
  infect (OUTSIDE, 0, &g, &degree);
  ptree->time = 0;
  node_host[0] = 0;
  gen_ind++;

  while (time < arguments.finish_time && rate_sum > 0)
    {
      step_count++;

#if DYNET_DEBUG > 0
      /* check that rate sum is being updated correctly */

      if (fabs (get_rate_sum () - rate_sum) > 1e-6)
	{
	  error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
			 "rate_sum - get_rate_sum ()= %g\n",
			 rate_sum - get_rate_sum ());

	}

      /* check that counts of susceptibles and infecteds are
       * correct when summed over all degree classes */

      int Scount, Icount, Scount2, Icount2;
      Scount = Icount = Scount2 = Icount2 = 0;
      for (i = 0; node_states[i] != 0; i++)
	{
	  if (node_states[i] == 's')
	    Scount++;
	  else if (node_states[i] == 'i')
	    Icount++;
	}
      for (i = 0; i <= max_degree; i++)
	{
	  Scount2 += S_k[i];
	  Icount2 += I_k[i];
	}
      if (Scount != Scount2 || Icount != Icount2)
	{
	  fprintf (stderr, "Debug: %s: %d: counts don't match\n",
		   __FILE__, __LINE__);
	  fprintf (stderr, "\tSc=%d Sc2=%d, Ic=%d, Ic2=%d\n", Scount, Scount2,
		   Icount, Icount2);
	}

      /* Check that counts of Susceptiable and infecteds by degree
       * class are being updated correctly. */

      int deg, *I_k_check, *S_k_check;
      I_k_check = (int *) calloc ((max_degree + 1), sizeof (int));
      if (!I_k_check)
	{
	  error_at_line (EXIT_FAILURE, errno, __FILE__, __LINE__,
			 "I_k_check");
	}
      S_k_check = (int *) calloc ((max_degree + 1), sizeof (int));
      if (!S_k_check)
	{
	  error_at_line (EXIT_FAILURE, errno, __FILE__, __LINE__,
			 "S_k_check");
	}
      for (i = 0; i < arguments.network_size; i++)
	{
	  if (node_states[i] == 's')
	    {
	      deg = VECTOR (degree)[i];
	      S_k_check[deg]++;
	    }
	  if (node_states[i] == 'i')
	    {
	      deg = VECTOR (degree)[i];
	      I_k_check[deg]++;
	    }
	}
      error_message_count = 0;
      for (i = 0; i <= max_degree; i++)
	{
	  if (S_k_check[i] != S_k[i])
	    {
	      error_at_line (0, 0, __FILE__, __LINE__,
			     "S_k[%d] = %d, should be %d",
			     i, S_k[i], S_k_check[i]);
	    }
	  if (I_k_check[i] != I_k[i])
	    {
	      error_at_line (0, 0, __FILE__, __LINE__,
			     "I_k[%d] = %d, should be %d",
			     i, I_k[i], I_k_check[i]);
	    }
	}
      if (error_message_count != 0)
	{
	  error (EXIT_FAILURE, 0, "%u errors found", error_message_count);
	}
      free (I_k_check);
      free (S_k_check);

      /* calculate p.I, i.e. fraction of arcs from susctpibles
       * going to infecteds */
      int num_S_nb, num_I_nb, sum_S_nb, sum_I_nb, sum_nb, susceptible_id;
      sum_nb = sum_S_nb = sum_I_nb = 0;
      for (i = 0; i < arguments.network_size; i++)
	{
	  if (node_states[i] == 's')
	    {
	      num_S_nb = num_I_nb = 0;
	      sum_nb += VECTOR (degree)[i];
	      susceptible_id = i;
	      neighbor_state_count (susceptible_id, &num_S_nb, &num_I_nb, &g);
	      sum_S_nb += num_S_nb;
	      sum_I_nb += num_I_nb;
	    }
	}
      p_S = (double) sum_S_nb / sum_nb;
      p_I = (double) sum_I_nb / sum_nb;
#endif
#if DYNET_DEBUG > 1

      /* Check for all events in the hash table, making sure the
       * hash table is consistent with the current state of the network. */
      error_message_count = 0;
      int infector_id;
      int event_count = 0;

      /*first check for and count one edge addition or deletion event that
       * can always happen */

      memset (&ev2, 0, sizeof (struct event));
      ev2.event_code = EDGE_ADD;
      HASH_FIND (hh, event_table, &ev2.event_code, keylen, ev1);
      if (ev1)
	{
	  event_count++;
	}
      else
	{
	  ev2.event_code = EDGE_DEL;
	  HASH_FIND (hh, event_table, &ev2.event_code, keylen, ev1);
	  if (ev1)
	    {
	      event_count++;
	    }
	}

      for (i = 0; i < arguments.network_size; i++)
	{
	  if (node_states[i] == 'i')
	    {
	      infector_id = i;
	      event_count += check_infector_events (infector_id, &g);
	    }
	}
      int hash_count = HASH_COUNT (event_table);
      if (error_message_count != 0)
	{
	  error (EXIT_FAILURE, 0, "%u events are missing in hash table "
		 "or are present but occuring at the wrong rate",
		 error_message_count);
	}
      if (hash_count != event_count)
	{
	  error (EXIT_FAILURE, 0, "%d extra events in hash table\n",
		 hash_count - event_count);
	}

#endif
#if DYNET_DEBUG_TEST == 1
      if (step_count == 10)
	{
	  /* set the rate of the first event in the table to 0 */
	  fprintf (stderr,
		   "Changing the rate of an event to test debugging tests\n");
	  event_table->rate = 0;
	}
#endif
#if DYNET_DEBUG_TEST == 2
      if (step_count == 70)
	{
	  /* add an extra event to the table */
	  fprintf (stderr, "Adding an event to test debugging tests\n");
	  ev1 = (struct event *) malloc (sizeof (struct event));
	  if (!ev1)
	    {
	      fprintf (stderr,
		       "Error: %s: %d: malloc failed\n", __FILE__, __LINE__);
	      return (1);
	    }
	  memset (ev1, 0, sizeof (struct event));
	  ev1->event_code = INFECT;
	  ev1->rate = 0;
	  HASH_ADD (hh, event_table, event_code, keylen, ev1);
	}
#endif
#if DYNET_DEBUG_TEST == 3
      if (step_count == 50)
	{
	  /* remove an event at beginning of the table */
	  fprintf (stderr, "Removing an event to test debugging tests\n");
	  delete_event (event_table);
	}
#endif
#if DYNET_DEBUG_TEST == 4
      if (step_count == 50)
	{
	  /* give two events rates with compensating erros */
	  fprintf (stderr, "Swapping event rates to test debugging tests\n");
	  ev = event_table;
	  ev->rate += 0.001;
	  ev = (struct event *) ev->hh.next;
	  ev->rate -= 0.001;
	}
#endif
#if DYNET_DEBUG_TEST == 5
      if (step_count == 50)
	{
	  /* change event code  */
	  fprintf (stderr,
		   "Changing event code of an event to test debugging tests\n");
	  ev = event_table;
	  ev->event_code = 'Z';
	}
#endif

      if (step_count > STEPMAX)
	{
	  fprintf (stderr, "Warning: %s: %d: reached STEPMAX, breaking\n",
		   __FILE__, __LINE__);
	  break;
	}

      /* determine time of next event */
      time += gsl_ran_exponential (rng, 1 / rate_sum);
      if (time > write_point)
	{
	  fprintf (of, "%9e   ", write_point);
	  write_point += arguments.interval;
	  for (i = 0; i <= max_degree; i++)
	    {
              int N_i = round(gsl_histogram_get (degree_dist, i));
              if (N_i)
                {
                  fprintf (of, "%7ld %7d %7d   ",  N_i, S_k[i], I_k[i]);
                }
	    }
          double incidence = get_infection_rate_sum();
          int num_infected = get_num_infected(max_degree);
          double gen_time = (double) num_infected / incidence;
/*          double icpm = get_approx_case_prod_mean (gen_time, num_infected);*/
          double icpm;
          double icpv;
          get_approx_case_prod_variance (arguments.trans_rate,
                                         gen_time,
                                         num_infected,
                                         &icpm, &icpv, &g);
          double vm_ratio = icpv / icpm;

          fprintf (of, "%9e  %12d  %9e  %9e  %9e  %9e\n",
                   incidence, num_infected, gen_time, icpv, icpm, vm_ratio );
/*	  fprintf (of, " %g %g\n", p_I, p_S); */
	}

      /* select one event to be next proportionally to it's rate */

      rand = rate_sum * gsl_rng_uniform (rng);
      cum_density = 0;
      for (ev = event_table; ev != NULL; ev = (struct event *) ev->hh.next)
	{
	  cum_density += ev->rate;
	  if (cum_density > rand)
	    {
	      break;
	    }
	}

      switch (ev->event_code)
	{
	case RECOVER:
	  (ptree + gen_ind)->time = time;
	  (ptree + gen_ind)->abv = node_host[ev->ego_id];
	  node_host[ev->ego_id] = gen_ind;
	  gen_ind++;
	  recover (ev->ego_id, &g, &degree);
	  break;
	case INFECT:
	  (ptree + gen_ind)->time = time;
	  if (gen_ind > 0)
	    {
	      (ptree + gen_ind)->abv = node_host[ev->ego_id];
	    }
	  node_host[ev->ego_id] = gen_ind;
	  node_host[ev->alter_id] = gen_ind;
	  gen_ind++;
	  infect (ev->ego_id, ev->alter_id, &g, &degree);
	  break;
	default:
	  fprintf (stderr, "Error: %s: %d: Illegal event code\n",
		   __FILE__, __LINE__);
	  ret = 1;
	}
      if (ret == 1)
	{
	  break;
	}
    }				/* end while() */
  delete_all_events ();
  close_file (of);
  free (node_states);
  free (S_k);
  free (I_k);

  FILE *f_tree = open_output ("tree");
  FILE *f_coal_times = open_output ("coal_times");
  FILE *f_samp_times = open_output ("samp_times");

  int *sampled_pnode_IDs, *all_pnode_IDs;
  int sample_size = round (arguments.sample_fraction * gen_ind);
  if (sample_size == 0)
    {
      sample_size = 1;
    }
  sampled_pnode_IDs = (int *) malloc (sample_size * sizeof (int));
  all_pnode_IDs = (int *) malloc (gen_ind * sizeof (int));
  if (!sampled_pnode_IDs || !all_pnode_IDs)
    {
      fprintf (stderr, "Error: %s: %d: Malloc of array failed\n",
	       __FILE__, __LINE__);
      return (1);
    }
  for (i = 0; i < gen_ind; i++)
    {
      all_pnode_IDs[i] = i;
      (ptree + i)->flag = 0;
    }
  gsl_ran_choose (rng, sampled_pnode_IDs, sample_size, all_pnode_IDs, gen_ind,
		  sizeof (int));
  for (i = 0; i < sample_size; i++)
    {
      int w = sampled_pnode_IDs[i];
      (ptree + w)->flag = 1;
      while (w != 0)
	{
	  w = (ptree + w)->abv;
	  (ptree + w)->flag = 1;
	}
    }
  prtree3 (ptree, gen_ind, &f_tree, &f_coal_times, &f_samp_times);
  gsl_rng_free (rng);
  free (all_pnode_IDs);
  free (sampled_pnode_IDs);

  free (ptree);
  free (initial_degree_seq);
  gsl_histogram_free (degree_dist);
  igraph_destroy (&g);
  close_file (f_tree);
  close_file (f_coal_times);
  close_file (f_samp_times);
  return 0;
}

int
infect (int infector_id, int infectee_id, igraph_t *g,
        igraph_vector_t *degree)
{
  igraph_vector_t neighbors;
  int neighbor_id;
  int infectee_id_deg;
  struct event *ev1;
  size_t i;

  /* update system dynamic variable */
  node_states[infectee_id] = 'i';
  infectee_id_deg = round (igraph_vector_e (degree, infectee_id));
  S_k[infectee_id_deg]--;
  I_k[infectee_id_deg]++;

  /* add event of infectee recovering from the table */
  ev1 = (struct event *) malloc (sizeof (struct event));
  if (!ev1)
    {
      fprintf (stderr, "Error: %s: %d: malloc failed\n", __FILE__, __LINE__);
      return (1);
    }
  memset (ev1, 0, sizeof (struct event));
  ev1->ego_id = infectee_id;
  ev1->rate = recov_rate;
  ev1->event_code = RECOVER;
  HASH_ADD (hh, event_table, event_code, keylen, ev1);
  rate_sum += recov_rate;

  /* add events of infectee infecting susceptible neighbors and
   * remove events of infected host getting infeced by infectious
   * neighbors */

  igraph_vector_init (&neighbors, 0);
  igraph_neighbors (g, &neighbors, (igraph_integer_t) infectee_id, IGRAPH_ALL);
  for (i = 0; i < igraph_vector_size (&neighbors); i++)
    {
      neighbor_id = VECTOR (neighbors)[i];
      switch (node_states[neighbor_id])
	{
	case 's':
	  ev1 = (struct event *) malloc (sizeof (struct event));
	  if (!ev1)
	    {
	      fprintf (stderr, "Error: %s: %d: malloc failed\n", __FILE__,
		       __LINE__);
	      return (1);
	    }
	  memset (ev1, 0, sizeof (struct event));
	  ev1->ego_id = infectee_id;
	  ev1->alter_id = neighbor_id;
	  ev1->rate = trans_rate;
	  ev1->event_code = INFECT;
	  HASH_ADD (hh, event_table, event_code, keylen, ev1);
	  rate_sum += trans_rate;
	  break;
	case 'i':
	  delete_event_by_key (INFECT, neighbor_id, infectee_id);
	  rate_sum -= trans_rate;
	  break;
	case 'r':
	  break;
	default:
	  fprintf (stderr,
		   "Error: %s: %d: invalid state for node\n",
		   __FILE__, __LINE__);
	  exit (1);
	}
    }
  igraph_vector_destroy (&neighbors);
  return 0;
}

int
neighbor_state_count (int susceptible_id, int *num_S_nb, int *num_I_nb,
                      igraph_t *g)
{

  int neighbor_id;
  igraph_vector_t neighbors;
  size_t i;

  *num_S_nb = 0;
  *num_I_nb = 0;

  /* check for event of infector infecting neighbors */

  igraph_vector_init (&neighbors, 0);
  igraph_neighbors (g, &neighbors, susceptible_id, IGRAPH_ALL);

  for (i = 0; i < igraph_vector_size (&neighbors); i++)
    {
      neighbor_id = VECTOR (neighbors)[i];
      switch (node_states[neighbor_id])
	{
	case 's':
	  ++*num_S_nb;
	  break;
	case 'i':
	  ++*num_I_nb;
	  break;
	case 'r':
	  break;
	default:
	  fprintf (stderr,
		   "Error: %s: %d: invalid state for node\n",
		   __FILE__, __LINE__);
	  exit (1);
	}
    }
  igraph_vector_destroy (&neighbors);
  return 0;

}

int
check_infector_events (int infector_id, igraph_t *g)
{

  int neighbor_id, event_count = 0;
  struct event *ev1;
  igraph_vector_t neighbors;
  size_t i;

  /* check for recovery event */

  check_event_by_key (RECOVER, infector_id, 0);
  event_count++;

  /* check for event of infector infecting neighbors */

  igraph_vector_init (&neighbors, 0);
  igraph_neighbors (g, &neighbors, infector_id, IGRAPH_ALL);

  for (i = 0; i < igraph_vector_size (&neighbors); i++)
    {
      neighbor_id = VECTOR (neighbors)[i];
      switch (node_states[neighbor_id])
	{
	case 's':
	  check_event_by_key (INFECT, infector_id, neighbor_id);
	  event_count++;
	  break;
	case 'i':
	  break;
	case 'r':
	  break;
	default:
	  fprintf (stderr,
		   "Error: %s: %d: invalid state for node\n",
		   __FILE__, __LINE__);
	  exit (1);
	}
    }
  igraph_vector_destroy (&neighbors);
  return event_count;
}

int
recover (int recoverer_id, igraph_t *g, igraph_vector_t *degree)
{

  int recoverer_id_deg, neighbor_id;
  struct event *ev1;
  igraph_vector_t neighbors;
  size_t i;

  /* update system dynamic variable */

  node_states[recoverer_id] = 'r';
  recoverer_id_deg = round (igraph_vector_e (degree, recoverer_id));
  I_k[recoverer_id_deg]--;

  /* remove event of recoverer recovering again */

  delete_event_by_key (RECOVER, recoverer_id, 0);
  rate_sum -= recov_rate;

  /* remove event of recoverer infecting neighbors */

  igraph_vector_init (&neighbors, 0);
  igraph_neighbors (g, &neighbors, recoverer_id, IGRAPH_ALL);
  for (i = 0; i < igraph_vector_size (&neighbors); i++)
    {
      neighbor_id = VECTOR (neighbors)[i];
      switch (node_states[neighbor_id])
	{
	case 's':
	  delete_event_by_key (INFECT, recoverer_id, neighbor_id);
	  rate_sum -= trans_rate;
	  break;
	case 'i':
	  break;
	case 'r':
	  break;
	default:
	  fprintf (stderr,
		   "Error: %s: %d: invalid state for node\n",
		   __FILE__, __LINE__);
	  exit (1);
	}
    }
  igraph_vector_destroy (&neighbors);
  return 0;
}

double
get_infection_rate_sum ()
{
  struct event *ev;
  double infection_rate_sum = 0;

  for (ev = event_table; ev != NULL; ev = (struct event *) ev->hh.next)
    {
      if ( ev->event_code == INFECT)
        {
          infection_rate_sum += ev->rate;
        }
    }
  return infection_rate_sum;
}

int
get_num_infected (int max_degree)
{
  int i;
  int num_infected = 0;

  for (i = 0; i <= max_degree; i++)
    {
      num_infected += I_k[i];
    }
  return num_infected;
}

double
get_approx_case_prod_mean (double gen_time, int num_infected)
{
  struct event *ev;
  double mean_sum = 0;

  /* probability of transmission over edge within GEN_TIME */
  double p;

  /* individual case production mean */
  double icpm;

  for (ev = event_table; ev != NULL; ev = (struct event *) ev->hh.next)
    {
      if ( ev->event_code == INFECT)
        {
          p = 1 - exp (-ev->rate * gen_time);
          mean_sum += p;
        }
    }
  icpm = mean_sum / num_infected;
  return icpm;
}

int
get_approx_case_prod_variance (double trans, double gen_time, int num_infected,
                               double *icpm, double *icpv, igraph_t *g)
{
  int num_S_nb;
  int num_I_nb;
  int i;

  double p;
  double np;

  double E_y = 0;
  double E_y_sq = 0;

  p = 1 - exp (-trans * gen_time );

  for (i = 0; node_states[i] != 0; i++)
    {
      if (node_states[i] == 'i')
        {
          neighbor_state_count (i, &num_S_nb, &num_I_nb, g);
          np = num_S_nb * p;
          E_y += np;
          E_y_sq += np * (np - p + 1);
        }
    }
  *icpm = (double) E_y / num_infected;
  *icpv = (double) E_y_sq / num_infected  - (*icpm) * (*icpm);
  return 0;
}

double
get_rate_sum ()
{
  struct event *ev;
  double rate_sum_check = 0;

  for (ev = event_table; ev != NULL; ev = (struct event *) ev->hh.next)
    {
      rate_sum_check += ev->rate;
    }
  return rate_sum_check;
}

void
prtree3 (struct node *ptree, int num_nodes,
         FILE** f_tree, FILE** f_coal_times, FILE** f_samp_times)
  /* print genealogy tree in Newick format */
  /* after Hudson's ms code */
{
  double t;
  int i, *descl, *descr, w, tmp;
  void parens (struct node *ptree, int *descl, int *descr, int noden,
               FILE** f_tree, FILE** f_coal_times, FILE** f_samp_times);

  descl = (int *) malloc ((unsigned) (num_nodes) * sizeof (int));
  descr = (int *) malloc ((unsigned) (num_nodes) * sizeof (int));
  for (i = 0; i < num_nodes; i++)
    descl[i] = descr[i] = -1;
  for (i = 1; i < num_nodes; i++)
    {
      if (descl[(ptree + i)->abv] == -1 && (ptree + i)->flag == 1)
	descl[(ptree + i)->abv] = i;
      else if ((ptree + i)->flag == 1)
	descr[(ptree + i)->abv] = i;
    }
  for (i = 1; i < num_nodes; i++)
    {
      if (descl[i] != -1 && descr[i] == -1)
	{
	  w = descl[i];
	  while (descr[w] == -1)
	    {
	      tmp = descl[w];
	      descl[w] = -1;
	      if (tmp == -1)
		{
		  break;
		}
	      w = tmp;
	    }
	  if (tmp != -1)
	    {
	      descl[i] = descl[w];
	      descr[i] = descr[w];
              (ptree + descl[i])->abv = i;
              (ptree + descr[i])->abv = i;
              (ptree + i)->time = (ptree + w)->time;
	      descl[w] = -1;
	      descr[w] = -1;
	    }
	  if (tmp == -1)
	    {
	      descl[i] = -1;
	      tmp = 0;
	      while (descl[tmp] != i && descr[tmp] != i)
		{
		  tmp++;
		}
	      if (descl[tmp] == i)
		{
		  descl[tmp] = w;
                  (ptree + w)->abv = tmp;
		}
	      if (descr[tmp] == i)
		{
		  descr[tmp] = w;
                  (ptree + w)->abv = tmp;
		}
	    }
	}
    }
  parens (ptree, descl, descr, 1, f_tree, f_coal_times, f_samp_times);
  fprintf (*f_coal_times, "\n");
  fprintf (*f_samp_times, "\n");

  free (descl);
  free (descr);
}

void
parens (struct node *ptree, int *descl, int *descr, int noden,
        FILE** f_tree, FILE** f_coal_times, FILE** f_samp_times)
  /* print genealogy tree in Newick format */
  /* after Hudson's ms code */
{
  double time;

  if (descl[noden] == -1)
    {
      time = (ptree + (ptree + noden)->abv)->time - (ptree + noden)->time;
      fprintf (*f_tree, "%d_%g:%g", noden + 1, (ptree + noden)->time, -1*time);
      fprintf (*f_samp_times, "%g ", (ptree + noden)->time);
    }
  else
    {
      fprintf (*f_coal_times, "%g ", (ptree + noden)->time);
      fprintf (*f_tree, "(");
      parens (ptree, descl, descr, descl[noden], f_tree,
              f_coal_times, f_samp_times);
      fprintf (*f_tree, ",");
      parens (ptree, descl, descr, descr[noden], f_tree,
              f_coal_times, f_samp_times);
      if ((ptree + noden)->abv == 0)
	{
	  fprintf (*f_tree, ");\n");
	}
      else
	{
	  time = (ptree + (ptree + noden)->abv)->time - (ptree + noden)->time;
	  fprintf (*f_tree, "):%g", -1*time);
	}
    }
}

int
get_degree_seq (int *resulting_degree_seq, int num_nodes,
                    double  par1, double par2, int net_type)
{
  int i, degree_sum;
  double mean_degree;
  const gsl_rng_type * T;
  gsl_rng *r;

  if (net_type == DYNET_PARETO)
    {
      gsl_rng_env_setup ();
      T = gsl_rng_default;
      r = gsl_rng_alloc (T);
      degree_sum = 0;
      for (i = 0; i < num_nodes; i++)
      {
        resulting_degree_seq[i] = gsl_ran_pareto (r, par1, par2);
        degree_sum += resulting_degree_seq[i];
      }
      if (degree_sum % 2 != 0)
        {
          fprintf (stderr,
                   "Warning: adjusting degree sequence so that it's even\n");
          resulting_degree_seq[i-1]++;
        }
      gsl_rng_free (r);
    }
  if (net_type == DYNET_EXP)
    {
      mean_degree = par1;
      gsl_rng_env_setup ();
      T = gsl_rng_default;
      r = gsl_rng_alloc (T);
      degree_sum = 0;
      for (i = 0; i < num_nodes; i++)
      {
        resulting_degree_seq[i] = round(gsl_ran_exponential (r, mean_degree));
        if (resulting_degree_seq[i] == 0)
          {
            resulting_degree_seq[i]++;
          }
        degree_sum += resulting_degree_seq[i];
      }
      if (degree_sum % 2 != 0)
        {
          fprintf (stderr,
                   "Warning: adjusting degree sequence so that it's even\n");
          resulting_degree_seq[i-1]++;
        }
      gsl_rng_free (r);
    }
  if (net_type == DYNET_REG)
    {
      mean_degree = par1;
      if (mean_degree - (int) mean_degree > 0 || (int) mean_degree *num_nodes % 2 !=0)
        {
          fprintf (stderr, "Error: %s: %d:\
 mean degree impossible for regular net\n", __FILE__, __LINE__);
          return 1;
        }
      for (i = 0; i <num_nodes; i++)
        {
          resulting_degree_seq[i] = (int) mean_degree;
        }
    }
  return 0;
}

int
sequence_to_graph(int *sequence, int num_nodes, igraph_t* g)
{

/* construct a graph from a degree sequence*/

  igraph_vector_t v;
  size_t i;

  igraph_vector_init(&v, 0);
  igraph_vector_resize(&v, num_nodes);

  for (i=0; i<num_nodes; i++)
    {
      VECTOR(v)[i]=sequence[i];
    }
  igraph_degree_sequence_game(g, &v, NULL, IGRAPH_DEGSEQ_VL);
/* TODO
 * figure out how to use these function so that initial graph can have multiple components and isolates
  igraph_degree_sequence_game(g, &v, NULL, IGRAPH_DEGSEQ_SIMPLE);
  igraph_simplify (g, 1, 1);
  */
  igraph_vector_destroy(&v);

  return 0;
}

FILE *
open_output (const char *filename)
     /* open filename for output; return NULL if problem */
     /* after p. 367, C A Reference Manual 5e */
{
  FILE *f;
  errno = 0;
  /* Functions below might choke on a NULL filename. */
  if (filename == NULL) filename = "";
  f = fopen (filename, "w");
  if (f == NULL)
    {
      error_at_line (EXIT_FAILURE, errno, __FILE__, __LINE__,
                     "Failed to open file %s", filename);
    }
  return f;
}

int
close_file (FILE *f)
     /* Close file f */
     /* after p. 367, C A Reference Manual 5e */
{
  int s = 0;
  if (f == NULL)
    {
      return 0;
    }
  errno = 0;
  s = fclose (f);
  if (s == EOF) perror ("Close failed");
  return s;
}

ext: c


text:
/*
VEGAS_TABLE Version 1.  5/27/16
report issues to Tim Waters,
waters@lanl.gov

**** Basic code description:
Input: 	An ascii table with entries
		on a uniform grid of (x,y) values
		and a (x,y) value for which an
    interpolated value is sought
Output: The interpolated value
     corresponding to that (x,y) pair

**** Detailed description:

>> Step 1. Instantiation:
The class constructor calls a method
to parse an ascii text file which contains
an (N+1)x(M+1) table:
The first row contains the M values of x.
The first column contains the N values of y.
The remainder of the table contains the NxM
values of z.

Here is a simple example of an input file:

4	1e0     1e1     1e2     1e3
1e4	1e-20	2e-20	3e-20	4e-20
1e5	5e-20	6e-20	7e-20	8e-20
1e6	1e-19	2e-19	3e-19	4e-19
1e7	5e-19	6e-19	7e-19	8e-19
1e8	1e-18	2e-18	3e-18	4e-18

In this example, N = 5 and M = 4.
M must be specified as the first element.
The remaining elements in row 1 are values
of the photoionization parameter, xi.
The left column below M are the temperature
values.  The remaining elements are net heating
rates that are functions of T and xi.

By default, the input file is limited to size
1000 x 1000.  To increase this size limit,
increase Nmax in the default constructor.

>> Step 2. Table initialization:
The file parsed above is read into a
table, but this is not the actual table
used.  The actual table is made upon calling
an initialization function, and optionally
array arguments can be passed specifying the
min/max of both T and xi. This may be useful
when combining multiple overlapping tables
and will allow for flexibility if tables
generated from non-uniform grids are needed,
although a table-search algorithm has not
been implemented.  The input table is then
reduced to a size that contains only the
values between these limits.

>> Step 3. Invocation:
The code returns the interpolated value
of whatever rate was in the input file.
The code uses the bicubic and bilinear
interpolation algorithms from the GSL
library.  By default, bicubic is used.

Importantly, the code assumes that a
nearly uniform grid was produced via
x = j*(xmax-xmin)/(M-1) + xmin
y = i*(ymax-ymin)/(N-1) + ymin
This allows avoiding an expensive
table search for the lookup operation.
For a uniform grid generated from the
above formulas, the table indices
for a desired value of (x,y) are
i_x = (M-1)*(x - xmin)/(xmax - xmin)
i_y = (N-1)*(y - ymin)/(ymax - ymin)
Since xstar may slightly modify the
grid locations (x,y) originally
specified, our lookup operation
consists of first guessing that the
lookup value indices are (i_x,i_y)
and then allowing for some wiggle
room by searching 1 index left, right,
above, and below upon finding that our
guess was incorrect.  The code will
terminate if this search fails.

**** Example usages:
1. Bicubic interpolation on a table that was
generated from a log-log spaced grid
VEGAS_LUT hc_xstar("input_table.tab","bicubic");
hc_xstar.initialize_table();
hc_rate = hc_xstar.get_rate(T,xi);

1a. If lin-lin spacing was used instead, use
VEGAS_LUT hc_xstar("input_table.tab","bicubic",false,false);

1c. To apply a bounding box, use instead
const double T_bbox[2] = {1e5,1e7};
const double xi_bbox[2] = {1e1,5e2};
hc_xstar.initialize_table(T_bbox,xi_bbox);
*/

#ifndef H_VEGAS_TABLE
#define H_VEGAS_TABLE

#include <iostream>		// std::cout, std::endl
#include <iomanip>      // std::setprecision
#include <fstream>		// std::ifstream
#include <string>
#include <vector>
#include <time.h>
#include <cmath>

#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

class VEGAS_LUT
{
public:
    //default constructor will parse the input file
    //log-log spacing for (T,xi) is assumed by default
	VEGAS_LUT(	const std::string a="VEGAS_table.tab",
				const std::string b="bicubic",
				bool c=true, bool d=true):
	filename(a),gsl_routine(b),logT(c),logxi(d)
	{
	 // allocate a (Nmax x Nmax) table for reading in the file
	  int Nmax = 1000;
	  input_table = new double*[Nmax];
	  for (int i = 0; i < Nmax; i++)
    	input_table[i] = new double[Nmax];
      xivals = new double[Nmax];
      Tvals = new double[Nmax];

      gsl_init();
      parse_file();
	}
	~VEGAS_LUT()  { free_memory(); }

	//user must call the initialize_table method
	void initialize_table(const double*, const double*);

	//user function that returns an interpolated value
	double get_rate(double, double);

	// actual table (T,xi) bounds:
	// hydro code needs to know Temp bounds when doing implicit root finding
	double y_min,y_max,x_min,x_max;

private:
	bool logT,logxi;					// determines log gridding
	double* xlim;						// input table xi bounds
	double* ylim;						// input table T bounds
	const double* T_bounds;				// user specified T bounds
	const double* xi_bounds;			// user specified xi bounds
	double** 	input_table; 			// pointer to input table
	double** 	table;					// pointer to actual table
	double* 	xivals; 				// xi-vals from input table
	double*		Tvals; 					// T-vals from input table
	double* 	xvals; 					// xi-vals for actual table
	double*		yvals; 					// T-vals for actual table
	unsigned int 	NT,	  // # rows of input table (# of Ts)
					Nxi,  // # cols of input table (# of xis)
					N,	  // # rows of actual table (# of Ts)
					M,	  // # cols of actual table (# of xis)
					i_T,  // T-offset between the input and actual tables
					i_xi; // xi-offset between the input and actual tables
	typedef std::vector<double> array_t;	//dynamic array data type
	std::string filename, gsl_routine;		//input strings

	//member functions for handling the input file
	void parse_file();
	bool read_value(std::ifstream&, const int, array_t&);
	void free_input_table();
	void free_memory();

	//other member functions
	void reduce_input_table();
	void make_table(unsigned int, unsigned int);
	void determine_box_indices(double, double, unsigned int&, unsigned int&);

	//GSL variables and member functions
	bool use_gsl_bicubic;
	const gsl_interp2d_type *gsl_scheme;
	typedef gsl_spline2d* spline_t;
	spline_t **spline_table;
    gsl_interp_accel *xacc, *yacc;
    size_t nx,ny;
    void gsl_init();
    void make_spline_table();
};

void VEGAS_LUT::gsl_init()
{
    if (gsl_routine == "bilinear")
    {
      gsl_scheme = gsl_interp2d_bilinear;
      use_gsl_bicubic = false;
      nx = 2;
      ny = 2;
    }
    else if (gsl_routine == "bicubic")
    {
      gsl_scheme = gsl_interp2d_bicubic;
      use_gsl_bicubic = true;
      nx = 4;
      ny = 4;
    }
    else
    {
      std::cout << "ERROR: GSL interpolation scheme " << gsl_routine
      << " not recognized. Choose from:\n"
      << "bilinear\n"
      << "bicubic\n"
      << std::endl;
      exit(0);
    }
}

bool VEGAS_LUT::read_value(std::ifstream& fin, const int j, array_t& input_data)
{
  double value;
  if (fin >> value)
  {
    input_data[j]=value;
    return true;
  }
  else return false;
}

void VEGAS_LUT::parse_file()
{
  std::ifstream fin(filename.c_str());
  if (!fin)
  {
    std::cout << "File " << filename << " does not exist!" << std::endl;
	exit(0);
  }
  bool file_is_open = false;
  bool end_of_line = false;

  // read in first value: Nxi
  array_t val1(1);
  file_is_open = read_value(fin,0,val1);
  Nxi = (unsigned int) val1[0];

  // read in remainder of first row: all xi-values
  array_t this_row(Nxi+1);
  int i=1,j;
  while(!end_of_line)
  {
    file_is_open = read_value(fin,i,this_row);
    xivals[i-1] = this_row[i];
    if( i!=0 && i%Nxi == 0 )
        end_of_line = true;
    else i++; //next row element
  }

  // read in 1st element of 2nd row: the 1st T value
  file_is_open = read_value(fin,0,this_row);
  Tvals[0] = this_row[0];

  // read in the rest of the file
  end_of_line = false;
  i = 1; j=0;
  while(file_is_open)
  {
    file_is_open = read_value(fin,i,this_row);

    if( i!=0 && i%Nxi == 0 )
    {
        end_of_line = true;
        i=0; // j++; //reset column, increment row
    }
    else i++; //next row element

    if(end_of_line)
    {
      // 1st element of this_row is T
      Tvals[j] = this_row[0];

      // remaining elements are rates
      for (int ii=0; ii<Nxi; ii++)
        //std::cout << "this_row[" << ii << "] = " << this_row[ii] << std::endl;
        input_table[j][ii] = this_row[ii+1];

      j++;
    }
    end_of_line = false;
  }

  NT = j; //number of temperature values
  fin.close();

  if (Nxi < 2 || NT < 2)
    {
    std::cout << "ERROR: Table dimensions must be at least 2x2!" << std::endl;
    exit(0);
    }
}

void VEGAS_LUT::make_table(unsigned int i0, unsigned int j0)
{

    // allocate memory
	table = new double*[N];
	for (int j = 0; j < N; j++)
    	table[j] = new double[M];

    // assign values
	for (int j = 0; j < N; j++)
	  for (int i = 0; i < M; i++)
	    table[j][i] = input_table[j0+j][i0+i];
}

void VEGAS_LUT::make_spline_table()
{

    // allocate memory for spline_table
	spline_table = new spline_t*[N];
	for (int j = 0; j < N; j++)
    	spline_table[j] = new spline_t[M];

    /* Each location in our table is to be thought of as a 4-corner box.
     * The corners have locations (xb[], yb[]) and function values (zb[])
     * and this data is collected ahead of time so that GSL only needs to
     * access it to carry out the interpolation.  For the case of bicubic
     * interpolation, GSL will also take numerical derivatives and store
     * all that data, so a lot of expense is saved by collecting this at
     * initialization. */
    double *xb = new double[nx];
	double *yb = new double[ny];
	double *zb = new double[nx*ny];

	/* Associations:
	T  = y(j), with j from (0,N-1)
	xi = x(i), with i from (0,M-1)
	*/
	for (int j=1; j < (N-2); j++)
      for (int i=1; i < (M-2); i++)
      {
      	// allocate storage for interpolation data
        gsl_spline2d *spline = gsl_spline2d_alloc(gsl_scheme, nx, ny);

        // evaluate position data (xb[],yb[]) and store function values (zb[]) to spline
	    for (int jb=0; jb < nx; jb++)
		  for (int ib=0; ib < ny; ib++)
		  {
		    if (use_gsl_bicubic)
		    {
		      xb[ib] = xvals[i+ib-1];
		      yb[jb] = yvals[j+jb-1];
			  gsl_spline2d_set(spline, zb, ib, jb, table[j+jb-1][i+ib-1]);
		    }
		    else
		    {
		      xb[ib] = xvals[i+ib];
		      yb[jb] = yvals[j+jb];
		      gsl_spline2d_set(spline, zb, ib, jb, table[j+jb][i+ib]);
		    }
		  }

          /* add position data
          * In the case of bicubic, this also internally calculates
          * derivative data and stores that as well. */
          gsl_spline2d_init(spline, xb, yb, zb, nx, ny);

          // store this interpolation data in our spline table
          spline_table[j][i] = spline;
      } //end loop

    /* The above loop excludes all boxes on table edges, which always require
     * doing bilinear interpolation, so add that data now */
    int j_edges[2] = {0, N-2};
    int i_edges[2] = {0, M-2};
    int i,j;

    for (int jj=0; jj < 2; jj++)
      for (int i=0; i < (M-1); i++)
      {
        j = j_edges[jj];

		gsl_spline2d *spline = gsl_spline2d_alloc(gsl_interp2d_bilinear, 2, 2);

		// evaluate position data (xb[],yb[]) and store function values (zb[]) to spline
		  for (int jb=0; jb < 2; jb++)
			for (int ib=0; ib < 2; ib++)
			{
				xb[ib] = xvals[i+ib];
				yb[jb] = yvals[j+jb];
				gsl_spline2d_set(spline, zb, ib, jb, table[j+jb][i+ib]);
			}

		/* add position data
        * In the case of bicubic, this also internally calculates
        * derivative data and stores that as well. */
        gsl_spline2d_init(spline, xb, yb, zb, 2, 2);

        // store this interpolation data in our spline table
        //std::cout<<  "(i,j) = (" << i << ", " << j << ")" << std::endl;
        spline_table[j][i] = spline;
      }

    for (int j=0; j < (N-1); j++)
      for (int ii=0; ii < 2; ii++)
      {
        i = i_edges[ii];

		gsl_spline2d *spline = gsl_spline2d_alloc(gsl_interp2d_bilinear, 2, 2);

		// evaluate position data (xb[],yb[]) and store function values (zb[]) to spline
		  for (int jb=0; jb < 2; jb++)
			for (int ib=0; ib < 2; ib++)
			{
				xb[ib] = xvals[i+ib];
				yb[jb] = yvals[j+jb];
				gsl_spline2d_set(spline, zb, ib, jb, table[j+jb][i+ib]);
			}

		/* add position data
        * In the case of bicubic, this also internally calculates
        * derivative data and stores that as well. */
        gsl_spline2d_init(spline, xb, yb, zb, 2, 2);

        // store this interpolation data in our spline table
        //std::cout<<  "(i,j) = (" << i << ", " << j << ")" << std::endl;
        spline_table[j][i] = spline;
      }

    // These I think are used when a lookup table must be searched
    // but since we know the lookup locations there is no point
    xacc = NULL; //gsl_interp_accel_alloc();
    yacc = NULL; //gsl_interp_accel_alloc();

}

void VEGAS_LUT::free_memory()
{
  for (int i = 0; i < N; i++)
  {
    delete[] table[i];
    delete[] spline_table[i];
  }

  //gsl_interp_accel_free(xacc);
  //gsl_interp_accel_free(yacc);

}

void VEGAS_LUT::free_input_table()
{
  for (int i = 0; i < N; i++)
    delete[] input_table[i];
}

void VEGAS_LUT::reduce_input_table()
{

    double Tmin,Tmax,ximin,ximax;
    Tmin = T_bounds[0];
    Tmax = T_bounds[1];
    ximin = xi_bounds[0];
    ximax = xi_bounds[1];

    if (Tmin < Tvals[0] || Tmax > Tvals[NT-1])
    {
    std::cout << "ERROR: Tmin or Tmax is not contained within file "
    << filename << "!"
    << std::endl;
    exit(0);
    }

    if (ximin < xivals[0] || ximax > xivals[Nxi-1])
    {
    std::cout << "ERROR: x_min or x_max is not contained within file "
    << filename << "!"
    << std::endl;
    exit(0);
    }

    // find the indices of the table corresponding to input min/max values
    int i = 0;
    while (Tvals[i] < Tmax && i < NT)
      i++;
    y_max = Tvals[i];  //so y_max is slightly greater than Tmax

    int ii = 0;
    while (Tvals[ii] <= Tmin && ii < i)
      ii++;
    y_min = Tvals[ii-1]; //so y_min is slightly less than Tmin

    // assign private variables
    i_T = ii-1; // start position in input table
    N = i - ii + 2; // # of T-vals in desired table

    i = 0;
    while (xivals[i] < ximax && i < Nxi)
      i++;
    x_max = xivals[i]; //so x_max is slightly greater than ximax

    ii = 0;
    while (xivals[ii] <= ximin && ii < i)
      ii++;
    x_min = xivals[ii-1]; //so x_min is slightly less than ximin

    // assign private variables
    i_xi = ii-1; // start position in input table
    M = i - ii + 2; // # of xi-vals in desired table
}

void VEGAS_LUT::initialize_table(const double * T_bbox = NULL, const double * xi_bbox = NULL)
{
	int i0,j0;

	if (T_bbox != NULL || xi_bbox != NULL)
	{
	  T_bounds = T_bbox;
	  xi_bounds = xi_bbox;
	  reduce_input_table(); // this sets i_xi, i_T, N, M, y_min, Tmax, x_min, x_max
	  i0 = i_xi;
	  j0 = i_T;
	}
	else // the actual table will be the same size as the input table
	{
	  y_min = Tvals[0];
	  y_max = Tvals[NT-1];
	  x_min = xivals[0];
	  x_max = xivals[Nxi-1];
	  N = NT;
	  M = Nxi;
	  i0 = 0;
	  j0 = 0;
	}

	// store the values of xi and T that will be used
    // they way they will be used (i.e. log or lin)
    xvals = new double[M];
    yvals = new double[N];
    double x,y;
    for (int i=0; i<M; i++)
    {
      if(logxi)	x = log10(xivals[i0+i]);
      else 		x = xivals[i0+i];
      xvals[i] = x;
    }
    for (int j=0; j<N; j++)
    {
      if (logT)	y = log10(Tvals[j0+j]);
      else 		y = Tvals[j0+j];
      yvals[j] = y;
    }

    // generate the actual table of H/C rates
	make_table(i0,j0);

	// alloc storage for max/min vals of T,xi
	xlim = new double[2]; // xi limits
  ylim = new double[2]; // T limits

    // set any log gridding: use limits of input table
	 if (logT) {
	   ylim[0] = log10(Tvals[0]);
	   ylim[1] = log10(Tvals[NT-1]); }
	  else {
	   ylim[0] = Tvals[0];
	   ylim[1] = Tvals[NT-1]; }

     if (logxi) {
       xlim[0] = log10(xivals[0]);
       xlim[1] = log10(xivals[Nxi-1]); }
     else {
       xlim[0] = xivals[0];
       xlim[1] = xivals[Nxi-1]; }

    // generate the spline table for GSL
	make_spline_table();

	// delete input_table
	free_input_table();

}

void VEGAS_LUT::determine_box_indices(double x, double y, unsigned int &i_x, unsigned int &i_y)
{
  unsigned int i,j;
  bool correct_guess = false;

  /* First guess the indices assuming the input table was
   * generated from a uniform grid: the actual xi and T values
   * that xstar produces data for may be slightly different
   * as xstar's stopping criteria algorithms may slightly adjust
   * the input values.  On rare occasions this guess can be off
   * by 1 index, so here we allow for that. */
  i = (Nxi-1) * (x - xlim[0]) / (xlim[1] - xlim[0]) - i_xi;
  j = (NT-1) * (y - ylim[0]) / (ylim[1] - ylim[0]) - i_T ;

  /* now we check our guess in x and if needed, adjust one box left/right */
  if (x > xvals[i] && x < xvals[i+1])  //then we guessed right
  {
    i_x = i;
    correct_guess = true;
  }

  if (!correct_guess)
  {
    if (i==0 || i==(M-1))
    {
      if (i==0 && x < xvals[2])
        i_x = 1;
      else if (i==(M-1) && x > xvals[M-2])
        i_x = M-2;
    }
    else if (x < xvals[i]) //then move to the left 1
      i_x = i-1;
    else if (x > xvals[i+1]) //then move to the right 1
      i_x = i+1;
    else // we give up
    {
      std::cout << "ERROR: search for the table xi-index failed!"
      			<< "\nDouble check inputs or try debugging "
      			<< "function determine_box_indices()"
      			<< std::endl;
      exit(0);
    }
  }

  /* now we check our guess in y and if needed, adjust one box up/down */
  correct_guess = false; // reset

  if (y > yvals[j] && y < yvals[j+1])  //then we guessed right
  {
    i_y = j;
    correct_guess = true;
  }

  if (!correct_guess)
  {
    if (j==0 || j==(N-1))
    {
      if (j==0 && y < yvals[2])
        i_y = 1;
      else if (j==(N-1) && y > yvals[N-2])
        i_y = N-2;
    }
    else if (y < yvals[j]) //then move down 1
      i_y = j-1;
    else if (y > yvals[j+1]) //then move up 1
      i_y = j+1;
    else // we give up
    {
      std::cout << "ERROR: search for the table T-index failed!"
      			<< "\nDouble check inputs or try debugging "
      			<< "function determine_box_indices()"
      			<< std::endl;
      exit(0);
    }
  }

}

double VEGAS_LUT::get_rate(double y_val, double x_val)
{
	unsigned int i_x,i_y;
	double x,y;

	/* terminate program if values exceed bounding box
    if (y_val < y_min || y_val > y_max || x_val < x_min || x_val > x_max)
    {
      std::cout
      << "\nFATAL ERROR: (y_val,x_val) = (" << y_val << ", " << x_val << ") "
      << "lies outside of table's bounding box!\n"
      << ">> Bounding Box:\n"
      << "(y_min,y_max) = (" << y_min << ", " << y_max << ")\n"
      << "(x_min,x_max) = (" << x_min << ", " << x_max << ")\n"
      << std::endl;
	  exit(0);
    }
    */

    /* use boundary rates for values outside bounding box */
    double eps = 1e-10;
    if (y_val < y_min) y_val = y_min*(1. + eps);
    if (y_val > y_max) y_val = y_max*(1. - eps);
    if (x_val < x_min) x_val = x_min*(1. + eps);
    if (x_val > x_max) x_val = x_max*(1. - eps);

	// apply any log scaling
	if (logT)  y = log10(y_val);
	else       y = y_val;
  if (logxi) x = log10(x_val);
  else       x = x_val;

	 determine_box_indices(x,y,i_x,i_y);

	return gsl_spline2d_eval(spline_table[i_y][i_x], x, y, xacc, yacc);
}

#endif //H_VEGAS_TABLE

ext: h

max_stars_repo_name: twaters/VegasTables

max_issues_repo_name: twaters/VegasTables

max_forks_repo_name: twaters/VegasTables


text:
#ifndef TAYLORSERIES_H
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <common/refcount.h>
#include "common/jexception.h"
#include <list>
#include <vector>

using namespace std;

double LogAdd(double x, double y);
double LogSub(double x, double y);
double logFactorial( unsigned int n );

#define NO_FILE_SYMBOL "NONE"
class nonamePdf {
 protected:
  vector<float>  _points; /* @note it must be positive and sorted. */
  int   _maxCoeffs;
  double **_coeffs;       /* [_points.size()][nCoeffs] */
  double **_coeffsLog;

 protected:
  void clear();

 public:
  nonamePdf();
  ~nonamePdf();
  bool loadCoeffDescFile( const String &coefDescfn );

 private:
  bool loadCoeffFiles( char *coeffn, char *logCoeffn, float a );
  int  loadCoeffFile(  char *coeffn, int preNPoints, int idx, bool logProb );
};

class gammaPdf : public nonamePdf {

 public:
  gammaPdf( int numberOfVariate = 2 );
  ~gammaPdf();
  double calcLog( double x, int N );
  double calcDerivative1( double x, int N );
  void   bi( int printLevel = 0);
  void   four( int printLevel = 0 );

  void   printCoeff();
 private:
  void allocate();
  int  indexOfCoeffArray( double x );
  bool interpolate( int printLevel );
};

typedef  refcount_ptr<nonamePdf> nonamePdfPtr;
typedef  Inherit<gammaPdf, nonamePdfPtr> gammaPdfPtr;

#endif


text:
/*
 * dgl.c -- Loesung einer DGL mit der GNU scientific library
 *
 * (c) 2016 Prof Dr Andreas Mueller, Hochschule Rapperswil
 */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

static long	evalcounter = 0;

int	f(double x, const double y[], double f[], void *params) {
	evalcounter++;
	double	omega = *(double *)params;
	f[0] = y[1];
	f[1] = -y[0] + sin(omega * x);
	return GSL_SUCCESS;
}

static inline double	sqr(double x) { return x * x; }

int	main(int argc, char *argv[]) {
	double	omega = 0.1;
	gsl_odeiv2_system	sys = { f, NULL, 2, &omega };
	gsl_odeiv2_driver	*driver
		= gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rkf45,
			1e-6, 1e-6, 0.0);
	double	x = 0.0;
	double	y[2] = { 0.0, 0.0 };
	long	lastcounter = evalcounter;
	for (int i = 1; i <= 10000; i++) {
		double	xnext = i;
		int	status = gsl_odeiv2_driver_apply(driver, &x, xnext, y);
		if (status != GSL_SUCCESS) {
			fprintf(stderr, "error: return value = %d\n", status);
		}
		double yexakt = (omega * sin(x) - sin(omega * x))
					/ (sqr(omega) - 1);
		double delta = y[0] - yexakt;
		long	evaldensity = evalcounter - lastcounter;
		lastcounter = evalcounter;
		if ((i < 10) ||
			((i < 100) && (0 == i % 10)) ||
			((i < 1000) && (0 == i % 100)) ||
			((i < 10000) && (0 == i % 1000)) ||
			((i < 100000) && (0 == i % 10000))) {
			printf("%5.0f& %12.8f& %12.8f& %12.8f&%ld\\\\\n",
				 x,    y[0],  yexakt, delta, evaldensity);
		}
	}
	return EXIT_SUCCESS;
}

max_stars_repo_name: MatthiasRubin/SeminarDGL

max_issues_repo_name: MatthiasRubin/SeminarDGL

max_forks_repo_name: MatthiasRubin/SeminarDGL


text:
#ifndef __MATH_VECTOROPS_INCLUDED
#define __MATH_VECTOROPS_INCLUDED

#include <cmath>
#include <limits>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
<<<<<<< vectorops.h
#include <gsl/gsl_cblas.h>
#include "math/specialfunc.h"
=======
#include <gsl/gsl_blas.h>
#include "specialfunc.h"
>>>>>>> 1.27

#ifndef M_PI
#define M_PI        3.14159265358979323846
#endif

#ifndef isnan
# define isnan(x) ((x) != (x))
#endif

/*
 * take the exponent of a vector
 *
 * If the exponent is infinite, then we replace the value with a
 * suitably large max_val
 */
void vexp(const gsl_vector* v,
          gsl_vector* exp_v,
          double max_val = std::numeric_limits<double>::infinity()) {
  assert(exp_v->size >= v->size);
  for (unsigned int ii = 0; ii < v->size; ++ii) {
    double val = exp(gsl_vector_get(v, ii));
    if (val == std::numeric_limits<double>::infinity() || val > max_val) {
      val = max_val;
    }
    gsl_vector_set(exp_v, ii, val);
  }
}

/* take the exponent of a matrix */
void mexp(const gsl_matrix* m, gsl_matrix* exp_m) {
  for (unsigned int ii = 0; ii < m->size1; ++ii) {
    for (unsigned int jj = 0; jj < m->size2; ++jj) {
      double val = exp(gsl_matrix_get(m, ii, jj));
      assert(!isnan(val));
      gsl_matrix_set(exp_m, ii, jj, val);
    }
  }
}

/* like vexp except that it also computes sum x log x */
double vexp_entropy(const gsl_vector* v, gsl_vector* exp_v) {
  double entropy = 0.0;
  for (unsigned int ii = 0; ii < v->size; ++ii) {
    double logval = gsl_vector_get(v, ii);
    double val = exp(logval);
    assert(!isnan(val));
    gsl_vector_set(exp_v, ii, val);
    if (val != 0) {
      entropy -= val * logval;
    }
  }
  return entropy;
}

double ventropy(const gsl_vector* v) {
  double entropy = 0.0;
  for (unsigned int ii = 0; ii < v->size; ++ii) {
    double val = gsl_vector_get(v, ii);
    if (val != 0) {
      entropy -= val * log(val);
    }
  }
  return entropy;
}

double lgamma(double x) {
    double x0,x2,xp,gl,gl0;
    int n=0,k=0;
    static double a[] = {
        8.333333333333333e-02,
       -2.777777777777778e-03,
        7.936507936507937e-04,
       -5.952380952380952e-04,
        8.417508417508418e-04,
       -1.917526917526918e-03,
        6.410256410256410e-03,
       -2.955065359477124e-02,
        1.796443723688307e-01,
       -1.39243221690590};

    x0 = x;
    if (x <= 0.0) return 1e308;
    else if ((x == 1.0) || (x == 2.0)) return 0.0;
    else if (x <= 7.0) {
        n = (int)(7-x);
        x0 = x+n;
    }
    x2 = 1.0/(x0*x0);
    xp = 2.0*M_PI;
    gl0 = a[9];
    for (k=8;k>=0;k--) {
        gl0 = gl0*x2 + a[k];
    }
    gl = gl0/x0+0.5*log(xp)+(x0-0.5)*log(x0)-x0;
    if (x <= 7.0) {
        for (k=1;k<=n;k++) {
            gl -= log(x0-1.0);
            x0 -= 1.0;
        }
    }
    return gl;
}

void mlog(const gsl_matrix* m, gsl_matrix* log_m) {
  for (unsigned int ii = 0; ii < m->size1; ++ii) {
    for (unsigned int jj = 0; jj < m->size2; ++jj) {
      gsl_matrix_set(log_m, ii, jj, log(gsl_matrix_get(m, ii, jj)));
    }
  }
}

void vlog(const gsl_vector* v, gsl_vector* log_v) {
  for (unsigned int ii = 0; ii < v->size; ++ii)
    gsl_vector_set(log_v, ii, log(gsl_vector_get(v, ii)));
}

void vlogit(const gsl_vector* v, gsl_vector* log_v) {
  for (unsigned int ii = 0; ii < v->size; ++ii) {
    double p = gsl_vector_get(v, ii);
    assert(p >= 0.0);
    assert(p <= 1.0);
    gsl_vector_set(log_v, ii, log(p / (1 - p)));
  }
}

void vsigmoid(const gsl_vector* v, gsl_vector* sig_v) {
  for (unsigned int ii = 0; ii < v->size; ++ii) {
    double p = gsl_vector_get(v, ii);
    gsl_vector_set(sig_v, ii, 1. / (1. + exp(-p)));
  }
}

double vlog_entropy(const gsl_vector* v, gsl_vector* log_v) {
  double entropy = 0;
  for (unsigned int ii = 0; ii < v->size; ++ii) {
    double val = gsl_vector_get(v, ii);
    entropy -= val * log(val);
    gsl_vector_set(log_v, ii, log(val));
  }
  return entropy;
}

double entropy(const gsl_vector* v) {
  double entropy = 0;
  for (unsigned int ii = 0; ii < v->size; ++ii) {
    double val = gsl_vector_get(v, ii);
    entropy -= val * log(val);
  }
  return entropy;
}

void vdigamma(const gsl_vector* v, gsl_vector* digamma_v) {
  for (unsigned int ii = 0; ii < v->size; ++ii)
    // gsl_sf_psi throws an error when its argument is 0, whereas digamma returns inf.
    //    gsl_vector_set(digamma_v, ii, gsl_sf_psi(gsl_vector_get(v, ii)));
    gsl_vector_set(digamma_v, ii, digamma(gsl_vector_get(v, ii)));
}

void vlgamma(const gsl_vector* v, gsl_vector* lgamma_v) {
  for (unsigned int ii = 0; ii < v->size; ++ii)
    gsl_vector_set(lgamma_v, ii, lgamma(gsl_vector_get(v, ii)));
}

double gsl_blas_dsum(const gsl_vector* v) {
  double sum = 0;
  for (unsigned int ii = 0; ii < v->size; ++ii) {
    sum += gsl_vector_get(v, ii);
  }
  return sum;
}

double gsl_blas_dsum(const gsl_matrix* v) {
  double sum = 0;
  for (unsigned int ii = 0; ii < v->size1; ++ii) {
    for (unsigned int jj = 0; jj < v->size2; ++jj) {
      sum += gsl_matrix_get(v, ii, jj);
    }
  }
  return sum;
}

double gsl_matrix_rowsum(const gsl_matrix* m, const int row) {
  double sum = 0;
  for(unsigned int i=0; i < m->size2; i++) {
    sum += gsl_matrix_get(m, row, i);
  }
  return sum;
}

double dot_product(const gsl_vector* a, const gsl_vector* b) {
  assert(a->size == b->size);
  double val = 0;
  for(unsigned i=0; i<a->size; i++) {
    val += gsl_vector_get(a, i) * gsl_vector_get(b, i);
  }
  return val;
}

void uniform(gsl_vector* v) {
	gsl_vector_set_all(v, 1.0 / (double)v->size);
}

double normalize(gsl_vector* v) {
  double sum = gsl_blas_dsum(v);
  gsl_blas_dscal(1 / sum, v);
  return sum;
}

/*
  This function takes as input a multinomial parameter vector and
  computes the "total" variance, i.e., the sum of the diagonal of the
  covariance matrix.

  If the multinomial parameter is unnormalized, then the variance of
  the normalized multinomial vector will be computed and then
  multiplied by the scale of the vector.
 */
double MultinomialTotalVariance(const gsl_vector* v) {
  double scale = gsl_blas_dsum(v);
  double variance = 0.0;
  for (size_t ii = 0; ii < v->size; ++ii) {
    double val = gsl_vector_get(v, ii) / scale;
    variance += val * (1. - val);
  }
  return variance * scale;
}

/*
  Computes covariance using the renormalization above and adds it to
  an existing matrix.
*/
void MultinomialCovariance(double alpha,
			   const gsl_vector* v,
			   gsl_matrix* m) {
  double scale = gsl_blas_dsum(v);
  gsl_blas_dger(-alpha / scale, v, v, m);
  gsl_vector_view diag = gsl_matrix_diagonal(m);
  gsl_blas_daxpy(alpha, v, &diag.vector);
}

double MatrixProductSum(const gsl_matrix* m1,
			const gsl_matrix* m2) {
  double val = 0;
  assert(m1->size1 == m2->size1);
  assert(m1->size2 == m2->size2);
  for (size_t ii = 0; ii < m1->size1; ++ii) {
    for (size_t jj = 0; jj < m2->size2; ++jj) {
      val += gsl_matrix_get(m1, ii, jj) *
	gsl_matrix_get(m2, ii, jj);
    }
  }
  return val;
}

double MatrixProductProductSum(const gsl_matrix* m1,
			       const gsl_matrix* m2,
			       const gsl_matrix* m3) {
  double val = 0;
  assert(m1->size1 == m2->size1);
  assert(m1->size2 == m2->size2);
  assert(m1->size1 == m3->size1);
  assert(m1->size2 == m3->size2);
  for (size_t ii = 0; ii < m1->size1; ++ii) {
    for (size_t jj = 0; jj < m2->size2; ++jj) {
      for (size_t kk = 0; kk < m3->size2; ++kk) {
	val += gsl_matrix_get(m1, ii, jj) *
	  gsl_matrix_get(m2, ii, jj) *
	  gsl_matrix_get(m3, ii, jj);
      }
    }
  }
  return val;
}

double SumLGamma(const gsl_vector* v) {
  double s = 0.0;
  for (size_t ii = 0; ii < v->size; ++ii) {
    s += lgamma(gsl_vector_get(v, ii));
  }
  return s;
}

void mtx_fprintf(const char* filename, const gsl_matrix * m)
{
    FILE* fileptr;
    fileptr = fopen(filename, "w");
    gsl_matrix_fprintf(fileptr, m, "%20.17e");
    fclose(fileptr);
}

void mtx_fscanf(const char* filename, gsl_matrix * m)
{
    FILE* fileptr;
    fileptr = fopen(filename, "r");
    gsl_matrix_fscanf(fileptr, m);
    fclose(fileptr);
}

double mtx_accum(const int i,
                 const int j,
                 const double contribution,
                 gsl_matrix* m) {

  double new_val = gsl_matrix_get(m, i, j) + contribution;
  gsl_matrix_set(m, i, j, new_val);
  return new_val;
}

void vct_fscanf(const char* filename, gsl_vector* v)
{
    FILE* fileptr;
    fileptr = fopen(filename, "r");
    gsl_vector_fscanf(fileptr, v);
    fclose(fileptr);
}

void vct_fprintf(const char* filename, gsl_vector* v)
{
    FILE* fileptr;
    fileptr = fopen(filename, "w");
    gsl_vector_fprintf(fileptr, v, "%20.17e");
    fclose(fileptr);
}

#endif

max_stars_repo_name: boomsbloom/dtm-fmri

max_issues_repo_name: boomsbloom/dtm-fmri

max_forks_repo_name: boomsbloom/dtm-fmri


text:
/* This file is part of the KDE project
   Copyright (C) 2006 Stefan Nikolaus <stefan.nikolaus@kdemail.net>

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

   This library 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
   Library General Public License for more details.

   You should have received a copy of the GNU Library General Public License
   along with this library; see the file COPYING.LIB.  If not, write to
   the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
   Boston, MA 02110-1301, USA.
*/

#ifndef KSPREAD_SOLVER
#define KSPREAD_SOLVER

#include <gsl/gsl_multimin.h>

#include <kxmlguiclient.h>

#include <Cell.h>

#include <QObject>
#include <QVariantList>

namespace Calligra
{
namespace Sheets
{
namespace Plugins
{

/**
 * \class Solver Function Optimizer
 * \author Stefan Nikolaus <stefan.nikolaus@kdemail.net>
 */
class Solver : public QObject, public KXMLGUIClient
{
    Q_OBJECT

public:
    struct Parameters {
        QList<Cell> cells;
    };

    /**
     * Constructor.
     */
    Solver(QObject* parent, const QVariantList& args);

    /**
     * Destructor.
     */
    ~Solver();

    double evaluate(const gsl_vector* vector, void* parameters);

protected Q_SLOTS:
    /**
     * Called when the Solver action is triggered.
     * Opens the dialog.
     */
    void showDialog();

    /**
     * This method does the real work.
     * Uses the parameters of the dialog to optimize the given function.
     */
    void optimize();

private:
    Q_DISABLE_COPY(Solver)

    class Private;
    Private * const d;
};

} // namespace Plugins
} // namespace Sheets
} // namespace Calligra

#endif

max_stars_repo_name: afarcat/QtSheetView

max_issues_repo_name: afarcat/QtSheetView

max_forks_repo_name: afarcat/QtSheetView


text:
#pragma once

#include "Cesium3DTiles/BoundingVolume.h"
#include "Cesium3DTiles/Library.h"
#include "Cesium3DTiles/Tile.h"
#include "Cesium3DTiles/TileContext.h"
#include "Cesium3DTiles/TileID.h"
#include "Cesium3DTiles/TileRefine.h"

#include <gsl/span>
#include <spdlog/fwd.h>

#include <cstddef>
#include <memory>

namespace Cesium3DTiles {
/**
 * @brief The information that is passed to a {@link TileContentLoader} to
 * create a {@link TileContentLoadResult}.
 *
 * For many types of tile content, only the `data` field is required. The other
 * members are used for content that can generate child tiles, like external
 * tilesets or composite tiles. These members are usually initialized from
 * the corresponding members of the {@link Tile} that the content belongs to.
 */
struct CESIUM3DTILES_API TileContentLoadInput {

  /**
   * @brief Creates a new, uninitialized instance for the given tile.
   *
   * The `data`, `contentType` and `url` will have default values,
   * and have to be initialized before this instance is passed to
   * one of the loader functions.
   *
   * @param pLogger_ The logger that will be used
   * @param tile_ The {@link Tile} that the content belongs to
   */
  TileContentLoadInput(
      const std::shared_ptr<spdlog::logger> pLogger_,
      const Tile& tile_)
      : pLogger(pLogger_),
        data(),
        contentType(),
        url(),
        tileID(tile_.getTileID()),
        tileBoundingVolume(tile_.getBoundingVolume()),
        tileContentBoundingVolume(tile_.getContentBoundingVolume()),
        tileRefine(tile_.getRefine()),
        tileGeometricError(tile_.getGeometricError()),
        tileTransform(tile_.getTransform()) {}

  /**
   * @brief Creates a new instance for the given tile.
   *
   * @param pLogger_ The logger that will be used
   * @param data_ The actual data that the tile content will be created from
   * @param contentType_ The content type, if the data was received via a
   * network response
   * @param url_ The URL that the data was loaded from
   * @param tile_ The {@link Tile} that the content belongs to
   */
  TileContentLoadInput(
      const std::shared_ptr<spdlog::logger> pLogger_,
      const gsl::span<const std::byte>& data_,
      const std::string& contentType_,
      const std::string& url_,
      const Tile& tile_)
      : pLogger(pLogger_),
        data(data_),
        contentType(contentType_),
        url(url_),
        tileID(tile_.getTileID()),
        tileBoundingVolume(tile_.getBoundingVolume()),
        tileContentBoundingVolume(tile_.getContentBoundingVolume()),
        tileRefine(tile_.getRefine()),
        tileGeometricError(tile_.getGeometricError()),
        tileTransform(tile_.getTransform()) {}

  /**
   * @brief Creates a new instance.
   *
   * For many types of tile content, only the `data` field is required. The
   * other parameters are used for content that can generate child tiles, like
   * external tilesets or composite tiles.
   *
   * @param pLogger_ The logger that will be used
   * @param data_ The actual data that the tile content will be created from
   * @param contentType_ The content type, if the data was received via a
   * network response
   * @param url_ The URL that the data was loaded from
   * @param context_ The {@link TileContext}
   * @param tileID_ The {@link TileID}
   * @param tileBoundingVolume_ The tile {@link BoundingVolume}
   * @param tileContentBoundingVolume_ The tile content {@link BoundingVolume}
   * @param tileRefine_ The {@link TileRefine} strategy
   * @param tileGeometricError_ The geometric error of the tile
   * @param tileTransform_ The tile transform
   */
  TileContentLoadInput(
      const std::shared_ptr<spdlog::logger> pLogger_,
      const gsl::span<const std::byte>& data_,
      const std::string& contentType_,
      const std::string& url_,
      const TileID& tileID_,
      const BoundingVolume& tileBoundingVolume_,
      const std::optional<BoundingVolume>& tileContentBoundingVolume_,
      TileRefine tileRefine_,
      double tileGeometricError_,
      const glm::dmat4& tileTransform_)
      : pLogger(pLogger_),
        data(data_),
        contentType(contentType_),
        url(url_),
        tileID(tileID_),
        tileBoundingVolume(tileBoundingVolume_),
        tileContentBoundingVolume(tileContentBoundingVolume_),
        tileRefine(tileRefine_),
        tileGeometricError(tileGeometricError_),
        tileTransform(tileTransform_) {}

  /**
   * @brief The logger that receives details of loading errors and warnings.
   */
  std::shared_ptr<spdlog::logger> pLogger;

  /**
   * @brief The raw input data.
   *
   * The {@link TileContentFactory} will try to determine the type of the
   * data using the first four bytes (i.e. the "magic header"). If this
   * does not succeed, it will try to determine the type based on the
   * `contentType` field.
   */
  gsl::span<const std::byte> data;

  /**
   * @brief The content type.
   *
   * If the data was obtained via a HTTP response, then this will be
   * the `Content-Type` of that response. The {@link TileContentFactory}
   * will try to interpret the data based on this content type.
   *
   * If the data was not directly obtained from an HTTP response, then
   * this may be the empty string.
   */
  std::string contentType;

  /**
   * @brief The source URL.
   */
  std::string url;

  /**
   * @brief The {@link TileID}.
   */
  TileID tileID;

  /**
   * @brief The tile {@link BoundingVolume}.
   */
  BoundingVolume tileBoundingVolume;

  /**
   * @brief Tile content {@link BoundingVolume}.
   */
  std::optional<BoundingVolume> tileContentBoundingVolume;

  /**
   * @brief The {@link TileRefine}.
   */
  TileRefine tileRefine;

  /**
   * @brief The geometric error.
   */
  double tileGeometricError;

  /**
   * @brief The tile transform
   */
  glm::dmat4 tileTransform;
};
} // namespace Cesium3DTiles

max_stars_repo_name: TJKoury/cesium-native

max_issues_repo_name: MAPSWorks/cesium-native

max_forks_repo_name: MAPSWorks/cesium-native


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_chebyshev.h>

/* For efficiency there are separate implementations of each of these
   functions */

double
gsl_cheb_eval (const gsl_cheb_series * cs, const double x)
{
  size_t i;
  double d1 = 0.0;
  double d2 = 0.0;

  double y = (2.0 * x - cs->a - cs->b) / (cs->b - cs->a);
  double y2 = 2.0 * y;

  for (i = cs->order; i >= 1; i--)
    {
      double temp = d1;
      d1 = y2 * d1 - d2 + cs->c[i];
      d2 = temp;
    }

  return y * d1 - d2 + 0.5 * cs->c[0];
}

double
gsl_cheb_eval_n (const gsl_cheb_series * cs, const size_t n, const double x)
{
  size_t i;
  double d1 = 0.0;
  double d2 = 0.0;

  size_t eval_order = GSL_MIN (n, cs->order);

  double y = (2.0 * x - cs->a - cs->b) / (cs->b - cs->a);
  double y2 = 2.0 * y;

  for (i = eval_order; i >= 1; i--)
    {
      double temp = d1;
      d1 = y2 * d1 - d2 + cs->c[i];
      d2 = temp;
    }

  return y * d1 - d2 + 0.5 * cs->c[0];
}

int
gsl_cheb_eval_err (const gsl_cheb_series * cs, const double x,
                   double *result, double *abserr)
{
  size_t i;
  double d1 = 0.0;
  double d2 = 0.0;

  double y = (2. * x - cs->a - cs->b) / (cs->b - cs->a);
  double y2 = 2.0 * y;

  double absc = 0.0;

  for (i = cs->order; i >= 1; i--)
    {
      double temp = d1;
      d1 = y2 * d1 - d2 + cs->c[i];
      d2 = temp;
    }

  *result = y * d1 - d2 + 0.5 * cs->c[0];

  /* Estimate cumulative numerical error */

  for (i = 0; i <= cs->order; i++)
    {
      absc += fabs(cs->c[i]);
    }

  /* Combine truncation error and numerical error */

  *abserr = fabs (cs->c[cs->order]) + absc * GSL_DBL_EPSILON;

  return GSL_SUCCESS;
}

int
gsl_cheb_eval_n_err (const gsl_cheb_series * cs,
                     const size_t n, const double x,
                     double *result, double *abserr)
{
  size_t i;
  double d1 = 0.0;
  double d2 = 0.0;

  double y = (2. * x - cs->a - cs->b) / (cs->b - cs->a);
  double y2 = 2.0 * y;

  double absc = 0.0;

  size_t eval_order = GSL_MIN (n, cs->order);

  for (i = eval_order; i >= 1; i--)
    {
      double temp = d1;
      d1 = y2 * d1 - d2 + cs->c[i];
      d2 = temp;
    }

  *result = y * d1 - d2 + 0.5 * cs->c[0];

  /* Estimate cumulative numerical error */

  for (i = 0; i <= eval_order; i++)
    {
      absc += fabs(cs->c[i]);
    }

  /* Combine truncation error and numerical error */

  *abserr = fabs (cs->c[eval_order]) + absc * GSL_DBL_EPSILON;

  return GSL_SUCCESS;
}

int
gsl_cheb_eval_mode_e (const gsl_cheb_series * cs,
                      const double x, gsl_mode_t mode,
                      double *result, double *abserr)
{
  size_t i;
  double d1 = 0.0;
  double d2 = 0.0;

  double y = (2. * x - cs->a - cs->b) / (cs->b - cs->a);
  double y2 = 2.0 * y;

  double absc = 0.0;

  size_t eval_order;

  if (GSL_MODE_PREC (mode) == GSL_PREC_DOUBLE)
    eval_order = cs->order;
  else
    eval_order = cs->order_sp;

  for (i = eval_order; i >= 1; i--)
    {
      double temp = d1;
      d1 = y2 * d1 - d2 + cs->c[i];
      d2 = temp;
    }

  *result = y * d1 - d2 + 0.5 * cs->c[0];

  /* Estimate cumulative numerical error */

  for (i = 0; i <= eval_order; i++)
    {
      absc += fabs(cs->c[i]);
    }

  /* Combine truncation error and numerical error */

  *abserr = fabs (cs->c[eval_order]) + absc * GSL_DBL_EPSILON;

  return GSL_SUCCESS;
}

double
gsl_cheb_eval_mode (const gsl_cheb_series * cs,
                    const double x, gsl_mode_t mode)
{
  double result, abserr;
  int status = gsl_cheb_eval_mode_e (cs, x, mode, &result, &abserr);

  if (status != GSL_SUCCESS)
    {
      GSL_ERROR_VAL("gsl_cheb_eval_mode", status, result);
    };

  return result;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified