text:
#pragma once
#include <gsl/gsl_vector.h>
#include "BooleanNodes.h"
#include "BooleanDAG.h"
#include <map>
#include "SymbolicEvaluator.h"
#include "ConflictGenerator.h"

#include <iostream>
#include "Util.h"

class SmartConflictGenerator: public ConflictGenerator {
	SymbolicEvaluator* eval;
	map<int, int>& imap;
	BooleanDAG* dag;
	set<int>& ignoredBoolNodes;
	set<int>& ctrlNodes;

	vector<set<int>> dependentInputs;
	vector<set<int>> dependentCtrls;

public:
	ConflictGenerator(SymbolicEvaluator* eval_, map<int, int>& imap_, BooleanDAG* dag_, set<int>& ignoredBoolNodes_, set<int>& ctrlNodes_): eval(eval_), imap(imap_), dag(dag_), ignoredBoolNodes(ignoredBoolNodes_), ctrlNodes(ctrlNodes_) {
		// process the dag and collect dependencies
		for (BooleanDAG::iterator node_it = dag->begin(); node_it != dag->end(); ++node_it) {
			bool_node* n = *node_it;
			set<int> inputs; set<int> ctrls;

			const vector<bool_node*>& parents = n->parents();
			for (int i = 0; i < parents.size(); i++) {
				bool_node* parent = parents[i];
				set<int>& parentInputs = dependentInputs[parent->id];
				set<int>& parentCtrls = dependentCtrls[parent->id];
				if (parent->getOtype() == OutType::BOOL) {
					if (ignoredBoolNodes.find(parent->id) == ignoredBoolNodes.end()) {
						inputs.insert(parent->id);
						//inputs.insert(parentInputs.begin(), parentInputs.end());
					}
				} else {
					inputs.insert(parentInputs.begin(), parentInputs.end());
					ctrls.insert(parentCtrls.begin(), parentCtrls.end());
				}

				if (ctrlNodes.find(parent->id) != ctrlNodes.end()) {
					ctrls.insert(parent->id);
				}
			}
			dependentInputs.push_back(inputs);
			dependentCtrls.push_back(ctrls);
		}

	}

	bool isConflict(set<int> influentialControls, set<int> myControls) {
		for (int i : influentialControls) {
			if (myControls.find(i) != myControls.end()) {
				return true;
			}
		}
		return false;
	}

	virtual vector<pair<int, int>> getConflicts(gsl_vector* state, vector<vector<int>>& allInputs, vector<int>& instanceIds, int rowid, int colid) {
		vector<pair<int, int>> conflicts;
		for (int i = 0; i < allInputs.size(); i++) {
			const map<int, int>& nodeValsMap = Util::getNodeToValMap(imap, allInputs[i]);
			eval->run(state, nodeValsMap);

			// First, find a failing node
			bool_node* failedNode = NULL;
			for (BooleanDAG::iterator node_it = dag->begin(); node_it != dag->end(); ++node_it) {
				bool_node* n = *node_it;
				if (n->type == bool_node::ASSERT) {
					if (!eval->check(n->mother, 1)) {
						failedNode = n;
						break;
					}
				}
				if (n->getOtype() == OutType::BOOL) {
					auto it = nodeValsMap.find(n->id);
					if (it != nodeValsMap.end()) {
						int val = it->second;
						if (!eval->check(n, val)) {
							failedNode = n;
							break;
						}
					}
				}
			}

			Assert(failedNode != NULL, "No conflict?");
			cout << failedNode->lprint() << endl;
			set<int> conflictNodes;
			// First, add all dependent inputs
			set<int>& depInputs = dependentInputs[failedNode->id];
			conflictNodes.insert(depInputs.begin(), depInputs.end());
			cout << "Dep inputs:" << endl;
			for (auto it = depInputs.begin(); it != depInputs.end(); it++) {
				cout << (*dag)[(*it)]->lprint() << endl;
			}

			set<int>& depCtrls = dependentCtrls[failedNode->id];
			cout << "Dep ctrls:" << endl;
			for (auto it = depCtrls.begin(); it != depCtrls.end(); it++) {
				cout << (*dag)[(*it)]->lprint() << endl;
			}
			// Next, collect all boolean nodes that are influenced by the controls that affect the failed node
			for (BooleanDAG::iterator node_it = dag->begin(); node_it != dag->end(); ++node_it) {
				bool_node* n = *node_it;
				if (n->type == bool_node::ASSERT) {
					if (isConflict(depCtrls, dependentCtrls[n->id])) {
						conflictNodes.insert(n->id);
					}
				}
				if (n->getOtype() == OutType::BOOL) {
					if (isConflict(depCtrls, dependentCtrls[n->id])) {
						conflictNodes.insert(n->id);
					}
				}
			}

			cout << "Final conflict nodes" << endl;

			//conflictNodes.insert(colid);

			// Process the conflicts
			for (int j = 0; j < imap.size(); j++) {
				if (allInputs[i][j] != EMPTY) {
					if (ignoredBoolNodes.find(imap[j]) == ignoredBoolNodes.end()) {
						if (imap[j] < 0 || conflictNodes.find(imap[j]) != conflictNodes.end()) {
							if (imap[j] >= 0) {
								cout << (*dag)[imap[j]]->lprint() << " " << allInputs[i][j] << endl;
							}
							conflicts.push_back(make_pair(instanceIds[i], j));
						}
					}
				}
			}
		}
		return conflicts;
	}

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


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

#ifndef __GSL_STATISTICS_LONG_DOUBLE_H__
#define __GSL_STATISTICS_LONG_DOUBLE_H__

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

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

__BEGIN_DECLS

GSL_EXPORT double gsl_stats_long_double_mean (const long double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_long_double_variance (const long double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_long_double_sd (const long double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_long_double_variance_with_fixed_mean (const long double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_long_double_sd_with_fixed_mean (const long double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_long_double_absdev (const long double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_long_double_skew (const long double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_long_double_kurtosis (const long double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_long_double_lag1_autocorrelation (const long double data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_long_double_covariance (const long double data1[], const size_t stride1,const long double data2[], const size_t stride2, const size_t n);

GSL_EXPORT double gsl_stats_long_double_variance_m (const long double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_long_double_sd_m (const long double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_long_double_absdev_m (const long double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_long_double_skew_m_sd (const long double data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_long_double_kurtosis_m_sd (const long double data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_long_double_lag1_autocorrelation_m (const long double data[], const size_t stride, const size_t n, const double mean);

GSL_EXPORT double gsl_stats_long_double_covariance_m (const long double data1[], const size_t stride1,const long double data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);

/* DEFINED FOR FLOATING POINT TYPES ONLY */

GSL_EXPORT double gsl_stats_long_double_wmean (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_long_double_wvariance (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_long_double_wsd (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_long_double_wvariance_with_fixed_mean (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_long_double_wsd_with_fixed_mean (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_long_double_wabsdev (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_long_double_wskew (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_long_double_wkurtosis (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_long_double_wvariance_m (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double wmean);
GSL_EXPORT double gsl_stats_long_double_wsd_m (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double wmean);
GSL_EXPORT double gsl_stats_long_double_wabsdev_m (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double wmean);
GSL_EXPORT double gsl_stats_long_double_wskew_m_sd (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double wmean, const double wsd);
GSL_EXPORT double gsl_stats_long_double_wkurtosis_m_sd (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double wmean, const double wsd);

/* END OF FLOATING POINT TYPES */

GSL_EXPORT double gsl_stats_long_double_pvariance (const long double data1[], const size_t stride1, const size_t n1, const long double data2[], const size_t stride2, const size_t n2);
GSL_EXPORT double gsl_stats_long_double_ttest (const long double data1[], const size_t stride1, const size_t n1, const long double data2[], const size_t stride2, const size_t n2);

GSL_EXPORT long double gsl_stats_long_double_max (const long double data[], const size_t stride, const size_t n);
GSL_EXPORT long double gsl_stats_long_double_min (const long double data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_long_double_minmax (long double * min, long double * max, const long double data[], const size_t stride, const size_t n);

GSL_EXPORT size_t gsl_stats_long_double_max_index (const long double data[], const size_t stride, const size_t n);
GSL_EXPORT size_t gsl_stats_long_double_min_index (const long double data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_long_double_minmax_index (size_t * min_index, size_t * max_index, const long double data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_long_double_median_from_sorted_data (const long double sorted_data[], const size_t stride, const size_t n) ;
GSL_EXPORT double gsl_stats_long_double_quantile_from_sorted_data (const long double sorted_data[], const size_t stride, const size_t n, const double f) ;

__END_DECLS

#endif /* __GSL_STATISTICS_LONG_DOUBLE_H__ */

ext: h

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <cblas.h>

int max(int x, int y){
   if(x > y) return x;
   else return y;
}
void PrintMatrix(float* pMatrix, const size_t nR, const size_t nC, const CBLAS_ORDER Order) {
    unsigned int i, j;
    if (Order == CblasRowMajor) {
        for (i = 0; i < nR; i++) {
            for (j = 0; j < nC; j++) {
                printf("%f \t ", pMatrix[i * nC + j]); // !!!
            }
            printf("\n"); // !!!
        }
    } else {
        for (i = 0; i < nR; i++) {
            for (j = 0; j < nC; j++) {
                printf("%f \t ", pMatrix[i + j* nR ]); // !!!
            }
            printf("\n"); // !!!
        }
    }
    printf("\n"); // !!!
}

int main(void)
{
    const int m = 4;
    const int n = 5;
    const int k = 1;

    float A[] = { 8, 4, 7, 3, 5, 1, 1, 3, 2, 1, 2, 3, 2, 0, 1, 1, 2, 3, 4, 1};
    float B[5] = { -1, 2, -1, 1, 2 };

    float alpha = 1.0, beta = 0.0;
    int lda, ldb, ldc;

    float * C = (float*) malloc(m * k * sizeof(float));
    for (int i = 0; i < m*k; i++) C[i] = 0.0;

    //working cblas_sgemv(CblasRowMajor, CblasNoTrans, m, n, alpha, A, n, B, k, beta, C, k);
    //lda = max(1,n); ldb = max(1,k); ldc = max(1,k);  // refer manual and switch n and k
    //cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m,k,n, alpha, A, lda, B, ldb, beta, C, ldc);

    // working cblas_sgemv(CblasColMajor, CblasNoTrans, m, n, alpha, A, m, B, k, beta, C, k);
    // working cblas_sgemv(CblasColMajor, CblasTrans, m, n, alpha, A, m, B, k, beta, C, k);
    lda = max(1,n); ldb = max(1,n); ldc = max(1,m);  // refer manual and switch n and k
    cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans, m,k,n, alpha, A, lda, B, ldb, beta, C, ldc);

    PrintMatrix(C, m, k, CblasRowMajor);

    free(C);

    return 0;
}

ext: c

max_stars_repo_name: ramcn/flappie-x-cpp

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: ramcn/flappie-x-cpp

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: ramcn/flappie-x-cpp

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


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

/* This is a maximally equidistributed combined, collision free
   Tausworthe generator, with a period ~2^{113}. The sequence is,

   x_n = (z1_n ^ z2_n ^ z3_n ^ z4_n)

   b = (((z1_n <<  6) ^ z1_n) >> 13)
   z1_{n+1} = (((z1_n & 4294967294) << 18) ^ b)
   b = (((z2_n <<  2) ^ z2_n) >> 27)
   z2_{n+1} = (((z2_n & 4294967288) <<  2) ^ b)
   b = (((z3_n << 13) ^ z3_n) >> 21)
   z3_{n+1} = (((z3_n & 4294967280) <<  7) ^ b)
   b = (((z4_n <<  3)  ^ z4_n) >> 12)
   z4_{n+1} = (((z4_n & 4294967168) << 13) ^ b)

   computed modulo 2^32. In the formulas above '^' means exclusive-or
   (C-notation), not exponentiation.
   The algorithm is for 32-bit integers, hence a bitmask is used to clear
   all but least significant 32 bits, after left shifts, to make the code
   work on architectures where integers are 64-bit.

   The generator is initialized with
   zi = (69069 * z{i+1}) MOD 2^32 where z0 is the seed provided
   During initialization a check is done to make sure that the initial seeds
   have a required number of their most significant bits set.
   After this, the state is passed through the RNG 10 times to ensure the
   state satisfies a recurrence relation.

   References:
   P. L'Ecuyer, "Tables of Maximally-Equidistributed Combined LFSR Generators",
   Mathematics of Computation, 68, 225 (1999), 261--269.
     http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme2.ps
   P. L'Ecuyer, "Maximally Equidistributed Combined Tausworthe Generators",
   Mathematics of Computation, 65, 213 (1996), 203--213.
     http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme.ps
   the online version of the latter contains corrections to the print version.
*/

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

#define LCG(n) ((69069UL * n) & 0xffffffffUL)
#define MASK 0xffffffffUL

static inline unsigned long int taus113_get (void *vstate);
static double taus113_get_double (void *vstate);
static void taus113_set (void *state, unsigned long int s);

typedef struct
{
  unsigned long int z1, z2, z3, z4;
}
taus113_state_t;

static inline unsigned long
taus113_get (void *vstate)
{
  taus113_state_t *state = (taus113_state_t *) vstate;
  unsigned long b1, b2, b3, b4;

  b1 = ((((state->z1 << 6UL) & MASK) ^ state->z1) >> 13UL);
  state->z1 = ((((state->z1 & 4294967294UL) << 18UL) & MASK) ^ b1);

  b2 = ((((state->z2 << 2UL) & MASK) ^ state->z2) >> 27UL);
  state->z2 = ((((state->z2 & 4294967288UL) << 2UL) & MASK) ^ b2);

  b3 = ((((state->z3 << 13UL) & MASK) ^ state->z3) >> 21UL);
  state->z3 = ((((state->z3 & 4294967280UL) << 7UL) & MASK) ^ b3);

  b4 = ((((state->z4 << 3UL) & MASK) ^ state->z4) >> 12UL);
  state->z4 = ((((state->z4 & 4294967168UL) << 13UL) & MASK) ^ b4);

  return (state->z1 ^ state->z2 ^ state->z3 ^ state->z4);

}

static double
taus113_get_double (void *vstate)
{
  return taus113_get (vstate) / 4294967296.0;
}

static void
taus113_set (void *vstate, unsigned long int s)
{
  taus113_state_t *state = (taus113_state_t *) vstate;

  if (!s)
    s = 1UL;                    /* default seed is 1 */

  state->z1 = LCG (s);
  if (state->z1 < 2UL)
    state->z1 += 2UL;
  state->z2 = LCG (state->z1);
  if (state->z2 < 8UL)
    state->z2 += 8UL;
  state->z3 = LCG (state->z2);
  if (state->z3 < 16UL)
    state->z3 += 16UL;
  state->z4 = LCG (state->z3);
  if (state->z4 < 128UL)
    state->z4 += 128UL;

  /* Calling RNG ten times to satify recurrence condition */
  taus113_get (state);
  taus113_get (state);
  taus113_get (state);
  taus113_get (state);
  taus113_get (state);
  taus113_get (state);
  taus113_get (state);
  taus113_get (state);
  taus113_get (state);
  taus113_get (state);

  return;
}

static const gsl_rng_type taus113_type = {
  "taus113",                    /* name */
  0xffffffffUL,                 /* RAND_MAX */
  0,                            /* RAND_MIN */
  sizeof (taus113_state_t),
  &taus113_set,
  &taus113_get,
  &taus113_get_double
};

const gsl_rng_type *gsl_rng_taus113 = &taus113_type;

/*  Rules for analytic calculations using GNU Emacs Calc:
    (used to find the values for the test program)

  [ LCG(n) := n * 69069 mod (2^32) ]

  [ b1(x) := rsh(xor(lsh(x, 6), x), 13),
  q1(x) := xor(lsh(and(x, 4294967294), 18), b1(x)),
  b2(x) := rsh(xor(lsh(x, 2), x), 27),
  q2(x) := xor(lsh(and(x, 4294967288), 2), b2(x)),
  b3(x) := rsh(xor(lsh(x, 13), x), 21),
  q3(x) := xor(lsh(and(x, 4294967280), 7), b3(x)),
  b4(x) := rsh(xor(lsh(x, 3), x), 12),
  q4(x) := xor(lsh(and(x, 4294967168), 13), b4(x))
  ]

  [ S([z1,z2,z3,z4]) := [q1(z1), q2(z2), q3(z3), q4(z4)] ]
*/

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
// tcc -lgsl -lglut -lGL -run gsl_opengl.c
// clang -lgsl -framework GLUT -framework OpenGL gsl_opengl.c
#include <gsl/gsl_qrng.h>
#include <GL/glut.h>

// #pragma comment(lib,"glut32")
// #pragma comment(lib,"gsl")
// #pragma comment(lib,"opengl32")
// #pragma comment(lib,"freeglut")

void glutInitialize()
{
    glClearColor(1,1,1,1);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glTranslatef(-1,-1,0);
    glScalef(2,2,1);
    glViewport(0,0,512,512);
    glPointSize(3.0);
}

void Display()
{
    glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
    const int n = 100;
    gsl_qrng* qrng = gsl_qrng_alloc(gsl_qrng_sobol,2);
    glBegin(GL_POINTS);
    glColor3f(0,0,1);
    for( int i=0; i<n; i++ )
    {
        double s[2];
        gsl_qrng_get(qrng,s);
        glVertex3f(s[0],s[1],0.0);
    }
    glEnd();
    gsl_qrng_free(qrng);
    glutSwapBuffers();
}

int main(int argc, char** argv)
{
    glutInit(&argc,argv);
    glutInitDisplayMode(GLUT_RGBA|GLUT_DOUBLE);
    glutInitWindowPosition((glutGet(GLUT_SCREEN_WIDTH)-512)/2,
                       (glutGet(GLUT_SCREEN_HEIGHT)-640)/2);
    glutInitWindowSize(512,512);
    glutCreateWindow("GSL Quasi Random Number Generator");
    glutInitialize();
    glutDisplayFunc(Display);
    glutMainLoop();
    return 0;
}


text:
/*
** compute edge densities
**
** G.Lohmann, Feb 2015
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_histogram.h>

#include "viaio/Vlib.h"
#include "viaio/VImage.h"
#include "viaio/mu.h"

#ifdef _OPENMP
#include <omp.h>
#endif /*_OPENMP*/

#define SQR(x) ((x) * (x))
#define ABS(x) ((x) > 0 ? (x) : -(x))

extern void GetRank(float *data,gsl_vector *v,gsl_permutation *perm,gsl_permutation *rank);
extern float EdgeCorr(const float *data1,const float *data2,int,int);

int GetAddr(VImage mapimage,int bi,int ri,int ci,int m,int k,int l)
{
  int b = bi+m;
  int r = ri+k;
  int c = ci+l;
  if (b < 0 || r < 0 || c < 0) return -1;
  if (b >= VImageNBands(mapimage) || r >= VImageNRows(mapimage) || c >= VImageNColumns(mapimage)) return -1;
  int ii = (long)VPixel(mapimage,b,r,c,VInteger);
  return ii;
}

int VEdgeNeigbours(size_t id,VImage map,VImage mapimage,int adjdef,int *x)
{
  int nslices = VImageNBands(mapimage);
  int nrows = VImageNRows(mapimage);
  int ncols = VImageNColumns(mapimage);

  int bi = VPixel(map,0,0,id,VShort);
  int ri = VPixel(map,0,1,id,VShort);
  int ci = VPixel(map,0,2,id,VShort);
  if (bi < 1 || bi >= nslices-2) return 0;
  if (ri < 1 || ri >= nrows-2) return 0;
  if (ci < 1 || ci >= ncols-2) return 0;

  int wn=1,rad2=0;
  if (adjdef == 3) {
    wn = 2;
    rad2 = 2*2;
  }

  int n=0;
  int k,l,m;
  for (m=-wn; m<=wn; m++) {
    for (k=-wn; k<=wn; k++) {
      for (l=-wn; l<=wn; l++) {
      	int jj = GetAddr(mapimage,bi,ri,ci,m,k,l);
      	if (jj < 0) continue;  /* outside of brain mask */

      	if (adjdef == 0) {     /* 6 adjacency */
      	  if (ABS(m)+ABS(k)+ABS(l) > 1) continue;
      	}
      	if (adjdef == 1) {     /* 18-adjacency */
      	  if (ABS(m) > 0 && ABS(k) > 0 && ABS(l) > 0) continue;
      	}
	if (adjdef == 3) {     /* sphere */
	  if (m*m + k*k + l*l > rad2) continue;
	}
      	x[n] = jj;
      	n++;
      }
    }
  }
  return n;
}

void VCheckImage(VImage src)
{
  VAttrList out_list = VCreateAttrList();
  VAppendAttr(out_list,"image",NULL,VImageRepn,src);
  FILE *out_file = fopen("test.v","w");
  VWriteFile (out_file, out_list);
}

size_t EdgeDensity(float *E,int *I,int *J,size_t nedges_estimated,
		   gsl_histogram *TedHist,gsl_matrix_float *SNR1,gsl_matrix_float *SNR2,int n1,int n2,
		   VImage roi,VImage map,VImage mapimage,int adjdef,float elength,float zthreshold,
		   int nperm,int numperm,int step,float noise_cutoff,int metric)
{
  size_t i;
  size_t nvox=SNR1->size1;
  int nt=SNR1->size2;
  size_t progress=0;
  int rad2 = (int)(elength*elength);
  float tiny=1.0e-6;
  gsl_set_error_handler_off ();

  int iisel=4;

  /* size of neighbourhood */
  int nadj=0;
  if (adjdef == 0) nadj = 7;
  else if (adjdef == 1) nadj = 19;
  else if (adjdef == 2) nadj = 27;
  else if (adjdef == 3) nadj = 33;
  else VError(" illegal adjdef %d",adjdef);
  size_t nedges=0;
  int minadj = 1;

  /* histogram init */
  size_t nbins = gsl_histogram_bins (TedHist);
  double hmax = gsl_histogram_max (TedHist);
  double hmin = gsl_histogram_min (TedHist);

#pragma omp parallel for shared(nedges,progress,TedHist,E,I,J) schedule(dynamic)
  for (i=0; i<nvox; i+=step) {
    if (i%1000 == 0) fprintf(stderr," %ld000\r",(long)progress++);

    int bi = (int)VPixel(map,0,0,i,VShort);
    int ri = (int)VPixel(map,0,1,i,VShort);
    int ci = (int)VPixel(map,0,2,i,VShort);

    int iselect=0;

    int roiflagi = 0;
    if (roi != NULL) {
      if (VGetPixel(roi,bi,ri,ci) > 0.5) roiflagi = 1;
    }

    const float *datax1 = gsl_matrix_float_const_ptr(SNR1,i,0);
    const float *datax2 = gsl_matrix_float_const_ptr(SNR2,i,0);

    int *x = (int *) VCalloc(nadj,sizeof(int));
    int *y = (int *) VCalloc(nadj,sizeof(int));

    gsl_histogram *tmphist = gsl_histogram_alloc (nbins);
    gsl_histogram_set_ranges_uniform (tmphist,hmin,hmax);
    gsl_histogram_reset(tmphist);

    int nadjx = VEdgeNeigbours(i,map,mapimage,adjdef,x);
    if (nadjx < minadj) continue;

    size_t j=0;
    for (j=0; j<i; j+=step) {
      int bj = (int)VPixel(map,0,0,j,VShort);
      int rj = (int)VPixel(map,0,1,j,VShort);
      int cj = (int)VPixel(map,0,2,j,VShort);
      int d = SQR(bi-bj) + SQR(ri-rj) + SQR(ci-cj);
      if (d < rad2) continue;

      int roiflagj = 0;
      if (roi != NULL) {
	if (VGetPixel(roi,bj,rj,cj) > 0.5) roiflagj = 1;
	if (roiflagi + roiflagj != 1) continue;
      }

      const float *datay1 = gsl_matrix_float_const_ptr(SNR1,j,0);
      const float *datay2 = gsl_matrix_float_const_ptr(SNR2,j,0);

      /* edge z-value */
      double z1 = EdgeCorr(datax1,datay1,nt,metric);
      double z2 = EdgeCorr(datax2,datay2,nt,metric);
      double z = (z1-z2);
      if (z < zthreshold) continue;

      int nadjy = VEdgeNeigbours(j,map,mapimage,adjdef,y);
      if (nadjy < minadj) continue;

      /* use every 'step' value, only when estimating noise cutoff, i.e. step > 1 */
      if (iselect%iisel == 0 && step > 1) {
	iselect = 0;
      }
      else if (iselect%iisel != 0 && step > 1) {
	iselect++;
	continue;
      }

      /* inspect local neighbourhood */
      int r,s;
      double kx=0.0,nx=0.0;
      for (s=0; s<nadjx; s++) {
      	size_t ii = (size_t) x[s];
        int bii = (int)VPixel(map,0,0,ii,VShort);
        int rii = (int)VPixel(map,0,1,ii,VShort);
        int cii = (int)VPixel(map,0,2,ii,VShort);

	const float *dataxx1 = gsl_matrix_float_const_ptr(SNR1,ii,0);
	const float *dataxx2 = gsl_matrix_float_const_ptr(SNR2,ii,0);

      	for (r=0; r<nadjy; r++) {
      	  size_t jj = (size_t) y[r];
          int bjj = (int)VPixel(map,0,0,jj,VShort);
          int rjj = (int)VPixel(map,0,1,jj,VShort);
          int cjj = (int)VPixel(map,0,2,jj,VShort);

          int d2 = SQR(bii-bjj) + SQR(rii-rjj) + SQR(cii-cjj);
          if (d2 < rad2) continue;

	  const float *datayy1 = gsl_matrix_float_const_ptr(SNR1,jj,0);
	  const float *datayy2 = gsl_matrix_float_const_ptr(SNR2,jj,0);

	  /* edge z-value */
	  double z1 = EdgeCorr(dataxx1,datayy1,nt,metric);
	  double z2 = EdgeCorr(dataxx2,datayy2,nt,metric);
	  double z = (z1-z2);
          if (ABS(z) < tiny) continue;
      	  if (z > zthreshold) kx++;
          nx++;
     	}
      }
      if (nx < 1.0) continue;
      double ted = kx/nx;

      if (ted < hmin) ted = hmin;
      if (ted > hmax-tiny) ted = hmax-tiny;
      gsl_histogram_increment (tmphist,ted);

      if (E == NULL) continue;
      if (nperm != numperm) continue;
      if (ted < noise_cutoff) continue;

      #pragma omp critical
      {
	if (nedges >= nedges_estimated) {
	  VError(" alloc, nedges_estimated:  %lu",nedges_estimated);
	}
	E[nedges] = ted;    /* task-based edge density */
	I[nedges] = i;      /* voxel address i */
	J[nedges] = j;      /* voxel address j */
	nedges++;
      }
    }
#pragma omp critical
    {
      gsl_histogram_add (TedHist,tmphist);
    }
    gsl_histogram_free (tmphist);
    VFree(x);
    VFree(y);
  }
  return nedges;
}

/* estimate number of truly needed edges */
size_t EstimateEdges(float *E,int *I,int *J,size_t old_estimate,
		     gsl_histogram *TedHist,gsl_matrix_float *SNR1,gsl_matrix_float *SNR2,int n1,int n2,
		     VImage roi,VImage map,VImage mapimage,int adjdef,float elength,float zthr,
		     float noise_cutoff,int metric)
{
  float tmp_cutoff = -1.0;
  int step = 2;
  int nperm=0,numperm=1;

  gsl_histogram_reset(TedHist);
  EdgeDensity(E,I,J,old_estimate,TedHist,SNR1,SNR2,n1,n2,roi,map,mapimage,(int)adjdef,elength,zthr,
	      nperm,numperm,step,tmp_cutoff,metric);

  double sd = gsl_histogram_sigma (TedHist);
  if (sd < 1.0e-6) VError(" stdev of TedHist: %f",sd);

  size_t i0 = 0;
  if (gsl_histogram_find (TedHist,(double)noise_cutoff,&i0) != GSL_SUCCESS)
    VError (" err finding noise cutoff index, %f",noise_cutoff);

  size_t hbins = gsl_histogram_bins(TedHist);
  gsl_histogram_pdf *cdf = gsl_histogram_pdf_alloc(hbins);
  gsl_histogram_pdf_init (cdf,TedHist);

  double qx = 1.0 - cdf->sum[i0];
  if (qx > 1) qx = 1.0;
  size_t new_estimate = (size_t) (qx*(double)old_estimate);

  gsl_histogram_reset(TedHist);
  gsl_histogram_pdf_free (cdf);
  return new_estimate;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2015 Microsoft Corporation. All rights reserved.
//
// This code is licensed under the MIT License (MIT).
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.
//
///////////////////////////////////////////////////////////////////////////////

#ifndef GSL_SPAN_H
#define GSL_SPAN_H

#include <gsl/gsl_assert> // for Expects
#include <gsl/gsl_byte>   // for byte
#include <gsl/gsl_util>   // for narrow_cast, narrow

#include <algorithm> // for lexicographical_compare
#include <array>     // for array
#include <cstddef>   // for ptrdiff_t, size_t, nullptr_t
#include <iterator>  // for reverse_iterator, distance, random_access_...
#include <limits>
#include <stdexcept>
#include <type_traits> // for enable_if_t, declval, is_convertible, inte...
#include <utility>
#include <memory> // for std::addressof

#ifdef _MSC_VER
#pragma warning(push)

// turn off some warnings that are noisy about our Expects statements
#pragma warning(disable : 4127) // conditional expression is constant
#pragma warning(disable : 4702) // unreachable code

// Turn MSVC /analyze rules that generate too much noise. TODO: fix in the tool.
#pragma warning(disable : 26495) // uninitalized member when constructor calls constructor
#pragma warning(disable : 26446) // parser bug does not allow attributes on some templates

#if _MSC_VER < 1910
#pragma push_macro("constexpr")
#define constexpr /*constexpr*/
#define GSL_USE_STATIC_CONSTEXPR_WORKAROUND

#endif // _MSC_VER < 1910
#endif // _MSC_VER

// See if we have enough C++17 power to use a static constexpr data member
// without needing an out-of-line definition
#if !(defined(__cplusplus) && (__cplusplus >= 201703L))
#define GSL_USE_STATIC_CONSTEXPR_WORKAROUND
#endif // !(defined(__cplusplus) && (__cplusplus >= 201703L))

// GCC 7 does not like the signed unsigned missmatch (size_t ptrdiff_t)
// While there is a conversion from signed to unsigned, it happens at
// compiletime, so the compiler wouldn't have to warn indiscriminently, but
// could check if the source value actually doesn't fit into the target type
// and only warn in those cases.
#if __GNUC__ > 6
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wsign-conversion"
#endif

namespace gsl
{

// [views.constants], constants
constexpr const std::ptrdiff_t dynamic_extent = -1;

template <class ElementType, std::ptrdiff_t Extent = dynamic_extent>
class span;

// implementation details
namespace details
{
    template <class T>
    struct is_span_oracle : std::false_type
    {
    };

    template <class ElementType, std::ptrdiff_t Extent>
    struct is_span_oracle<gsl::span<ElementType, Extent>> : std::true_type
    {
    };

    template <class T>
    struct is_span : public is_span_oracle<std::remove_cv_t<T>>
    {
    };

    template <class T>
    struct is_std_array_oracle : std::false_type
    {
    };

    template <class ElementType, std::size_t Extent>
    struct is_std_array_oracle<std::array<ElementType, Extent>> : std::true_type
    {
    };

    template <class T>
    struct is_std_array : public is_std_array_oracle<std::remove_cv_t<T>>
    {
    };

    template <std::ptrdiff_t From, std::ptrdiff_t To>
    struct is_allowed_extent_conversion
        : public std::integral_constant<bool, From == To || From == gsl::dynamic_extent ||
                                                  To == gsl::dynamic_extent>
    {
    };

    template <class From, class To>
    struct is_allowed_element_type_conversion
        : public std::integral_constant<bool, std::is_convertible<From (*)[], To (*)[]>::value>
    {
    };

    template <class Span, bool IsConst>
    class span_iterator
    {
        using element_type_ = typename Span::element_type;

    public:
#ifdef _MSC_VER
        // Tell Microsoft standard library that span_iterators are checked.
        using _Unchecked_type = typename Span::pointer;
#endif

        using iterator_category = std::random_access_iterator_tag;
        using value_type = std::remove_cv_t<element_type_>;
        using difference_type = typename Span::index_type;

        using reference = std::conditional_t<IsConst, const element_type_, element_type_>&;
        using pointer = std::add_pointer_t<reference>;

        span_iterator() = default;

        constexpr span_iterator(const Span* span, typename Span::index_type idx) noexcept
            : span_(span), index_(idx)
        {}

        friend span_iterator<Span, true>;
        template <bool B, std::enable_if_t<!B && IsConst>* = nullptr>
        constexpr span_iterator(const span_iterator<Span, B>& other) noexcept
            : span_iterator(other.span_, other.index_)
        {}

        GSL_SUPPRESS(bounds.1) // NO-FORMAT: attribute
        constexpr reference operator*() const
        {
            Expects(index_ != span_->size());
            return *(span_->data() + index_);
        }

        constexpr pointer operator->() const
        {
            Expects(index_ != span_->size());
            return span_->data() + index_;
        }

        constexpr span_iterator& operator++()
        {
            Expects(0 <= index_ && index_ != span_->size());
            ++index_;
            return *this;
        }

        constexpr span_iterator operator++(int)
        {
            auto ret = *this;
            ++(*this);
            return ret;
        }

        constexpr span_iterator& operator--()
        {
            Expects(index_ != 0 && index_ <= span_->size());
            --index_;
            return *this;
        }

        constexpr span_iterator operator--(int)
        {
            auto ret = *this;
            --(*this);
            return ret;
        }

        constexpr span_iterator operator+(difference_type n) const
        {
            auto ret = *this;
            return ret += n;
        }

        friend constexpr span_iterator operator+(difference_type n, span_iterator const& rhs)
        {
            return rhs + n;
        }

        constexpr span_iterator& operator+=(difference_type n)
        {
            Expects((index_ + n) >= 0 && (index_ + n) <= span_->size());
            index_ += n;
            return *this;
        }

        constexpr span_iterator operator-(difference_type n) const
        {
            auto ret = *this;
            return ret -= n;
        }

        constexpr span_iterator& operator-=(difference_type n) { return *this += -n; }

        constexpr difference_type operator-(span_iterator rhs) const
        {
            Expects(span_ == rhs.span_);
            return index_ - rhs.index_;
        }

        constexpr reference operator[](difference_type n) const { return *(*this + n); }

        constexpr friend bool operator==(span_iterator lhs, span_iterator rhs) noexcept
        {
            return lhs.span_ == rhs.span_ && lhs.index_ == rhs.index_;
        }

        constexpr friend bool operator!=(span_iterator lhs, span_iterator rhs) noexcept
        {
            return !(lhs == rhs);
        }

        constexpr friend bool operator<(span_iterator lhs, span_iterator rhs) noexcept
        {
            return lhs.index_ < rhs.index_;
        }

        constexpr friend bool operator<=(span_iterator lhs, span_iterator rhs) noexcept
        {
            return !(rhs < lhs);
        }

        constexpr friend bool operator>(span_iterator lhs, span_iterator rhs) noexcept
        {
            return rhs < lhs;
        }

        constexpr friend bool operator>=(span_iterator lhs, span_iterator rhs) noexcept
        {
            return !(rhs > lhs);
        }

#ifdef _MSC_VER
        // MSVC++ iterator debugging support; allows STL algorithms in 15.8+
        // to unwrap span_iterator to a pointer type after a range check in STL
        // algorithm calls
        friend constexpr void _Verify_range(span_iterator lhs, span_iterator rhs) noexcept
        { // test that [lhs, rhs) forms a valid range inside an STL algorithm
            Expects(lhs.span_ == rhs.span_        // range spans have to match
                    && lhs.index_ <= rhs.index_); // range must not be transposed
        }

        constexpr void _Verify_offset(const difference_type n) const noexcept
        { // test that the iterator *this + n is a valid range in an STL
            // algorithm call
            Expects((index_ + n) >= 0 && (index_ + n) <= span_->size());
        }

        GSL_SUPPRESS(bounds.1) // NO-FORMAT: attribute
        constexpr pointer _Unwrapped() const noexcept
        { // after seeking *this to a high water mark, or using one of the
            // _Verify_xxx functions above, unwrap this span_iterator to a raw
            // pointer
            return span_->data() + index_;
        }

        // Tell the STL that span_iterator should not be unwrapped if it can't
        // validate in advance, even in release / optimized builds:
#if defined(GSL_USE_STATIC_CONSTEXPR_WORKAROUND)
        static constexpr const bool _Unwrap_when_unverified = false;
#else
        static constexpr bool _Unwrap_when_unverified = false;
#endif
        GSL_SUPPRESS(con.3) // NO-FORMAT: attribute // TODO: false positive
        constexpr void _Seek_to(const pointer p) noexcept
        { // adjust the position of *this to previously verified location p
            // after _Unwrapped
            index_ = p - span_->data();
        }
#endif

    protected:
        const Span* span_ = nullptr;
        std::ptrdiff_t index_ = 0;
    };

    template <std::ptrdiff_t Ext>
    class extent_type
    {
    public:
        using index_type = std::ptrdiff_t;

        static_assert(Ext >= 0, "A fixed-size span must be >= 0 in size.");

        constexpr extent_type() noexcept {}

        template <index_type Other>
        constexpr extent_type(extent_type<Other> ext)
        {
            static_assert(Other == Ext || Other == dynamic_extent,
                          "Mismatch between fixed-size extent and size of initializing data.");
            Expects(ext.size() == Ext);
        }

        constexpr extent_type(index_type size) { Expects(size == Ext); }

        constexpr index_type size() const noexcept { return Ext; }
    };

    template <>
    class extent_type<dynamic_extent>
    {
    public:
        using index_type = std::ptrdiff_t;

        template <index_type Other>
        explicit constexpr extent_type(extent_type<Other> ext) : size_(ext.size())
        {}

        explicit constexpr extent_type(index_type size) : size_(size) { Expects(size >= 0); }

        constexpr index_type size() const noexcept { return size_; }

    private:
        index_type size_;
    };

    template <class ElementType, std::ptrdiff_t Extent, std::ptrdiff_t Offset, std::ptrdiff_t Count>
    struct calculate_subspan_type
    {
        using type = span<ElementType, Count != dynamic_extent
                                           ? Count
                                           : (Extent != dynamic_extent ? Extent - Offset : Extent)>;
    };
} // namespace details

// [span], class template span
template <class ElementType, std::ptrdiff_t Extent>
class span
{
public:
    // constants and types
    using element_type = ElementType;
    using value_type = std::remove_cv_t<ElementType>;
    using index_type = std::ptrdiff_t;
    using pointer = element_type*;
    using reference = element_type&;

    using iterator = details::span_iterator<span<ElementType, Extent>, false>;
    using const_iterator = details::span_iterator<span<ElementType, Extent>, true>;
    using reverse_iterator = std::reverse_iterator<iterator>;
    using const_reverse_iterator = std::reverse_iterator<const_iterator>;

    using size_type = index_type;

#if defined(GSL_USE_STATIC_CONSTEXPR_WORKAROUND)
    static constexpr const index_type extent{Extent};
#else
    static constexpr index_type extent{Extent};
#endif

    // [span.cons], span constructors, copy, assignment, and destructor
    template <bool Dependent = false,
              // "Dependent" is needed to make "std::enable_if_t<Dependent || Extent <= 0>" SFINAE,
              // since "std::enable_if_t<Extent <= 0>" is ill-formed when Extent is greater than 0.
              class = std::enable_if_t<(Dependent || Extent <= 0)>>
    constexpr span() noexcept : storage_(nullptr, details::extent_type<0>())
    {}

    constexpr span(pointer ptr, index_type count) : storage_(ptr, count) {}

    constexpr span(pointer firstElem, pointer lastElem)
        : storage_(firstElem, std::distance(firstElem, lastElem))
    {}

    template <std::size_t N>
    constexpr span(element_type (&arr)[N]) noexcept
        : storage_(KnownNotNull{std::addressof(arr[0])}, details::extent_type<N>())
    {}

    template <std::size_t N, class = std::enable_if_t<(N > 0)>>
    constexpr span(std::array<std::remove_const_t<element_type>, N>& arr) noexcept
        : storage_(KnownNotNull{arr.data()}, details::extent_type<N>())
    {
    }

    constexpr span(std::array<std::remove_const_t<element_type>, 0>&) noexcept
        : storage_(static_cast<pointer>(nullptr), details::extent_type<0>())
    {
    }

    template <std::size_t N, class = std::enable_if_t<(N > 0)>>
    constexpr span(const std::array<std::remove_const_t<element_type>, N>& arr) noexcept
        : storage_(KnownNotNull{arr.data()}, details::extent_type<N>())
    {
    }

    constexpr span(const std::array<std::remove_const_t<element_type>, 0>&) noexcept
        : storage_(static_cast<pointer>(nullptr), details::extent_type<0>())
    {
    }

    // NB: the SFINAE here uses .data() as a incomplete/imperfect proxy for the requirement
    // on Container to be a contiguous sequence container.
    template <class Container,
              class = std::enable_if_t<
                  !details::is_span<Container>::value && !details::is_std_array<Container>::value &&
                  std::is_convertible<typename Container::pointer, pointer>::value &&
                  std::is_convertible<typename Container::pointer,
                                      decltype(std::declval<Container>().data())>::value>>
    constexpr span(Container& cont) : span(cont.data(), narrow<index_type>(cont.size()))
    {}

    template <class Container,
              class = std::enable_if_t<
                  std::is_const<element_type>::value && !details::is_span<Container>::value &&
                  std::is_convertible<typename Container::pointer, pointer>::value &&
                  std::is_convertible<typename Container::pointer,
                                      decltype(std::declval<Container>().data())>::value>>
    constexpr span(const Container& cont) : span(cont.data(), narrow<index_type>(cont.size()))
    {}

    constexpr span(const span& other) noexcept = default;

    template <
        class OtherElementType, std::ptrdiff_t OtherExtent,
        class = std::enable_if_t<
            details::is_allowed_extent_conversion<OtherExtent, Extent>::value &&
            details::is_allowed_element_type_conversion<OtherElementType, element_type>::value>>
    constexpr span(const span<OtherElementType, OtherExtent>& other)
        : storage_(other.data(), details::extent_type<OtherExtent>(other.size()))
    {}

    ~span() noexcept = default;
    constexpr span& operator=(const span& other) noexcept = default;

    // [span.sub], span subviews
    template <std::ptrdiff_t Count>
    constexpr span<element_type, Count> first() const
    {
        Expects(Count >= 0 && Count <= size());
        return {data(), Count};
    }

    template <std::ptrdiff_t Count>
    GSL_SUPPRESS(bounds.1) // NO-FORMAT: attribute
    constexpr span<element_type, Count> last() const
    {
        Expects(Count >= 0 && size() - Count >= 0);
        return {data() + (size() - Count), Count};
    }

    template <std::ptrdiff_t Offset, std::ptrdiff_t Count = dynamic_extent>
    GSL_SUPPRESS(bounds.1) // NO-FORMAT: attribute
    constexpr auto subspan() const ->
        typename details::calculate_subspan_type<ElementType, Extent, Offset, Count>::type
    {
        Expects((Offset >= 0 && size() - Offset >= 0) &&
                (Count == dynamic_extent || (Count >= 0 && Offset + Count <= size())));

        return {data() + Offset, Count == dynamic_extent ? size() - Offset : Count};
    }

    constexpr span<element_type, dynamic_extent> first(index_type count) const
    {
        Expects(count >= 0 && count <= size());
        return {data(), count};
    }

    constexpr span<element_type, dynamic_extent> last(index_type count) const
    {
        return make_subspan(size() - count, dynamic_extent, subspan_selector<Extent>{});
    }

    constexpr span<element_type, dynamic_extent> subspan(index_type offset,
                                                         index_type count = dynamic_extent) const
    {
        return make_subspan(offset, count, subspan_selector<Extent>{});
    }

    // [span.obs], span observers
    constexpr index_type size() const noexcept { return storage_.size(); }
    constexpr index_type size_bytes() const noexcept
    {
        return size() * narrow_cast<index_type>(sizeof(element_type));
    }
    constexpr bool empty() const noexcept { return size() == 0; }

    // [span.elem], span element access
    GSL_SUPPRESS(bounds.1) // NO-FORMAT: attribute
    constexpr reference operator[](index_type idx) const
    {
        Expects(CheckRange(idx, storage_.size()));
        return data()[idx];
    }

    constexpr reference at(index_type idx) const { return this->operator[](idx); }
    constexpr reference operator()(index_type idx) const { return this->operator[](idx); }
    constexpr pointer data() const noexcept { return storage_.data(); }

    // [span.iter], span iterator support
    constexpr iterator begin() const noexcept { return {this, 0}; }
    constexpr iterator end() const noexcept { return {this, size()}; }

    constexpr const_iterator cbegin() const noexcept { return {this, 0}; }
    constexpr const_iterator cend() const noexcept { return {this, size()}; }

    constexpr reverse_iterator rbegin() const noexcept { return reverse_iterator{end()}; }
    constexpr reverse_iterator rend() const noexcept { return reverse_iterator{begin()}; }

    constexpr const_reverse_iterator crbegin() const noexcept
    {
        return const_reverse_iterator{cend()};
    }
    constexpr const_reverse_iterator crend() const noexcept
    {
        return const_reverse_iterator{cbegin()};
    }

#ifdef _MSC_VER
    // Tell MSVC how to unwrap spans in range-based-for
    constexpr pointer _Unchecked_begin() const noexcept { return data(); }
    constexpr pointer _Unchecked_end() const noexcept
    {
        GSL_SUPPRESS(bounds.1) // NO-FORMAT: attribute
        return data() + size();
    }
#endif // _MSC_VER

private:
    static bool CheckRange(index_type idx, index_type size)
    {
        // Optimization:
        //
        // idx >= 0 && idx < size
        // =>
        // static_cast<size_t>(idx) < static_cast<size_t>(size)
        //
        // because size >=0 by span construction, and negative idx will
        // wrap around to a value always greater than size when casted.

        // check if we have enough space to wrap around
        if (sizeof(index_type) <= sizeof(size_t))
        {
            return narrow_cast<size_t>(idx) < narrow_cast<size_t>(size);
        }
        else
        {
            return idx >= 0 && idx < size;
        }
    }

    // Needed to remove unnecessary null check in subspans
    struct KnownNotNull
    {
        pointer p;
    };

    // this implementation detail class lets us take advantage of the
    // empty base class optimization to pay for only storage of a single
    // pointer in the case of fixed-size spans
    template <class ExtentType>
    class storage_type : public ExtentType
    {
    public:
        // KnownNotNull parameter is needed to remove unnecessary null check
        // in subspans and constructors from arrays
        template <class OtherExtentType>
        constexpr storage_type(KnownNotNull data, OtherExtentType ext)
            : ExtentType(ext), data_(data.p)
        {
            Expects(ExtentType::size() >= 0);
        }

        template <class OtherExtentType>
        constexpr storage_type(pointer data, OtherExtentType ext) : ExtentType(ext), data_(data)
        {
            Expects(ExtentType::size() >= 0);
            Expects(data || ExtentType::size() == 0);
        }

        constexpr pointer data() const noexcept { return data_; }

    private:
        pointer data_;
    };

    storage_type<details::extent_type<Extent>> storage_;

    // The rest is needed to remove unnecessary null check
    // in subspans and constructors from arrays
    constexpr span(KnownNotNull ptr, index_type count) : storage_(ptr, count) {}

    template <std::ptrdiff_t CallerExtent>
    class subspan_selector
    {
    };

    template <std::ptrdiff_t CallerExtent>
    span<element_type, dynamic_extent> make_subspan(index_type offset, index_type count,
                                                    subspan_selector<CallerExtent>) const
    {
        const span<element_type, dynamic_extent> tmp(*this);
        return tmp.subspan(offset, count);
    }

    GSL_SUPPRESS(bounds.1) // NO-FORMAT: attribute
    span<element_type, dynamic_extent> make_subspan(index_type offset, index_type count,
                                                    subspan_selector<dynamic_extent>) const
    {
        Expects(offset >= 0 && size() - offset >= 0);

        if (count == dynamic_extent) { return {KnownNotNull{data() + offset}, size() - offset}; }

        Expects(count >= 0 && size() - offset >= count);
        return {KnownNotNull{data() + offset}, count};
    }
};

#if defined(GSL_USE_STATIC_CONSTEXPR_WORKAROUND)
template <class ElementType, std::ptrdiff_t Extent>
constexpr const typename span<ElementType, Extent>::index_type span<ElementType, Extent>::extent;
#endif

// [span.comparison], span comparison operators
template <class ElementType, std::ptrdiff_t FirstExtent, std::ptrdiff_t SecondExtent>
constexpr bool operator==(span<ElementType, FirstExtent> l, span<ElementType, SecondExtent> r)
{
    return std::equal(l.begin(), l.end(), r.begin(), r.end());
}

template <class ElementType, std::ptrdiff_t Extent>
constexpr bool operator!=(span<ElementType, Extent> l, span<ElementType, Extent> r)
{
    return !(l == r);
}

template <class ElementType, std::ptrdiff_t Extent>
constexpr bool operator<(span<ElementType, Extent> l, span<ElementType, Extent> r)
{
    return std::lexicographical_compare(l.begin(), l.end(), r.begin(), r.end());
}

template <class ElementType, std::ptrdiff_t Extent>
constexpr bool operator<=(span<ElementType, Extent> l, span<ElementType, Extent> r)
{
    return !(l > r);
}

template <class ElementType, std::ptrdiff_t Extent>
constexpr bool operator>(span<ElementType, Extent> l, span<ElementType, Extent> r)
{
    return r < l;
}

template <class ElementType, std::ptrdiff_t Extent>
constexpr bool operator>=(span<ElementType, Extent> l, span<ElementType, Extent> r)
{
    return !(l < r);
}

namespace details
{
    // if we only supported compilers with good constexpr support then
    // this pair of classes could collapse down to a constexpr function

    // we should use a narrow_cast<> to go to std::size_t, but older compilers may not see it as
    // constexpr
    // and so will fail compilation of the template
    template <class ElementType, std::ptrdiff_t Extent>
    struct calculate_byte_size
        : std::integral_constant<std::ptrdiff_t,
                                 static_cast<std::ptrdiff_t>(sizeof(ElementType) *
                                                             static_cast<std::size_t>(Extent))>
    {
    };

    template <class ElementType>
    struct calculate_byte_size<ElementType, dynamic_extent>
        : std::integral_constant<std::ptrdiff_t, dynamic_extent>
    {
    };
} // namespace details

// [span.objectrep], views of object representation
template <class ElementType, std::ptrdiff_t Extent>
span<const byte, details::calculate_byte_size<ElementType, Extent>::value>
as_bytes(span<ElementType, Extent> s) noexcept
{
    GSL_SUPPRESS(type.1) // NO-FORMAT: attribute
    return {reinterpret_cast<const byte*>(s.data()), s.size_bytes()};
}

template <class ElementType, std::ptrdiff_t Extent,
          class = std::enable_if_t<!std::is_const<ElementType>::value>>
span<byte, details::calculate_byte_size<ElementType, Extent>::value>
as_writeable_bytes(span<ElementType, Extent> s) noexcept
{
    GSL_SUPPRESS(type.1) // NO-FORMAT: attribute
    return {reinterpret_cast<byte*>(s.data()), s.size_bytes()};
}

//
// make_span() - Utility functions for creating spans
//
template <class ElementType>
constexpr span<ElementType> make_span(ElementType* ptr,
                                      typename span<ElementType>::index_type count)
{
    return span<ElementType>(ptr, count);
}

template <class ElementType>
constexpr span<ElementType> make_span(ElementType* firstElem, ElementType* lastElem)
{
    return span<ElementType>(firstElem, lastElem);
}

template <class ElementType, std::size_t N>
constexpr span<ElementType, N> make_span(ElementType (&arr)[N]) noexcept
{
    return span<ElementType, N>(arr);
}

template <class Container>
constexpr span<typename Container::value_type> make_span(Container& cont)
{
    return span<typename Container::value_type>(cont);
}

template <class Container>
constexpr span<const typename Container::value_type> make_span(const Container& cont)
{
    return span<const typename Container::value_type>(cont);
}

template <class Ptr>
constexpr span<typename Ptr::element_type> make_span(Ptr& cont, std::ptrdiff_t count)
{
    return span<typename Ptr::element_type>(cont, count);
}

template <class Ptr>
constexpr span<typename Ptr::element_type> make_span(Ptr& cont)
{
    return span<typename Ptr::element_type>(cont);
}

// Specialization of gsl::at for span
template <class ElementType, std::ptrdiff_t Extent>
constexpr ElementType& at(span<ElementType, Extent> s, index i)
{
    // No bounds checking here because it is done in span::operator[] called below
    return s[i];
}

} // namespace gsl

#ifdef _MSC_VER
#if _MSC_VER < 1910
#undef constexpr
#pragma pop_macro("constexpr")

#endif // _MSC_VER < 1910

#pragma warning(pop)
#endif // _MSC_VER

#if __GNUC__ > 6
#pragma GCC diagnostic pop
#endif // __GNUC__ > 6

#endif // GSL_SPAN_H

max_stars_repo_name: Redchards/CVRP

max_issues_repo_name: Redchards/CVRP

max_forks_repo_name: Redchards/CVRP


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_core_test_Primitive_test_util_h_
#define SQ_INCLUDE_GUARD_core_test_Primitive_test_util_h_

#include "core/Primitive.h"
#include "core/typeutil.h"

#include <gsl/gsl>
#include <string_view>

namespace sq::test {

SQ_ND Primitive to_primitive(PrimitiveString &&v);
SQ_ND Primitive to_primitive(const PrimitiveString &v);
SQ_ND Primitive to_primitive(std::string_view v);
SQ_ND Primitive to_primitive(gsl::czstring<> v);
SQ_ND Primitive to_primitive(PrimitiveInt v);
SQ_ND Primitive to_primitive(int v);
SQ_ND Primitive to_primitive(PrimitiveFloat v);
SQ_ND Primitive to_primitive(PrimitiveBool v);
SQ_ND Primitive to_primitive(PrimitiveNull v);

} // namespace sq::test

#endif // SQ_INCLUDE_GUARD_core_test_Primitive_test_util_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_spline.h>
// #include "mex.h"

#ifndef verbose
#define verbose 0
#endif
/*
   $ ./fwhm > interp.dat
   $ graph -T ps < interp.dat > interp.ps
   */

struct dbg_data {
  double xm; // location of maxima
  double ym; // value at maxima
  double xleft;
  double xright;
};

FILE * logFile;

// my_f and my_f_params defines the interpolation function that the
// root-finding functions are run on

struct my_f_params {
  gsl_spline * spline;
  gsl_interp_accel * acc;
  double  offset; // shifts the function
};

double my_f(double , void * );
int findmin(double , double , gsl_spline *, gsl_interp_accel * , size_t , double * , double * );
static void createGaussian(double * , double * , size_t , double );
int fwhm(double * , double * , size_t , double * );

double my_f(double x, void * p)
{
  struct my_f_params * params = (struct my_f_params*) p;
#if verbose > 0
  printf("f(%f) + %f=  ", x, params->offset); fflush(stdout);
#endif

  double y =  -gsl_spline_eval(params->spline, x, params->acc) + params->offset;
#if verbose > 0
  printf("%f\n", y); fflush(stdout);
#endif

  return(y);
}

int findmin(double a, double b,
    gsl_spline *spline, gsl_interp_accel * acc,
    size_t N,
    double * xm, double * ym)
{

  // a, b: range of function
  int iter = 0, max_iter = 100;
  double m = 0; // expected location of minima

  int status;
  const gsl_min_fminimizer_type *T;
  gsl_min_fminimizer *s;

  struct my_f_params f_params;
  f_params.spline = spline;
  f_params.acc = acc;
  f_params.offset = 0;

  gsl_function F;
  F.function = &my_f;
  F.params = (void *) &f_params;

  T = gsl_min_fminimizer_brent;
  s = gsl_min_fminimizer_alloc (T);
  if(gsl_min_fminimizer_set (s, &F, m, a, b) == GSL_EINVAL)
  {
    // If no minima is enclosed ... i.e. the value at (a+b)/2 isn't lower than at a and b
    gsl_min_fminimizer_free(s);
    return 1;
  }

#if verbose > 0
  fflush(stdout);
  printf("Interval: [%f, %f]\n", a, b);
  double m_expected = -10;
  printf ("method: '%s'\n", gsl_min_fminimizer_name (s));
  printf ("%5s [%9s, %9s] %9s %10s %9s\n", "iter", "lower", "upper", "min", "err", "err(est)");
  printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n", iter, a, b, m, m - m_expected, b - a);
#endif

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

    a = gsl_min_fminimizer_x_lower (s);
    b = gsl_min_fminimizer_x_upper (s);

    status = gsl_min_test_interval (a, b, 0.001, 0.0);

#if verbose > 0
    float m = gsl_min_fminimizer_x_minimum (s);
    if (status == GSL_SUCCESS) {
      printf ("Converged:\n");
    } else { printf("No convergence %d\n", status); }

    printf ("%5d [%.7f, %.7f] " "%.7f %+.7f %.7f\n", iter, a, b, m, m - m_expected, b - a);
#endif

  } while (status == GSL_CONTINUE && iter < max_iter);

  xm[0] = (a+b)/2;
  ym[0] = my_f(xm[0], &f_params);

  gsl_min_fminimizer_free(s);

#if verbose > 0
  printf("Status: %d\n", status);
#endif

  return status;
}

int fwhm(double * x, double * y, size_t N, double * w)
{

  int useLog = 0;

  if(verbose)
  {
    useLog = 1;
  }

  if(useLog)
  {
    logFile = fopen("/tmp/fwhmlog", "w");
  }

#ifndef debug
gsl_set_error_handler_off();
#endif

  int N2 = (N-1)/2; // even number

  // Set up interpolation
  gsl_interp_accel *acc = gsl_interp_accel_alloc();  // Neccessary?
  gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, N);
  gsl_spline_init(spline, x, y, N);

  // See that the correct value is given at x=0
#if verbose > 0
  printf("y[%f] = %f\n", x[N2], gsl_spline_eval(spline, x[N2], acc));
  fflush(stdout);
#endif

  // Find position of max, i.e., centre in [xmin, xmax]
  double xmin = -2.5;
  double xmax = 2.501;

  double xm = 10e99;
  double ym = 10e99;
#if verbose > 0
  printf("calling findmin\n"); fflush(stdout);
#endif
  if( findmin(xmin, xmax, spline, acc, N, &xm, &ym) )
  {
#if verbose > 0
    printf("Could not find the position of the maxima!\n");
#endif
    gsl_spline_free (spline);
    gsl_interp_accel_free (acc);

    if(useLog)
    {
      fclose(logFile);
    }
    return 1;
  }

  if(useLog)
  {
    fprintf(logFile, "xm: %f, ym: %f\n", xm, ym);
  }

#if verbose > 0
  printf("-> Minima: f(%f) = %f\n", xm, ym);
#endif

  // Determine background
  double bg = (y[0]+y[N-1])/2;
#if verbose > 0
  printf("Background level: %f\n", bg);
#endif

  if(useLog)
  {
    fprintf(logFile, "bg: %f\n", bg);
  }

  if(bg <= ym)
  {
    fprintf(logFile, "stopping, bg <= tm\n");
    if(useLog)
    {
      fclose(logFile);
    }
    gsl_spline_free (spline);
    gsl_interp_accel_free (acc);

    return 1;
  }

  // Find intersection .5*(max-bg) for each side
  // ROOT FINDING

  const gsl_root_fsolver_type * Tsolve = gsl_root_fsolver_bisection;
  gsl_root_fsolver * rsolve = gsl_root_fsolver_alloc(Tsolve);

  // Search for left and right intersections

  double intersections[2];
  intersections[0] = 0; intersections[1] = 0;

  for(int dire = 0; dire<2; dire++)
  {
    // Search in range [x_lo, x_hi]
    double x_lo;
    double x_hi;

    if(dire==0)
    {
      x_lo = -N2;
      x_hi = xm;
    }

    if(dire==1)
    {
      x_lo = xm;
      x_hi = N2;
    }

    struct my_f_params f_params;
    f_params.spline = spline;
    f_params.acc = acc;
    f_params.offset = (bg+(-ym-bg)/2);

    if(useLog)
    {
      fprintf(logFile, "x_lo: %f x_hi: %f\n", x_lo, x_hi);
      fprintf(logFile, "offset: %f, bg: %f, ym: %f\n", f_params.offset, bg, ym);
    }

    gsl_function F;
    F.function = &my_f;
    F.params = (void *) &f_params;

    /*
     * Function: int gsl_root_fsolver_set (gsl_root_fsolver * s, gsl_function * f, double x_lower, double x_upper)
     * This function initializes, or reinitializes, an existing solver s to use the function f and the initial
     * search interval [x_lower, x_upper]
     */

    int status =  gsl_root_fsolver_set(rsolve, &F, x_lo, x_hi);

    size_t iter = 0;
    size_t max_iter = 10;

#if verbose > 0
    double r_expected = (x_lo+x_hi)/2;
    printf ("using %s method\n", gsl_root_fsolver_name (rsolve));
    printf ("%5s [%9s, %9s] %9s %10s %9s\n", "iter", "lower", "upper", "root", "err", "err(est)");
#endif

    do {
      iter++;
      status = gsl_root_fsolver_iterate (rsolve);
      // update bounds
      x_lo = gsl_root_fsolver_x_lower (rsolve);
      x_hi = gsl_root_fsolver_x_upper (rsolve);
      // see if converged
      status = gsl_root_test_interval (x_lo, x_hi, 0, 0.001);
#if verbose > 0
      double r = gsl_root_fsolver_root (rsolve);
      if (status == GSL_SUCCESS)
        printf ("Converged:\n");
      printf ("%5lu [%.7f, %.7f] %.7f %+.7f %.7f\n", iter, x_lo, x_hi,
          r, r - r_expected,
          x_hi - x_lo);
#endif
    }
    while (status == GSL_CONTINUE && iter < max_iter);

#if verbose > 0
  printf("Found intersection at %f\n", (x_lo+x_hi)/2);
#endif
  intersections[dire] = (x_lo+x_hi)/2;

  if(useLog)
  {
    fprintf(logFile, "intersection %d: %f\n", dire, (x_lo+x_hi)/2);
  }
  }

  gsl_root_fsolver_free(rsolve);

  // Clean up
  gsl_spline_free (spline);
  gsl_interp_accel_free (acc);

  // fwhm: right-left positions
  w[0] = intersections[1] - intersections[0];

  if(useLog)
  {
    fclose(logFile);
  }

  return 0; // Ok
}

static void createGaussian(
    double * x,
    double * y,
    size_t N,
    double x0)
{
  /* Create a gaussian shaped test signal, y
   * over the domain x
   * x0 is the offset from center
   * */
  int N2 = (N-1)/2; // middle element

  for(int kk=0; kk<N; kk++)
  {
    double p = 2*((double) kk - N2)/N; // [-1,1];
    p = (p*3-x0);
    y[kk] =  exp(-p*p);
    x[kk] = kk-N2; // Set x to pixel position from centre
#if verbose>0
    printf("%f %f\n", x[kk], y[kk]);
#endif
  }
}

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

  int N = 11; // Size of test signal
  double * y = malloc(N*sizeof(double));
  double * x = malloc(N*sizeof(double));
  memset(x, 0, N*sizeof(double));
  memset(y, 0, N*sizeof(double));
  createGaussian(x,y, N, 0);

  double w = -1; // output

  if(1){
    for(double delta = -1; delta <=1; delta +=0.1)
    {

      createGaussian(x,y, N, delta);
      if(fwhm(x, y, N, &w))
      {
        printf("Error: Could not calculate fwhm!\n");
      }

      printf("offset: % .2f pixels, fwhm: %f pixels\n", delta, w);
    }
  }

  for(int kk = 0; kk<N; kk++)
  {
    y[kk] = 0;
  }

  if(1){
  printf("constant y\n");
  if(fwhm(x, y, N, &w))
  {
    printf("# Error!\n");
    free(y); free(x);
    return(1);
  }
  }

  free(y);
  free(x);
  return(0);
}

max_stars_repo_name: elgw/dotter

max_issues_repo_name: elgw/dotter

max_forks_repo_name: elgw/dotter


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>

#include "ccl.h"

/* BCM correction */
// See Schneider & Teyssier (2015) for details of the model.
double ccl_bcm_model_fka(ccl_cosmology * cosmo, double k, double a, int *status) {
  double fkz;
  double b0;
  double bfunc, bfunc4;
  double kg;
  double gf,scomp;
  double kh;
  double z;

  z = 1./a - 1.;
  kh = k / cosmo->params.h;
  b0 = 0.105*cosmo->params.bcm_log10Mc - 1.27;
  bfunc = b0 / (1. + pow(z/2.3, 2.5));
  bfunc4 = (1-bfunc) * (1-bfunc) * (1-bfunc) * (1-bfunc);
  kg = 0.7 * bfunc4 * pow(cosmo->params.bcm_etab, -1.6);
  gf = bfunc / (1 + pow(kh/kg, 3.)) + 1. - bfunc; //k in h/Mpc
  scomp = 1 + (kh / cosmo->params.bcm_ks) * (kh / cosmo->params.bcm_ks); //k in h/Mpc
  fkz = gf * scomp;
  return fkz;
}

void ccl_bcm_correct(ccl_cosmology *cosmo, ccl_f2d_t *psp, int *status)
{
  size_t nk, na;
  double *x, *z, *y2d=NULL;

  //Find lk array
  if(psp->fk != NULL) {
    nk = psp->fk->size;
    x = psp->fk->x;
  }
  else {
    nk = psp->fka->interp_object.xsize;
    x = psp->fka->xarr;
  }

  //Find a array
  if(psp->fa != NULL) {
    na = psp->fa->size;
    z = psp->fa->x;
  }
  else {
    na = psp->fka->interp_object.ysize;
    z = psp->fka->yarr;
  }

  //Allocate pka array
  y2d = malloc(nk * na * sizeof(double));
  if (y2d == NULL) {
    *status = CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo,
                                     "ccl_bcm.c: ccl_bcm_correct(): "
                                     "memory allocation\n");
  }

  if (*status == 0) {
    for (int j = 0; j<na; j++) {
      for (int i=0; i<nk; i++) {
        if (*status == 0) {
          double pk = ccl_f2d_t_eval(psp, x[i], z[j], cosmo, status);
          double fbcm = ccl_bcm_model_fka(cosmo, exp(x[i]), z[j], status);
          if(psp->is_log)
            y2d[j*nk + i] = log(pk*fbcm);
          else
            y2d[j*nk + i] = pk*fbcm;
        }
      }
    }
  }

  if (*status == 0) {
    gsl_spline2d *fka = gsl_spline2d_alloc(gsl_interp2d_bicubic, nk, na);

    if (fka == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_bcm.c: ccl_bcm_correct(): "
                                       "memory allocation\n");
    }
    if(*status == 0) {
      int spstatus = gsl_spline2d_init(fka, x, z, y2d, nk, na);
      if(spstatus) {
        *status = CCL_ERROR_MEMORY;
        ccl_cosmology_set_status_message(cosmo,
                                         "ccl_bcm.c: ccl_bcm_correct(): "
                                         "Error initializing spline\n");
      }
    }
    if(*status == 0) {
      if(psp->fa != NULL)
        gsl_spline_free(psp->fa);
      if(psp->fk != NULL)
        gsl_spline_free(psp->fk);
      if(psp->fka != NULL)
        gsl_spline2d_free(psp->fka);
      psp->fka = fka;
      psp->is_factorizable = 0;
      psp->is_k_constant = 0;
      psp->is_a_constant = 0;
    }
    else
      gsl_spline2d_free(fka);
  }

  free(y2d);
}

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL


text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

	1. Redistributions of source code must retain the above copyright notice,
		this list of conditions and the following disclaimer.

	2. Redistributions in binary form must reproduce the above copyright notice,
		this list of conditions and the following disclaimer in the
		documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file rk_4_2.c
 * \brief Source file to optimize Runge-Kutta 4 steps 2nd order methods.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <string.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#include "config.h"
#include "utils.h"
#include "optimize.h"
#include "rk.h"
#include "rk_4_2.h"

#define DEBUG_RK_4_2 0          ///< macro to debug.

/**
 * Function to obtain the coefficients of a 4 steps 2nd order Runge-Kutta
 * method.
 */
int
rk_tb_4_2 (Optimize * optimize) ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_4_2
  fprintf (stderr, "rk_tb_4_2: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t4 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b31 (tb) = r[4];
  b32 (tb) = r[5];
  b41 (tb) = r[6];
  b42 (tb) = r[7];
  b43 (tb) = (0.5L - b41 (tb) * t1 (tb) - b42 (tb) * t2 (tb)) / t3 (tb);
  if (isnan (b43 (tb)))
    return 0;
  rk_b_4 (tb);
#if DEBUG_RK_4_2
  fprintf (stderr, "rk_tb_4_2: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 4 steps 2nd order, 3rd order in
 * equations depending only in time, Runge-Kutta method.
 */
int
rk_tb_4_2t (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_4_2
  fprintf (stderr, "rk_tb_4_2t: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t4 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b31 (tb) = r[4];
  b32 (tb) = r[5];
  b41 (tb) = r[6];
  b42 (tb) = (1.L / 3.L - 0.5L * t3 (tb)
              - b41 (tb) * t1 (tb) * (t1 (tb) - t3 (tb)))
    / (t2 (tb) * (t2 (tb) - t3 (tb)));
  if (isnan (b42 (tb)))
    return 0;
  b43 (tb) = (0.5L - b41 (tb) * t1 (tb) - b42 (tb) * t2 (tb)) / t3 (tb);
  if (isnan (b43 (tb)))
    return 0;
  rk_b_4 (tb);
#if DEBUG_RK_4_2
  fprintf (stderr, "rk_tb_4_2t: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 4 steps 1st-2nd order Runge-Kutta
 * pair.
 */
int
rk_tb_4_2p (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_4_2
  fprintf (stderr, "rk_tb_4_2p: start\n");
#endif
  if (!rk_tb_4_2 (optimize))
    return 0;
  tb = optimize->coefficient;
  e41 (tb) = e42 (tb) = 0.L;
  rk_e_4 (tb);
#if DEBUG_RK_4_2
  rk_print_e (optimize, "rk_tb_4_2p", stderr);
  fprintf (stderr, "rk_tb_4_2p: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 4 steps 1st-2nd order, 1st-3rd order
 * in equations depending only in time, Runge-Kutta method.
 */
int
rk_tb_4_2tp (Optimize * optimize)       ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_4_2
  fprintf (stderr, "rk_tb_4_2tp: start\n");
#endif
  if (!rk_tb_4_2t (optimize))
    return 0;
  tb = optimize->coefficient;
  e41 (tb) = e42 (tb) = 0.L;
  rk_e_4 (tb);
#if DEBUG_RK_4_2
  rk_print_tb (optimize, "rk_tb_4_2tp", stderr);
  fprintf (stderr, "rk_tb_4_2pt: end\n");
#endif
  return 1;
}

/**
 * Function to calculate the objective function of a 4 steps 2nd order
 * Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_4_2 (RK * rk)   ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_4_2
  fprintf (stderr, "rk_objective_tb_4_2: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b43 (tb) < 0.L)
    o += b43 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), t3 (tb))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_4_2
  fprintf (stderr, "rk_objective_tb_4_2: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_4_2: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 4 steps 2nd order, 3rd
 * order in equations depending only in time, Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_4_2t (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_4_2
  fprintf (stderr, "rk_objective_tb_4_2t: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b42 (tb) < 0.L)
    o += b42 (tb);
  if (b43 (tb) < 0.L)
    o += b43 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), t3 (tb))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_4_2
  fprintf (stderr, "rk_objective_tb_4_2t: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_4_2t: end\n");
#endif
  return o;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


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

#include <fastpm/libfastpm.h>
#include <fastpm/logging.h>

#define LENGTH_FERMI_DIRAC_TABLE 4000 //sets the length of the table on which CDF
                                      //will be evaluated
#define MAX_FERMI_DIRAC          20.0 //Upper limit of F-D distribution in units of
                                      // p/T

#define NEXT(n, i)  (((n) + (i)/(n)) >> 1) // Needed for Healpix routine

//needed for healpix
unsigned int isqrt(int number) {
    unsigned int n  = 1;
    unsigned int n1 = NEXT(n, number);

    while(abs(n1 - n) > 1) {
        n  = n1;
        n1 = NEXT(n, number);
    }
    while(n1*n1 > number)
        n1--;
    return n1;
}

//Converts from pixel number to unit vector (needed for healpix)
void pix2vec (int pix, double *vec, int n_side)
{
    double z, phi;
    int nside_ = n_side;
    long ncap_=nside_*(nside_-1)*2;
    long npix_=12*nside_*nside_;
    double fact2_  = 4./npix_;
    if (pix<ncap_) /* North Polar cap */
    {
        int iring = (int)(0.5*(1+isqrt(1+2*pix))); /* counted from North pole */
        int iphi  = (pix+1) - 2*iring*(iring-1);

        z = 1.0 - (iring*iring)*fact2_;
        phi = (iphi-0.5) * 0.5*M_PI/iring;
    }
    else if (pix<(npix_-ncap_)) /* Equatorial region */
    {
        double fact1_  = (nside_<<1)*fact2_;
        int ip  = pix - ncap_;
        int iring = ip/(4*nside_) + nside_; /* counted from North pole */
        int iphi  = ip%(4*nside_) + 1;
        /* 1 if iring+nside is odd, 1/2 otherwise */
        double fodd = ((iring+nside_)&1) ? 1 : 0.5;

        int nl2 = 2*nside_;
        z = (nl2-iring)*fact1_;
        phi = (iphi-fodd) * M_PI/nl2;
    }
    else /* South Polar cap */
    {
        int ip = npix_ - pix;
        int iring = (int)(0.5*(1+isqrt(2*ip-1))); /* counted from South pole */
        int iphi  = 4*iring + 1 - (ip - 2*iring*(iring-1));

        z = -1.0 + (iring*iring)*fact2_;
        phi = (iphi-0.5) * 0.5*M_PI/iring;
    }

    double v[3];
    v[0] = sin(acos(z))*cos(phi);
    v[1] = sin(acos(z))*sin(phi);
    v[2] = z;

    /* rotate the vector to break degeneracy with grid axes */
    vec[0] = 0.5*v[0] -0.5*v[1] +0.70710678*v[2];
    vec[1] = 0.85355339*v[0] +0.14644661*v[1] -0.5*v[2];
    vec[2] = 0.14644661*v[0] +0.85355339*v[1] +0.5*v[2];
}

/* make params struct for functions we will integrate with gsl */
typedef struct kernel_params {
    double m[3];
    int n;
} kernel_params;

double fermi_dirac_kernel_vol_1(double x)
{
    /* Functional form of FD distribution for 1 particle
       over the full volume (i.e. no *x^2) */
    return 1. / (exp(x) + 1);
}

double fermi_dirac_kernel_vol(double x, void * p)
{
    /* FD kernel in the case of multiple particles
       over the full volume (i.e. no *x^2) */
    struct kernel_params * params = (struct kernel_params *)p;

    double result = fermi_dirac_kernel_vol_1(x);
    // only loop from i=1. i=0 term has been included in result already (r=1)
    for(int i = 1; i < params->n; i ++) {
        double r = params->m[i] / params->m[0];
        result += (r*r*r*r) * fermi_dirac_kernel_vol_1(x * r);
    }
    return result;
}

/* now define the 1D kernel and moments */
double fermi_dirac_kernel(double x, void * p)
{
    return x*x * fermi_dirac_kernel_vol(x, p);
}

/* Samples the low velocity end more effectively */
double low_vel_kernel(double x, void * p)
{
    return x * fermi_dirac_kernel_vol(x, p);
}

/* Needed for calculatin the velocity dispersion */
double fermi_dirac_dispersion(double x, void * p)
{
     return x*x*x*x * fermi_dirac_kernel_vol(x, p);
}

void
divide_fd(double *vel_table, double *mass,
          kernel_params * params,
          int n_shells,
          int lvk)
{
    double fermi_dirac_vel_ncdm[LENGTH_FERMI_DIRAC_TABLE];
    double fermi_dirac_cdf_ncdm[LENGTH_FERMI_DIRAC_TABLE]; //stores CDF

    int i,j;
    double v_bin,u,bin_average,bin_mass;
    double vel_edges[n_shells];
    double result,error;
    gsl_integration_workspace * w
        = gsl_integration_workspace_alloc (1000);

    gsl_function F,G,H;
    if(lvk)
        F.function = &low_vel_kernel;
    else
        F.function = &fermi_dirac_kernel;
    G.function = &fermi_dirac_kernel;
    H.function = &fermi_dirac_dispersion;

    F.params = (void*) params;
    G.params = (void*) params;
    H.params = (void*) params;

    // Initialize the CDF table
    for(i = 0; i < LENGTH_FERMI_DIRAC_TABLE; i++){
        fermi_dirac_vel_ncdm[i] = MAX_FERMI_DIRAC * i / (LENGTH_FERMI_DIRAC_TABLE - 1.0);
        gsl_integration_qags (&F, 0, fermi_dirac_vel_ncdm[i], 0,1e-7 , 1000,
			  w, &result, &error);
        fermi_dirac_cdf_ncdm[i] = result;
    }

    //Normalize to 1
    for(i = 0; i < LENGTH_FERMI_DIRAC_TABLE; i++)
        fermi_dirac_cdf_ncdm[i] /= fermi_dirac_cdf_ncdm[LENGTH_FERMI_DIRAC_TABLE - 1];

    //Define the edges of the velocity bins (lower edge held to 0)
    for(i=0;i<n_shells;i++){
        v_bin = (i+1)/(n_shells*1.0);
        j=0;
        while(j < LENGTH_FERMI_DIRAC_TABLE - 2)
            if(v_bin > fermi_dirac_cdf_ncdm[j + 1])
                j++;
            else
                break;

        u = (v_bin - fermi_dirac_cdf_ncdm[j]) / (fermi_dirac_cdf_ncdm[j + 1] - fermi_dirac_cdf_ncdm[j]);

        vel_edges[i] = fermi_dirac_vel_ncdm[j] * (1 - u) + fermi_dirac_vel_ncdm[j + 1] * u;
    }

    //Get the bin velocities and bin masses
    double total_mass;
    gsl_integration_qags (&G, 0, fermi_dirac_vel_ncdm[LENGTH_FERMI_DIRAC_TABLE-1], 0, 1e-7, 1000,
			  w, &result, &error);
    total_mass = result;
    for(i=0;i<n_shells;i++){
        //Special case for lowest bin - left edge has to be zero
        if (i==0){
            gsl_integration_qags (&H, 0,vel_edges[i], 0, 1e-7, 1000,
			  w, &result, &error);
            bin_average = result;
            gsl_integration_qags (&G, 0, vel_edges[i], 0, 1e-7, 1000,
			  w, &result, &error);
            bin_average /= result;

            bin_mass = result/total_mass;
        }
        else{
            gsl_integration_qags (&H, vel_edges[i-1], vel_edges[i], 0, 1e-7, 1000,
			  w, &result, &error);
            bin_average = result;
            gsl_integration_qags (&G, vel_edges[i-1], vel_edges[i], 0, 1e-7, 1000,
			  w, &result, &error);
            bin_average /= result;

            bin_mass = result/total_mass;
        }
        vel_table[i] = sqrt(bin_average);
        mass[i] = bin_mass;
    }

}

void divide_sphere_healpix(double *vec_table, int n_side)  //so is the direction totally fixed? could we rotate the divisions relative to north pole (i.e. relative to the box)?
{
    int i,k;
    double v_sq[3];
    for(k=0;k<3;k++)
        v_sq[k] = 0.;

    // Get unit vectors for all pixels
    for(i=0;i<12*n_side*n_side;i++){
        pix2vec(i,&vec_table[i*3],n_side);
        for(k=0;k<3;k++)
            v_sq[k] += vec_table[i*3+k]*vec_table[i*3+k];
    }
    // Isotropize the velocity dispersion - important for low n_side    [should we remove this if for high n_side???]
    for(k=0;k<3;k++){
        v_sq[k] /= 12*n_side*n_side;
        v_sq[k] /= 1./3.; // Set each direction to have 1/3 of total dispersion
    }
    for(i=0;i<12*n_side*n_side;i++)
        for(k=0;k<3;k++)
            vec_table[i*3+k] /= sqrt(v_sq[k]);
}

void divide_sphere_fibonacci(double *vec_table, int n_side)
{
    int n_fibonacci = n_side;
    double lat, lon;
    int i, j;
    //int N_tot = 2 * n_fibonacci + 1;

    for (i=-n_fibonacci;i<n_fibonacci+1;i++){
        lat = asin(2.0*i/(2.0*n_fibonacci+1));
        lon = 2.0*M_PI*i*2.0/(1.0+sqrt(5.0));
        j = i + n_fibonacci;
        vec_table[j*3+0] = cos(lat)*cos(lon);
        vec_table[j*3+1] = cos(lat)*sin(lon);
        vec_table[j*3+2] = sin(lat);
    }
}

void divide_sphere(FastPMncdmSphereScheme ncdm_sphere_scheme,
                   double *vec_table, int n_side)
{
    switch (ncdm_sphere_scheme){
        case FASTPM_NCDM_SPHERE_HEALPIX:
            divide_sphere_healpix(vec_table, n_side);
        break;
        case FASTPM_NCDM_SPHERE_FIBONACCI:
            divide_sphere_fibonacci(vec_table, n_side);
        break;
        default:
            fastpm_raise(-1, "Wrong ncdm sphere scheme.\n");
    }
}

static void _fastpm_ncdm_init_fill(FastPMncdmInitData* nid);

//create and fill the table
FastPMncdmInitData*
fastpm_ncdm_init_create(
    double BoxSize,
    FastPMCosmology * c,
    double z,
    int n_shells,
    int n_side,
    int lvk,
    FastPMncdmSphereScheme ncdm_sphere_scheme)
{
    FastPMncdmInitData* nid = malloc(sizeof(nid[0]));

    nid->BoxSize = BoxSize;
    nid->m_ncdm_sum = 0;
    for(int i = 0; i < c->N_ncdm; i ++) {
        nid->m_ncdm[i] = c->m_ncdm[i];
        nid->m_ncdm_sum += c->m_ncdm[i];
    }
    nid->n_ncdm = c->N_ncdm;
    nid->Omega_ncdm = c->Omega_ncdm;
    nid->z = z;
    fastpm_info("ncdm reference redshift = %g\n", z);
    nid->n_shells = n_shells;
    nid->n_side = n_side;
    nid->lvk = lvk;
    nid->ncdm_sphere_scheme = ncdm_sphere_scheme;

    /* this is the total number of velocity vectors produced */
    switch (ncdm_sphere_scheme){
        case FASTPM_NCDM_SPHERE_HEALPIX:
            nid->n_sphere = 12 * n_side * n_side;
        break;
        case FASTPM_NCDM_SPHERE_FIBONACCI:
            nid->n_sphere = 2 * n_side + 1;
        break;
        default:
            fastpm_raise(-1, "Wrong ncdm sphere scheme.\n");
    }

    nid->n_split = n_shells * nid->n_sphere;

    /* recall vel is a pointer to a 3 element array
    so lets make into multidim array of dim (n_split x 3) */
    //for now store all (3) neutrinos in the same nid table.
    nid->vel = calloc(nid->n_split, sizeof(double[3]));
    nid->mass = calloc(nid->n_split, sizeof(double));   //a 1d array to store all the masses

    _fastpm_ncdm_init_fill(nid);

    return nid;
}

//free the table
void
fastpm_ncdm_init_free(FastPMncdmInitData* nid)
{
    free(nid->vel);
    free(nid->mass);
    free(nid);
}

//append the table
static void
_fastpm_ncdm_init_fill(FastPMncdmInitData* nid)
{
    int n_ncdm = nid->n_ncdm;
    int n_shells = nid->n_shells;
    int n_side = nid->n_side;
    int n_sphere = nid->n_sphere;
    int lvk = nid->lvk;

    double *vel_table, *vec_table, *masstab;
    vel_table = malloc(sizeof(double)*n_shells);
    /* masstab is the distribution integrated per shell, sums to 1 */
    masstab = calloc(n_shells,sizeof(double));
    vec_table = malloc(sizeof(double)*n_sphere*3);

    /* define the kernel params for dividing fd */
    kernel_params* params = malloc(sizeof(params[0]));
    for(int i = 0; i < n_ncdm; i ++) {
        params->m[i] = nid->m_ncdm[i];
    }
    params->n = n_ncdm;

    divide_fd(vel_table, masstab, params, n_shells, lvk);
    divide_sphere(nid->ncdm_sphere_scheme, vec_table, n_side);

    /* velocity_conversion_factor converts the dimless exponent
       in the FD distribution to velocity in FastPM units:
       conjugate momentum (a^2 xdot, where x is comoving dist).
       kTc = 50.3 eV/c^2 km/s */
    double m0 = nid->m_ncdm[0];
    double velocity_conversion_factor = 50.3 / m0 / HubbleConstant;

    int i, j, d;
    int s = 0;                          //row num  (s for split index)
    for(i = 0; i < n_sphere; i ++){
        for(j = 0; j < n_shells; j ++){
            /* define mass st sum over split gives 1 */
            nid->mass[s] = masstab[j] / n_sphere;
            for(d = 0; d < 3; d ++){
                nid->vel[s][d] = vel_table[j] * vec_table[i*3 + d] * velocity_conversion_factor;
            }
            s++;
        }
    }
    free(vel_table);
    free(masstab);
    free(vec_table);
}

void
fastpm_split_ncdm(FastPMncdmInitData * nid,
                  FastPMStore * src,
                  FastPMStore * dest,
                  MPI_Comm comm)
{
    /*
    Takes store src and splits according to the ncdm
    init data, and uses this to populate an dest store.
    src is fully populated after calling the setup_lpt
    func, but dest has only been init'd.
    */

    dest->np = src->np * nid->n_split;    // FIXME: ? Is proc imbalance an issue that needs to be accounted for?
    if(dest->np > dest->np_upper) {
        fastpm_raise(-1, "exceeding limit on the number limit of particles for the ncdm species \n");
    }

    size_t np_total = fastpm_store_get_np_total(src, comm);     //np_total is num of UNsplit ncdms

    /* avg mass of UNsplit ncdm site */
    double M0 = nid->Omega_ncdm * FASTPM_CRITICAL_DENSITY * pow(nid->BoxSize, 3) / np_total;
    double Mestate = M0 / nid->n_ncdm;
    fastpm_info("average mass of a ncdm particle is %g\n", Mestate);

    /* copy and amend meta-data */
    memmove(&dest->meta, &src->meta, sizeof(src->meta));

    dest->meta.M0 = 0.; /* will have a mass per particle */

    /* Compute the factor used to displace ncdm so that the expanded spheres almost touch.
       This assumes the largest speed particle is in the final entry of vel table */
    int n_shells = nid->n_shells;
    int n_sphere = nid->n_sphere;
    double vthm_max = 0;
    for (int d = 0; d < 3; d ++)
        vthm_max += nid->vel[n_sphere*n_shells-1][d]*nid->vel[n_sphere*n_shells-1][d];
    vthm_max = sqrt(vthm_max);
    double disp_factor = 0.5 * nid->BoxSize / nid->n_ncdm / vthm_max * (n_shells - 1) / n_shells;

    ptrdiff_t i, s, d;
    ptrdiff_t r = 0;    //row index of dest
    int c;
    double vthm;
    for(i = 0; i < src->np; i ++) {    //loop thru each src particle
        for(s = 0; s < nid->n_split; s ++) {  //loop thru split velocities
            /* copy all cols el by el */
            for(c = 0; c < 32; c ++) {
                if (!dest->columns[c] || !src->columns[c]) continue;

                size_t elsize = dest->_column_info[c].elsize;
                memcpy(dest->columns[c] + r * elsize,
                       src->columns[c] + i * elsize,
                       elsize);
            }

            /* give id, mass and add thm vel */
            dest->id[r] = s * src->meta._q_size + src->id[i];
            dest->mass[r] = nid->mass[s] * M0;    // ensures sum over split gives M0

            for(d = 0; d < 3; d ++){
                vthm = nid->vel[s][d];
                dest->v[r][d] = vthm;
                dest->x[r][d] += vthm * disp_factor;
                if (dest->q) dest->q[r][d] += vthm * disp_factor;
            }
            r ++;
        }
    }
}

max_stars_repo_name: sbird/FastPMRunner

max_issues_repo_name: sbird/FastPMRunner

max_forks_repo_name: sbird/FastPMRunner


text:
#ifndef ALE_INCLUDE_ALE_H
#define ALE_INCLUDE_ALE_H

#include <nlohmann/json.hpp>
#include <string>
#include <vector>

#include <gsl/gsl_interp.h>

#include <nlohmann/json.hpp>
using json = nlohmann::json;

namespace ale {

  /// Interpolation enum for defining different methods of interpolation
  enum interpolation {
    /// Interpolate using linear interpolation
    linear,
    /// Interpolate using spline interpolation
    spline
  };

  /**
   *@brief Get the position of the spacecraft at a given time based on a set of coordinates, and their associated times
   *@param coords A vector of double vectors of coordinates
   *@param times A double vector of times
   *@param time Time to observe the spacecraft's position at
   *@param interp Interpolation type
   *@return A vector double of the spacecrafts position
   */
  std::vector<double> getPosition(std::vector<std::vector<double>> coords,
                                  std::vector<double> times,
                                  double time, interpolation interp);
  /**
   *@brief Get the velocity of the spacecraft at a given time based on a set of coordinates, and their associated times
   *@param coords A vector of double vectors of coordinates
   *@param times A double vector of times
   *@param time Time to observe the spacecraft's velocity at
   *@param interp Interpolation type
   *@return A vector double of the spacecrafts velocity
   */
  std::vector<double> getVelocity(std::vector<std::vector<double>> coords,
                                  std::vector<double> times,
                                  double time, const interpolation interp);
  /**
   *@brief Get the position of the spacecraft at a given time based on a derived function from a set of coeffcients
   *@param coeffs A vector of double vector of coeffcients
   *@param time Time to observe the spacecraft's position at
   *@return A vector double of the spacecrafts position
   */
  std::vector<double> getPosition(std::vector<std::vector<double>> coeffs, double time);

  /**
   *@brief Get the velocity of the spacecraft at a given time based on a derived function from a set of coeffcients
   *@param coeffs A vector of double vector of coeffcients
   *@param time Time to observe the spacecraft's velocity at
   *@return A vector double of the spacecrafts velocity
   */
  std::vector<double> getVelocity(std::vector<std::vector<double>> coeffs, double time);

  /**
   *@brief Get the rotation of the spacecraft at a given time based on a set of rotations, and their associated times
   *@param rotations A vector of double vector of rotations
   *@param times A double vector of times
   *@param time Time to observe the spacecraft's rotation at
   *@param interp Interpolation type
   *@return A vector double of the spacecrafts rotation
   */
  std::vector<double> getRotation(std::vector<std::vector<double>> rotations,
                                  std::vector<double> times,
                                  double time, interpolation interp);

  /**
   *@brief Get the angular velocity of the spacecraft at a given time based on a set of rotations, and their associated times
   *@param rotations A vector of double vector of rotations
   *@param times A double vector of times
   *@param time Time to observe the spacecraft's angular velocity at
   *@param interp Interpolation type
   *@return A vector double of the spacecrafts angular velocity
   */
  std::vector<double> getAngularVelocity(std::vector<std::vector<double>> rotations,
                                         std::vector<double> times,
                                         double time, interpolation interp);

   /**
    *@brief Get the rotation of the spacecraft at a given time based on a derived function from a set of coeffcients
    *@param coeffs A vector of double vector of coeffcients
    *@param time Time to observe the spacecraft's rotation at
    *@return A vector double of the spacecrafts rotation
    */
  std::vector<double> getRotation(std::vector<std::vector<double>> coeffs, double time);

  /**
   *@brief Get the angular velocity of the spacecraft at a given time based on a derived function from a set of coeffcients
   *@param coeffs A vector of double vector of coeffcients
   *@param time Time to observe the spacecraft's angular velocity at
   *@return A vector double of the spacecrafts angular velocity
   */
  std::vector<double> getAngularVelocity(std::vector<std::vector<double>> coeffs, double time);

  /**
   *@brief Generates a derivatives in respect to time from a polynomial constructed using the given coeffcients, time, and derivation number
   *@param coeffs A double vector of coefficients can be any number of coefficients
   *@param time Time to use when deriving
   *@param d The order of the derivative to generate (Currently supports 0, 1, and 2)
   *@return Evalutation of the given polynomial as a double
   */
  double evaluatePolynomial(std::vector<double> coeffs, double time, int d);

  /**
   *@brief Interpolates the spacecrafts position along a path generated from a set of points,
                times, and a time of observation
   *@param points A double vector of points
   *@param times A double vector of times
   *@param time A double to use as the time of observation
   *@param interp An interpolation enum dictating what type of interpolation to use
   *@param d The order of the derivative to generate when interpolating
                      (Currently supports 0, 1, and 2)
   *@return
   */
  double interpolate(std::vector<double> points, std::vector<double> times, double time, interpolation interp, int d);
  std::string loads(std::string filename, std::string props="", std::string formatter="usgscsm", bool verbose=true);

  json load(std::string filename, std::string props="", std::string formatter="usgscsm", bool verbose=true);

}

#endif // ALE_H

max_stars_repo_name: kaitlyndlee/ale

max_issues_repo_name: kaitlyndlee/ale

max_forks_repo_name: kaitlyndlee/ale


text:
#include <gsl/gsl_rng.h>

double GrowthFactor(double astart, double aend);
void   print_spec(void);
int    FatalError(int errnum);
void   displacement_fields(void);
void   initialize_ffts(void);
void   set_units(void);
void   assemble_particles(void);
void   free_ffts(void);
double fnl(double x);

int find_files(char *fname);

void   assemble_grid(void);
void   read_power_table(void);
double periodic_wrap(double x);

double PowerSpec(double kmag);
double PowerSpec_Efstathiou(double k);
double PowerSpec_EH(double k);
double PowerSpec_Tabulated(double k);
double PowerSpec_DM_2ndSpecies(double k);
double PowerSpec_Tabulated2nd(double k);

void   initialize_powerspectrum(void);
double GrowthFactor(double astart, double aend);
double growth(double a);
double growth_int(double, void *param);
double qromb(double (*func)(double), double a, double b);
double sigma2_int(double k, void *param);
double TopHatSigma2(double R);
double F_Omega(double a);

void   combine_particle_data(void);
int    compare_logk(const void *a, const void *b);

void  write_particle_data(void);
void  read_parameterfile(char *fname);
void  read_glass(char *fname);

double tk_eh(double k);

size_t my_fread(void *ptr, size_t size, size_t nmemb, FILE * stream);
size_t my_fwrite(void *ptr, size_t size, size_t nmemb, FILE * stream);

void save_local_data(void);
void add_WDM_thermal_speeds(float *vel);

int compare_type(const void *a, const void *b);

double get_fermi_dirac_vel_nu(void);
void fermi_dirac_init_nu(void);
void add_NU_thermal_speeds(float *vel);

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp