text:
#ifndef MATH_FNS_H
#define MATH_FNS_H

#ifdef USE_OPENBLAS
extern "C"
{
#include <cblas.h>
}
#endif

struct QuantizationParams
{
	float scale;
	uint8_t zero_point;
};

template<typename Dtype>
void cpu_gemm(bool transA, bool transB, uint64_t M, uint64_t N, uint64_t K, Dtype alpha, Dtype* A, Dtype* B, Dtype beta, Dtype* C);

template<typename Dtype>
void cpu_axpy(uint64_t N, Dtype alpha, Dtype* X, Dtype* Y, Dtype* C); // c = alpha * x + y

template<typename Dtype>
void cpu_axpby(uint64_t N, Dtype alpha, Dtype* X_ptr, Dtype beta, Dtype* Y_ptr, Dtype* C_ptr); // c = alpha * x + beta * y

template<typename Dtype>
void cpu_mul(uint64_t N, Dtype* A_ptr, Dtype* B_ptr, Dtype* C_ptr); // c = a * b

template<typename Dtype>
void cpu_mul(uint64_t N, Dtype alpha, Dtype* A_ptr, Dtype* B_ptr, Dtype beta, Dtype* C_ptr); // c = beta * c + alpha * a * b

template<typename Dtype>
void cpu_mul(uint64_t N, Dtype alpha, Dtype* A_ptr, Dtype* B_ptr); // b = alpha * a

template<typename Dtype>
void cpu_add(uint64_t N, Dtype alpha, Dtype* A_ptr, Dtype* B_ptr); // b = alpha + a

template<typename Dtype>
void cpu_div(uint64_t N, Dtype* A_ptr, Dtype* B_ptr, Dtype* C_ptr); // c += a / b

template<typename Dtype>
void cpu_div(uint64_t N, Dtype alpha, Dtype* A_ptr, Dtype* B_ptr, Dtype beta, Dtype* C_ptr); // c = (beta * c) +  (alpha * a / b)

template<typename Dtype>
void cpu_div_back(uint64_t N, Dtype* A_ptr, Dtype* B_ptr, Dtype* C_ptr, Dtype* D_ptr); // special function for processing div_backward during backpropagation (calculates d += a * (-b) / (c * c))

template<typename Dtype>
void cpu_sig(uint64_t N, const Dtype* A_ptr, Dtype* B_ptr); // b = sig(a)

template<typename Dtype>
void cpu_tanh(uint64_t N, const Dtype* A_ptr, Dtype* B_ptr); // b = tanh(a)

template<typename Dtype>
void cpu_powx(uint64_t N, const Dtype* A_ptr, Dtype x, Dtype* B_ptr); // b = pow(a,x)

template<typename Dtype>
void cpu_copy(uint64_t N, Dtype* A_ptr, Dtype* B_ptr); // b = a

template<typename Dtype>
void cpu_max(const Dtype* src, Dtype* dst, uint64_t* indices, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void cpu_max_backward(Dtype* dst, const Dtype* src, const uint64_t* indices, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void cpu_sum(const Dtype* src, Dtype* dst, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void cpu_sum_backward(Dtype* dst, const Dtype* src, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void cpu_mean(const Dtype* src, Dtype* dst, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void cpu_mean_backward(Dtype* dst, const Dtype* src, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void cpu_var(const Dtype* src, Dtype* dst, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void cpu_var_backward(Dtype* dst, const Dtype* src, const Dtype* op1, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void cpu_std(const Dtype* src, Dtype* dst, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride, bool sample_mode = true);

template<typename Dtype>
void cpu_std_backward(Dtype* dst, const Dtype* src, const Dtype* op1, const Dtype* std, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride, bool sample_mode = true);

template<typename Dtype>
void cpu_dropout(Dtype* dst, Dtype* src, unsigned int* mask, unsigned int threshold, Dtype scale, uint64_t len);

template<typename Dtype>
void cpu_sig_backward(Dtype* bottom, const Dtype* top, const Dtype* middle, const uint64_t numels);

template<typename Dtype>
void cpu_tanh_backward(Dtype* bottom, const Dtype* top, const Dtype* middle, const uint64_t numels);

void cpu_transpose(float* src, float* dst, int dim_1, int dim_2,
	int stride_src_dim_1, int stride_src_dim_1_minus_1, int stride_src_dim_2, int stride_src_dim_2_minus_1,
	int stride_trn_dim_1, int stride_trn_dim_1_minus_1, int stride_trn_dim_2, int stride_trn_dim_2_minus_1, uint64_t numels);

void quantized_matmul(bool traspose_A, bool traspose_B, uint64_t M, uint64_t N, uint64_t K, uint8_t alpha, uint8_t* A, uint8_t* B, uint8_t beta, uint8_t* C, QuantizationParams* qparms, int* bias, int* workspace); // workspace must be at least (M + N) * sizeof(int) bytes

//-------------------------------------CUDA functions----------------------------------------------------------------

template<typename Dtype>
void gpu_sum(Dtype* A, Dtype* B, Dtype* C, uint64_t height_A, uint64_t width_A, uint64_t height_B, uint64_t width_B);

template<typename Dtype>
void gpu_sub(Dtype* A, Dtype* B, Dtype* C, uint64_t height_A, uint64_t width_A, uint64_t height_B, uint64_t width_B);

// C =  A * B
template<typename Dtype>
void gpu_mul(uint64_t N, Dtype* A_ptr, Dtype* B_ptr, Dtype* C_ptr);

// C = beta * C + A * B
// supports broadcast semantics
// use this when C dims >= both A and B dims (avoid atomics)
template<typename Dtype>
void gpu_mul(Dtype* A, Dtype* B, Dtype* C, uint64_t height_A, uint64_t width_A, uint64_t height_B, uint64_t width_B, Dtype beta);

// C = C + A * B
// supports broadcast semantics
// use this when C dims < either A or B dims (uses atomics, but no support for uint8_t)
void gpu_mul(float* A, float* B, float* C, uint64_t height_A, uint64_t width_A, uint64_t height_B, uint64_t width_B, uint64_t height_C, uint64_t width_C);

template<typename Dtype>
void gpu_mul(uint64_t N, Dtype alpha, Dtype* A, Dtype* B, Dtype beta, Dtype* C); // c = beta * c + alpha * a * b

template<typename Dtype>
void gpu_mul(uint64_t N, Dtype alpha, Dtype* A_ptr, Dtype* B_ptr); // b = alpha * a

// C =  A / B
// supports broadcast semantics
// use this when C dims >= both A and B dims (avoid atomics)
template<typename Dtype>
void gpu_div(Dtype* A, Dtype* B, Dtype* C, uint64_t height_A, uint64_t width_A, uint64_t height_B, uint64_t width_B);

// C = C + A / B
// use this when C numels == A numels == B numels (no broadcast therefore fastest)
template<typename Dtype>
void gpu_div(uint64_t N, Dtype* A, Dtype* B, Dtype* C);

// C = A / B
// use this when C numels == A numels == B numels (no broadcast therefore fastest)
//template<typename Dtype>
//void gpu_div2(uint64_t N, Dtype* A, Dtype* B, Dtype* C);

// C = C + A / B
// supports broadcast semantics
// use this when C dims < either A or B dims (uses atomics, but no support for uint8_t)
void gpu_div(float* A, float* B, float* C, uint64_t height_A, uint64_t width_A, uint64_t height_B, uint64_t width_B, uint64_t height_C, uint64_t width_C);

// special function for processing div_backward during backpropagation (calculates d += a * (-b) / (c * c))
// supports broadcast semantics
// use this when D dims < A or B or C dims (uses atomics, but no support for uint8_t)
//template<typename Dtype>
void gpu_div_back(float* A, float* B, float* C, float* D, uint64_t height_A, uint64_t width_A, uint64_t height_B, uint64_t width_B, uint64_t height_C, uint64_t width_C, uint64_t height_D, uint64_t width_D);

template<typename Dtype>
void gpu_sum(Dtype* data, Dtype* sum, uint64_t len);

template<typename Dtype>
void gpu_nll(Dtype* input, Dtype* target, Dtype* loss, uint64_t len, uint64_t batches);

template<typename Dtype>
void gpu_scalar_mul(Dtype* A, Dtype* B, Dtype scalar, uint64_t len); // B = scaler * A

// C += alpha * A
// supports broadcast semantics
// use this when C dims < A dims (uses atomics, but no support for uint8_t)
void gpu_scalar_mul(float alpha, float* A, float* C, uint64_t height_A, uint64_t width_A, uint64_t height_C, uint64_t width_C);

template<typename Dtype>
void gpu_fill(Dtype* memory, uint64_t size, Dtype value);

template<typename Dtype>
void gpu_fill(Dtype* memory, uint64_t len, Dtype* value);

void gpu_sgd_step(float* weight_ptr, float* weight_grad_ptr, float* velocity_ptr, int64_t numels, float mo, float wd, float lr);

template<typename Dtype>
void gpu_max(const Dtype* src, Dtype* dst, uint64_t* indices, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void gpu_max_backward(Dtype* dst, const Dtype* src, const uint64_t* indices, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void gpu_powx(uint64_t N, const Dtype* A, Dtype x, Dtype* B); // b = pow(a,x)

template<typename Dtype>
void gpu_exp(Dtype* dst, const Dtype* src, const uint64_t numels);

template<typename Dtype>
void gpu_log(Dtype* dst, const Dtype* src, const uint64_t numels);

template<typename Dtype>
void gpu_sig(Dtype* dst, const Dtype* src, const uint64_t numels);

template<typename Dtype>
void gpu_sig_backward(Dtype* bottom, const Dtype* top, const Dtype* middle, const uint64_t numels);

template<typename Dtype>
void gpu_tanh(Dtype* dst, const Dtype* src, const uint64_t numels);

template<typename Dtype>
void gpu_tanh_backward(Dtype* bottom, const Dtype* top, const Dtype* middle, const uint64_t numels);

template<typename Dtype>
void gpu_sum(const Dtype* src, Dtype* dst, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void gpu_sum_backward(Dtype* dst, const Dtype* src, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void gpu_add(uint64_t N, Dtype alpha, Dtype* A_ptr, Dtype* B_ptr); // b = alpha + a

template<typename Dtype>
void gpu_axpy(uint64_t N, Dtype alpha, Dtype* X_ptr, Dtype* Y_ptr, Dtype* C_ptr);

template<typename Dtype>
void gpu_axpby(uint64_t N, Dtype alpha, Dtype* X_ptr, Dtype beta, Dtype* Y_ptr, Dtype* C_ptr); // c = alpha * x + beta * y

template<typename Dtype>
void gpu_relu(Dtype* dst, Dtype* src, uint64_t len);

template<typename Dtype>
void gpu_relu_backward(Dtype* bottom, const Dtype* top, const Dtype* middle, const uint64_t len);

template<typename Dtype>
void gpu_dropout(Dtype* dst, Dtype* src, unsigned int* mask, unsigned int threshold, Dtype scale, uint64_t len);

template<typename Dtype>
void gpu_transpose(Dtype* src, Dtype* dst, int dim_1, int dim_2,
	int stride_src_dim_1, int stride_src_dim_1_minus_1, int stride_src_dim_2, int stride_src_dim_2_minus_1,
	int stride_trn_dim_1, int stride_trn_dim_1_minus_1, int stride_trn_dim_2, int stride_trn_dim_2_minus_1, uint64_t numels);

template<typename Dtype>
void gpu_cat(Dtype* dest, Dtype* op1, Dtype* op2, uint64_t dest_stride_1, uint64_t dest_stride_2, uint64_t op1_stride, uint64_t op2_stride_1, uint64_t op2_stride_2, uint64_t dim_offset, uint64_t op1_numels, uint64_t op2_numels);

template<typename Dtype>
void gpu_cat_backward(Dtype* dest, Dtype* src, uint64_t dest_stride, uint64_t src_stride, uint64_t dest_numels);

template<typename Dtype>
void gpu_cat_backward(Dtype* dest, Dtype* src, uint64_t dest_stride_1, uint64_t dest_stride_2, uint64_t src_stride_1, uint64_t src_stride_2, uint64_t dim_offset, uint64_t op1_numels, uint64_t dest_numels);

template<typename Dtype>
void gpu_embedding(Dtype* dst, Dtype* wts, int* indices, uint64_t numels, uint64_t indices_per_batch, unsigned int embedding_dim);

template<typename Dtype>
void gpu_embedding_backward(Dtype* dst, Dtype* wts, int* indices, uint64_t numels, uint64_t indices_per_batch, unsigned int embedding_dim);

template<typename Dtype>
void gpu_sqrt(Dtype* dst, const Dtype* src, const uint64_t numels);

template<typename Dtype>
void gpu_sqrt_backward(Dtype* bottom, const Dtype* top, const Dtype* middle, const uint64_t numels);

template<typename Dtype>
void gpu_mean(const Dtype* src, Dtype* dst, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void gpu_mean_backward(Dtype* dst, const Dtype* src, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void gpu_var(const Dtype* src, Dtype* dst, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void gpu_var_backward(Dtype* dst, const Dtype* src, const Dtype* op1, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride);

template<typename Dtype>
void gpu_std(const Dtype* src, Dtype* dst, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride, bool sample_mode = true);

template<typename Dtype>
void gpu_std_backward(Dtype* dst, const Dtype* src, const Dtype* op1, const Dtype* std, const uint64_t numels, const uint64_t ratio, const uint64_t dim_size, const uint64_t stride, bool sample_mode = true);

#endif // MATH_FNS_H

max_stars_repo_name: adeobootpin/light-tensor

max_issues_repo_name: adeobootpin/light-tensor

max_forks_repo_name: adeobootpin/light-tensor


text:
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

#include "core_allvars.h"
#include "core_proto.h"
#include "UVmag/UVmag.h"

// Local Proto-Types //

int32_t init_delayedSN(void);
int32_t init_nionlookup(void);
int32_t init_metalcooling(void);
void read_snap_list(void);
void set_units(void);
int32_t init_fesc(void);
int32_t determine_fescMH_constants(void);

// External Functions //

void sage_init(void)
{
  int32_t i, status;

  count_gal = 0;
  random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);
  gsl_rng_set(random_generator, 42);	 // start-up seed

  set_units();
  srand((unsigned) time(NULL));

  read_snap_list();
  for(i = 0; i < MAXSNAPS; i++)
  {
    ZZ[i] = 1 / AA[i] - 1;
    Age[i] = time_to_present(ZZ[i]);
  }

  a0 = 1.0 / (1.0 + Reionization_z0);
  ar = 1.0 / (1.0 + Reionization_zr);

  read_cooling_functions();

  status = init_fesc();
  if (status != EXIT_SUCCESS)
  {
    ABORT(0);
  }

  count_onehalo = 0;

  zeromass_count = 0;
  suppression_count = 0;
  previous_tree = 0;

  smallest_mass = 100000;
  lowmass_halo = 0;

  outside_box = 0;
  inside_box = 0;

  count_Mvir = 0;
  count_Len = 0;

  if (IMF == 0)
  {
    // Salpeter IMF //
    IMF_norm =0.1706;
    IMF_slope = -2.35;
    Eta_SNII = 7.432e-3; //Number fraction of stars that will end their lives as type II supernovae.
    m_SNII = 0.144; // Mass fraction of stars that will end their lives as type II supernovae.

  } else if (IMF == 1)
  {
    // Chabrier IMF //
    IMF_norm = 0.23638;
    IMF_slope = -2.3;
    Eta_SNII = 1.1893e-2; //Number fraction of stars that will end their lives as type II supernovae.
    m_SNII = 0.23638; // Mass fraction of stars that will end their lives as type II supernovae.

  }

  // Tiamat Parameters //

  V_energy = 70.0;
  alpha_energy = 0.5;
  beta_energy = 10.0;

  V_mass = 70.0;
  alpha_mass = 6.0;
  beta_mass = 10.0;

  epsilon_mass_max = 30.0;

  if (IRA == 0)
  {
    status = init_delayedSN();
    if (status != EXIT_SUCCESS)
    {
      ABORT(EXIT_FAILURE);
    }
  }

  if (PhotonPrescription == 1)
  {
    status = init_nionlookup();
    if (status != EXIT_SUCCESS)
    {
      ABORT(EXIT_FAILURE);
    }
  }

  if (calcUVmag == 1)
  {
    status = init_UVlookup();
    if (status != EXIT_SUCCESS)
    {
      ABORT(EXIT_FAILURE);
    }
  }

  if ((PhotonPrescription == 1) || (calcUVmag == 1))
  {
    StellarTracking_Len = STELLAR_TRACKING_TIME / TimeResolutionStellar;
  }

  mergedgal_mallocs = 0;
  gal_mallocs = 0 ;

  mergedgal_frees = 0;
  gal_frees = 0;

  Ngamma_HI_Total = 0.0;

}

// Local Functions //

int32_t init_delayedSN(void)
{

  // We keep the past 50 Myrs of star formation stored for each galaxy.
  // Based on the TimeResolution specified need to find out how many elements this corresponds to.

  if(TimeResolutionSN > 50)
  {
    fprintf(stderr, "The selected time resolution for SN feedback (TimeResolutionSN) is set too high (%d Myr).  Using TimeResolutionSN > 50Myr is the same as using the instantaneous recycling approximation; set 'IRA' to 1 instead!\n", TimeResolutionSN);
    ABORT(EXIT_FAILURE);
  } else if(TimeResolutionSN > 35)
  {
    fprintf(stderr, "Your selected time resolution for SN feedback (TimeResolutionSN) is quite high (%d Myr).  Beyond 50Myr the instantaneous recycling approximation is valid hence with your value it would likely be correct to set 'IRA' to 1.\n", TimeResolutionSN);
  } else
  {
    Time_SFH = 0;
    SN_Array_Len = 0;
    while(Time_SFH < 50)
    {
      Time_SFH += TimeResolutionSN;
      ++SN_Array_Len;
    }

  }

  // For the delayed SN scheme the number of SN for each time step depends on the IMF.
  // As we will need to calculate this for every substep we create look-up tables for efficiency.

  // First we need to know given a timestep, what is the minimum mass star that will go supernova?

  double a = 0.7473; // Fits from Portinari et al. (1998).
  double b = -2.6979;
  double c = -4.7659;
  double d = 0.5934;

  coreburning_tbins_low = 2; // Any time below 2Myr will result in stars > 120Myr to explode. This is beyond the limits of the IMF.
  coreburning_tbins_high = 45; // Any time above 45Myr will result in stars > 8Myr to explode. Can just return 8.0 in this case.
  coreburning_tbins_delta = 0.00001;

  N_tbins = (coreburning_tbins_high - coreburning_tbins_low) / (coreburning_tbins_delta);
  int32_t bin_idx;

  coreburning_times = malloc(sizeof(*(coreburning_times)) * N_tbins);
  if (coreburning_times == NULL)
  {
    fprintf(stderr, "Could not allocate memory for the grid to contains the coreburning times.\n");
    return EXIT_FAILURE;
  }

  for (bin_idx = 0; bin_idx < N_tbins; ++bin_idx)
  {
    double t = coreburning_tbins_low + ((double)bin_idx * coreburning_tbins_delta);
    coreburning_times[bin_idx] = exp10(a/log10(t) + b * exp(c/log10(t)) + d);
  }

  // Now need to know given a mass range (defined by the timestep) what is the number and mass fraction of SN that go supernova.

  m_IMFbins_low = 8.0; // The IMF range is from 8.0 to 120.0 Msun.
  m_IMFbins_high = 120.0;
  m_IMFbins_delta = 0.00001;

  N_massbins = (m_IMFbins_high - m_IMFbins_low) / (m_IMFbins_delta);

  IMF_massgrid_eta = malloc(sizeof(*(IMF_massgrid_eta)) * N_massbins);
  if (IMF_massgrid_eta == NULL)
  {
    fprintf(stderr, "Could not allocate memory for the grid to contain the IMF eta masses for delayed supernova.\n");
    return EXIT_FAILURE;
  }

  IMF_massgrid_m = malloc(sizeof(*(IMF_massgrid_m)) * N_massbins);
  if (IMF_massgrid_m == NULL)
  {
    fprintf(stderr, "Could not allocate memory for the grid to contain the IMF m masses for delayed supernova.\n");
    return EXIT_FAILURE;
  }

  for (bin_idx = 0; bin_idx < N_massbins; ++bin_idx)
  {
    IMF_massgrid_eta[bin_idx] = pow(m_IMFbins_low + ((double)bin_idx * m_IMFbins_delta), IMF_slope + 1.0);
    IMF_massgrid_m[bin_idx] = pow(m_IMFbins_low + ((double)bin_idx * m_IMFbins_delta), IMF_slope + 2.0);
  }

  return EXIT_SUCCESS;

}

int32_t init_nionlookup(void)
{

  // For PhotonPrescription == 1 we wish to explicitly track the stellar ages of a galaxy.
  // Then we determine the number of ionizing photons using the age of the stellar population.

  // Using STARBURST99 it was determined that the number of ionizing photons emitted from an instantaneous starburst depends on the
  // mass of stars formed and the time since the starburst.  Furthermore, the number of ionizing photons scales linearly with the mass of stars formed.
  // That is, a starburst that forms 7.0e10Msun worth of stars will emit 10x as many photons as a starburst that forms 6.0e10Msun worth of stars.

  // So if we read in a table that contains the number of ionizing photons emitted from a starburst for a 6.0e10Msun episode, then we can scale our values to this lookup table
  // using log10 Ngamma(Msun, t) = (log10 M* - 6.0) + log10 Ngamma(6.0, t).

#define MAXBINS 10000

  char buf[MAX_STRING_LEN], fname[MAX_STRING_LEN];
  FILE *niontable;
  int32_t i = 0, num_lines = 0;
  float t, HI, HI_L, HeI, HeI_L, HeII, HeII_L, L;

  snprintf(fname, MAX_STRING_LEN - 1, ROOT_DIR "/extra/nion_table.txt");
  niontable = fopen(fname, "r");
  if (niontable == NULL)
  {
    fprintf(stderr, "Could not open file %s\n", fname);
    return EXIT_FAILURE;
  }

  stars_tbins = calloc(MAXBINS, sizeof(*(stars_tbins)));
  if (stars_tbins == NULL)
  {
    fprintf(stderr, "Could not allocate memory for the time bins for the tracking of stellar populations.\n");
    return EXIT_FAILURE;
  }

  stars_Ngamma = calloc(MAXBINS, sizeof(*(stars_Ngamma)));
  if (stars_Ngamma == NULL)
  {
    fprintf(stderr, "Could not allocate memory for the Ngamma HI bins for the tracking of stellar populations.\n");
    return EXIT_FAILURE;
  }

  // The first 8 lines are all header information. So move past them.
  while (i < 8)
  {
    fgets(buf, MAX_STRING_LEN, niontable);
    ++i;
  }

  while (fscanf(niontable, "%f %f %f %f %f %f %f %f", &t, &HI, &HI_L, &HeI, &HeI_L, &HeII, &HeII_L, &L) == 8)
  {

    stars_tbins[num_lines] = t;
    stars_Ngamma[num_lines] = HI;

    ++num_lines;
    if (num_lines == MAXBINS - 1)
    {
      fprintf(stderr, "Exceeding the maximum bins for the tracking of stellar populations.\n");
      return EXIT_FAILURE;
    }
  }
  fclose(niontable);

  // Check that the Nion lookup table had enough datapoints to cover the time we're tracking the ages for.
  if (stars_tbins[num_lines - 1] / 1.0e6 < STELLAR_TRACKING_TIME)
  {
    fprintf(stderr, "The final time specified in the Nion lookup table is %.4f Myr. However we specified to track stellar ages over %d Myr.\n", stars_tbins[num_lines - 1] / 1.0e6, STELLAR_TRACKING_TIME);
    fprintf(stderr, "Either update the Nion lookup table or reduce the value of STELLAR_TRACKING_TIME in `core_allvars.h`.\n");
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;

#undef MAXBINS
}

void set_units(void)
{

  UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s;
  UnitTime_in_Megayears = UnitTime_in_s / SEC_PER_MEGAYEAR;
  UnitDensity_in_cgs = UnitMass_in_g / pow(UnitLength_in_cm, 3);
  UnitPressure_in_cgs = UnitMass_in_g / UnitLength_in_cm / pow(UnitTime_in_s, 2);
  UnitCoolingRate_in_cgs = UnitPressure_in_cgs / UnitTime_in_s;
  UnitEnergy_in_cgs = UnitMass_in_g * pow(UnitLength_in_cm, 2) / pow(UnitTime_in_s, 2);

  EnergySNcode = EnergySN / UnitEnergy_in_cgs * Hubble_h;

  sage_G = GRAVITY / pow(UnitLength_in_cm, 3) * UnitMass_in_g * pow(UnitTime_in_s, 2);

  // convert some physical input parameters to internal units
  sage_Hubble = HUBBLE * UnitTime_in_s;

  // compute a few quantitites
  RhoCrit = 3 * sage_Hubble * sage_Hubble / (8 * M_PI * sage_G);
}

void read_snap_list(void)
{
  FILE *fd;
  char fname[1000];

  sprintf(fname, "%s", FileWithSnapList);

  if(!(fd = fopen(fname, "r")))
  {
    printf("can't read output list in file '%s'\n", fname);
    ABORT(0);
  }

  Snaplistlen = 0;
  do
  {
    if(fscanf(fd, " %lg ", &AA[Snaplistlen]) == 1)
      Snaplistlen++;
    else
      break;
  }
  while(Snaplistlen < MAXSNAPS);

  fclose(fd);
}

double time_to_present(double z)
{
#define WORKSIZE 1000
  gsl_function F;
  gsl_integration_workspace *workspace;
  double time, result, abserr;

  workspace = gsl_integration_workspace_alloc(WORKSIZE);
  F.function = &integrand_time_to_present;

  gsl_integration_qag(&F, 1.0 / (z + 1), 1.0, 1.0 / sage_Hubble,
    1.0e-8, WORKSIZE, GSL_INTEG_GAUSS21, workspace, &result, &abserr);

  time = 1 / sage_Hubble * result;

  gsl_integration_workspace_free(workspace);

  // return time to present as a function of redshift
  return time;
}

double integrand_time_to_present(double a, void *param)
{
  (void)(param); // Avoids triggering unused-parameter warning.
  return 1 / sqrt(Omega / a + (1 - Omega - OmegaLambda) + OmegaLambda * a * a);
}

int32_t init_fesc(void)
{
  switch(fescPrescription)
  {

    case 0:
      if (ThisTask == 0)
      {
        printf("\n\nUsing a constant escape fraction of %.4f\n", beta);
      }
      break;

    case 1:
      if (ThisTask == 0)
      {
        printf("\n\nUsing an fesc prescription that scales with the fraction of ejected mass in the galaxy.\nThis takes the form A*fej + B with A = %.4e and B = %.4e\n", alpha, beta);
      }
      break;

    case 2:
      if (ThisTask == 0)
      {
        printf("\n\nUsing an fesc prescription that depends upon quasar activity.\n");
        printf("\nFor a galaxy that had a quasar event within %.2f dynamical times go the escape fraction will be %.2f.  Otherwise it will have a constant value of %.2f\n", N_dyntime, quasar_boosted, quasar_baseline);
      }
      break;

    case 3:
      if (ThisTask == 0)
      {
        printf("\n\nUsing Anne's functional form for an escape fraction that increases for descreasing halo mass.\n");
        printf("MH_low = %.4e\tMH_high = %.4e\tfesc_low = %.4f\tfesc_high =  %.4f.\n", MH_low, MH_high, fesc_low, fesc_high);
      }
      XASSERT(fesc_low > fesc_high, "Input file contain fesc_low = %.2f and fesc_high = %.2f. For this prescription (fescPrescription == 3), we require fesc_low > fesc_high\n", fesc_low, fesc_high);

      break;

    case 4:
      if (ThisTask == 0)
      {
        printf("\n\nUsing Anne's functional form for an escape fraction that increases for increasing halo mass.\n");
        printf("MH_low = %.4e\tMH_high = %.4e\tfesc_low = %.4f\tfesc_high =  %.4f.\n", MH_low, MH_high, fesc_low, fesc_high);
      }
      XASSERT(fesc_low < fesc_high, "Input file contain fesc_low = %.2f and fesc_high = %.2f. For this prescription (fescPrescription == 4), we require fesc_low < fesc_high\n", fesc_low, fesc_high);
      break;

    case 5:
      if (ThisTask == 0)
      {
        printf("\n\nUsing an fesc prescription that depends upon the SFR of the galaxy.\n");
        printf("This takes the function form of a logistic function, fesc = delta / (1.0 + exp(-alpha*(log10(SFR)-beta))), with range [0.0, delta].  `alpha` controls the steepness of the curve and beta defines the value of log10(SFR) gives fesc = delta/2.\n");
        printf("alpha = %.4f\tbeta = %.4f\tdelta = %.4f\n", alpha, beta, delta);
      }
      break;

    default:
      printf("\n\nOnly escape fraction prescriptions 0, 1, 2, 3, 4 or 5 are permitted.\n");
      return EXIT_FAILURE;
  }
  return EXIT_SUCCESS;
}

/*
For some escape fraction prescriptions, the functional form is a power law with the constants calculated using two fixed points.
This function determines these constants based on the fixed points specified.
Refer to the .ini file or `determine_fesc()` function for full details on each `fescPrescription`.

Parameters
----------
None.  All variables used/adjusted are global.

Returns
----------
None. All variables adjusted are global
.
Pointer Updates
----------
None.

Units
----------
The Halo Masses used to specify the fixed points (MH_low and MH_high) are in Msun.
*/

int32_t determine_fescMH_constants(void)
{

  double A, B, log_A;

  switch(fescPrescription)
  {
    case 3:
      log_A = (log10(fesc_high) - (log10(fesc_low)*log10(MH_high)/log10(MH_low))) * pow(1 - (log10(MH_high) / log10(MH_low)), -1);
      B = (log10(fesc_low) - log_A) / log10(MH_low);
      A = pow(10, log_A);

      fescMH_alpha = A;
      fescMH_beta = B;

      printf("Fixing the points (%.4e, %.2f) and (%.4e, %.2f)\n", MH_low, fesc_low, MH_high, fesc_high);
      printf("This gives a power law with constants A = %.4e, B = %.4e\n", fescMH_alpha, fescMH_beta);

      break;

  }

  return EXIT_SUCCESS;

}

ext: c

max_stars_repo_name: jacobseiler/rsage

max_issues_repo_name: jacobseiler/rsage

max_forks_repo_name: jacobseiler/rsage


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

/* Author:  G. Jungman */

/* Miscellaneous implementations of use
 * for evaluation of hypergeometric functions.
 */
#ifndef _HYPERG_H_
#define _HYPERG_H_

#include <gsl/gsl_sf_result.h>

/* Direct implementation of 1F1 series.
 */
int
gsl_sf_hyperg_1F1_series_e(const double a, const double b, const double x, gsl_sf_result * result);

/* Implementation of the 1F1 related to the
 * incomplete gamma function: 1F1(1,b,x), b >= 1.
 */
int
gsl_sf_hyperg_1F1_1_e(double b, double x, gsl_sf_result * result);

/* 1F1(1,b,x) for integer b >= 1
 */
int
gsl_sf_hyperg_1F1_1_int_e(int b, double x, gsl_sf_result * result);

/* Implementation of large b asymptotic.
 * [Bateman v. I, 6.13.3 (18)]
 * [Luke, The Special Functions and Their Approximations v. I, p. 129, 4.8 (4)]
 *
 * a^2 << b, |x|/|b| < 1 - delta
 */
int
gsl_sf_hyperg_1F1_large_b_e(const double a, const double b, const double x, gsl_sf_result * result);

/* Implementation of large b asymptotic.
 *
 * Assumes a > 0 is small, x > 0, and |x|<|b|.
 */
int
gsl_sf_hyperg_U_large_b_e(const double a, const double b, const double x,
                             gsl_sf_result * result,
                             double * ln_multiplier
                             );

/* Implementation of 2F0 asymptotic series.
 */
int
gsl_sf_hyperg_2F0_series_e(const double a, const double b, const double x, int n_trunc,
                              gsl_sf_result * result);

#endif /* !_HYPERG_H_ */

ext: h

max_stars_repo_name: vinej/sml

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_qrng.h>
#include <gsl/gsl_errno.h>

/* Setup */
gsl_qrng *mgsl_qrng_alloc(int type, unsigned int d)
{
  const gsl_qrng_type *types[] = { gsl_qrng_niederreiter_2, gsl_qrng_sobol, gsl_qrng_halton, gsl_qrng_reversehalton };
  const gsl_qrng_type *T;
  gsl_qrng *q;
  T = types[type];
  q = gsl_qrng_alloc(T, d);
  return q;
}

max_stars_repo_name: frithnanth/raku-Math-Libgsl-QuasiRandom

max_issues_repo_name: frithnanth/raku-Math-Libgsl-QuasiRandom

max_forks_repo_name: frithnanth/raku-Math-Libgsl-QuasiRandom


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C code for example generation of waveforms with the EOBNRv2HM reduced order model,
 * processing through the Fourier-domain LISA response and SNR calculation.
 *
 */

#define _XOPEN_SOURCE 500

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

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

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

#include "constants.h"
#include "struct.h"
#include "EOBNRv2HMROMstruct.h"
#include "EOBNRv2HMROM.h"
#include "wip.h"
#include "LISAgeometry.h"
#include "LISAFDresponse.h"

/* Parameters for the generation of a ROM waveform (in the form of a list of modes) */
/* All parameters are to be given in SI units! */
typedef struct tagROMParams {
  int nbmode;                /* Number of modes to generate (starting with the 22) - defaults to 1 (22 mode only) */
  double tRef;               /* shift in time with respect to the 22-fit-removed waveform */
  double phiRef;             /* phase at fRef */
  double fRef;               /* reference frequency */
  double m1;                 /* mass of companion 1 */
  double m2;                 /* mass of companion 2 */
  double distance;           /* distance of source */
  char outname[256];         /* file to which output should be written */
} ROMParams;

/* Parameters for the generation of a ROM waveform (in the form of a list of modes) */
/* All parameters are in SI units in the internals */
/* Angle definitions are taken from the Krolak&al paper gr-qc/0401108 */
typedef struct tagLISAParams {
  int nbmode;                /* Number of modes to generate (starting with the 22) - defaults to 1 (22 mode only) */
  double tRef;               /* shift in time with respect to the 22-fit-removed waveform */
  double phiRef;             /* phase at fRef */
  double fRef;               /* reference frequency */
  double m1;                 /* mass of companion 1 */
  double m2;                 /* mass of companion 2 */
  double distance;           /* distance of source */
  double inclination;        /* inclination of L relative to line of sight */
  double lambda;             /* First angle for the position in the sky of the source */
  double beta;               /* Second angle for the position in the sky of the source */
  double psi;                /* Polarization angle */
  char outname[256];         /* file to which output should be written */
} LISAParams;

/* Parse command line and return a newly allocated ROMParams object
 * Masses are input in solar masses and distances in Mpc - converted in SI for the internals */
static ROMParams* parse_args_ROM(ssize_t argc, char **argv) {
    ssize_t i;
    ROMParams* params;
    params = (ROMParams*) malloc(sizeof(ROMParams));
    memset(params, 0, sizeof(ROMParams));

    /* Set default values to the arguments */
    params->nbmode = 1;
    params->tRef = 0.;
    params->phiRef = 0.;
    params->fRef = 0.;
    params->m1 = 1. * 1e6 * MSUN_SI;
    params->m2 = 1. * 1e6 * MSUN_SI;
    params->distance = 1. * 1e9 * PC_SI;

    /* consume command line */
    for (i = 1; i < argc; ++i) {
        if (strcmp(argv[i], "--nbmode") == 0) {
            params->nbmode = atof(argv[++i]);
        } else if (strcmp(argv[i], "--tRef") == 0) {
            params->tRef = atof(argv[++i]);
        } else if (strcmp(argv[i], "--phiRef") == 0) {
            params->phiRef = atof(argv[++i]);
        } else if (strcmp(argv[i], "--fRef") == 0) {
            params->fRef = atof(argv[++i]);
        } else if (strcmp(argv[i], "--m1") == 0) {
            params->m1 = atof(argv[++i]) * MSUN_SI;
        } else if (strcmp(argv[i], "--m2") == 0) {
            params->m2 = atof(argv[++i]) * MSUN_SI;
        } else if (strcmp(argv[i], "--distance") == 0) {
            params->distance = atof(argv[++i]) * 1e6 * PC_SI;
        } else if (strcmp(argv[i], "--outname") == 0) {
            strncpy(params->outname, argv[++i], 256);
        } else {
            printf("Error: invalid option: %s\n", argv[i]);
            goto fail;
        }
    }

    return params;

    fail:
    free(params);
    exit(1);
}

/* Parse command line and return a newly allocated LISAParams object
 * Masses are input in solar masses and distances in Mpc - converted in SI for the internals */
static LISAParams* parse_args_LISA(ssize_t argc, char **argv) {
    ssize_t i;
    LISAParams* params;
    params = (LISAParams*) malloc(sizeof(LISAParams));
    memset(params, 0, sizeof(LISAParams));

    /* Set default values to the arguments */
    params->nbmode = 1;
    params->tRef = 0.;
    params->phiRef = 0.;
    params->fRef = 0.;
    params->m1 = 1. * 1e6 * MSUN_SI;
    params->m2 = 1. * 1e6 * MSUN_SI;
    params->distance = 1. * 1e9 * PC_SI;
    params->inclination = 0.;
    params->lambda = 0.;
    params->beta = 0.;
    params->psi = 0.;
    sprintf(params->outname, "");

    /* consume command line */
    for (i = 1; i < argc; ++i) {
        if (strcmp(argv[i], "--nbmode") == 0) {
            params->nbmode = atof(argv[++i]);
        } else if (strcmp(argv[i], "--tRef") == 0) {
            params->tRef = atof(argv[++i]);
        } else if (strcmp(argv[i], "--phiRef") == 0) {
            params->phiRef = atof(argv[++i]);
        } else if (strcmp(argv[i], "--fRef") == 0) {
            params->fRef = atof(argv[++i]);
        } else if (strcmp(argv[i], "--m1") == 0) {
            params->m1 = atof(argv[++i]) * MSUN_SI;
        } else if (strcmp(argv[i], "--m2") == 0) {
            params->m2 = atof(argv[++i]) * MSUN_SI;
        } else if (strcmp(argv[i], "--distance") == 0) {
            params->distance = atof(argv[++i]) * 1e6 * PC_SI;
        } else if (strcmp(argv[i], "--inclination") == 0) {
            params->inclination = atof(argv[++i]);
        } else if (strcmp(argv[i], "--lambda") == 0) {
            params->lambda = atof(argv[++i]);
        } else if (strcmp(argv[i], "--beta") == 0) {
            params->beta = atof(argv[++i]);
        } else if (strcmp(argv[i], "--psi") == 0) {
            params->psi = atof(argv[++i]);
        } else if (strcmp(argv[i], "--outname") == 0) {
            strncpy(params->outname, argv[++i], 256);
        } else {
            printf("Error: invalid option: %s\n", argv[i]);
            goto fail;
        }
    }

    return params;

    fail:
    free(params);
    exit(1);
}

/* Function to output to a file an AmpPhaseFrequencySeries; each mode will be output to a separate file using this function */
static int Write_CAmpPhaseFrequencySeries(FILE* f, CAmpPhaseFrequencySeries* freqseries) {
    gsl_vector* freq = freqseries->freq;
    gsl_vector* amp_real = freqseries->amp_real;
    gsl_vector* amp_imag = freqseries->amp_imag;
    gsl_vector* phase = freqseries->phase;

    int len = (int) freq->size;
    /*Here, we could add a check on the length of the gsl_vectors*/

    fprintf(f, "# f amp_re amp_im phase\n");
    for (int i=0; i<len; i++) {
      fprintf(f, "%.16e %.16e %.16e %.16e\n", gsl_vector_get(freq, i), gsl_vector_get(amp_real, i), gsl_vector_get(amp_imag, i), gsl_vector_get(phase, i));
    }

    return 0;
}

/*
 * Swhitenoise
 * a simple noise function
 */
double Swhitenoise(double x){
  return 1;
};

/*
 * main
 */
int main (int argc , char **argv) {
    FILE* f;
    ListmodesCAmpPhaseFrequencySeries* listROM = NULL;
    ListmodesCAmpPhaseFrequencySeries* listA = NULL;
    ListmodesCAmpPhaseFrequencySeries* listE = NULL;
    ListmodesCAmpPhaseFrequencySeries* listT = NULL;
    //ROMParams* params;
    LISAParams* params;

    /* parse commandline */
    //params = parse_args_ROM(argc, argv);
    params = parse_args_LISA(argc, argv);

    /* Generate the waveform with the ROM */
    SimEOBNRv2HMROM(&listROM, params->nbmode, params->tRef, params->phiRef, params->fRef, params->m1, params->m2, params->distance);
    /* Process the waveform through the LISA response */
    LISASimFDResponseTDI(&listROM, &listA, &listE, &listT, params->inclination, params->lambda, params->beta, params->psi);

    ListmodesCAmpPhaseFrequencySeries* listelementA;
    ListmodesCAmpPhaseFrequencySeries* listelementE;
    ListmodesCAmpPhaseFrequencySeries* listelementT;
    int l,m;
    int status = 0;
    /* Loop over the modes */
    for(int i=0; i<params->nbmode; i++){
      l = listmode[i][0];
      m = listmode[i][1];
      listelementA = ListmodesCAmpPhaseFrequencySeries_GetMode(listA, l, m);
      listelementE = ListmodesCAmpPhaseFrequencySeries_GetMode(listE, l, m);
      listelementT = ListmodesCAmpPhaseFrequencySeries_GetMode(listT, l, m);
      /* Write files - suffix _A,E,T_lm.dat imposed to each file name, attached to the string given by outname */
      /* If outname is still the default empty string, we do not output and skip this stage */
      if(!(strcmp(params->outname, "") == 0)){
	char *filenameA = malloc(strlen(params->outname)+64);
	char *filenameE = malloc(strlen(params->outname)+64);
	char *filenameT = malloc(strlen(params->outname)+64);
	sprintf(filenameA, "%s%s%d%d%s", params->outname, "_A_", l, m, ".dat");
	sprintf(filenameE, "%s%s%d%d%s", params->outname, "_E_", l, m, ".dat");
	sprintf(filenameT, "%s%s%d%d%s", params->outname, "_T_", l, m, ".dat");
	f = fopen(filenameA, "w");
	status |= Write_CAmpPhaseFrequencySeries(f, listelementA->freqseries);
	fclose(f);
	f = fopen(filenameE, "w");
	status |= Write_CAmpPhaseFrequencySeries(f, listelementE->freqseries);
	fclose(f);
	f = fopen(filenameT, "w");
	status |= Write_CAmpPhaseFrequencySeries(f, listelementT->freqseries);
	fclose(f);
	if (status) goto fail;
      }

      /* Example SNR calculation, taking the A observable */
      double *f1=listelementA->freqseries->freq->data;
      int n1=listelementA->freqseries->freq->size;
      double *s1Ar=listelementA->freqseries->amp_real->data;
      double *s1Ai=listelementA->freqseries->amp_imag->data;
      double *s1p=listelementA->freqseries->phase->data;
      double *f2=listelementA->freqseries->freq->data;
      int n2=listelementA->freqseries->freq->size;
      double *s2Ar=listelementA->freqseries->amp_real->data;
      double *s2Ai=listelementA->freqseries->amp_imag->data;
      double *s2p=listelementA->freqseries->phase->data;

      printf("n1: %d\n", n1);
      printf("s1Ar[299]: %g\n", s1Ar[299]);

      double start=((double)clock())/CLOCKS_PER_SEC;
      //JGB: Set these somehwere else where it makes sense.  f_min or f_max <= 0 means use intersection of signal domains.
      double f_min=-1.0;
      double f_max=-1.0;
      double rho2= wip_phase (f1, n1, f2, n2, s1Ar, s1Ai, s1p, s2Ar, s2Ai, s2p, Swhitenoise, 1.0, f_min, f_max);
      double end=((double)clock())/CLOCKS_PER_SEC;
      printf( "SNR2 = %g,   SNR time = %g\n", rho2, end-start);
      printf( "SNR = %g,   SNR time = %g\n",sqrt(rho2), end-start);

    }

    /* clean up */
    free(params);
    ListmodesCAmpPhaseFrequencySeries_Destroy(listROM);
    ListmodesCAmpPhaseFrequencySeries_Destroy(listA);
    ListmodesCAmpPhaseFrequencySeries_Destroy(listE);
    ListmodesCAmpPhaseFrequencySeries_Destroy(listT);
    return 0;

    fail:
    free(params);
    ListmodesCAmpPhaseFrequencySeries_Destroy(listROM);
    ListmodesCAmpPhaseFrequencySeries_Destroy(listA);
    ListmodesCAmpPhaseFrequencySeries_Destroy(listE);
    ListmodesCAmpPhaseFrequencySeries_Destroy(listT);
    return 1;
}

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
/*
 * Partition.h
 * Provides class prototype that simplifies partition vectors.
 *
   Copyright 2017 Dariusz Kuchta

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

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

#pragma once

#include <vector>
#include <gsl.h>

namespace spectre::unsupervised
{
/// <summary>
/// Class storing simplified representation of partition vectors.
/// </summary>
class Partition
{
public:

    /// <summary>
    /// Default constructor. Deleted so no dataless instance can be created.
    /// </summary>
    Partition() = delete;

    /// <summary>
    /// Constructor taking partition data and computing its simplified version.
    /// </summary>
    /// <param name="partition">Input partition data.</param>
    explicit Partition(gsl::span<unsigned int> partition);

    /// <summary>
    /// Getter for simplified partition data.
    /// </summary>
    /// <returns>Stored simplified partition.</returns>
    const std::vector<unsigned int>& Get() const;

    /// <summary>
    /// Method for comparing specified partitions.
    /// </summary>
    /// <param name="lhs">The first partition.</param>
    /// <param name="rhs">The second partition.</param>
    /// <param name="tolerance">The tolerance rate of mismatch.</param>
    static bool Compare(const Partition &lhs, const Partition &rhs, double tolerance = 0);

    /// <summary>
    /// Equality operator. Compares with tolerance = 0.
    /// </summary>
    /// <param name="other">The second partition.</param>
    bool operator==(const Partition &other) const;

    /// <summary>
    /// Inequality operator. Compares with tolerance = 0.
    /// </summary>
    /// <param name="other">The second partition.</param>
    bool operator!=(const Partition &other) const;

private:
    std::vector<unsigned int> m_Partition;
};
}

max_stars_repo_name: spectre-team/native-algorithms

max_issues_repo_name: spectre-team/native-algorithms

max_forks_repo_name: spectre-team/native-algorithms


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

/* goldensection.c -- goldensection minimum finding algorithm */

#include <config.h>

#include <stddef.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <float.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_min.h>

#include "min.h"

typedef struct
  {
    double dummy;
  }
goldensection_state_t;

static int goldensection_init (void * vstate, gsl_function * f, double x_minimum, double f_minimum, double x_lower, double f_lower, double x_upper, double f_upper);
static int goldensection_iterate (void * vstate, gsl_function * f, double * x_minimum, double * f_minimum, double * x_lower, double * f_lower, double * x_upper, double * f_upper);

static int
goldensection_init (void * vstate, gsl_function * f, double x_minimum, double f_minimum, double x_lower, double f_lower, double x_upper, double f_upper)
{
  goldensection_state_t * state = (goldensection_state_t *) vstate;

  /* no initialization required, prevent warnings about unused variables */

  state = 0;
  f = 0;
  x_minimum = 0;
  f_minimum = 0;
  x_lower = 0;
  f_lower = 0;
  x_upper = 0;
  f_upper = 0;

  return GSL_SUCCESS;
}

static int
goldensection_iterate (void * vstate, gsl_function * f, double * x_minimum, double * f_minimum, double * x_lower, double * f_lower, double * x_upper, double * f_upper)
{
  goldensection_state_t * state = (goldensection_state_t *) vstate;

  const double x_center = *x_minimum ;
  const double x_left = *x_lower ;
  const double x_right = *x_upper ;

  const double f_min = *f_minimum;

  const double golden = 0.3819660; /* golden = (3 - sqrt(5))/2 */

  const double w_lower = (x_center - x_left);
  const double w_upper = (x_right - x_center);

  double x_new, f_new;

  state = 0 ; /* avoid warning about unused parameters */

  x_new = x_center + golden * ((w_upper > w_lower) ? w_upper : -w_lower) ;

  SAFE_FUNC_CALL (f, x_new, &f_new);

  if (f_new < f_min)
    {
      *x_minimum = x_new ;
      *f_minimum = f_new ;
      return GSL_SUCCESS;
    }
  else if (x_new < x_center && f_new > f_min)
    {
      *x_lower = x_new ;
      *f_lower = f_new ;
      return GSL_SUCCESS;
    }
  else if (x_new > x_center && f_new > f_min)
    {
      *x_upper = x_new ;
      *f_upper = f_new ;
      return GSL_SUCCESS;
    }
  else
    {
      return GSL_FAILURE;
    }
}

static const gsl_min_fminimizer_type goldensection_type =
{"goldensection",                               /* name */
 sizeof (goldensection_state_t),
 &goldensection_init,
 &goldensection_iterate};

const gsl_min_fminimizer_type  * gsl_min_fminimizer_goldensection = &goldensection_type;

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef DIRICHLETSAMPLE_H
#define DIRICHLETSAMPLE_H

#include <stdio.h>
#include <iostream>
#include <vector>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <sys/time.h>
using namespace std;

std::vector<double> DirSample(std::vector<double> alpha, int k, int N);

#endif

max_stars_repo_name: KathrynLaing/CPNLearning

max_issues_repo_name: KathrynLaing/CPNLearning

max_forks_repo_name: KathrynLaing/CPNLearning


text:
/** @file matrix.h
 *  @author T J Atherton
 *
 *  @brief Veneer class over the objectmatrix type that interfaces with blas and lapack
 */

#ifndef matrix_h
#define matrix_h

#include <stdio.h>
#include "veneer.h"
/** Use Apple's Accelerate library for LAPACK and BLAS */
#ifdef __APPLE__
#ifdef MORPHO_LINALG_USE_ACCELERATE
#include <Accelerate/Accelerate.h>
#define MATRIX_LAPACK_PRESENT
#endif
#endif

/** Otherwise, use LAPACKE */
#ifndef MATRIX_LAPACK_PRESENT
#include <cblas.h>
#include <lapacke.h>
#define MORPHO_LINALG_USE_LAPACKE
#define MATRIX_LAPACK_PRESENT
#endif

/* -------------------------------------------------------
 * Matrix objects
 * ------------------------------------------------------- */

extern objecttype objectmatrixtype;
#define OBJECT_MATRIX objectmatrixtype

/** Matrices are a purely numerical collection type oriented toward linear algebra.
    Elements are stored in column-major format, i.e.
        [ 1 2 ]
        [ 3 4 ]
    is stored ( 1, 3, 2, 4 ) in memory. This is for compatibility with standard linear algebra packages */

typedef struct {
    object obj;
    unsigned int nrows;
    unsigned int ncols;
    double *elements;
    double matrixdata[];
} objectmatrix;

/** Tests whether an object is a matrix */
#define MORPHO_ISMATRIX(val) object_istype(val, OBJECT_MATRIX)

/** Gets the object as an matrix */
#define MORPHO_GETMATRIX(val)   ((objectmatrix *) MORPHO_GETOBJECT(val))

/** Creates a matrix object */
objectmatrix *object_newmatrix(unsigned int nrows, unsigned int ncols, bool zero);

/** Creates a new matrix from an array */
objectmatrix *object_matrixfromarray(objectarray *array);

/** Creates a new matrix from an existing matrix */
objectmatrix *object_clonematrix(objectmatrix *array);

/** @brief Use to create static matrices on the C stack
    @details Intended for small matrices; Caller needs to supply a double array of size nr*nc. */
#define MORPHO_STATICMATRIX(darray, nr, nc)      { .obj.type=OBJECT_MATRIX, .obj.status=OBJECT_ISUNMANAGED, .obj.next=NULL, .elements=darray, .nrows=nr, .ncols=nc }

/** Macro to decide if a matrix is 'small' or 'large' and hence static or dynamic allocation should be used. */
#define MATRIX_ISSMALL(m) (m->nrows*m->ncols<MORPHO_MAXIMUMSTACKALLOC)

/* -------------------------------------------------------
 * Matrix class
 * ------------------------------------------------------- */

#define MATRIX_CLASSNAME "Matrix"

#define MATRIX_TRANSPOSE_METHOD "transpose"
#define MATRIX_TRACE_METHOD "trace"
#define MATRIX_INNER_METHOD "inner"
#define MATRIX_DET_METHOD "det"
#define MATRIX_EIGENVALUES_METHOD "eigenvalues"
#define MATRIX_EIGENSYSTEM_METHOD "eigensystem"
#define MATRIX_NORM_METHOD "norm"
#define MATRIX_GETCOLUMN_METHOD "column"
#define MATRIX_SETCOLUMN_METHOD "setcolumn"

#define MATRIX_DIMENSIONS_METHOD "dimensions"

#define MATRIX_INDICESOUTSIDEBOUNDS       "MtrxBnds"
#define MATRIX_INDICESOUTSIDEBOUNDS_MSG   "Matrix index out of bounds."

#define MATRIX_INVLDINDICES               "MtrxInvldIndx"
#define MATRIX_INVLDINDICES_MSG           "Matrix indices must be integers."

#define MATRIX_INVLDNUMINDICES            "MtrxInvldNumIndx"
#define MATRIX_INVLDNUMINDICES_MSG        "Matrix expects two arguments for indexing."

#define MATRIX_CONSTRUCTOR                "MtrxCns"
#define MATRIX_CONSTRUCTOR_MSG            "Matrix() constructor should be called either with dimensions or an array, list or matrix initializer."

#define MATRIX_INVLDARRAYINIT             "MtrxInvldInit"
#define MATRIX_INVLDARRAYINIT_MSG         "Invalid initializer passed to Matrix()."

#define MATRIX_ARITHARGS                  "MtrxInvldArg"
#define MATRIX_ARITHARGS_MSG              "Matrix arithmetic methods expect a matrix or number as their argument."

#define MATRIX_INCOMPATIBLEMATRICES       "MtrxIncmptbl"
#define MATRIX_INCOMPATIBLEMATRICES_MSG   "Matrices have incompatible shape."

#define MATRIX_SINGULAR                   "MtrxSnglr"
#define MATRIX_SINGULAR_MSG               "Matrix is singular."

#define MATRIX_NOTSQ                      "MtrxNtSq"
#define MATRIX_NOTSQ_MSG                  "Matrix is not square."

#define MATRIX_SETCOLARGS                 "MtrxStClArgs"
#define MATRIX_SETCOLARGS_MSG             "Method setcolumn expects an integer column index and a column matrix as arguments."

/* -------------------------------------------------------
 * Matrix errors
 * ------------------------------------------------------- */

typedef enum { MATRIX_OK, MATRIX_INCMPTBLDIM, MATRIX_SING, MATRIX_INVLD, MATRIX_BNDS, MATRIX_NSQ, MATRIX_ALLOC } objectmatrixerror;

/* -------------------------------------------------------
 * Matrix interface
 * ------------------------------------------------------- */

bool matrix_getarraydimensions(objectarray *array, unsigned int dim[], unsigned int maxdim, unsigned int *ndim);
value matrix_getarrayelement(objectarray *array, unsigned int ndim, unsigned int *indx);

bool matrix_getlistdimensions(objectlist *list, unsigned int dim[], unsigned int maxdim, unsigned int *ndim);
bool matrix_getlistelement(objectlist *list, unsigned int ndim, unsigned int *indx, value *val);

bool matrix_setelement(objectmatrix *matrix, unsigned int row, unsigned int col, double value);
bool matrix_getelement(objectmatrix *matrix, unsigned int row, unsigned int col, double *value);

bool matrix_getcolumn(objectmatrix *matrix, unsigned int col, double **v);
bool matrix_setcolumn(objectmatrix *matrix, unsigned int col, double *v);
bool matrix_addtocolumn(objectmatrix *m, unsigned int col, double alpha, double *v);

objectmatrixerror matrix_copy(objectmatrix *a, objectmatrix *out);
objectmatrixerror matrix_add(objectmatrix *a, objectmatrix *b, objectmatrix *out);
objectmatrixerror matrix_accumulate(objectmatrix *a, double lambda, objectmatrix *b);
objectmatrixerror matrix_sub(objectmatrix *a, objectmatrix *b, objectmatrix *out);
objectmatrixerror matrix_mul(objectmatrix *a, objectmatrix *b, objectmatrix *out);
objectmatrixerror matrix_inner(objectmatrix *a, objectmatrix *b, double *out);
objectmatrixerror matrix_divs(objectmatrix *a, objectmatrix *b, objectmatrix *out);
objectmatrixerror matrix_divl(objectmatrix *a, objectmatrix *b, objectmatrix *out);
objectmatrixerror matrix_inverse(objectmatrix *a, objectmatrix *out);
objectmatrixerror matrix_transpose(objectmatrix *a, objectmatrix *out);
objectmatrixerror matrix_trace(objectmatrix *a, double *out);
objectmatrixerror matrix_scale(objectmatrix *a, double scale);
objectmatrixerror matrix_identity(objectmatrix *a);
double matrix_sum(objectmatrix *a);
//objectmatrixerror matrix_det(objectmatrix *a, double *out);
//objectmatrixerror matrix_eigensystem(objectmatrix *a, double *val, objectmatrix *vec);

void matrix_print(objectmatrix *m);

void matrix_initialize(void);

#endif /* matrix_h */

max_stars_repo_name: mattsep/morpho

max_issues_repo_name: mattsep/morpho

max_forks_repo_name: mattsep/morpho


text:
#pragma once

#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#else
#include "CustomSolver.h"
#endif
#include "ValueGrad.h"
#include "BooleanNodes.h"
#include "BooleanDAG.h"
#include "NodeVisitor.h"
#include "VarStore.h"
#include <map>
#include <set>
#include "SymbolicEvaluator.h"
#include "Interface.h"

#include <iostream>

using namespace std;

class Path {
public:
	map<int, int> nodeToVal;

	Path() { }
	void add(map<int, int>& initMap) {
		nodeToVal.insert(initMap.begin(), initMap.end());
	}
	void addCond(int nodeid, int val) {
		Assert(nodeToVal.find(nodeid) == nodeToVal.end(), "Node already has value");
		nodeToVal[nodeid] = val;
	}

	string toString() {
		stringstream s;
		for (auto it = nodeToVal.begin(); it != nodeToVal.end(); it++) {
			s << "(" << it->first << "," << it->second << ");";
		}
		return s.str();
	}

	int getVal(int nodeid) {
		if (nodeToVal.find(nodeid) == nodeToVal.end()) {
			return -1;
		} else {
			return nodeToVal[nodeid];
		}
	}

	void empty() {
		nodeToVal.clear();
	}

	int length() {
		return nodeToVal.size();
	}

	int lowestId() {
		if (nodeToVal.size() == 0) {
			return -1;
		}
		return nodeToVal.begin()->first;
	}

	static bool isCompatible(Path* p1, Path* p2) {
		int size1 = p1->nodeToVal.size();
		int size2 = p2->nodeToVal.size();
		if (size1 == 0 || size2 == 0) {
			return true;
		}
		Path* pl;
		Path* ps;
		if (size1 > size2) {
			pl = p1;
			ps = p2;
		} else {
			pl = p2;
			ps = p1;
		}
		for (auto it = pl->nodeToVal.begin(); it != pl->nodeToVal.end(); it++) {
			int nodeid = it->first;
			int val = it->second;
			int val1 = ps->getVal(nodeid);
			if (val1 != -1 && val != val1) {
				return false;
			}
		}
		return true;
	}

	static void combinePaths(Path* p1, Path* p2, Path* p) {
		// assumes that paths are compatible
		p->nodeToVal.clear();
		p->add(p1->nodeToVal);
		p->add(p2->nodeToVal);
	}

	static void combinePath(Path* p1, Path* p) {
		p->nodeToVal.clear();
		p->add(p1->nodeToVal);
	}

	static void combinePaths(vector<Path*>& paths, Path* p2, Path* p) {
		p->nodeToVal.clear();
		p->add(p2->nodeToVal);
		// TODO: need to add common conditions in paths to p
	}

	static Path* pIntersect(Path* p1, Path* p2) {
		// assumes that the paths are compatible
		Path* p = new Path();
		p->add(p1->nodeToVal);
		p->add(p2->nodeToVal);
		return p;
	}

	// only one cond  should differ
	static int isConjugate(Path* p1, Path* p2) {
		//cout << "Conjugate check: " << p1->toString() << " | " << p2->toString() << endl;
		if (p1->length() != p2->length()) {
			return -1;
		}
		int conjIdx = -1;
		int ctr = 0;
		for (auto it = p1->nodeToVal.begin(); it != p1->nodeToVal.end(); it++) {
			if (p2->getVal(it->first) == -1) {
				return -1;
			}
			if (p2->getVal(it->first) != it->second) {
				if (conjIdx == -1) {
					conjIdx = ctr;
				} else {
					return -1;
				}
			}
			ctr++;
		}
		//cout << conjIdx << endl;
		return conjIdx;
	}

	static Path* pUnion(Path* p1, Path* p2) {
		// assumes that the paths are conjugates
		Path* p = new Path();
		for (auto it = p1->nodeToVal.begin(); it != p1->nodeToVal.end(); it++) {
			if (p2->getVal(it->first) == it->second) {
				p->addCond(it->first, it->second);
			}

		}
		return p;

	}

	static bool isSubset(Path* p1, Path* p2) {
		if (p1->length() == 0) return true;
		for (auto it = p1->nodeToVal.begin(); it != p1->nodeToVal.end(); it++) {
			if (p2->getVal(it->first) != it->second) {
				return false;
			}
		}
		return true;

	}

	// no common condition
	static bool isDisjoint(Path* p1, Path* p2) {
		if (p1->length() == 0 || p2->length() == 0) return false;
		for (auto it = p1->nodeToVal.begin(); it != p1->nodeToVal.end(); it++) {
			if (p2->getVal(it->first) != -1) {
				return false;
			}
		}
		return true;
	}
};

class MergeNodesList {
	vector<tuple<int, int, int, int>> nodesToMerge;

public:
	void add(int node1, int idx1, int node2, int idx2) {
		nodesToMerge.push_back(make_tuple(node1, idx1, node2, idx2));
	}

	string print() {
		stringstream s;
		for (int i = 0; i < nodesToMerge.size(); i++) {
			auto& t = nodesToMerge[i];
			s << "(" << get<0>(t) << "," << get<1>(t) << "," << get<2>(t) << "," << get<3>(t) << ") ";
		}
		return s.str();
	}

	int size() {
		return nodesToMerge.size();
	}

	tuple<int, int, int, int>& getTuple(int idx) {
		return nodesToMerge[idx];
	}
};

class KLocalityAutoDiff: public NodeVisitor, public SymbolicEvaluator
{
	BooleanDAG& bdag;
	map<string, int>& floatCtrls; // Maps float ctrl names to indices within grad vector
	int nctrls; // number of float ctrls
	gsl_vector* ctrls; // ctrl values
	vector<vector<ValueGrad*>> values; // Keeps track of values along with gradients for each node
	vector<vector<DistanceGrad*>> distances; // Keeps track of distance metric for boolean nodes
	vector<vector<Path*>> paths;
	vector<int> sizes;
    Interface* inputValues;

    set<string> pathStrings;

    int MAX_REGIONS = 16 + 1; // 0th idx is for final merging

    vector<vector<MergeNodesList*>> mergeNodes;

public:
    int DEFAULT_INP = -32;
    KLocalityAutoDiff(BooleanDAG& bdag_p, map<string, int>& floatCtrls_p);

	~KLocalityAutoDiff(void);

	virtual void visit( SRC_node& node );
	virtual void visit( DST_node& node );
	virtual void visit( CTRL_node& node );
	virtual void visit( PLUS_node& node );
	virtual void visit( TIMES_node& node );
	virtual void visit( ARRACC_node& node );
	virtual void visit( DIV_node& node );
	virtual void visit( MOD_node& node );
	virtual void visit( NEG_node& node );
	virtual void visit( CONST_node& node );
	virtual void visit( LT_node& node );
	virtual void visit( EQ_node& node );
	virtual void visit( AND_node& node );
	virtual void visit( OR_node& node );
	virtual void visit( NOT_node& node );
	virtual void visit( ARRASS_node& node );
	virtual void visit( UFUN_node& node );
	virtual void visit( TUPLE_R_node& node );
	virtual void visit( ASSERT_node& node );

    virtual void setInputs(Interface* inputValues_p);

	virtual void run(const gsl_vector* ctrls_p, const set<int>& nodesSubset) {
		Assert(false, "Not yet supported");
	}
    virtual void run(const gsl_vector* ctrls_p);

    virtual double getErrorOnConstraint(int nodeid, gsl_vector* grad);
    double getSqrtError(bool_node* node, gsl_vector* grad);
    double getAssertError(bool_node* node, gsl_vector* grad);
    double getBoolCtrlError(bool_node* node, gsl_vector* grad);
    double getBoolExprError(bool_node* node, gsl_vector* grad);

    virtual double getErrorOnConstraint(int nodeid);
    double getSqrtError(bool_node* node);
    double getAssertError(bool_node* node);
    double getBoolCtrlError(bool_node* node);
    double getBoolExprError(bool_node* node);

    virtual double getErrorForAsserts(const set<int>& assertIds, gsl_vector* grad) {
    	Assert(false, "TODO");
    }

	virtual double getErrorForAssert(int assertId, gsl_vector* grad) {
		Assert(false, "TODO");
	}

	void setvalue(bool_node& bn, int idx, ValueGrad* v) {
		values[bn.id][idx] = v;
	}

	ValueGrad* v(bool_node& bn, int idx) {
		ValueGrad* val = values[bn.id][idx];
		if (val == NULL) {
			gsl_vector* g = gsl_vector_alloc(nctrls);
			GradUtil::default_grad(g);
			val = new ValueGrad(0, g);
			setvalue(bn, idx, val);
		}
		return val;
	}

	ValueGrad* v(bool_node* bn, int idx) {
		return v(*bn, idx);
	}

	void setdistance(bool_node& bn, int idx, DistanceGrad* d) {
		distances[bn.id][idx] = d;
	}

	DistanceGrad* d(bool_node& bn, int idx) {
		DistanceGrad* dist = distances[bn.id][idx];
		if (dist == NULL) {
			gsl_vector* g = gsl_vector_alloc(nctrls);
			GradUtil::default_grad(g);
			dist = new DistanceGrad(0, g);
			setdistance(bn, idx, dist);
		}
		return dist;
	}

	DistanceGrad* d(bool_node* bn, int idx) {
		return d(*bn, idx);
	}

	int size(bool_node& bn) {
		return sizes[bn.id];
	}

	int size(bool_node* bn) {
		return size(*bn);
	}

	void setsize(bool_node& bn, int sz) {
		sizes[bn.id] = sz;
	}

	void setsize(bool_node* n, int sz) {
		setsize(*n, sz);
	}

	void setpath(bool_node& bn, int idx, Path* p) {
		paths[bn.id][idx] = p;
	}

	Path* path(bool_node& bn, int idx) {
		Path* p = paths[bn.id][idx];
		if (p == NULL) {
			p = new Path();
			setpath(bn, idx, p);
		}
		return p;
	}

	Path* path(bool_node* bn, int idx) {
		return path(*bn, idx);
	}

	virtual double dist(int nid) {
		ValueGrad* val = v(bdag[nid], 0);
		return val->getVal();
	}

	virtual double dist(int nid, gsl_vector* grad) {
		ValueGrad* val = v(bdag[nid], 0);
		gsl_vector_memcpy(grad, val->getGrad());
		return val->getVal();
	}

	virtual void print() {
		/*for (int i = 0; i < bdag.size(); i++) {
			if (bdag[i]->type == bool_node::ASSERT) {
				DistanceGrad* dist = d(bdag[i]->mother);
				if (dist->set) {
				double gmag = gsl_blas_dnrm2(dist->grad);
				if (dist->dist < 0.1) {
					cout << bdag[i]->mother->lprint() << endl;
					cout << dist->printFull() << endl;
				}
				}
			}
		}
		return;*/
		for (int i = 0; i < bdag.size(); i++) {
			cout << bdag[i]->lprint() << " ";
			if (bdag[i]->getOtype() == OutType::FLOAT) {
				if (v(bdag[i], 0)->set) {
					cout << v(bdag[i], 0)->print() << endl;
				} else {
					cout << "UNSET" << endl;
				}
			} else {
				if (d(bdag[i], 0)->set) {
					cout << d(bdag[i], 0)->print() << endl;
				} else {
					cout << "UNSET" << endl;
				}
			}
		}
	}
	virtual void printFull() {
		for (int i = 0; i < bdag.size(); i++) {
			cout << bdag[i]->lprint() << endl;
			if (bdag[i]->getOtype() == OutType::FLOAT) {
				if (v(bdag[i],0)->set) {
					cout << v(bdag[i],0)->printFull() << endl;
				} else {
					cout << "UNSET" <<endl;
				}
			} else {
				if(d(bdag[i],0)->set) {
					cout << d(bdag[i],0)->printFull() << endl;
				} else {
					cout << "UNSET" << endl;
				}
			}
			if (bdag[i]->type == bool_node::PLUS) {
				cout << v(bdag[i]->mother(),0)->getVal() << " + " << v(bdag[i]->father(), 0)->getVal() << endl;
			}
		}
	}

	bool isFloat(bool_node& bn) {
		return (bn.getOtype() == OutType::FLOAT);
	}

	bool isFloat(bool_node* bn) {
		return (bn->getOtype() == OutType::FLOAT);
	}

	int getInputValue(bool_node& bn) {
		if (inputValues->hasValue(bn.id)) {
			int val = inputValues->getValue(bn.id);
			return val;
		} else {
			return DEFAULT_INP;
		}
	}

	int getInputValue(bool_node* bn) {
		return getInputValue(*bn);
	}

	void doUfun(UFUN_node& node, ValueGrad* mval, ValueGrad* val);

	void copyNodes(bool_node& node, bool_node* m);
	double merge(bool_node* n, gsl_vector* grad);

	void mergePaths(vector<tuple<double, vector<tuple<int, int, int>>, Path*>>& mergedPaths);
	void getMergeNodesBinop(bool_node* n);
	void getMergeNodesUnop(bool_node* n, bool_node* m);
	void getMergeNodesArracc(bool_node* node);
	void getMergeNodes();

	void combineDistance(bool_node& node, int nodeid1, int idx1, int nodeid2, int idx2, double& dist, gsl_vector* dg_grad);
	void doPair(bool_node& node, int nodeid1, int idx1, int nodeid2, int idx2, double& val, gsl_vector* val_grad);
	void pairNodes(bool_node& node);
	void normalize(bool_node* n);

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
#include <gsl/gsl_rng.h>

struct ising
{
    int Lx;                     // number of spins in x direction
    int Ly;                     // number of spins in y direction
    int N;                      // total number of spins
    int steps;                  // steps so far
    int **s;                    // the spins
    double J;                   // ferromagnetic coupling
    double T;                   // temperature
    double H;                   // magnetic field
    gsl_rng *r;                 // random number generator
    unsigned long seed;         // rng seed
    double w[25][3];            // Boltzmann factors
    int **cluster;              // spin cluster labels
    double addprobability;      // probability to add like spin
};


text:
#ifndef _WAVELET_H_
#define _WAVELET_H_

#include "gsl/gsl_wavelet.h"
//#include <gsl/gsl_wavelet2d.h>

// N-to-M conversion helper functions
template <typename T> T *** arrayto3D(T * data, size_t x, size_t y, size_t z);
template <typename T> T ** arrayto2D(T * data, size_t x, size_t y);
template <typename T> T * arrayto1D(T *** data, size_t x, size_t y, size_t z);
template <typename T> T * arrayto1D(T ** data, size_t x, size_t y);

// Helper functions to compute size of wavelet hierarchy at a level
static int binary_logn (const size_t n);
int binary_upto_logn (const size_t n);
int binary_upto_value (const size_t n, const int log_n);

// Debug output of discrete wavelet signal
void gsl_wavelet_print (const gsl_wavelet * w);

// Native stepping functions that apply signal to data
template <typename T> static void dwt_step (const gsl_wavelet * w, T *a, size_t stride, size_t n, gsl_wavelet_direction dir, gsl_wavelet_workspace * work);
template <typename T> static void dwt_step_proto(const gsl_wavelet * w, T *a, size_t stride, size_t n, gsl_wavelet_direction dir, gsl_wavelet_workspace * work);

// 1D-support
template <typename T> int wavelet_transform(const gsl_wavelet * w, T *&data, size_t stride, size_t n, gsl_wavelet_direction dir, gsl_wavelet_workspace * work);

// 2D-support
template <typename T> int wavelet2d_transform(const gsl_wavelet * w, T *data, size_t tda, size_t size1, size_t size2, gsl_wavelet_direction dir, gsl_wavelet_workspace * work);
template <typename T> int wavelet2d_nstransform(const gsl_wavelet * w, T *data, size_t tda, size_t size1, size_t size2, gsl_wavelet_direction dir, gsl_wavelet_workspace * work);

// 3D-support versions of (GSL) wavelet transform
template <typename T> int gsl_wavelet3d_transform (const gsl_wavelet * w, T *data, size_t tda, size_t size1, size_t size2,size_t size3, gsl_wavelet_direction dir, gsl_wavelet_workspace * work);
template <typename T> int gsl_wavelet3d_nstransform (const gsl_wavelet * w, T *data, size_t tda, size_t size1, size_t size2, size_t size3, gsl_wavelet_direction dir, gsl_wavelet_workspace * work);
template <typename T> int gsl_wavelet3d_nstransform (const gsl_wavelet * w, T *data, size_t tda, size_t size1, size_t size2, size_t size3, gsl_wavelet_direction dir, gsl_wavelet_workspace * work, int levels);
template <typename T> int gsl_wavelet3d_nstransform_proto(const gsl_wavelet * w, T *data, size_t tda, size_t size1, size_t size2, size_t size3, gsl_wavelet_direction dir, gsl_wavelet_workspace * work);

// Functions to compute dimensions for each coefficient block (quadrant) in a wavelet hierarchy
int *** getCoeffDims (int dimx, int dimy, int dimz);
int *** getCoeffDims (int dimx, int dimy);

// Set/Extract functions used in conjunction with above to extract coefficient subregions
template <typename T> T * getCutout3d(T * in, int * codims, int * dims);
template <typename T> T * getCutout2d(T * in, int * codims, int * dims);
template <typename T> void setCutout3d(T * in, T * cutout, int * codims, int * dims);
template <typename T> void setCutout2d(T * in, T * cutout, int * codims, int * dims);

template <typename T> int grow(int val, T *** data, size_t x, size_t y, size_t z);
template <typename T> int shrink(int val, T *** data, size_t x, size_t y, size_t z);

#endif

#include "wavelet.inl"

max_stars_repo_name: jpulidojr/VizAly-SNWPAC

max_issues_repo_name: jpulidojr/VizAly-SNWPAC

max_forks_repo_name: jpulidojr/VizAly-SNWPAC


text:
/*
 * Monte Carlo
 */

#include <gsl/gsl_rng.h>

double mc_pi(long m, gsl_rng *r);

double mc_pi(long m, gsl_rng *r)
{

  unsigned long seed = 1UL;
  double p[NP], mean, sd;
  int i;

/* allocate random number generator */
   r = gsl_rng_alloc (gsl_rng_taus2);

/* seed the random number generator */
   gsl_rng_set (r, seed);

  for (i = 0; i < NP; i++)
  {
  p[i] = gsl_rng_uniform (r);
  }

  gsl_rng_free (r);

  return 0;

   // for (n = 0; n < m; n++) {
      /* generate random coordinates */
     // x = gsl_rng_uniform(r);
      // y = gsl_rng_uniform(r);

      /* count when the point lands in the circle */
      // if ((x*x + y*y) <= 1.0) {
      //  count += 1;
  //    }
   // }

   // return(4.0 * ((double)count)/((double)m));

}


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

#ifndef __GSL_VECTOR_COMPLEX_FLOAT_H__
#define __GSL_VECTOR_COMPLEX_FLOAT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_float.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_block_complex_float.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  float *data;
  gsl_block_complex_float *block;
  int owner;
} gsl_vector_complex_float;

typedef struct
{
  gsl_vector_complex_float vector;
} _gsl_vector_complex_float_view;

typedef _gsl_vector_complex_float_view gsl_vector_complex_float_view;

typedef struct
{
  gsl_vector_complex_float vector;
} _gsl_vector_complex_float_const_view;

typedef const _gsl_vector_complex_float_const_view gsl_vector_complex_float_const_view;

/* Allocation */

GSL_FUN gsl_vector_complex_float *gsl_vector_complex_float_alloc (const size_t n);
GSL_FUN gsl_vector_complex_float *gsl_vector_complex_float_calloc (const size_t n);

GSL_FUN gsl_vector_complex_float *
gsl_vector_complex_float_alloc_from_block (gsl_block_complex_float * b,
                                           const size_t offset,
                                           const size_t n,
                                           const size_t stride);

GSL_FUN gsl_vector_complex_float *
gsl_vector_complex_float_alloc_from_vector (gsl_vector_complex_float * v,
                                             const size_t offset,
                                             const size_t n,
                                             const size_t stride);

GSL_FUN void gsl_vector_complex_float_free (gsl_vector_complex_float * v);

/* Views */

GSL_FUN _gsl_vector_complex_float_view
gsl_vector_complex_float_view_array (float *base,
                                     size_t n);

GSL_FUN _gsl_vector_complex_float_view
gsl_vector_complex_float_view_array_with_stride (float *base,
                                                 size_t stride,
                                                 size_t n);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_view_array (const float *base,
                                           size_t n);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_view_array_with_stride (const float *base,
                                                       size_t stride,
                                                       size_t n);

GSL_FUN _gsl_vector_complex_float_view
gsl_vector_complex_float_subvector (gsl_vector_complex_float *base,
                                         size_t i,
                                         size_t n);

GSL_FUN _gsl_vector_complex_float_view
gsl_vector_complex_float_subvector_with_stride (gsl_vector_complex_float *v,
                                                size_t i,
                                                size_t stride,
                                                size_t n);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_subvector (const gsl_vector_complex_float *base,
                                               size_t i,
                                               size_t n);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_subvector_with_stride (const gsl_vector_complex_float *v,
                                                      size_t i,
                                                      size_t stride,
                                                      size_t n);

GSL_FUN _gsl_vector_float_view
gsl_vector_complex_float_real (gsl_vector_complex_float *v);

GSL_FUN _gsl_vector_float_view
gsl_vector_complex_float_imag (gsl_vector_complex_float *v);

GSL_FUN _gsl_vector_float_const_view
gsl_vector_complex_float_const_real (const gsl_vector_complex_float *v);

GSL_FUN _gsl_vector_float_const_view
gsl_vector_complex_float_const_imag (const gsl_vector_complex_float *v);

/* Operations */

GSL_FUN void gsl_vector_complex_float_set_zero (gsl_vector_complex_float * v);
GSL_FUN void gsl_vector_complex_float_set_all (gsl_vector_complex_float * v,
                                       gsl_complex_float z);
GSL_FUN int gsl_vector_complex_float_set_basis (gsl_vector_complex_float * v, size_t i);

GSL_FUN int gsl_vector_complex_float_fread (FILE * stream,
                                    gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_fwrite (FILE * stream,
                                     const gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_fscanf (FILE * stream,
                                     gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_fprintf (FILE * stream,
                                      const gsl_vector_complex_float * v,
                                      const char *format);

GSL_FUN int gsl_vector_complex_float_memcpy (gsl_vector_complex_float * dest, const gsl_vector_complex_float * src);

GSL_FUN int gsl_vector_complex_float_reverse (gsl_vector_complex_float * v);

GSL_FUN int gsl_vector_complex_float_swap (gsl_vector_complex_float * v, gsl_vector_complex_float * w);
GSL_FUN int gsl_vector_complex_float_swap_elements (gsl_vector_complex_float * v, const size_t i, const size_t j);

GSL_FUN int gsl_vector_complex_float_equal (const gsl_vector_complex_float * u,
                                    const gsl_vector_complex_float * v);

GSL_FUN int gsl_vector_complex_float_isnull (const gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_ispos (const gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_isneg (const gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_isnonneg (const gsl_vector_complex_float * v);

GSL_FUN int gsl_vector_complex_float_add (gsl_vector_complex_float * a, const gsl_vector_complex_float * b);
GSL_FUN int gsl_vector_complex_float_sub (gsl_vector_complex_float * a, const gsl_vector_complex_float * b);
GSL_FUN int gsl_vector_complex_float_mul (gsl_vector_complex_float * a, const gsl_vector_complex_float * b);
GSL_FUN int gsl_vector_complex_float_div (gsl_vector_complex_float * a, const gsl_vector_complex_float * b);
GSL_FUN int gsl_vector_complex_float_scale (gsl_vector_complex_float * a, const gsl_complex_float x);
GSL_FUN int gsl_vector_complex_float_add_constant (gsl_vector_complex_float * a, const gsl_complex_float x);
GSL_FUN int gsl_vector_complex_float_axpby (const gsl_complex_float alpha, const gsl_vector_complex_float * x, const gsl_complex_float beta, gsl_vector_complex_float * y);

GSL_FUN INLINE_DECL gsl_complex_float gsl_vector_complex_float_get (const gsl_vector_complex_float * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_complex_float_set (gsl_vector_complex_float * v, const size_t i, gsl_complex_float z);
GSL_FUN INLINE_DECL gsl_complex_float *gsl_vector_complex_float_ptr (gsl_vector_complex_float * v, const size_t i);
GSL_FUN INLINE_DECL const gsl_complex_float *gsl_vector_complex_float_const_ptr (const gsl_vector_complex_float * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
gsl_complex_float
gsl_vector_complex_float_get (const gsl_vector_complex_float * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      gsl_complex_float zero = {{0, 0}};
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, zero);
    }
#endif
  return *GSL_COMPLEX_FLOAT_AT (v, i);
}

INLINE_FUN
void
gsl_vector_complex_float_set (gsl_vector_complex_float * v,
                              const size_t i, gsl_complex_float z)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  *GSL_COMPLEX_FLOAT_AT (v, i) = z;
}

INLINE_FUN
gsl_complex_float *
gsl_vector_complex_float_ptr (gsl_vector_complex_float * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_FLOAT_AT (v, i);
}

INLINE_FUN
const gsl_complex_float *
gsl_vector_complex_float_const_ptr (const gsl_vector_complex_float * v,
                                    const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_FLOAT_AT (v, i);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_COMPLEX_FLOAT_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA