text:
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include <cconfigspace.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

#define NUM_SAMPLES 10000

static void test_create_normal_distribution() {
	ccs_distribution_t      distrib = NULL;
	ccs_result_t            err = CCS_SUCCESS;
	int32_t                 refcount;
	ccs_object_type_t       otype;
	ccs_distribution_type_t dtype;
	ccs_scale_type_t        stype;
	ccs_numeric_type_t      data_type;
	ccs_numeric_t           quantization;
	ccs_float_t             mu, sigma;
	ccs_interval_t          interval;

	err = ccs_create_normal_distribution(
		CCS_NUM_FLOAT,
		1.0,
		2.0,
		CCS_LINEAR,
		CCSF(0.0),
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_object_get_type(distrib, &otype);
	assert( err == CCS_SUCCESS );
	assert( otype == CCS_DISTRIBUTION );

	err = ccs_distribution_get_type(distrib, &dtype);
	assert( err == CCS_SUCCESS );
	assert( dtype == CCS_NORMAL );

	err = ccs_distribution_get_data_types(distrib, &data_type);
	assert( err == CCS_SUCCESS );
	assert( data_type == CCS_NUM_FLOAT );

	err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_FLOAT );
	assert( interval.lower.f == -CCS_INFINITY );
	assert( interval.lower_included == CCS_FALSE );
	assert( interval.upper.f == CCS_INFINITY );
	assert( interval.upper_included == CCS_FALSE );

	err = ccs_normal_distribution_get_parameters(distrib, &mu, &sigma, &stype, &quantization);
	assert( err == CCS_SUCCESS );
	assert( mu == 1.0 );
	assert( sigma == 2.0 );
	assert( stype == CCS_LINEAR );
	assert( quantization.f == 0.0 );

	err = ccs_object_get_refcount(distrib, &refcount);
	assert( err == CCS_SUCCESS );
	assert( refcount == 1 );

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
}

static void test_create_normal_distribution_errors() {
	ccs_distribution_t      distrib = NULL;
	ccs_result_t            err = CCS_SUCCESS;

	// check wrong data_type
	err = ccs_create_normal_distribution(
		(ccs_numeric_type_t)CCS_STRING,
		1.0,
		2.0,
		CCS_LINEAR,
		CCSF(0.0),
		&distrib);
	assert( err == -CCS_INVALID_TYPE );

	// check wrong data_type
	err = ccs_create_normal_distribution(
		CCS_NUM_FLOAT,
		1.0,
		2.0,
		(ccs_scale_type_t)0xdeadbeef,
		CCSF(0.0),
		&distrib);
	assert( err == -CCS_INVALID_SCALE );

	// check wrong quantization
	err = ccs_create_normal_distribution(
		CCS_NUM_FLOAT,
		1.0,
		2.0,
		CCS_LINEAR,
		CCSF(-1.0),
		&distrib);
	assert( err == -CCS_INVALID_VALUE );

	// check wrong pointer
	err = ccs_create_normal_distribution(
		CCS_NUM_FLOAT,
		1.0,
		2.0,
		CCS_LINEAR,
		CCSF(0.0),
		NULL);
	assert( err == -CCS_INVALID_VALUE );

}

static void to_float(ccs_numeric_t * values, size_t length) {
	for (size_t i = 0; i < length; i++) {
		values[i].f = values[i].i;
	}
}

static void to_log(ccs_numeric_t * values, size_t length) {
	for (size_t i = 0; i < length; i++) {
		values[i].f = log(values[i].f);
	}
}

static void test_normal_distribution_int() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	const ccs_float_t  mu = 1;
	const ccs_float_t  sigma = 2;
	ccs_numeric_t      samples[NUM_SAMPLES];
	double             mean, sig;
	ccs_interval_t     interval;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_normal_distribution(
		CCS_NUM_INTEGER,
		mu,
		sigma,
		CCS_LINEAR,
		CCSI(0),
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_INTEGER );
	assert( interval.lower.i == CCS_INT_MIN );
	assert( interval.lower_included == CCS_TRUE );
	assert( interval.upper.i == CCS_INT_MAX );
	assert( interval.upper_included == CCS_TRUE );

	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	to_float(samples, num_samples);
	mean = gsl_stats_mean((double*)samples, 1, num_samples);
	assert( mean < mu + 0.1 );
	assert( mean > mu - 0.1 );
	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, mu);
	assert( sig < sigma + 0.1 );
	assert( sig > sigma - 0.1 );

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_normal_distribution_float() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	const ccs_float_t  mu = 1;
	const ccs_float_t  sigma = 2;
	ccs_numeric_t      samples[NUM_SAMPLES];
	double             mean, sig;
	ccs_interval_t     interval;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_normal_distribution(
		CCS_NUM_FLOAT,
		mu,
		sigma,
		CCS_LINEAR,
		CCSF(0.0),
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_FLOAT );
	assert( interval.lower.f == -CCS_INFINITY );
	assert( interval.lower_included == CCS_FALSE );
	assert( interval.upper.f == CCS_INFINITY );
	assert( interval.upper_included == CCS_FALSE );

	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	mean = gsl_stats_mean((double*)samples, 1, num_samples);
	assert( mean < mu + 0.1 );
	assert( mean > mu - 0.1 );
	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, mu);
	assert( sig < sigma + 0.1 );
	assert( sig > sigma - 0.1 );

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_normal_distribution_int_log() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	const ccs_float_t  mu = 1;
	const ccs_float_t  sigma = 2;
	ccs_numeric_t      samples[NUM_SAMPLES];
	double             mean, sig;
	double             tmean, tsigma, alpha, zee, pdfa;
	ccs_interval_t     interval;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_normal_distribution(
		CCS_NUM_INTEGER,
		mu,
		sigma,
		CCS_LOGARITHMIC,
		CCSI(0),
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_INTEGER );
	assert( interval.lower.i == 1 );
	assert( interval.lower_included == CCS_TRUE );
	assert( interval.upper.i == CCS_INT_MAX );
	assert( interval.upper_included == CCS_TRUE );

	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	to_float(samples, num_samples);
	to_log(samples, num_samples);
	mean = gsl_stats_mean((double*)samples, 1, num_samples);
	// cutoff at 0.0 to have exp(v) >= 1
	// see https://en.wikipedia.org/wiki/Truncated_normal_distribution
	alpha = (log(0.5) - mu)/sigma;
	zee = (1.0 - gsl_cdf_ugaussian_P(alpha));
	pdfa = gsl_ran_ugaussian_pdf(alpha);
	tmean = mu + pdfa * sigma / zee;
	assert( mean < tmean + 0.1 );
	assert( mean > tmean - 0.1 );
	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, tmean);
	tsigma = sqrt( sigma * sigma * ( 1.0 + alpha * pdfa / zee - ( pdfa * pdfa )/( zee * zee ) ) );
	assert( sig < tsigma + 0.1 );
	assert( sig > tsigma - 1.1 );

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_normal_distribution_float_log() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	const ccs_float_t  mu = 1;
	const ccs_float_t  sigma = 2;
	ccs_numeric_t      samples[NUM_SAMPLES];
	double             mean, sig;
	ccs_interval_t     interval;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_normal_distribution(
		CCS_NUM_FLOAT,
		mu,
		sigma,
		CCS_LOGARITHMIC,
		CCSF(0.0),
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_FLOAT );
	assert( interval.lower.f == 0.0 );
	assert( interval.lower_included == CCS_FALSE );
	assert( interval.upper.f == CCS_INFINITY );
	assert( interval.upper_included == CCS_FALSE );

	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	to_log(samples, num_samples);
	mean = gsl_stats_mean((double*)samples, 1, num_samples);
	assert( mean < mu + 0.1 );
	assert( mean > mu - 0.1 );
	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, mu);
	assert( sig < sigma + 0.1 );
	assert( sig > sigma - 0.1 );

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_normal_distribution_int_quantize() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	const ccs_float_t  mu = 1;
	const ccs_float_t  sigma = 2;
	const ccs_int_t    q = 2L;
	ccs_numeric_t      samples[NUM_SAMPLES];
	double             mean, sig;
	ccs_interval_t     interval;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_normal_distribution(
		CCS_NUM_INTEGER,
		mu,
		sigma,
		CCS_LINEAR,
		CCSI(q),
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_INTEGER );
	assert( interval.lower.i == (CCS_INT_MIN/q)*q );
	assert( interval.lower_included == CCS_TRUE );
	assert( interval.upper.i == (CCS_INT_MAX/q)*q );
	assert( interval.upper_included == CCS_TRUE );

	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	to_float(samples, num_samples);
	mean = gsl_stats_mean((double*)samples, 1, num_samples);
	assert( mean < mu + 0.1 );
	assert( mean > mu - 0.1 );
	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, mu);
	assert( sig < sigma + 0.1 );
	assert( sig > sigma - 0.1 );

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_normal_distribution_float_quantize() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	const ccs_float_t  mu = 1;
	const ccs_float_t  sigma = 2;
	ccs_numeric_t      samples[NUM_SAMPLES];
	double             mean, sig;
	ccs_interval_t     interval;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_normal_distribution(
		CCS_NUM_FLOAT,
		mu,
		sigma,
		CCS_LINEAR,
		CCSF(0.2),
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_FLOAT );
	assert( interval.lower.f == -CCS_INFINITY );
	assert( interval.lower_included == CCS_FALSE );
	assert( interval.upper.f == CCS_INFINITY );
	assert( interval.upper_included == CCS_FALSE );

	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	mean = gsl_stats_mean((double*)samples, 1, num_samples);
	assert( mean < mu + 0.1 );
	assert( mean > mu - 0.1 );
	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, mu);
	assert( sig < sigma + 0.1 );
	assert( sig > sigma - 0.1 );

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_normal_distribution_int_log_quantize() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	const ccs_float_t  mu = 3;
	const ccs_float_t  sigma = 2;
	const ccs_int_t    quantize = 2L;
	ccs_numeric_t      samples[NUM_SAMPLES];
	double             mean, sig;
	double             tmean, tsigma, alpha, zee, pdfa;
	ccs_interval_t     interval;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_normal_distribution(
		CCS_NUM_INTEGER,
		mu,
		sigma,
		CCS_LOGARITHMIC,
		CCSI(quantize),
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_INTEGER );
	assert( interval.lower.i == quantize );
	assert( interval.lower_included == CCS_TRUE );
	assert( interval.upper.i == (CCS_INT_MAX/quantize)*quantize );
	assert( interval.upper_included == CCS_TRUE );

	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	to_float(samples, num_samples);
	to_log(samples, num_samples);
	mean = gsl_stats_mean((double*)samples, 1, num_samples);
	// cutoff at 0.0 to have exp(v) >= 1
	// see https://en.wikipedia.org/wiki/Truncated_normal_distribution
	alpha = (log(0.5*quantize) - mu)/sigma;
	zee = (1.0 - gsl_cdf_ugaussian_P(alpha));
	pdfa = gsl_ran_ugaussian_pdf(alpha);
	tmean = mu + pdfa * sigma / zee;
	assert( mean < tmean + 0.1 );
	assert( mean > tmean - 0.1 );
	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, tmean);
	tsigma = sqrt( sigma * sigma * ( 1.0 + alpha * pdfa / zee - ( pdfa * pdfa )/( zee * zee ) ) );
	assert( sig < tsigma + 0.1 );
	assert( sig > tsigma - 1.1 );

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_normal_distribution_float_log_quantize() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	const ccs_float_t  mu = 3;
	const ccs_float_t  sigma = 2;
	const ccs_float_t  quantization = 2.0;
	ccs_numeric_t      samples[NUM_SAMPLES];
	double             mean, sig;
	double             tmean, tsigma, alpha, zee, pdfa;
	ccs_interval_t     interval;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_normal_distribution(
		CCS_NUM_FLOAT,
		mu,
		sigma,
		CCS_LOGARITHMIC,
		CCSF(quantization),
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_FLOAT );
	assert( interval.lower.f == quantization );
	assert( interval.lower_included == CCS_TRUE );
	assert( interval.upper.f == CCS_INFINITY );
	assert( interval.upper_included == CCS_FALSE );

	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	to_log(samples, num_samples);
	mean = gsl_stats_mean((double*)samples, 1, num_samples);
	// cutoff at log(quantization/2.0) to have quantized(exp(v)) >= quantization
	// see https://en.wikipedia.org/wiki/Truncated_normal_distribution
	alpha = (log(quantization*0.5) - mu)/sigma;
	zee = (1.0 - gsl_cdf_ugaussian_P(alpha));
	pdfa = gsl_ran_ugaussian_pdf(alpha);
	tmean = mu + pdfa * sigma / zee;
	assert( mean < tmean + 0.1 );
	assert( mean > tmean - 0.1 );
	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, tmean);
	tsigma = sqrt( sigma * sigma * ( 1.0 + alpha * pdfa / zee - ( pdfa * pdfa )/( zee * zee ) ) );
	assert( sig < tsigma + 0.1 );
	assert( sig > tsigma - 0.1 );

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_normal_distribution_strided_samples() {
	ccs_distribution_t distrib1 = NULL;
	ccs_distribution_t distrib2 = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	const ccs_float_t  mu1 = 1;
	const ccs_float_t  sigma1 = 2;
	const ccs_float_t  mu2 = 0;
	const ccs_float_t  sigma2 = 2;
	ccs_numeric_t      samples[NUM_SAMPLES*2];
	double             mean, sig;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_normal_distribution(
		CCS_NUM_FLOAT,
		mu1,
		sigma1,
		CCS_LINEAR,
		CCSF(0.0),
		&distrib1);
	assert( err == CCS_SUCCESS );

	err = ccs_create_normal_distribution(
		CCS_NUM_FLOAT,
		mu2,
		sigma2,
		CCS_LINEAR,
		CCSF(0.0),
		&distrib2);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_strided_samples(distrib1, rng, num_samples, 2, samples);
	assert( err == CCS_SUCCESS );
	err = ccs_distribution_strided_samples(distrib2, rng, num_samples, 2, &(samples[0]) + 1);
	assert( err == CCS_SUCCESS );

	mean = gsl_stats_mean((double*)samples, 2, num_samples);
	assert( mean < mu1 + 0.1 );
	assert( mean > mu1 - 0.1 );
	sig    = gsl_stats_sd_m((double*)samples, 2, num_samples, mu1);
	assert( sig < sigma1 + 0.1 );
	assert( sig > sigma1 - 0.1 );

	mean = gsl_stats_mean((double*)samples + 1, 2, num_samples);
	assert( mean < mu2 + 0.1 );
	assert( mean > mu2 - 0.1 );
	sig    = gsl_stats_sd_m((double*)samples + 1, 2, num_samples, mu2);
	assert( sig < sigma2 + 0.1 );
	assert( sig > sigma2 - 0.1 );

	err = ccs_release_object(distrib1);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(distrib2);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

static void test_normal_distribution_soa_samples() {
	ccs_distribution_t distrib = NULL;
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_samples = NUM_SAMPLES;
	const ccs_float_t  mu = 1;
	const ccs_float_t  sigma = 2;
	ccs_numeric_t      samples[NUM_SAMPLES];
	double             mean, sig;
	ccs_interval_t     interval;
	ccs_numeric_t     *p_samples;

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );
	err = ccs_create_normal_distribution(
		CCS_NUM_FLOAT,
		mu,
		sigma,
		CCS_LINEAR,
		CCSF(0.0),
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_FLOAT );
	assert( interval.lower.f == -CCS_INFINITY );
	assert( interval.lower_included == CCS_FALSE );
	assert( interval.upper.f == CCS_INFINITY );
	assert( interval.upper_included == CCS_FALSE );

	p_samples = &(samples[0]);
	err = ccs_distribution_soa_samples(distrib, rng, num_samples, &p_samples);
	assert( err == CCS_SUCCESS );

	mean = gsl_stats_mean((double*)samples, 1, num_samples);
	assert( mean < mu + 0.1 );
	assert( mean > mu - 0.1 );
	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, mu);
	assert( sig < sigma + 0.1 );
	assert( sig > sigma - 0.1 );

	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

int main() {
	ccs_init();
	test_create_normal_distribution();
	test_create_normal_distribution_errors();
	test_normal_distribution_int();
	test_normal_distribution_int_log();
	test_normal_distribution_int_quantize();
	test_normal_distribution_int_log_quantize();
	test_normal_distribution_float();
	test_normal_distribution_float_log();
	test_normal_distribution_float_quantize();
	test_normal_distribution_float_log_quantize();
	test_normal_distribution_strided_samples();
	test_normal_distribution_soa_samples();
	ccs_fini();
	return 0;
}

max_stars_repo_name: deephyper/CCS

max_issues_repo_name: deephyper/CCS

max_forks_repo_name: deephyper/CCS


text:
//#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp.h>

static double fun(double x)
{
  return x + cos(x*x);
}

int main (void)
{
  const double a = 1.0;
  const double b = 10.0;
  const int steps = 10;
  double xi, yi, x[100], y[100], dx;
  FILE *input, *output;
  int i;

  input = fopen("wartosci.txt", "w");
  output = fopen("gsl_interp_steffen.txt", "w");

  dx = (b-a) / (double) steps;

  for (i = 0; i <= steps; ++i) {
	  x[i] = a + (double)i * dx + 0.5 * sin((double)i * dx);
      y[i] = fun(x[i]);
      fprintf (input,"%g %g\n", x[i], y[i]);
  }

  {
    gsl_interp_accel *acc
      = gsl_interp_accel_alloc ();

    gsl_spline *spline
      = gsl_spline_alloc (gsl_interp_steffen, steps + 1);

    gsl_spline_init (spline, x, y, steps + 1);

    for (xi = a; xi <= b; xi += 0.01) {
        yi = gsl_spline_eval (spline, xi, acc);
        fprintf (output,"%g %g\n", xi, yi);
	}

    gsl_spline_free (spline);
    gsl_interp_accel_free(acc);
  }
  return 0;
}

//Porownanie roznych interpolacji
//Opisanie komend, screeny

//https://stackoverflow.com/questions/10792015/gnuplot-plotting-multiple-line-graphs
//https://stackoverflow.com/questions/40296851/gnuplot-log-plot-y-axis
//http://lowrank.net/gnuplot/plot3d2-e.html
//http://gnuplot-surprising.blogspot.com/2012/05/how-to-pick-out-maximum-and-minimum.html
//http://gnuplot.sourceforge.net/docs_4.2/node237.html
//http://gnuplot.sourceforge.net/demo/stats.html

// stats "dane2.dat" u 3 nooutput
// set label 1 "Maximun" at STATS_pos_max, STATS_max offset 1,-0.5
// set label 2 "Minimun" at STATS_pos_min, STATS_min offset 1,0.5
// splot "dane2.dat" w p pt 3 lc rgb"#ff0000" notitle, \
//      STATS_min w l lc rgb"#00ffff" notitle, \
//      STATS_max w l lc rgb"#00ffff" notitle

stats "dane2.dat" u 1:3 nooutput
a = STATS_pos_max_y
stats "dane2.dat" u 2:3 nooutput
b = STATS_pos_max_y
stats "dane2.dat" u 3 nooutput
c = STATS_max
set arrow 1 from -5, -2.1, c+5 to a,b,c fill
splot "dane2.dat" w p pt 3 lc rgb"#ff0000" notitle, \
     STATS_min w l lc rgb"#00ffff" notitle, \
     STATS_max w l lc rgb"#00ffff" notitle

n=100
max=3.
min=-3.
width=(max-min)/n
hist(x,width)=width*floor(x/width)+width/2.0
set xrange [min:max]
set yrange[-5:5]
set offset graph 0.05,0.05,0.05,0.0
set xtics min,(max-min)/5,max
set boxwidth width*0.82
set tics out nomirror
set xlabel "x"
set ylabel "Frequency"
plot 2*cos(x*sin(x)) smooth freq w boxes lc rgb"blue", \
    sin (x**5) lt rgb"#00ff00", \
    3*sin(x) lt rgb"red", \
    "fun1.txt" with yerrorbars lt rgb"red"

ext: c


text:
#ifndef _PARRSB_H_
#define _PARRSB_H_

#include <gslib.h>

typedef struct {
  /* General options */
  int global_partitioner; // -1 - None, 0 - RSB, 1 - RCB, 2 - RIB (Default: 0)
  int local_partitioner;  // -1 - None, 0 - RSB, 1 - RCB, 2 - RIB (Default: -1)
  int debug_level;        // 0, 1, 2, .. etc (Default: 0)
  int print_timing_info;  // 0 or 1 (Default: 0)

  /* RSB specific */
  int rsb_algo;         // 0 - Lanczos, 1 - MG (Default: 0)
  int rsb_prepartition; // 0 - None, 1 - RCB , 2 - RIB (Default: 1)
  int rsb_grammian;     // 0 or 1 (Default: 1)
} parRSB_options;

extern parRSB_options parrsb_default_options;

#define fparRSB_partMesh FORTRAN_UNPREFIXED(fparrsb_partmesh, FPARRSB_PARTMESH)
void fparRSB_partMesh(int *part, int *seq, long long *vtx, double *coord,
                      int *nel, int *nve, int *options, int *comm, int *err);

int parRSB_partMesh(int *part, int *seq, long long *vtx, double *coord, int nel,
                    int nv, parRSB_options *options, MPI_Comm comm);

#define fparRSB_findConnectivity                                               \
  FORTRAN_UNPREFIXED(fparrsb_findconnectivity, FPARRSB_FINDCONNECTIVITY)
void fparRSB_findConnectivity(long long *vertexId, double *coord, int *nel,
                              int *nDim, long long *periodicInfo,
                              int *nPeriodicFaces, double *tol, MPI_Fint *fcomm,
                              int *verbose, int *err);

int parRSB_findConnectivity(long long *vertexid, double *coord, int nel,
                            int nDim, long long *periodicInfo,
                            int nPeriodicFaces, double tol, MPI_Comm comm,
                            int verbose);

#endif

ext: h

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


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

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_slauum - Computes the product U * U' or L' * L, where the triangular
 *  factor U or L is stored in the upper or lower triangular part of
 *  the array A.
 *
 *  If UPLO = 'U' or 'u' then the upper triangle of the result is stored,
 *  overwriting the factor U in A.
 *  If UPLO = 'L' or 'l' then the lower triangle of the result is stored,
 *  overwriting the factor L in A.
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          = PlasmaUpper: Upper triangle of A is stored;
 *          = PlasmaLower: Lower triangle of A is stored.
 *
 * @param[in] N
 *          The order of the triangular factor U or L.  N >= 0.
 *
 * @param[in,out] A
 *          On entry, the triangular factor U or L.
 *          On exit, if UPLO = 'U', the upper triangle of A is
 *          overwritten with the upper triangle of the product U * U';
 *          if UPLO = 'L', the lower triangle of A is overwritten with
 *          the lower triangle of the product L' * L.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,N).
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slauum = PCORE_slauum
#define CORE_slauum PCORE_slauum
#endif
void CORE_slauum(PLASMA_enum uplo, int N, float *A, int LDA)
{
    LAPACKE_slauum_work(LAPACK_COL_MAJOR, lapack_const(uplo), N, A, LDA );
}


text:
#include "ccv_nnc.h"
#include "ccv_nnc_easy.h"
#include "ccv_nnc_internal.h"
#include "ccv_internal.h"
#include "3rdparty/khash/khash.h"
#ifdef HAVE_GSL
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#else
#include "3rdparty/sfmt/SFMT.h"
#endif

KHASH_MAP_INIT_INT64(ctx, ccv_array_t*)

struct ccv_cnnp_dataframe_s {
	int row_count;
	int column_size;
	int* shuffled_idx;
#ifdef HAVE_GSL
	gsl_rng* rng;
#else
	sfmt_t sfmt;
#endif
	khash_t(ctx)* data_ctx; // The stream context based cache for data entity of columns. This helps us to avoid allocations when iterate through data.
	ccv_array_t* derived_column_data;
	ccv_cnnp_column_data_t column_data[1];
};

typedef struct {
	int stream_type;
	int column_idx_size;
	int* column_idxs;
	ccv_cnnp_column_data_enum_f data_enum;
	ccv_cnnp_column_data_deinit_f data_deinit;
	void* context;
	ccv_cnnp_column_data_context_deinit_f context_deinit;
	ccv_cnnp_column_data_map_f map;
} ccv_cnnp_derived_column_data_t;

ccv_cnnp_dataframe_t* ccv_cnnp_dataframe_new(const ccv_cnnp_column_data_t* const column_data, const int column_size, const int row_count)
{
	assert(column_size >= 0);
	ccv_cnnp_dataframe_t* const dataframe = (ccv_cnnp_dataframe_t*)cccalloc(1, sizeof(ccv_cnnp_dataframe_t) + sizeof(ccv_cnnp_column_data_t) * (column_size - 1));
	dataframe->row_count = row_count;
	dataframe->column_size = column_size;
	dataframe->data_ctx = kh_init(ctx);
	if (column_size > 0)
		memcpy(dataframe->column_data, column_data, sizeof(ccv_cnnp_column_data_t) * column_size);
	return dataframe;
}

void ccv_cnnp_dataframe_shuffle(ccv_cnnp_dataframe_t* const dataframe)
{
	assert(dataframe->row_count);
	int i;
	if (!dataframe->shuffled_idx)
	{
		dataframe->shuffled_idx = (int*)ccmalloc(sizeof(int) * dataframe->row_count);
		for (i = 0; i < dataframe->row_count; i++)
			dataframe->shuffled_idx[i] = i;
#ifdef HAVE_GSL
		assert(!dataframe->rng);
		gsl_rng_env_setup();
		dataframe->rng = gsl_rng_alloc(gsl_rng_default);
		gsl_rng_set(dataframe->rng, (unsigned long int)dataframe);
#else
		sfmt_init_gen_rand(&dataframe->sfmt, (uint32_t)dataframe);
#endif
	}
#ifdef HAVE_GSL
	gsl_ran_shuffle(dataframe->rng, dataframe->shuffled_idx, dataframe->row_count, sizeof(int));
#else
	sfmt_genrand_shuffle(&dataframe->sfmt, dataframe->shuffled_idx, dataframe->row_count, sizeof(int));
#endif
}

int ccv_cnnp_dataframe_row_count(ccv_cnnp_dataframe_t* const dataframe)
{
	return dataframe->row_count;
}

int ccv_cnnp_dataframe_add(ccv_cnnp_dataframe_t* const dataframe, ccv_cnnp_column_data_enum_f data_enum, const int stream_type, ccv_cnnp_column_data_deinit_f data_deinit, void* const context, ccv_cnnp_column_data_context_deinit_f context_deinit)
{
	if (!dataframe->derived_column_data)
		dataframe->derived_column_data = ccv_array_new(sizeof(ccv_cnnp_derived_column_data_t), 1, 0);
	ccv_cnnp_derived_column_data_t column_data = {
		.stream_type = stream_type,
		.data_enum = data_enum,
		.data_deinit = data_deinit,
		.context = context,
		.context_deinit = context_deinit,
	};
	ccv_array_push(dataframe->derived_column_data, &column_data);
	return dataframe->column_size + dataframe->derived_column_data->rnum - 1;
}

int ccv_cnnp_dataframe_map(ccv_cnnp_dataframe_t* const dataframe, ccv_cnnp_column_data_map_f map, const int stream_type, ccv_cnnp_column_data_deinit_f data_deinit, const int* const column_idxs, const int column_idx_size, void* const context, ccv_cnnp_column_data_context_deinit_f context_deinit)
{
	assert(column_idx_size > 0);
	if (!dataframe->derived_column_data)
		dataframe->derived_column_data = ccv_array_new(sizeof(ccv_cnnp_derived_column_data_t), 1, 0);
	const int column_size = dataframe->column_size + dataframe->derived_column_data->rnum;
	int i;
	for (i = 0; i < column_idx_size; i++)
		{ assert(column_idxs[i] < column_size); }
	ccv_cnnp_derived_column_data_t column_data = {
		.stream_type = stream_type,
		.column_idx_size = column_idx_size,
		.column_idxs = (int*)ccmalloc(sizeof(int) * column_idx_size),
		.map = map,
		.data_deinit = data_deinit,
		.context = context,
		.context_deinit = context_deinit,
	};
	memcpy(column_data.column_idxs, column_idxs, sizeof(int) * column_idx_size);
	ccv_array_push(dataframe->derived_column_data, &column_data);
	return dataframe->column_size + dataframe->derived_column_data->rnum - 1;
}

void* ccv_cnnp_dataframe_column_context(const ccv_cnnp_dataframe_t* const dataframe, const int column_idx)
{
	assert(column_idx >= 0);
	const int column_size = dataframe->column_size + (dataframe->derived_column_data ? dataframe->derived_column_data->rnum : 0);
	assert(column_idx < column_size);
	if (column_idx < dataframe->column_size)
		return dataframe->column_data[column_idx].context;
	assert(dataframe->derived_column_data);
	ccv_cnnp_derived_column_data_t* const derived_column_data = (ccv_cnnp_derived_column_data_t*)ccv_array_get(dataframe->derived_column_data, column_idx - dataframe->column_size);
	return derived_column_data->context;
}

typedef struct {
	int flag; // Mark this as cached or not.
	uint64_t ctx; // The stream context.
	void* data;
} ccv_cnnp_dataframe_data_item_t;

typedef struct {
	ccv_nnc_stream_context_t* stream_context;
	ccv_nnc_stream_signal_t* signal;
} ccv_cnnp_dataframe_column_ctx_t;

KHASH_MAP_INIT_INT64(iter_ctx, ccv_cnnp_dataframe_column_ctx_t*)

struct ccv_cnnp_dataframe_iter_s {
	int idx;
	int prefetch_head;
	int prefetch_tail;
	int column_idx_size;
	int fetched_size; // The size of fetched data.
	ccv_cnnp_dataframe_t* dataframe;
	void**** derived_data; // This is ridiculous, but it is true.
	void** fetched_data; // The cache to store fetched data.
	khash_t(iter_ctx)* column_ctx; // Column context specific to a stream context. The key will be a parent stream context and value will be child stream context + signal.
	ccv_array_t* prefetches; // The prefetch contents.
	int* column_idxs;
	ccv_cnnp_dataframe_data_item_t cached_data[1]; // The data cached when deriving data.
};

#define INDEX_DATA(iter) ((int*)((iter)->fetched_data))
#define FETCHED_DATA(iter, idx) ((iter)->fetched_data + ((idx) + 1) * (iter)->fetched_size)

ccv_cnnp_dataframe_iter_t* ccv_cnnp_dataframe_iter_new(ccv_cnnp_dataframe_t* const dataframe, const int* const column_idxs, const int column_idx_size)
{
	assert(column_idx_size > 0);
	const int column_size = dataframe->column_size + (dataframe->derived_column_data ? dataframe->derived_column_data->rnum : 0);
	int i;
	for (i = 0; i < column_idx_size; i++)
		{ assert(column_idxs[i] < column_size); }
	ccv_cnnp_dataframe_iter_t* const iter = (ccv_cnnp_dataframe_iter_t*)cccalloc(1, sizeof(ccv_cnnp_dataframe_iter_t) + sizeof(ccv_cnnp_dataframe_data_item_t) * column_size + sizeof(void*) * (column_idx_size - 1) + sizeof(int) * column_idx_size);
	iter->dataframe = dataframe;
	iter->prefetch_tail = -1;
	iter->column_idx_size = column_idx_size;
	iter->column_idxs = (int*)(iter->cached_data + column_size);
	memcpy(iter->column_idxs, column_idxs, sizeof(int) * column_idx_size);
	// Preallocate fetched data.
	iter->fetched_size = 1;
	iter->fetched_data = (void**)ccmalloc(sizeof(void*) * (column_size + 1));
	return iter;
}

static void _ccv_cnnp_dataframe_enqueue_data(ccv_cnnp_dataframe_t* const dataframe, void* const data, const int column_idx, const uint64_t ctx)
{
	if (!data)
		return;
	khash_t(ctx)* const data_ctx = dataframe->data_ctx;
	int ret = 0;
	khiter_t k = kh_put(ctx, data_ctx, ctx, &ret);
	assert(ret >= 0);
	const int column_size = dataframe->column_size + (dataframe->derived_column_data ? dataframe->derived_column_data->rnum : 0);
	assert(column_idx < column_size);
	// If ret == 0, the key already exist, we can get the columns directly, otherwise, create and assign back.
	ccv_array_t* const columns = (ret == 0) ? kh_val(data_ctx, k) : ccv_array_new(sizeof(ccv_array_t*), column_size, 0);
	if (ret != 0)
		kh_val(data_ctx, k) = columns;
	if (columns->rnum < column_size)
		ccv_array_resize(columns, column_size);
	ccv_array_t* column = *(ccv_array_t**)ccv_array_get(columns, column_idx);
	if (!column)
	{
		column = ccv_array_new(sizeof(void*), 1, 0);
		*(ccv_array_t**)ccv_array_get(columns, column_idx) = column;
	}
	ccv_array_push(column, &data);
}

static void* _ccv_cnnp_dataframe_dequeue_data(ccv_cnnp_dataframe_t* const dataframe, const int column_idx, ccv_nnc_stream_context_t* const stream_context)
{
	const uint64_t ctx = (uint64_t)(uintptr_t)stream_context;
	khash_t(ctx)* const data_ctx = dataframe->data_ctx;
	khiter_t k = kh_get(ctx, data_ctx, ctx);
	if (k == kh_end(data_ctx))
		return 0;
	ccv_array_t* const columns = kh_val(data_ctx, k);
	if (column_idx >= columns->rnum)
		return 0;
	ccv_array_t* const column = *(ccv_array_t**)ccv_array_get(columns, column_idx);
	if (!column || column->rnum == 0)
		return 0;
	void* const data = *(void**)ccv_array_get(column, column->rnum - 1);
	--column->rnum;
	return data;
}

static ccv_cnnp_dataframe_column_ctx_t _ccv_cnnp_child_column_ctx_for_stream_type(ccv_cnnp_dataframe_t* const dataframe, ccv_cnnp_dataframe_iter_t* const iter, const int column_idx, ccv_nnc_stream_context_t* const stream_context, const int stream_type)
{
	ccv_cnnp_dataframe_column_ctx_t child_ctx = {
		.stream_context = stream_context,
	};
	if (stream_context && ccv_nnc_stream_context_type(stream_context) != stream_type && stream_type != 0)
	{
		if (!iter->column_ctx)
			iter->column_ctx = kh_init(iter_ctx);
		khash_t(iter_ctx)* const column_ctx = iter->column_ctx;
		int ret = 0;
		khiter_t k = kh_put(iter_ctx, column_ctx, (uint64_t)(uintptr_t)stream_context, &ret);
		assert(ret >= 0);
		const int column_size = dataframe->column_size + (dataframe->derived_column_data ? dataframe->derived_column_data->rnum : 0);
		ccv_cnnp_dataframe_column_ctx_t* const ctx = (ret == 0) ? kh_val(column_ctx, k) : cccalloc(column_size, sizeof(ccv_cnnp_dataframe_column_ctx_t));
		if (ret != 0)
			kh_val(column_ctx, k) = ctx;
		if (!ctx[column_idx].stream_context)
			ctx[column_idx].stream_context = ccv_nnc_stream_context_new(stream_type);
		if (!ctx[column_idx].signal)
			ctx[column_idx].signal = ccv_nnc_stream_signal_new(stream_type);
		child_ctx = ctx[column_idx];
	}
	return child_ctx;
}

static void _ccv_cnnp_dataframe_column_data(ccv_cnnp_dataframe_t* const dataframe, ccv_cnnp_dataframe_iter_t* const iter, ccv_cnnp_dataframe_data_item_t* const cached_data, void** const fetched_data, const int* const row_idxs, const int row_size, const int column_idx, const int cached_step, ccv_nnc_stream_context_t* const stream_context)
{
	int i;
	if (cached_data[column_idx * cached_step].flag)
	{
		for (i = 1; i < row_size; i++)
			{ assert(cached_data[i + column_idx * cached_step].flag); }
		for (i = 0; i < row_size; i++)
			fetched_data[i] = cached_data[i + column_idx * cached_step].data;
		return;
	} else {
		for (i = 1; i < row_size; i++)
			{ assert(!cached_data[i + column_idx * cached_step].flag); }
		for (i = 0; i < row_size; i++)
			fetched_data[i] = _ccv_cnnp_dataframe_dequeue_data(dataframe, column_idx, stream_context);
	}
	if (column_idx >= dataframe->column_size)
	{
		assert(dataframe->derived_column_data);
		const int derived_column_idx = column_idx - dataframe->column_size;
		const ccv_cnnp_derived_column_data_t* const derived_column_data = (ccv_cnnp_derived_column_data_t*)ccv_array_get(dataframe->derived_column_data, derived_column_idx);
		ccv_cnnp_dataframe_column_ctx_t child_ctx = _ccv_cnnp_child_column_ctx_for_stream_type(dataframe, iter, column_idx, stream_context, derived_column_data->stream_type);
		const int column_idx_size = derived_column_data->column_idx_size;
		if (derived_column_data->map)
		{
			int i;
			if (!iter->derived_data)
				iter->derived_data = (void****)cccalloc(dataframe->derived_column_data->rnum, sizeof(void***));
			if (!iter->derived_data[derived_column_idx])
				iter->derived_data[derived_column_idx] = (void***)cccalloc(derived_column_data->column_idx_size, sizeof(void**));
			void*** const derived_data = iter->derived_data[derived_column_idx];
			for (i = 0; i < column_idx_size; i++)
			{
				derived_data[i] = FETCHED_DATA(iter, derived_column_data->column_idxs[i]);
				_ccv_cnnp_dataframe_column_data(dataframe, iter, cached_data, derived_data[i], row_idxs, row_size, derived_column_data->column_idxs[i], cached_step, stream_context);
			}
			derived_column_data->map(derived_data, derived_column_data->column_idx_size, row_size, fetched_data, derived_column_data->context, child_ctx.stream_context);
		} else
			derived_column_data->data_enum(column_idx, row_idxs, row_size, fetched_data, derived_column_data->context, child_ctx.stream_context);
		if (child_ctx.stream_context != stream_context)
		{
			ccv_nnc_stream_context_emit_signal(child_ctx.stream_context, child_ctx.signal);
			ccv_nnc_stream_context_wait_signal(stream_context, child_ctx.signal);
		}
	} else {
		const ccv_cnnp_column_data_t* const column_data = dataframe->column_data + column_idx;
		ccv_cnnp_dataframe_column_ctx_t child_ctx = _ccv_cnnp_child_column_ctx_for_stream_type(dataframe, iter, column_idx, stream_context, column_data->stream_type);
		column_data->data_enum(column_idx, row_idxs, row_size, fetched_data, column_data->context, child_ctx.stream_context);
		if (child_ctx.stream_context != stream_context)
		{
			ccv_nnc_stream_context_emit_signal(child_ctx.stream_context, child_ctx.signal);
			ccv_nnc_stream_context_wait_signal(stream_context, child_ctx.signal);
		}
	}
	for (i = 0; i < row_size; i++)
	{
		cached_data[i + column_idx * cached_step].flag = 1;
		cached_data[i + column_idx * cached_step].ctx = (uint64_t)(uintptr_t)stream_context;
		cached_data[i + column_idx * cached_step].data = fetched_data[i];
	}
}

int ccv_cnnp_dataframe_iter_next(ccv_cnnp_dataframe_iter_t* const iter, void** const data_ref, const int column_idx_size, ccv_nnc_stream_context_t* const stream_context)
{
	ccv_cnnp_dataframe_t* const dataframe = iter->dataframe;
	assert(column_idx_size <= iter->column_idx_size);
	const int column_size = dataframe->column_size + (dataframe->derived_column_data ? dataframe->derived_column_data->rnum : 0);
	int i;
	// Push existing data back to reusable state (note, these may not be reused immediately because they may be on a different stream context).
	for (i = 0; i < column_size; i++)
		if (iter->cached_data[i].flag)
		{
			_ccv_cnnp_dataframe_enqueue_data(dataframe, iter->cached_data[i].data, i, iter->cached_data[i].ctx);
			iter->cached_data[i].flag = 0;
			iter->cached_data[i].data = 0;
			iter->cached_data[i].ctx = 0;
		}
	const int idx = iter->idx;
	if (idx == dataframe->row_count)
		return -1;
	if (iter->prefetch_tail != -1) // If there is something in prefetch log.
	{
		ccv_array_t* const prefetches = iter->prefetches;
		assert(prefetches);
		const int lines = prefetches->rnum / column_size;
		if (iter->prefetch_head == iter->prefetch_tail) // Only one item.
			iter->prefetch_tail = -1;
		ccv_cnnp_dataframe_data_item_t* const cached_data = (ccv_cnnp_dataframe_data_item_t*)ccv_array_get(iter->prefetches, iter->prefetch_head);
		for (i = 0; i < column_size; i++)
		{
			if (!cached_data[i * lines].flag)
				continue;
			if (cached_data[i * lines].ctx == (uint64_t)(uintptr_t)stream_context) // If match existing stream context.
				iter->cached_data[i] = cached_data[i * lines];
			else // Recycle
				_ccv_cnnp_dataframe_enqueue_data(dataframe, cached_data[i * lines].data, i, cached_data[i * lines].ctx);
		}
		++iter->prefetch_head;
		assert(prefetches->rnum % column_size == 0);
		if (iter->prefetch_head >= lines)
			iter->prefetch_head = 0;
	}
	for (i = 0; i < column_idx_size; i++)
		_ccv_cnnp_dataframe_column_data(dataframe, iter, iter->cached_data, data_ref + i, dataframe->shuffled_idx ? dataframe->shuffled_idx + idx : &idx, 1, iter->column_idxs[i], 1, stream_context);
	++iter->idx;
	return 0;
}

static void _ccv_cnnp_null_prefetches(ccv_cnnp_dataframe_iter_t* const iter)
{
	ccv_cnnp_dataframe_t* const dataframe = iter->dataframe;
	assert(dataframe);
	int i, j;
	const int column_size = dataframe->column_size + (dataframe->derived_column_data ? dataframe->derived_column_data->rnum : 0);
	if (iter->prefetch_head <= iter->prefetch_tail)
	{
		assert(iter->prefetches);
		const int lines = iter->prefetches->rnum / column_size;
		for (i = iter->prefetch_head; i <= iter->prefetch_tail; i++)
		{
			ccv_cnnp_dataframe_data_item_t* const cached_data = ccv_array_get(iter->prefetches, i);
			for (j = 0; j < column_size; j++)
				if (cached_data[j * lines].flag)
					_ccv_cnnp_dataframe_enqueue_data(dataframe, cached_data[j * lines].data, j, cached_data[j * lines].ctx);
		}
	} else if (iter->prefetch_tail >= 0) { // -1 means no item.
		assert(iter->prefetches);
		const int lines = iter->prefetches->rnum / column_size;
		for (i = iter->prefetch_head; i < lines; i++)
		{
			ccv_cnnp_dataframe_data_item_t* const cached_data = ccv_array_get(iter->prefetches, i);
			for (j = 0; j < column_size; j++)
				if (cached_data[j * lines].flag)
					_ccv_cnnp_dataframe_enqueue_data(dataframe, cached_data[j * lines].data, j, cached_data[j * lines].ctx);
		}
		for (i = 0; i <= iter->prefetch_tail; i++)
		{
			ccv_cnnp_dataframe_data_item_t* const cached_data = ccv_array_get(iter->prefetches, i);
			for (j = 0; j < column_size; j++)
				if (cached_data[j * lines].flag)
					_ccv_cnnp_dataframe_enqueue_data(dataframe, cached_data[j * lines].data, j, cached_data[j * lines].ctx);
		}
	}
	iter->prefetch_head = 0;
	iter->prefetch_tail = -1;
}

static void _ccv_cnnp_prefetch_cached_data(ccv_cnnp_dataframe_iter_t* const iter, ccv_cnnp_dataframe_data_item_t* const cached_data, const int idx, const int max_to_prefetch, ccv_nnc_stream_context_t* const stream_context)
{
	ccv_cnnp_dataframe_t* const dataframe = iter->dataframe;
	assert(dataframe);
	const int column_size = dataframe->column_size + (dataframe->derived_column_data ? dataframe->derived_column_data->rnum : 0);
	assert(iter->prefetches);
	const int lines = iter->prefetches->rnum / column_size;
	int i, j;
	// Reset
	for (i = 0; i < column_size; i++)
		for (j = 0; j < max_to_prefetch; j++)
		{
			cached_data[j + i * lines].flag = 0;
			cached_data[j + i * lines].data = 0;
			cached_data[j + i * lines].ctx = 0;
		}
	if (iter->fetched_size < max_to_prefetch)
	{
		iter->fetched_data = ccrealloc(iter->fetched_data, sizeof(void*) * max_to_prefetch * (column_size + 1));
		iter->fetched_size = max_to_prefetch;
	}
	if (dataframe->shuffled_idx)
		for (i = 0; i < iter->column_idx_size; i++)
			_ccv_cnnp_dataframe_column_data(dataframe, iter, cached_data, FETCHED_DATA(iter, iter->column_idxs[i]), dataframe->shuffled_idx + idx, max_to_prefetch, iter->column_idxs[i], lines, stream_context);
	else {
		for (i = 0; i < max_to_prefetch; i++)
			INDEX_DATA(iter)[i] = idx + i;
		for (i = 0; i < iter->column_idx_size; i++)
			_ccv_cnnp_dataframe_column_data(dataframe, iter, cached_data, FETCHED_DATA(iter, iter->column_idxs[i]), INDEX_DATA(iter), max_to_prefetch, iter->column_idxs[i], lines, stream_context);
	}
}

int ccv_cnnp_dataframe_iter_prefetch(ccv_cnnp_dataframe_iter_t* const iter, const int prefetch_count, ccv_nnc_stream_context_t* const stream_context)
{
	ccv_cnnp_dataframe_t* const dataframe = iter->dataframe;
	assert(dataframe);
	const int column_size = dataframe->column_size + (dataframe->derived_column_data ? dataframe->derived_column_data->rnum : 0);
	int i, j;
	assert(iter->idx <= dataframe->row_count);
	int lines, next, max_to_prefetch;
	if (iter->prefetch_tail == -1)
	{
		if (iter->idx == dataframe->row_count)
			return -1; // Cannot be done.
		max_to_prefetch = ccv_min(dataframe->row_count - iter->idx, prefetch_count);
		if (!iter->prefetches)
		{
			iter->prefetches = ccv_array_new(sizeof(ccv_cnnp_dataframe_data_item_t), max_to_prefetch * column_size, 0);
			ccv_array_resize(iter->prefetches, max_to_prefetch * column_size);
		}
		iter->prefetch_tail = iter->prefetch_head = 0; // Advance!
		next = iter->idx;
		lines = iter->prefetches->rnum / column_size;
		// Reset to enough space.
		if (lines < max_to_prefetch)
		{
			ccv_array_resize(iter->prefetches, max_to_prefetch * column_size);
			lines = max_to_prefetch;
		}
	} else {
		assert(iter->prefetches);
		ccv_array_t* const prefetches = iter->prefetches;
		assert(prefetches->rnum % column_size == 0);
		lines = prefetches->rnum / column_size;
		const int prefetched = iter->prefetch_tail >= iter->prefetch_head ? iter->prefetch_tail - iter->prefetch_head + 1: lines - iter->prefetch_head + iter->prefetch_tail + 1;
		if (iter->idx + prefetched == dataframe->row_count) // Nothing to prefetch.
			return -1;
		max_to_prefetch = ccv_min(dataframe->row_count - (iter->idx + prefetched), prefetch_count);
		// Not enough space, need to resize.
		if (prefetched + max_to_prefetch > lines)
		{
			const int new_lines = prefetched + max_to_prefetch;
			ccv_array_resize(prefetches, new_lines * column_size);
			// These are overlap moves, have to make sure start from the end and move it up to the beginning.
			if (iter->prefetch_head > iter->prefetch_tail)
			{
				const int offset = new_lines - lines;
				for (i = column_size - 1; i >= 0; i--)
				{
					for (j = lines - 1; j >= iter->prefetch_head; j--)
						*(ccv_cnnp_dataframe_data_item_t*)ccv_array_get(prefetches, j + offset + i * new_lines) = *(ccv_cnnp_dataframe_data_item_t*)ccv_array_get(prefetches, j + i * lines);
					for (j = iter->prefetch_tail; j >= 0; j--)
						*(ccv_cnnp_dataframe_data_item_t*)ccv_array_get(prefetches, j + i * new_lines) = *(ccv_cnnp_dataframe_data_item_t*)ccv_array_get(prefetches, j + i * lines);
				}
				iter->prefetch_head += offset;
			} else {
				for (i = column_size - 1; i >= 0; i--)
					for (j = iter->prefetch_tail; j >= iter->prefetch_head; j--)
						*(ccv_cnnp_dataframe_data_item_t*)ccv_array_get(prefetches, j + i * new_lines) = *(ccv_cnnp_dataframe_data_item_t*)ccv_array_get(prefetches, j + i * lines);
			}
			lines = new_lines;
		}
		++iter->prefetch_tail; // Move to the next ready tail.
		if (iter->prefetch_tail >= lines)
			iter->prefetch_tail = 0;
		next = iter->idx + prefetched;
	}
	ccv_array_t* const prefetches = iter->prefetches;
	ccv_cnnp_dataframe_data_item_t* const cached_data = (ccv_cnnp_dataframe_data_item_t*)ccv_array_get(prefetches, iter->prefetch_tail);
	// If the tail is before the head, we must have enough space for the max_to_prefetch
	if (iter->prefetch_tail < iter->prefetch_head)
	{
		assert(iter->prefetch_tail + max_to_prefetch - 1 < iter->prefetch_head);
		_ccv_cnnp_prefetch_cached_data(iter, cached_data, next, max_to_prefetch, stream_context);
		iter->prefetch_tail += max_to_prefetch - 1;
	} else {
		// First, fetch to the end.
		const int fetch_to_end = ccv_min(max_to_prefetch, lines - iter->prefetch_tail);
		_ccv_cnnp_prefetch_cached_data(iter, cached_data, next, fetch_to_end, stream_context);
		if (fetch_to_end == max_to_prefetch)
			iter->prefetch_tail += fetch_to_end - 1;
		else {
			// Need to fetch more.
			ccv_cnnp_dataframe_data_item_t* const more_data = (ccv_cnnp_dataframe_data_item_t*)ccv_array_get(prefetches, 0);
			assert(max_to_prefetch > fetch_to_end);
			_ccv_cnnp_prefetch_cached_data(iter, more_data, next + fetch_to_end, max_to_prefetch - fetch_to_end, stream_context);
			iter->prefetch_tail = max_to_prefetch - fetch_to_end - 1;
		}
	}
	return 0;
}

int ccv_cnnp_dataframe_iter_set_cursor(ccv_cnnp_dataframe_iter_t* const iter, const int idx)
{
	ccv_cnnp_dataframe_t* const dataframe = iter->dataframe;
	assert(dataframe);
	if (idx >= dataframe->row_count)
		return -1;
	if (idx == iter->idx)
		return 0;
	iter->idx = idx;
	_ccv_cnnp_null_prefetches(iter);
	return 0;
}

void ccv_cnnp_dataframe_iter_free(ccv_cnnp_dataframe_iter_t* const iter)
{
	ccv_cnnp_dataframe_t* const dataframe = iter->dataframe;
	const int column_size = dataframe->column_size + (dataframe->derived_column_data ? dataframe->derived_column_data->rnum : 0);
	int i;
	// Push existing data back to reusable state (note, these may not be reused immediately because they may be on a different stream context).
	for (i = 0; i < column_size; i++)
		if (iter->cached_data[i].flag)
			_ccv_cnnp_dataframe_enqueue_data(dataframe, iter->cached_data[i].data, i, iter->cached_data[i].ctx);
	// Push prefetches back to reusable state.
	_ccv_cnnp_null_prefetches(iter);
	if (iter->prefetches)
		ccv_array_free(iter->prefetches);
	if (iter->derived_data)
	{
		assert(dataframe->derived_column_data);
		for (i = 0; i < dataframe->derived_column_data->rnum; i++)
			if (iter->derived_data[i])
				ccfree(iter->derived_data[i]);
		ccfree(iter->derived_data);
	}
	ccfree(iter->fetched_data);
	if (iter->column_ctx)
	{
		khash_t(iter_ctx)* const column_ctx = iter->column_ctx;
		khiter_t k;
		for (k = kh_begin(column_ctx); k != kh_end(column_ctx); ++k)
		{
			if (!kh_exist(column_ctx, k))
				continue;
			ccv_cnnp_dataframe_column_ctx_t* const ctx = kh_val(column_ctx, k);
			for (i = 0; i < column_size; i++)
			{
				if (ctx[i].stream_context)
					ccv_nnc_stream_context_free(ctx[i].stream_context);
				if (ctx[i].signal)
					ccv_nnc_stream_signal_free(ctx[i].signal);
			}
		}
		kh_destroy(iter_ctx, column_ctx);
	}
	ccfree(iter);
}

void ccv_cnnp_dataframe_free(ccv_cnnp_dataframe_t* const dataframe)
{
	int i, j;
	khash_t(ctx)* const data_ctx = dataframe->data_ctx;
	khiter_t k;
	const int column_size = dataframe->column_size + (dataframe->derived_column_data ? dataframe->derived_column_data->rnum : 0);
	for (k = kh_begin(data_ctx); k != kh_end(data_ctx); ++k)
	{
		if (!kh_exist(data_ctx, k))
			continue;
		ccv_array_t* const columns = kh_val(data_ctx, k);
		assert(columns->rnum <= column_size);
		for (i = 0; i < columns->rnum; i++)
		{
			ccv_array_t* const column = *(ccv_array_t**)ccv_array_get(columns, i);
			void* context;
			ccv_cnnp_column_data_deinit_f data_deinit;
			if (i < dataframe->column_size)
			{
				data_deinit = dataframe->column_data[i].data_deinit;
				context = dataframe->column_data[i].context;
			} else {
				assert(dataframe->derived_column_data);
				ccv_cnnp_derived_column_data_t* const derived_column_data = (ccv_cnnp_derived_column_data_t*)ccv_array_get(dataframe->derived_column_data, i - dataframe->column_size);
				data_deinit = derived_column_data->data_deinit;
				context = derived_column_data->context;
			}
			if (data_deinit)
				for (j = 0; j < column->rnum; j++)
					data_deinit(*(void**)ccv_array_get(column, j), context);
			ccv_array_free(column);
		}
		ccv_array_free(columns);
	}
	kh_destroy(ctx, data_ctx);
	if (dataframe->derived_column_data)
	{
		for (i = 0; i < dataframe->derived_column_data->rnum; i++)
		{
			ccv_cnnp_derived_column_data_t* const derived_column_data = (ccv_cnnp_derived_column_data_t*)ccv_array_get(dataframe->derived_column_data, i);
			if (derived_column_data->context_deinit)
				derived_column_data->context_deinit(derived_column_data->context);
			ccfree(derived_column_data->column_idxs);
		}
		ccv_array_free(dataframe->derived_column_data);
	}
	for (i = 0; i < dataframe->column_size; i++)
		if (dataframe->column_data[i].context_deinit)
			dataframe->column_data[i].context_deinit(dataframe->column_data[i].context);
	if (dataframe->shuffled_idx)
		ccfree(dataframe->shuffled_idx);
#ifdef HAVE_GSL
	if (dataframe->rng)
		gsl_rng_free(dataframe->rng);
#endif
	ccfree(dataframe);
}


text:
/*! /file global.h
 *  /brief Declarations of global variables and functions. */

#ifndef GLOBAL_H
#define GLOBAL_H

#ifdef COOLING_CPU
#include <gsl/gsl_spline.h>
#include <gsl/gsl_spline2d.h>
#endif

#if PRECISION==1
#ifndef TYPEDEF_DEFINED_REAL
typedef float Real;
#endif
#endif
#if PRECISION==2
#ifndef TYPEDEF_DEFINED_REAL
typedef double Real;
#endif
#endif

#define MAXLEN 140
#define TINY_NUMBER 1.0e-20
#define PI 3.141592653589793
#define MP 1.672622e-24 // mass of proton, grams
#define KB 1.380658e-16 // boltzmann constant, cgs
//#define GN 6.67259e-8 // gravitational constant, cgs
#define GN 4.49451e-18 // gravitational constant, kpc^3 / M_sun / kyr^2

#define MYR 31.536e12 //Myears in secs
#define KPC 3.086e16 // kpc in km
#define G_COSMO 4.300927161e-06; // gravitational constant, kpc km^2 s^-2 Msun^-1
#define MSUN_CGS 1.98847e33; //Msun in gr
#define KPC_CGS 3.086e21;  //kpc in cm
#define KM_CGS 1e5; //km in cm

#define TIME_UNIT 3.15569e10 // 1 kyr in s
#define LENGTH_UNIT 3.08567758e21 // 1 kpc in cm
#define MASS_UNIT 1.98855e33 // 1 solar mass in grams
#define DENSITY_UNIT (MASS_UNIT/(LENGTH_UNIT*LENGTH_UNIT*LENGTH_UNIT))
#define VELOCITY_UNIT (LENGTH_UNIT/TIME_UNIT)
#define ENERGY_UNIT (DENSITY_UNIT*VELOCITY_UNIT*VELOCITY_UNIT)
#define PRESSURE_UNIT (DENSITY_UNIT*VELOCITY_UNIT*VELOCITY_UNIT)
#define SP_ENERGY_UNIT (VELOCITY_UNIT*VELOCITY_UNIT)

#define LOG_FILE_NAME "run_output.log"

//Conserved Floor Values
#define TEMP_FLOOR 1e-3
#define DENS_FLOOR 1e-5

//Parameter for Enzo dual Energy Condition
#define DE_ETA_1 0.001 //Ratio of U to E for wich  Inetrnal Energy is used to compute the Pressure
#define DE_ETA_2 0.035 //Ratio of U to max(E_local) used to select wich Internal Energy is used for the update.

// Maximum time step for cosmological simulations
#define MAX_DELTA_A 0.001
#define MAX_EXPANSION_RATE 0.01  // Limit delta(a)/a

#ifdef COOLING_GRACKLE
  #ifdef GRACKLE_METALS
  #define NSCALARS 7
  #else
  #define NSCALARS 6
  #endif // GRACKLE_METALS
#else
#ifdef SCALAR
// Set Number of scalar fields when not using grackle
#define NSCALARS 1
#endif//SCALAR
#endif//COOLING_GRACKLE

// Inital Chemistry fractions
#define INITIAL_FRACTION_HI        0.75984603480
#define INITIAL_FRACTION_HII       1.53965115054e-4
#define INITIAL_FRACTION_HEI       0.23999999997
#define INITIAL_FRACTION_HEII      9.59999999903e-15
#define INITIAL_FRACTION_HEIII     9.59999999903e-18
#define INITIAL_FRACTION_ELECTRON  1.53965115054e-4
#define INITIAL_FRACTION_METAL     1.00000000000e-10

//Default Gravity Compiler Flags
#define GRAVITY_LONG_INTS
#define COUPLE_GRAVITATIONAL_WORK

//Default Particles Compiler Flags
#define PARTICLES_LONG_INTS
#define PARTICLES_KDK

#ifdef GRAVITY
#ifdef GRAVITY_5_POINTS_GRADIENT
#ifdef PARTICLES
#define N_GHOST_POTENTIAL 3 // 3 ghost cells are needed for 5 point gradient, ( one is for the CIC interpolation of the potential )
#else
#define N_GHOST_POTENTIAL 2 // 2 ghost cells are needed for 5 point gradient
#endif //PARTICLES

#else
#ifdef PARTICLES
#define N_GHOST_POTENTIAL 2 // 2 ghost cells are needed for 3 point gradient, ( one is for the CIC interpolation of the potential )
#else
#define N_GHOST_POTENTIAL 1 // 1 ghost cells are needed for 3 point gradient
#endif //PARTICLES
#endif //GRAVITY_5_POINTS_GRADIENT

#ifdef GRAVITY_LONG_INTS
typedef long int grav_int_t;
#else
typedef int grav_int_t;
#endif//GRAVITY_LONG_INTS
#endif

#ifdef PARTICLES
#ifdef PARTICLES_LONG_INTS
typedef long int part_int_t;
#else
typedef int part_int_t;
#endif//PARTICLES_LONG_INTS

#include <vector>
typedef std::vector<Real> real_vector_t;
typedef std::vector<part_int_t> int_vector_t;
#ifdef MPI_CHOLLA
// Constants for the inital size of the buffers for particles transfer
// and the number of data transfered for each particle
extern int N_PARTICLES_TRANSFER;
extern int N_DATA_PER_PARTICLE_TRANSFER;
#endif//MPI_CHOLLA

#ifdef AVERAGE_SLOW_CELLS
#define SLOW_FACTOR 10
#endif//AVERAGE_SLOW_CELLS

#endif//PARTICLES

#define SIGN(a) ( ((a) < 0.) ? -1. : 1. )

/* Global variables */
extern Real gama; // Ratio of specific heats
extern Real C_cfl; // CFL number (0 - 0.5)
extern Real t_comm;
extern Real t_other;

#ifdef COOLING_CPU
extern gsl_interp_accel *acc;
extern gsl_interp_accel *xacc;
extern gsl_interp_accel *yacc;
extern gsl_spline *highT_C_spline;
extern gsl_spline2d *lowT_C_spline;
extern gsl_spline2d *lowT_H_spline;
#endif
#ifdef COOLING_GPU
extern float *cooling_table;
extern float *heating_table;
#endif

/*! \fn void Set_Gammas(Real gamma_in)
 *  \brief Set gamma values for Riemann solver. */
extern void Set_Gammas(Real gamma_in);

/*! \fn double get_time(void)
 *  \brief Returns the current clock time. */
extern double get_time(void);

/*! \fn int sgn
 *  \brief Mathematical sign function. Returns sign of x. */
extern int sgn(Real x);

#ifndef CUDA
/*! \fn Real calc_eta(Real cW[], Real gamma)
 *  \brief Calculate the eta value for the H correction. */
extern Real calc_eta(Real cW[], Real gamma);
#endif

struct parameters
{
  int nx;
  int ny;
  int nz;
  double tout;
  double outstep;
  int n_steps_output;
  Real gamma;
  char init[MAXLEN];
  int nfile;
  int outstep_hydro;
  int outstep_particle;
  int outstep_projection;
  int outstep_rotated_projection;
  int outstep_slice;
  Real xmin;
  Real ymin;
  Real zmin;
  Real xlen;
  Real ylen;
  Real zlen;
  int xl_bcnd;
  int xu_bcnd;
  int yl_bcnd;
  int yu_bcnd;
  int zl_bcnd;
  int zu_bcnd;
#ifdef   MPI_CHOLLA
  int xlg_bcnd;
  int xug_bcnd;
  int ylg_bcnd;
  int yug_bcnd;
  int zlg_bcnd;
  int zug_bcnd;
#endif /*MPI_CHOLLA*/
  char custom_bcnd[MAXLEN];
  char outdir[MAXLEN];
  char indir[MAXLEN]; //Folder to load Initial conditions from
  Real rho;
  Real vx;
  Real vy;
  Real vz;
  Real P;
  Real A;
  Real rho_l;
  Real v_l;
  Real P_l;
  Real rho_r;
  Real v_r;
  Real P_r;
  Real diaph;
#ifdef ROTATED_PROJECTION
  int nxr;
  int nzr;
  Real delta;
  Real theta;
  Real phi;
  Real Lx;
  Real Lz;
  int n_delta;
  Real ddelta_dt;
  int flag_delta;
#endif /*ROTATED_PROJECTION*/
#ifdef COSMOLOGY
  Real H0;
  Real Omega_M;
  Real Omega_L;
  Real Omega_b;
  Real Init_redshift;
  Real End_redshift;
  char scale_outputs_file[MAXLEN]; //File for the scale_factor output values for cosmological simulations
#endif //COSMOLOGY
#ifdef TILED_INITIAL_CONDITIONS
  Real tile_length;
#endif //TILED_INITIAL_CONDITIONS

#ifdef SET_MPI_GRID
  // Set the MPI Processes grid [n_proc_x, n_proc_y, n_proc_z]
  int n_proc_x;
  int n_proc_y;
  int n_proc_z;
#endif
  int bc_potential_type;
#ifdef COOLING_GRACKLE
  char UVB_rates_file[MAXLEN];
#endif
#ifdef ANALYSIS
  char analysis_scale_outputs_file[MAXLEN]; //File for the scale_factor output values for cosmological simulations {{}}
  char analysisdir[MAXLEN];
  int lya_skewers_stride;
  Real lya_Pk_d_log_k;
  #ifdef OUTPUT_SKEWERS
  char skewersdir[MAXLEN];
  #endif
#endif
};

/*! \fn void parse_params(char *param_file, struct parameters * parms);
 *  \brief Reads the parameters in the given file into a structure. */
extern void parse_params (char *param_file, struct parameters * parms);

/*! \fn int is_param_valid(char *name);
 * \brief Verifies that a param is valid (even if not needed).  Avoids "warnings" in output. */
extern int is_param_valid(const char *name);

#endif //GLOBAL_H

max_stars_repo_name: bvillasen/cholla

max_issues_repo_name: bvillasen/cholla

max_forks_repo_name: bvillasen/cholla


text:
#pragma once
#include <gsl/gsl>
#include <map>
#include "halley/text/halleystring.h"
#include "halley/maths/vector4.h"

namespace Halley
{
	class Path;
	class Image;
	struct ImageData;

	class AsepriteReader
	{
	public:
		static std::vector<ImageData> importAseprite(String baseName, gsl::span<const gsl::byte> fileData, bool trim);
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>

#include "allvars.h"
#include "proto.h"

/*! \file conduction.c
*  \brief Computes conduction bases on an implicit diffusion solver
*
*/

#ifdef CONDUCTION

#define MAX_COND_ITER 25
#define COND_ITER_ACCURACY 1.0e-4

static struct conductiondata_in
{
  MyDouble Pos[3];
  MyFloat Hsml;
  MyFloat Density;
  MyFloat Kappa;
  int NodeList[NODELISTLENGTH];
}
 *ConductionDataIn, *ConductionDataGet;

static struct conductiondata_out
{
  MyFloat Out;
  MyFloat Sum;
}
 *ConductionDataResult, *ConductionDataOut;

static double *Energy, *EnergyOld;
static double *Residual, *DVec, *QVec;
static double *Kappa;

/* we will use the conjugate gradient method to compute a solution
   of the implicitly formulate diffusion equation */
/* Note: the conduction equation we solve is really formulated with u instead of T, i.e.
   the factor (gamma-1)*mu*mp/k_B that converts from T to u is implicitely absorbed in a
   redefinition of kappa */

void conduction(void)
{
  int i, iter;
  double delta0, delta1, alpha, beta, a3inv, dt, rel_change, loc_max_rel_change, glob_max_rel_change;
  double sumnew, sumold, sumtransfer, sumnew_tot, sumold_tot, sumtransfer_tot;

  if(ThisTask == 0)
    {
      printf("Start thermal conduction...\n");
      fflush(stdout);
    }

  Energy = (double *) mymalloc(N_gas * sizeof(double));
  EnergyOld = (double *) mymalloc(N_gas * sizeof(double));
  Residual = (double *) mymalloc(N_gas * sizeof(double));
  DVec = (double *) mymalloc(N_gas * sizeof(double));
  QVec = (double *) mymalloc(N_gas * sizeof(double));

  Kappa = (double *) mymalloc(N_gas * sizeof(double));

  if(All.ComovingIntegrationOn)
    a3inv = 1 / (All.Time * All.Time * All.Time);
  else
    a3inv = 1.0;

  dt = (All.Conduction_Ti_endstep - All.Conduction_Ti_begstep) * All.Timebase_interval;

  if(All.ComovingIntegrationOn)
    dt *= All.Time / hubble_function(All.Time);

  if(ThisTask == 0)
    {
      printf("dt=%g\n", dt);
    }

  /* First, let's compute the thermal energies per unit mass and
     conductivities for all particles */

  for(i = 0; i < N_gas; i++)
    {
      if(P[i].Type == 0)
	{
	  /* this gives the thermal energy per unit mass for particle i */
	  Energy[i] = EnergyOld[i] = SphP[i].Entropy *
	    pow(SphP[i].d.Density * a3inv, GAMMA_MINUS1) / GAMMA_MINUS1;

#ifdef CONDUCTION_CONSTANT
	  Kappa[i] = All.ConductionCoeff;
#else
	  Kappa[i] = All.ConductionCoeff * pow(EnergyOld[i], 2.5);
#ifdef CONDUCTION_SATURATION
	  electron_free_path = All.ElectronFreePathFactor * Energy[i] * Energy[i] / (SphP[i].Density * a3inv);
	  temp_scale_length =
	    atime * fabs(smoothentr) / sqrt(gradentr[0] * gradentr[0] +
					    gradentr[1] * gradentr[1] + gradentr[2] * gradentr[2]);
	  Kappa[i] /= (1 + 4.2 * electron_free_path / temp_scale_length);
#endif
#endif

#ifdef SFR
	  if(SphP[i].Density * a3inv >= All.PhysDensThresh)
	    Kappa[i] = 0;
#endif

	  /* we'll factor the timestep into the conductivities, for simplicity */
	  Kappa[i] *= dt;
	}
    }

  /* Let's start the Conjugate Gradient Algorithm */

  /* Initialization */

  conduction_matrix_multiply(EnergyOld, Residual);

  for(i = 0; i < N_gas; i++)
    {
      if(P[i].Type == 0)
	{
	  Residual[i] = EnergyOld[i] - Residual[i];
	  DVec[i] = Residual[i];
	}
    }

  delta1 = conduction_vector_multiply(Residual, Residual);
  delta0 = delta1;

  iter = 0;			/* iteration counter */
  glob_max_rel_change = 1 + COND_ITER_ACCURACY;	/* to make sure that we enter the iteration */

  while(iter < MAX_COND_ITER && glob_max_rel_change > COND_ITER_ACCURACY && delta1 > 0)
    {
      conduction_matrix_multiply(DVec, QVec);

      alpha = delta1 / conduction_vector_multiply(DVec, QVec);

      for(i = 0, loc_max_rel_change = 0; i < N_gas; i++)
	{
	  Energy[i] += alpha * DVec[i];
	  Residual[i] -= alpha * QVec[i];

	  rel_change = alpha * DVec[i] / Energy[i];
	  if(loc_max_rel_change < rel_change)
	    loc_max_rel_change = rel_change;
	}

      delta0 = delta1;
      delta1 = conduction_vector_multiply(Residual, Residual);

      beta = delta1 / delta0;

      for(i = 0; i < N_gas; i++)
	DVec[i] = Residual[i] + beta * DVec[i];

      iter++;

      MPI_Allreduce(&loc_max_rel_change, &glob_max_rel_change, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);

      if(ThisTask == 0)
	{
	  printf("conduction iter=%d  delta1=%g delta1/delta0=%g  max-rel-change=%g\n",
		 iter, delta1, delta1 / delta0, glob_max_rel_change);
	  fflush(stdout);
	}
    }

  /* Now we have the solution vector in Energy[] */
  /* assign it to the entropies, and update the pressure */

  for(i = 0, sumnew = sumold = sumtransfer = 0; i < N_gas; i++)
    {
      if(P[i].Type == 0)
	{
	  sumnew += P[i].Mass * Energy[i];
	  sumold += P[i].Mass * EnergyOld[i];
	  sumtransfer += P[i].Mass * fabs(Energy[i] - EnergyOld[i]);

	  SphP[i].Entropy = Energy[i] / pow(SphP[i].d.Density * a3inv, GAMMA_MINUS1) * GAMMA_MINUS1;
	}
    }

  MPI_Allreduce(&sumnew, &sumnew_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&sumold, &sumold_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&sumtransfer, &sumtransfer_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

  if(ThisTask == 0)
    {
      printf("\nconduction finished. energy_before=%g energy_after=%g  rel-change=%g rel-transfer=%g\n\n",
	     sumold_tot, sumnew_tot, (sumnew_tot - sumold_tot) / sumold_tot, sumtransfer_tot / sumold_tot);
      fflush(stdout);
    }

  myfree(Kappa);
  myfree(QVec);
  myfree(DVec);
  myfree(Residual);
  myfree(EnergyOld);
  myfree(Energy);
}

double conduction_vector_multiply(double *a, double *b)
{
  int i;
  double sum, sumall;

  for(i = 0, sum = 0; i < N_gas; i++)
    if(P[i].Type == 0)
      sum += a[i] * b[i];

  MPI_Allreduce(&sum, &sumall, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

  return sumall;
}

void conduction_matrix_multiply(double *in, double *out)
{
  int i, j, k, ngrp, ndone, ndone_flag, dummy;
  int sendTask, recvTask, nexport, nimport, place;
  double *sum;

  /* allocate buffers to arrange communication */

  Ngblist = (int *) mymalloc(NumPart * sizeof(int));

  All.BunchSize =
    (int) ((All.BufferSize * 1024 * 1024) / (sizeof(struct data_index) + sizeof(struct data_nodelist) +
					     sizeof(struct conductiondata_in) +
					     sizeof(struct conductiondata_out) +
					     sizemax(sizeof(struct conductiondata_in),
						     sizeof(struct conductiondata_out))));
  DataIndexTable = (struct data_index *) mymalloc(All.BunchSize * sizeof(struct data_index));
  DataNodeList = (struct data_nodelist *) mymalloc(All.BunchSize * sizeof(struct data_nodelist));

  sum = (double *) mymalloc(N_gas * sizeof(double));

  i = 0;			/* need to go over all gas particles */

  do
    {
      for(j = 0; j < NTask; j++)
	{
	  Send_count[j] = 0;
	  Exportflag[j] = -1;
	}

      /* do local particles and prepare export list */
      for(nexport = 0; i < N_gas; i++)
	if(P[i].Type == 0)
	  {
	    if(conduction_evaluate(i, 0, in, out, sum, &nexport, Send_count) < 0)
	      break;
	  }

#ifdef MYSORT
      mysort_dataindex(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);
#else
      qsort(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);
#endif

      MPI_Allgather(Send_count, NTask, MPI_INT, Sendcount_matrix, NTask, MPI_INT, MPI_COMM_WORLD);

      for(j = 0, nimport = 0, Recv_offset[0] = 0, Send_offset[0] = 0; j < NTask; j++)
	{
	  Recv_count[j] = Sendcount_matrix[j * NTask + ThisTask];
	  nimport += Recv_count[j];

	  if(j > 0)
	    {
	      Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
	      Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
	    }
	}

      ConductionDataGet = (struct conductiondata_in *) mymalloc(nimport * sizeof(struct conductiondata_in));
      ConductionDataIn = (struct conductiondata_in *) mymalloc(nexport * sizeof(struct conductiondata_in));

      /* prepare particle data for export */

      for(j = 0; j < nexport; j++)
	{
	  place = DataIndexTable[j].Index;

	  for(k = 0; k < 3; k++)
	    ConductionDataIn[j].Pos[k] = P[place].Pos[k];

	  ConductionDataIn[j].Hsml = PPP[place].Hsml;
	  ConductionDataIn[j].Density = SphP[place].d.Density;
	  ConductionDataIn[j].Kappa = Kappa[place];

	  memcpy(ConductionDataIn[j].NodeList,
		 DataNodeList[DataIndexTable[j].IndexGet].NodeList, NODELISTLENGTH * sizeof(int));
	}

      /* exchange particle data */
      for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	{
	  sendTask = ThisTask;
	  recvTask = ThisTask ^ ngrp;

	  if(recvTask < NTask)
	    {
	      if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		{
		  /* get the particles */
		  MPI_Sendrecv(&ConductionDataIn[Send_offset[recvTask]],
			       Send_count[recvTask] * sizeof(struct conductiondata_in), MPI_BYTE,
			       recvTask, TAG_HYDRO_A,
			       &ConductionDataGet[Recv_offset[recvTask]],
			       Recv_count[recvTask] * sizeof(struct conductiondata_in), MPI_BYTE,
			       recvTask, TAG_HYDRO_A, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		}
	    }
	}

      myfree(ConductionDataIn);
      ConductionDataResult =
	(struct conductiondata_out *) mymalloc(nimport * sizeof(struct conductiondata_out));
      ConductionDataOut = (struct conductiondata_out *) mymalloc(nexport * sizeof(struct conductiondata_out));

      /* now do the particles that were sent to us */
      for(j = 0; j < nimport; j++)
	conduction_evaluate(j, 1, in, out, sum, &dummy, &dummy);

      if(i >= N_gas)
	ndone_flag = 1;
      else
	ndone_flag = 0;

      MPI_Allreduce(&ndone_flag, &ndone, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

      /* get the result */
      for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	{
	  sendTask = ThisTask;
	  recvTask = ThisTask ^ ngrp;
	  if(recvTask < NTask)
	    {
	      if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		{
		  /* send the results */
		  MPI_Sendrecv(&ConductionDataResult[Recv_offset[recvTask]],
			       Recv_count[recvTask] * sizeof(struct conductiondata_out),
			       MPI_BYTE, recvTask, TAG_HYDRO_B,
			       &ConductionDataOut[Send_offset[recvTask]],
			       Send_count[recvTask] * sizeof(struct conductiondata_out),
			       MPI_BYTE, recvTask, TAG_HYDRO_B, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		}
	    }
	}

      /* add the result to the local particles */
      for(j = 0; j < nexport; j++)
	{
	  place = DataIndexTable[j].Index;

	  out[place] += ConductionDataOut[j].Out;
	  sum[place] += ConductionDataOut[j].Sum;
	}

      myfree(ConductionDataOut);
      myfree(ConductionDataResult);
      myfree(ConductionDataGet);
    }
  while(ndone < NTask);

  /* do final operations  */

  for(i = 0; i < N_gas; i++)
    if(P[i].Type == 0)
      {
	out[i] += in[i] * (1 + sum[i]);
      }

  myfree(sum);

  myfree(DataNodeList);
  myfree(DataIndexTable);
  myfree(Ngblist);
}

int conduction_evaluate(int target, int mode, double *in, double *out, double *sum,
			int *nexport, int *nsend_local)
{
  int startnode, numngb, listindex = 0;
  int j, n;
  MyDouble *pos;
  MyFloat h_i, rho;
  double dx, dy, dz;
  double h_i2, hinv_i, hinv4_i, hinv_j, hinv4_j;
  double dwk_i, h_j, dwk_j, dwk;
  double r, r2, u, Kappa_i, kappa_mean, w, out_sum, w_sum;

  if(mode == 0)
    {
      pos = P[target].Pos;
      h_i = PPP[target].Hsml;
      rho = SphP[target].d.Density;
      Kappa_i = Kappa[target];
    }
  else
    {
      pos = ConductionDataGet[target].Pos;
      h_i = ConductionDataGet[target].Hsml;
      rho = ConductionDataGet[target].Density;
      Kappa_i = ConductionDataGet[target].Kappa;
    }

  h_i2 = h_i * h_i;
  hinv_i = 1.0 / h_i;
#ifndef  TWODIMS
  hinv4_i = hinv_i * hinv_i * hinv_i * hinv_i;
#else
  hinv4_i = hinv_i * hinv_i * hinv_i / boxSize_Z;
#endif

  /* initialize variables before SPH loop is started */
  out_sum = 0;
  w_sum = 0;

  /* Now start the actual SPH computation for this particle */

  if(mode == 0)
    {
      startnode = All.MaxPart;	/* root node */
    }
  else
    {
      startnode = ConductionDataGet[target].NodeList[0];
      startnode = Nodes[startnode].u.d.nextnode;	/* open it */
    }

  while(startnode >= 0)
    {
      while(startnode >= 0)
	{
	  numngb = ngb_treefind_pairs(pos, h_i, target, &startnode, mode, nexport, nsend_local);

	  if(numngb < 0)
	    return -1;

	  for(n = 0; n < numngb; n++)
	    {
	      j = Ngblist[n];

	      dx = pos[0] - P[j].Pos[0];
	      dy = pos[1] - P[j].Pos[1];
	      dz = pos[2] - P[j].Pos[2];

#ifdef PERIODIC			/*  now find the closest image in the given box size  */
	      if(dx > boxHalf_X)
		dx -= boxSize_X;
	      if(dx < -boxHalf_X)
		dx += boxSize_X;
	      if(dy > boxHalf_Y)
		dy -= boxSize_Y;
	      if(dy < -boxHalf_Y)
		dy += boxSize_Y;
	      if(dz > boxHalf_Z)
		dz -= boxSize_Z;
	      if(dz < -boxHalf_Z)
		dz += boxSize_Z;
#endif
	      r2 = dx * dx + dy * dy + dz * dz;

	      h_j = PPP[j].Hsml;
	      if(r2 < h_i2 || r2 < h_j * h_j)
		{
		  r = sqrt(r2);
		  if(r > 0)
		    {
		      if(r2 < h_i2)
			{
			  u = r * hinv_i;
			  if(u < 0.5)
			    dwk_i = hinv4_i * u * (KERNEL_COEFF_3 * u - KERNEL_COEFF_4);
			  else
			    dwk_i = hinv4_i * KERNEL_COEFF_6 * (1.0 - u) * (1.0 - u);
			}
		      else
			dwk_i = 0;

		      if(r2 < h_j * h_j)
			{
			  hinv_j = 1.0 / h_j;
#ifndef  TWODIMS
			  hinv4_j = hinv_j * hinv_j * hinv_j * hinv_j;
#else
			  hinv4_j = hinv_j * hinv_j * hinv_j / boxSize_Z;
#endif
			  u = r * hinv_j;
			  if(u < 0.5)
			    dwk_j = hinv4_j * u * (KERNEL_COEFF_3 * u - KERNEL_COEFF_4);
			  else
			    dwk_j = hinv4_j * KERNEL_COEFF_6 * (1.0 - u) * (1.0 - u);
			}
		      else
			dwk_j = 0;

		      /* conduction equation kernel */
		      if((Kappa_i + Kappa[j]) > 0)
			{
			  kappa_mean = 2 * (Kappa_i * Kappa[j]) / (Kappa_i + Kappa[j]);
			  dwk = 0.5 * (dwk_i + dwk_j);

			  w = 2.0 * P[j].Mass / (rho * SphP[j].d.Density) * kappa_mean * (-dwk) / r;

			  out_sum += (-w * in[j]);
			  w_sum += w;
			}
		    }
		}
	    }
	}

      if(mode == 1)
	{
	  listindex++;
	  if(listindex < NODELISTLENGTH)
	    {
	      startnode = ConductionDataGet[target].NodeList[listindex];
	      if(startnode >= 0)
		startnode = Nodes[startnode].u.d.nextnode;	/* open it */
	    }
	}
    }

  /* Now collect the result at the right place */
  if(mode == 0)
    {
      out[target] = out_sum;
      sum[target] = w_sum;
    }
  else
    {
      ConductionDataResult[target].Out = out_sum;
      ConductionDataResult[target].Sum = w_sum;
    }

  return 0;
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
#pragma once

#include <petsc.h>

// TODO: include this in a different header?
#ifdef PETSC_USE_64BIT_INDICES
  #define builtin_parity __builtin_parityl
#else
  #define builtin_parity __builtin_parity
#endif

/* define a struct & enum to hold subspace information */
// TODO: more descriptive name?
typedef enum _subspace_type
{
  FULL,
  PARITY,
  AUTO
} subspace_type;

typedef struct _subspaces_t
{
  subspace_type left_type;
  subspace_type right_type;
  void *left_data;
  void *right_data;
} subspaces_t;

/***** FULL *****/

typedef struct _data_Full
{
  PetscInt L;
} data_Full;

static inline PetscErrorCode CopySubspaceData_Full(data_Full** out_p, const data_Full* in) {
  PetscErrorCode ierr;
  ierr = PetscMalloc1(1, out_p);CHKERRQ(ierr);
  ierr = PetscMemcpy(*out_p, in, sizeof(data_Full));CHKERRQ(ierr);
  return ierr;
}

static inline PetscErrorCode DestroySubspaceData_Full(data_Full* data) {
  PetscErrorCode ierr;
  ierr = PetscFree(data);CHKERRQ(ierr);
  return ierr;
}

static inline PetscInt Dim_Full(const data_Full* data) {
  return (PetscInt)1 << data->L;
}

static inline PetscInt S2I_Full(PetscInt state, const data_Full* data) {
  return state;
}

static inline PetscInt S2I_nocheck_Full(PetscInt state, const data_Full* data) {
  return state;
}

static inline PetscInt I2S_Full(PetscInt idx, const data_Full* data) {
  return idx;
}

static inline void S2I_Full_array(int n, const data_Full* data, const PetscInt* states, PetscInt* idxs) {
  PetscMemcpy(idxs, states, n*sizeof(PetscInt));
}

static inline void I2S_Full_array(int n, const data_Full* data, const PetscInt* idxs, PetscInt* states) {
  PetscMemcpy(states, idxs, n*sizeof(PetscInt));
}

/***** PARITY *****/

typedef struct _data_Parity
{
  PetscInt L;
  PetscInt space;
} data_Parity;

static inline PetscErrorCode CopySubspaceData_Parity(data_Parity** out_p, const data_Parity* in) {
  PetscErrorCode ierr;
  ierr = PetscMalloc1(1, out_p);CHKERRQ(ierr);
  ierr = PetscMemcpy(*out_p, in, sizeof(data_Parity));CHKERRQ(ierr);
  return ierr;
}

static inline PetscErrorCode DestroySubspaceData_Parity(data_Parity* data) {
  PetscErrorCode ierr;
  ierr = PetscFree(data);CHKERRQ(ierr);
  return ierr;
}

static inline PetscInt Dim_Parity(const data_Parity* data) {
  return (PetscInt)1 << (data->L-1);
}

static inline PetscInt S2I_Parity(PetscInt state, const data_Parity* data) {
  if (builtin_parity(state) == data->space) {
    return state>>1;
  }
  else {
    return (PetscInt)(-1);
  }
}

static inline PetscInt S2I_nocheck_Parity(PetscInt state, const data_Parity* data) {
  return state>>1;
}

static inline PetscInt I2S_Parity(PetscInt idx, const data_Parity* data) {
  return (idx<<1) | (builtin_parity(idx) ^ data->space);
}

static inline void S2I_Parity_array(int n, const data_Parity* data, const PetscInt* states, PetscInt* idxs) {
  PetscInt i;
  for (i = 0; i < n; ++i) {
    idxs[i] = S2I_Parity(states[i], data);
  }
}

static inline void I2S_Parity_array(int n, const data_Parity* data, const PetscInt* idxs, PetscInt* states) {
  PetscInt i;
  for (i = 0; i < n; ++i) {
    states[i] = I2S_Parity(idxs[i], data);
  }
}

/***** AUTO *****/

typedef struct _data_Auto
{
  PetscInt L;
  PetscInt dim;
  PetscInt* state_map;
  PetscInt* rmap_indices;
  PetscInt* rmap_states;
} data_Auto;

static inline PetscErrorCode CopySubspaceData_Auto(data_Auto** out_p, const data_Auto* in) {
  PetscErrorCode ierr;
  ierr = PetscMalloc1(1, out_p);CHKERRQ(ierr);
  ierr = PetscMemcpy(*out_p, in, sizeof(data_Auto));CHKERRQ(ierr);

  ierr = PetscMalloc1(in->dim, &((*out_p)->state_map));CHKERRQ(ierr);
  ierr = PetscMemcpy((*out_p)->state_map, in->state_map, in->dim*sizeof(PetscInt));CHKERRQ(ierr);

  ierr = PetscMalloc1(in->dim, &((*out_p)->rmap_indices));CHKERRQ(ierr);
  ierr = PetscMemcpy((*out_p)->rmap_indices, in->rmap_indices, in->dim*sizeof(PetscInt));CHKERRQ(ierr);

  ierr = PetscMalloc1(in->dim, &((*out_p)->rmap_states));CHKERRQ(ierr);
  ierr = PetscMemcpy((*out_p)->rmap_states, in->rmap_states, in->dim*sizeof(PetscInt));CHKERRQ(ierr);

  return ierr;
}

static inline PetscErrorCode DestroySubspaceData_Auto(data_Auto* data) {
  PetscErrorCode ierr;
  ierr = PetscFree(data->state_map);CHKERRQ(ierr);
  ierr = PetscFree(data->rmap_indices);CHKERRQ(ierr);
  ierr = PetscFree(data->rmap_states);CHKERRQ(ierr);
  ierr = PetscFree(data);CHKERRQ(ierr);
  return ierr;
}

static inline PetscInt Dim_Auto(const data_Auto* data) {
  return data->dim;
}

static inline PetscInt S2I_Auto(PetscInt state, const data_Auto* data) {
  /* do a binary search on rmap_states */
  PetscInt left, right, mid;
  left = 0;
  right = data->dim;
  while (left <= right) {
    mid = left + (right-left)/2;
    if (data->rmap_states[mid] == state) {
      return data->rmap_indices[mid];
    }

    if (data->rmap_states[mid] < state) {
      left = mid + 1;
    }
    else {
      right = mid - 1;
    }
  }
  /* element was not in the array */
  return -1;
}

static inline PetscInt I2S_Auto(PetscInt idx, const data_Auto* data) {
  return data->state_map[idx];
}

static inline void S2I_Auto_array(int n, const data_Auto* data, const PetscInt* states, PetscInt* idxs) {
  PetscInt i;
  for (i = 0; i < n; ++i) {
    idxs[i] = S2I_Auto(states[i], data);
  }
}

static inline void I2S_Auto_array(int n, const data_Auto* data, const PetscInt* idxs, PetscInt* states) {
  PetscInt i;
  for (i = 0; i < n; ++i) {
    states[i] = I2S_Auto(idxs[i], data);
  }
}

max_stars_repo_name: GregDMeyer/dynamite

max_issues_repo_name: GregDMeyer/dynamite

max_forks_repo_name: GregDMeyer/dynamite


text:
#ifndef L_1_LASSO_C_H
#define L_1_LASSO_C_H

#include "Matrix.h"
#include "PDCD.h"
#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <math.h>

//This class solves problem of the form lambda3\|Ax- b\|_1+ g(x) by SMART_CD;
// where g(x)=\frac{lambda2}{2}\|x\|_2+lambda1\|x\|_1.
// phi_i(x)= 0, h_i(x)= lambda3\|x- b_i\|_1.

template<typename L, typename D>
class L_1_Lasso_c: public PDCD<L, D>
{
private:

D lambda1;

D lambda2;

D lambda3;

Matrix<L,D> my_M;

D val_lambda_f;

protected:

public:

  L_1_Lasso_c(const char* Matrix_file,D val_lambda1, D val_lambda2, D val_lambda3)
  :PDCD<L,D>(),my_M(Matrix_file)
  {
    lambda1=val_lambda1;
    lambda2=val_lambda2;
    lambda3=val_lambda3;
    val_lambda_f=0;
    this->mu_f=0;
    this->mu_g=lambda2;
  }

  L get_n(){return my_M.nfeatures;}
  L get_m(){return my_M.nsamples;}

  inline D gradient_of_phi_j(D x1, L i){
      return 0;
  }

  inline D value_of_g_i(D x, L i){
        return lambda2*x*x/2+lambda1*fabs(x);
      }

  inline D value_of_phi_j(D x1, L i){return 0;}

  inline D value_of_h_j(D x, L i){
  	return lambda3*fabs(x- my_M.b[i]);
  }

  inline D value_of_h_star_j(D x, L j){
  	if (x<= lambda3 && x>= -lambda3){
  	return my_M.b[j]*x;
  		}
  	else{
  		cout<< "error in h*(x)"<< endl;
  		return std::numeric_limits<double>::max();
	  }
   }

  inline D prox_of_g_i(D x1,D x2,D x3, L i){
        return compute_one_step(1.0/x2, x1, x3)-x3;
      }

      D compute_one_step(D tau, D u, D x){
        D new_x;
        if(x>tau*(lambda1+u))
         new_x=(x-tau*(lambda1+u))/(1+lambda2*tau);
        else if(x<tau*(u-lambda1))
         new_x=(x-tau*(u-lambda1))/(1+lambda2*tau);
        else
         new_x=0;
        return new_x;
      }

// compute argmin_x h^*(x)+ x2/2*(x- x1)*(x- x1)
  inline D prox_of_h_star_j(D x1, D x2, L j){
  	if (x1- my_M.b[j]/x2> lambda3){
  		return lambda3;
   }
	  else if (x1- my_M.b[j]/x2< -lambda3){
	  	return -lambda3;
	  }
	  else{
	  	return x1- my_M.b[j]/x2;
	  }
	  }

  inline D value_of_phistar_i(D x,L i) {return 0;}

  inline void set_matrix_M(){
    this->data_M=my_M;
  }

  inline void set_matrix_A(){
    this->data_A.nsamples=0;
    this->data_A.nfeatures=this->data_M.nfeatures;
    this->data_A.nnz=0;
    this->data_A.ptr.resize(1,0);
    this->data_A.ptr_t.resize(this->data_M.nfeatures+1,0);
  }

  void SMART_CD_solver(D beta_0, vector<D> & x0,vector<D> & y0,L val_tau, L max_nb_outer, L p_N, string filename1, D time){
    this->PDCD_solver(beta_0, x0,y0, val_tau, max_nb_outer,  p_N, val_lambda_f, filename1, time);
}

};

#endif


text:
/****************************************************************
  gsl.h

  Defines wrappers for gsl functions.

  Language: C++

  Mark A. Caprio
  University of Notre Dame

  - 3/26/17 (mac): Created, based on code from spncci utilities (aem).

****************************************************************/

#ifndef MCUTILS_GSL_H_
#define MCUTILS_GSL_H_

#include <gsl/gsl_sf_gamma.h>

namespace mcutils {

  // GSL factorial wrappers
  //
  // Use more efficient int implementation for int arguments, else use
  // generic gamma function for double arguments.

  inline double Factorial(int x)
  {
    return gsl_sf_fact(x);
  }

  inline double Factorial(double x)
  {
    return gsl_sf_gamma(x+1);
  }

  // GSL combinatorial function wrappers

  inline int Choose(int x, int y)
  {
    int choose=0;
    if ((x>=y)&&(y>=0))
      {
        gsl_sf_result result;
        gsl_sf_choose_e(x,y,&result);
        choose=result.val;
      }
    return choose;
  }

} // namespace

#endif

max_stars_repo_name: nd-nuclear-theory/mcutils

max_issues_repo_name: nd-nuclear-theory/mcutils

max_forks_repo_name: nd-nuclear-theory/mcutils


text:
#pragma once

#include <cmath>

#include <gsl/gsl_math.h>

bool to_fraction(double num, long& numerator, long& denominator){
    double sign = num < 0 ? -1 : 1;
    double g = std::abs(num);
    long a = 0, b = 1, c = 1, d = 0;
    long s;
    int iter = 0;
    numerator = 0, denominator = 1;
    while (iter++ < 1e6){
        s = (long) std::floor(g);
        numerator = a + s*c;
        denominator = b + s*d;
        a = c;
        b = d;
        c = numerator;
        d = denominator;
        g = 1.0/(g-s);
        if(gsl_fcmp(sign*numerator, num*denominator, 1e-13) == 0){
            numerator *= sign;
            return true;
        }
    }
    return false;
}

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
#include <math.h>

#include <gsl/gsl_blas.h>

#include "qdm.h"

#define TAU_INITIAL_GUESS 0.5
#define TAU_EPS 0.001
#define TAU_ITERATION_MAX 1000

#define RESET_SIZE 30
const double RESET_VALUES[RESET_SIZE] = {
  0.010000000000000000,
  0.010100000000000000,
  0.010200000000000001,
  0.010300000000000000,
  0.010400000000000000,
  0.010500000000000001,
  0.010600000000000000,
  0.010699999999999999,
  0.010800000000000001,
  0.010900000000000000,
  0.010999999999999999,
  0.010999999999999999,
  0.150714285714285717,
  0.290428571428571425,
  0.430142857142857160,
  0.569857142857142840,
  0.709571428571428520,
  0.849285714285714310,
  0.988999999999999990,
  0.988999999999999990,
  0.989099999999999979,
  0.989199999999999968,
  0.989299999999999957,
  0.989399999999999946,
  0.989500000000000046,
  0.989600000000000035,
  0.989700000000000024,
  0.989800000000000013,
  0.989900000000000002,
  0.989999999999999991,
};

/* Helper: Compute ispline by mmm dot product. */
static
int
ispline_mmm(
    double *result,

    double tau,
    size_t spline_df,
    const gsl_vector *knots,
    const gsl_vector *mmm
)
{
  int status = 0;

  size_t m = (knots->size - spline_df) + 1;
  gsl_vector *ispline = gsl_vector_alloc(m);

  qdm_ispline_vector(ispline, tau, spline_df, knots);
  status = gsl_blas_ddot(ispline, mmm, result);
  if (status != 0) {
    goto cleanup;
  }

cleanup:
  gsl_vector_free(ispline);

  return status;
}

/* Helper: Compute mspline by mmm dot product. */
static
int
mspline_mmm(
    double *result,

    double tau,
    size_t spline_df,
    const gsl_vector *knots,
    const gsl_vector *mmm
)
{
  int status = 0;

  size_t m = (knots->size - spline_df) + 1;
  gsl_vector *mspline = gsl_vector_alloc(m);

  qdm_mspline_vector(mspline, tau, spline_df, knots);
  status = gsl_blas_ddot(mspline, mmm, result);
  if (status != 0) {
    goto cleanup;
  }

cleanup:
  gsl_vector_free(mspline);

  return status;
}

int
qdm_find_tau(
    double *result,

    double v,
    size_t spline_df,
    const gsl_vector *knots,
    const gsl_vector *mmm
)
{
  int status = 0;

  double tau = TAU_INITIAL_GUESS;
  double qi_u = 0;
  double qm_u = 0;

  size_t reset_i = 0;
  for (size_t i = 0; i < TAU_ITERATION_MAX; i++) {
    if (reset_i >= RESET_SIZE) {
      tau = TAU_INITIAL_GUESS;

      break;
    }

    /* Calculate position... */
    status = ispline_mmm(&qi_u, tau, spline_df, knots, mmm);
    if (status != 0) {
      goto cleanup;
    }

    /* Calculate slope... */
    status = mspline_mmm(&qm_u, tau, spline_df, knots, mmm);
    if (status != 0) {
      goto cleanup;
    }

    /* Check if the update is within our desired interval [0, 1]. */
    double update = tau - (qi_u - v) / qm_u;
    if (update < 0 || update > 1) {
      tau = RESET_VALUES[reset_i];
      reset_i++;

      continue;
    } else {
      tau = update;
    }

    if (fabs(qi_u - v) <= TAU_EPS) {
      break;
    }
  }

  *result = tau;

cleanup:

  return status;
}

int
qdm_logl(
    double *log_likelihood,
    double *tau,

    double v,
    size_t spline_df,
    const gsl_vector *knots,
    const gsl_vector *mmm,

    double tau_low,
    double tau_high,
    double xi_low,
    double xi_high
)
{
  int status = 0;

  double low_threshold = 0;
  double high_threshold = 0;

  status = ispline_mmm(&low_threshold, tau_low, spline_df, knots, mmm);
  if (status != 0) {
    goto cleanup;
  }

  status = ispline_mmm(&high_threshold, tau_high, spline_df, knots, mmm);
  if (status != 0) {
    goto cleanup;
  }

  if (v <= low_threshold) {
    double sigma_low = 0;

    status = mspline_mmm(&sigma_low, tau_low, spline_df, knots, mmm);
    if (status != 0) {
      goto cleanup;
    }

    sigma_low *= tau_low;

    double z = low_threshold - v;

    *log_likelihood = log(tau_low / sigma_low * pow(1 + (xi_low * z) / sigma_low, -1 / xi_low - 1));
    *tau = (1 - (1 - pow(1 + xi_low * z / sigma_low, -1 / xi_low))) * tau_low;
  } else if (v >= high_threshold) {
    double sigma_high = 0;

    status = mspline_mmm(&sigma_high, tau_high, spline_df, knots, mmm);
    if (status != 0) {
      goto cleanup;
    }

    sigma_high *= (1 - tau_high);

    double z = v - high_threshold;

    *log_likelihood = log((1 - tau_high) / sigma_high * pow(1 + (xi_high * z) / sigma_high, -1 / xi_high - 1));
    *tau = (1 - pow(1 + xi_high * z / sigma_high, -1 / xi_high)) * (1 - tau_high) + tau_high;
  } else {
    double d = 0;

    status = qdm_find_tau(tau, v, spline_df, knots, mmm);
    if (status != 0) {
      goto cleanup;
    }

    status = mspline_mmm(&d, *tau, spline_df, knots, mmm);
    if (status != 0) {
      goto cleanup;
    }

    *log_likelihood = log(1 / d);
  }

cleanup:

  return status;
}

int
qdm_logl_2(
    double *log_likelihood,
    double *tau,

    double x,
    double y,

    double tau_low,
    double tau_high,

    double xi_low,
    double xi_high,

    size_t spline_df,

    const gsl_matrix *theta,
    const gsl_vector *knots
)
{
  int status = 0;

  double xi_data[2] = {1, x};
  gsl_vector_view xi = gsl_vector_view_array(xi_data, 2);

  double mmm_data[theta->size2];
  gsl_vector_view mmm = gsl_vector_view_array(mmm_data, theta->size2);

  status = gsl_blas_dgemv(CblasTrans, 1, theta, &xi.vector, 0, &mmm.vector);
  if (status != 0) {
    return status;
  }

  status = qdm_logl(
      log_likelihood,
      tau,

      y,
      spline_df,
      knots,
      &mmm.vector,

      tau_low,
      tau_high,
      xi_low,
      xi_high
  );
  if (status != 0) {
    return status;
  }

  return status;
}

void
qdm_logl_3(
    double *log_likelihood,
    double *tau,

    double x,
    double y,

    const qdm_tau *t,
    const gsl_vector *xi,

    const gsl_matrix *theta
)
{
  double xi_intercept_data[2] = {1, x};
  gsl_vector_view xi_intercept = gsl_vector_view_array(xi_intercept_data, 2);

  double mmm_data[theta->size2];
  gsl_vector_view mmm = gsl_vector_view_array(mmm_data, theta->size2);

  gsl_blas_dgemv(CblasTrans, 1, theta, &xi_intercept.vector, 0, &mmm.vector);

  double low_threshold = qdm_tau_ispline_mmm(t, t->low, &mmm.vector);
  double high_threshold = qdm_tau_ispline_mmm(t, t->high, &mmm.vector);

  /*
  fprintf(stderr, "x    : %f\n", x);
  fprintf(stderr, "y    : %f\n", y);

  fprintf(stderr, "tl   : %f\n", t->low);
  fprintf(stderr, "th   : %f\n", t->high);

  fprintf(stderr, "low  : %f\n", low_threshold);
  fprintf(stderr, "high : %f\n", high_threshold);

  qdm_matrix_csv_fwrite(stderr, theta);
  */

  double xi_low = gsl_vector_get(xi, 0);
  double xi_high = gsl_vector_get(xi, 1);

  if (y <= low_threshold) {
    double sigma_low = qdm_tau_mspline_mmm(t, t->low, &mmm.vector) * t->low;
    double z = low_threshold - y;

    *log_likelihood = log(t->low / sigma_low * pow(1 + (xi_low * z) / sigma_low, -1 / xi_low - 1));
    *tau = (1 - (1 - pow(1 + xi_low * z / sigma_low, -1 / xi_low))) * t->low;
  } else if (y >= high_threshold) {
    double sigma_high = qdm_tau_mspline_mmm(t, t->high, &mmm.vector) * (1 - t->high);
    double z = y - high_threshold;

    *log_likelihood = log((1 - t->high) / sigma_high * pow(1 + (xi_high * z) / sigma_high, -1 / xi_high - 1));
    *tau = (1 - pow(1 + xi_high * z / sigma_high, -1 / xi_high)) * (1 - t->high) + t->high;
  } else {
    *tau = qdm_tau_find(t, y, &mmm.vector);
    *log_likelihood = log(1 / qdm_tau_mspline_mmm(t, *tau, &mmm.vector));
  }

  /*
  fprintf(stderr, "ll   : %f\n", *log_likelihood);
  fprintf(stderr, "tau  : %f\n", *tau);
  */
}

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
/*
 *
 * Copyright (c) Kresimir Fresl 2002
 *
 * Distributed under the Boost Software License, Version 1.0.
 * (See accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * Author acknowledges the support of the Faculty of Civil Engineering,
 * University of Zagreb, Croatia.
 *
 */

#ifndef BOOST_NUMERIC_BINDINGS_BLAS_DETAIL_CBLAS_H
#define BOOST_NUMERIC_BINDINGS_BLAS_DETAIL_CBLAS_H

//
// MKL-specific CBLAS include
//
#if defined BOOST_NUMERIC_BINDINGS_BLAS_MKL

extern "C" {
#include <mkl_cblas.h>
//#include <mkl_service.h>
//
// mkl_types.h defines P4 macro which breaks MPL, undefine it here.
//
#undef P4
}

//
// Default CBLAS include
//
#else

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

#endif
#endif

max_stars_repo_name: siconos/siconos-deb

max_issues_repo_name: siconos/siconos-deb

max_forks_repo_name: siconos/siconos-deb


text:
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

int normalizeVector( gsl_vector *vector )
{
	double min = gsl_vector_min( vector );
	double norm = gsl_vector_max( vector ) - min;
	gsl_vector_add_constant( vector, -min );
	gsl_vector_scale( vector, 1/norm );
	return 0;
}

max_stars_repo_name: jakeinater/spectralILU

max_issues_repo_name: jakeinater/spectralILU

max_forks_repo_name: jakeinater/spectralILU