text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <time.h>
#define USAGE "./mcmc_loktavolterra.x n_steps n_burn"

void rg4th(float *resultado, float x, float y, float a, float b, float c, float d, float h);
void x_prime(float *valor, float x, float y, float alpha, float beta);
void y_prime(float *valor, float x, float y, float gamma, float delta);

void my_model(float *mine_presas,float *mine_depredadores,float *tiempo,float ini_presas, float ini_depredadores,float a, float b, float c, float d, int n_puntos);
void likelihood(float *chi,float *presas,float *depredadores,float *mine_presas,float *mine_depredadores,int n_puntos);
void mcmc(int n_steps,float *a,float *b,float *c,float *d,float *l,float *tiempo,float *presas,float *depredadores,int n_puntos,float *mine_presas,float *mine_depredadores, float ini_presas, float ini_depredadores);
void min(float *mine_presas,float *mine_depredadores,float ini_presas,float ini_depredadores,float *tiempo,float *a,float *b,float *c,float *d,float *l,int n_steps,int n_puntos);

void leer(float *tiempo, float *presas, float *depredadores, int n_puntos);
float *reserva(int n_puntos);
void print_array(float *a, float *b, float *c, float *d , float *l, int n_steps, int n_burn);

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

	srand48(time(NULL));

	//Valores que entran
	int n_steps=atof(argv[1]);
	int n_burn=atof(argv[2]);
	//Inicializacion de los arrays que guardan los datos
	int n_puntos= 96;
	float *tiempo;
	float *presas;
	float *depredadores;
	tiempo=reserva(n_puntos);
	presas=reserva(n_puntos);
	depredadores=reserva(n_puntos);

	leer(tiempo, presas, depredadores, n_puntos);
	//Condiciones iniciales para resolver las ecuaciones diferenciales
	float ini_presas, ini_depredadores;
	ini_presas=presas[0];
	ini_depredadores=depredadores[0];

	//Array con los resultados del ajuste
	float *mine_presas;
	mine_presas=reserva(n_puntos);
	float *mine_depredadores;
	mine_depredadores=reserva(n_puntos);

	//Arrays que guardan los datos de las variables
	float *a; //alpha
	float *b; //beta
	float *c; //gamma
	float *d; //delta
	float *l;
	a=reserva(n_steps);
	b=reserva(n_steps);
	c=reserva(n_steps);
	d=reserva(n_steps);
	l=reserva(n_steps);
	//Primer paso en la cadena
	a[0]=drand48();
	b[0]=drand48();
	c[0]=drand48();
	d[0]=drand48();
	float chi;
	my_model(mine_presas,mine_depredadores,tiempo,ini_presas,ini_depredadores,a[0],b[0],c[0],d[0],n_puntos);
	likelihood(&chi,presas,depredadores,mine_presas,mine_depredadores,n_puntos);
	l[0]=chi;

	mcmc(n_steps,a,b,c,d,l,tiempo,presas,depredadores,n_puntos,mine_presas,mine_depredadores,ini_presas,ini_depredadores);
	min(mine_presas,mine_depredadores,ini_presas,ini_depredadores,tiempo,a,b,c,d,l,n_steps,n_puntos);

	print_array(a,b,c,d,l,n_steps,n_burn);

	free(tiempo);
	free(presas);
	free(mine_presas);
	free(depredadores);
	free(mine_depredadores);
	free(a);
	free(b);
	free(c);
	free(d);

	return 0;
}

//OBTENCION DEL MEJOR AJUSTE
void min(float *mine_presas,float *mine_depredadores,float ini_presas,float ini_depredadores,float *tiempo,float *a,float *b,float *c,float *d,float *l,int n_steps,int n_puntos){
	int i;
	float temp, a2, b2, c2, d2;
	temp=l[0];
	a2=a[0];
	b2=b[0];
	c2=c[0];
	d2=d[0];
	for(i=1;i<n_steps;i++){
		if(temp > l[i] ){
			temp=l[i];
			a2=a[i];
			b2=b[i];
			c2=c[i];
			d2=d[i];
		}
	}
	my_model(mine_presas,mine_depredadores,tiempo,ini_presas,ini_depredadores,a2,b2,c2,d2,n_puntos);
}

//ITERACIONES
void mcmc(int n_steps,float *a,float *b,float *c,float *d,float *l,float *tiempo,float *presas,float *depredadores,int n_puntos,float *mine_presas,float *mine_depredadores, float ini_presas, float ini_depredadores){
	int i;
	float a2,b2,c2,d2,chi,chi2,alpha,beta;
	const gsl_rng_type * T;
	gsl_rng * r;
	gsl_rng_env_setup();
	T = gsl_rng_default;
	r = gsl_rng_alloc (T);

    for(i=0;i<(n_steps-1);i++){
		a2 = a[i] + gsl_ran_gaussian(r, 0.1);
		b2 = b[i] + gsl_ran_gaussian(r, 0.1);
		c2 = c[i] + gsl_ran_gaussian(r, 0.1);
		d2 = c[i] + gsl_ran_gaussian(r, 0.1);

		float *presas_init;
		float *presas_prime;
		presas_init=reserva(n_puntos);
		presas_prime=reserva(n_puntos);
		float *depredadores_init;
		float *depredadores_prime;
		depredadores_init=reserva(n_puntos);
		depredadores_prime=reserva(n_puntos);

		my_model(presas_init,depredadores_init,tiempo,ini_presas,ini_depredadores,a[i],b[i],c[i],d[i],n_puntos);
		likelihood(&chi,presas,depredadores,presas_init,depredadores_init,n_puntos);

		my_model(presas_prime,depredadores_prime,tiempo,ini_presas,ini_depredadores,a2,b2,c2,d2,n_puntos);
		likelihood(&chi2,presas,depredadores,presas_prime,depredadores_prime,n_puntos);

		alpha = -chi2 + chi;
		if(alpha >= 0.0){
			a[i+1]=a2;
			b[i+1]=b2;
			c[i+1]=c2;
			d[i+1]=d2;
			l[i+1]=chi2;
		}
		else{
			beta = drand48();
			if(alpha < -6000.0 ){
				a[i+1]=a[i];
				b[i+1]=b[i];
				c[i+1]=c[i];
				d[i+1]=d[i];
				l[i+1]=chi;
			}
			else if(beta <= exp(alpha)){
				a[i+1]=a2;
				b[i+1]=b2;
				c[i+1]=c2;
				d[i+1]=d2;
				l[i+1]=chi2;
			}
			else{
				a[i+1]=a[i];
				b[i+1]=b[i];
				c[i+1]=c[i];
				d[i+1]=d[i];
				l[i+1]=chi;
			}
		}

		free(presas_init);
		free(presas_prime);
		free(depredadores_init);
		free(depredadores_prime);
	}

  gsl_rng_free (r);
}

//IMPRESION DE LOS ARRAYS
void print_array(float *a, float *b, float *c, float *d , float *l, int n_steps, int n_burn){
  int i;
  for(i=n_burn;i<n_steps;i++){
    printf("%f %f %f %f %f\n", a[i], b[i], c[i] , d[i], l[i]);
  }
}

//DEFINICION DEL CALCULO DEL AJUSTE
void my_model(float *mine_presas,float *mine_depredadores,float *tiempo,float ini_presas, float ini_depredadores,float a, float b, float c, float d, int n_puntos){

	//Condiciones iniciales
	mine_presas[0] = ini_presas;
	mine_depredadores[0] = ini_depredadores;

	//Recorrido R-K
	int i;
	for(i=1;i<n_puntos;i++){
		float *r;
		float h;
		h = tiempo[i]-tiempo[i-1];
		r=reserva(2);
		rg4th(r, mine_presas[i-1], mine_depredadores[i-1], a, b, c, d, h);
		mine_presas[i]=r[0];
		mine_depredadores[i]=r[1];
		free(r);
	}
}

//RUNGE-KUTTA DE CUARTO ORDEN
void rg4th(float *resultado, float x, float y, float a, float b, float c, float d, float h){
    float k1_x,k1_y;
    x_prime(&k1_x,x,y,a,b);
    y_prime(&k1_y,x,y,c,d);

    //First step
    float x1,y1;
    x1 = x + (h/2.0) * k1_x;
    y1 = y + (h/2.0) * k1_y;

	float k2_x,k2_y;
    x_prime(&k2_x,x1,y1,a,b);
    y_prime(&k2_y,x1,y1,c,d);

    //Second step
    float x2,y2;
    x2 = x + (h/2.0) * k2_x;
    y2 = y + (h/2.0) * k2_y;

	float k3_x,k3_y;
	x_prime(&k3_x,x2,y2,a,b);
    y_prime(&k3_y,x2,y2,c,d);

    //Third step
    float x3,y3;
    x3 = x + h * k3_x;
    y3 = y + h * k3_y;

	float k4_x,k4_y;
	x_prime(&k4_x,x3,y3,a,b);
    y_prime(&k4_y,x3,y3,c,d);

    //Fourth step
    float average_x, average_y;
    average_x = (1.0/6.0)*(k1_x + 2.0*k2_x + 2.0*k3_x + k4_x);
    average_y = (1.0/6.0)*(k1_y + 2.0*k2_y + 2.0*k3_y + k4_y);

    resultado[0] = x + h * average_x;
    resultado[1] = y + h * average_y;
}

//ECUACION DIFERENCIAL DE PRESAS
void x_prime(float *valor, float x, float y, float alpha, float beta){
	*valor= x * (alpha - beta * y);
}

//ECUACION DIFERENCIAL DE DEPREDADORES
void y_prime(float *valor, float x, float y, float gamma, float delta){
	*valor = -1 * y * (gamma - delta * x);
}

//RESERVA DEL AJUSTE
float *reserva(int n_puntos){
  float *array;
  int i;
  if(!(array = malloc(n_puntos * sizeof(float)))){
    printf("Problema en reserva\n");
    exit(1);
  }
  for(i=0;i<n_puntos;i++){
    array[i] = 0.0;
  }
  return array;
}

//CALCULO DE LA VARIACION DEL AJUSTE
void likelihood(float *chi,float *presas,float *depredadores,float *mine_presas,float *mine_depredadores,int n_puntos){
	int i;
	float temp1=0;
	float temp2=0;
	for(i=0;i<n_puntos;i++){
		temp1=temp1+pow(presas[i]-mine_presas[i],2.0);
		temp2=temp2+pow(depredadores[i]-mine_depredadores[i],2.0);
	}
	*chi=(1.0/2.0)*(temp1+temp2);
}

//LECTURA DEL ARCHIVO
void leer(float *tiempo, float *presas, float *depredadores,int n_puntos){

	FILE *in;
	int i;
	in=fopen("archivo.dat","r");
	if(!in){
    printf("problems opening the file %s\n","archivo.dat");
    exit(1);
	}
	for(i=0;i<(n_puntos);i++){

		  fscanf(in, "%f %f %f\n", &tiempo[i], &presas[i], &depredadores[i]);

	}
	fclose(in);
}


text:
#include "als.h"
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>

#define MAX_NAME_LEN 10

/* Define global variables */
global_info info;

feature_t *movieMatrix; // feature matrix of movies
feature_t *userMatrix;  // feature matrix of users

// store rating matrix in a compressed way
int userNum;
int movieNum;
int ratingNum;

int* userStartIdx;
int* movieId;
double* movieRating;

int* movieStartIdx;
int* userId;
double* userRating;

struct movie* movie_hashtable; // to record the number of rating for each movie
struct user* user_hashtable;   // to record the number of rating for each user

void computePredictionRMSE(gsl_matrix * M, gsl_matrix * U, gsl_matrix * R);

//---------------------------------------------
//------- helper function for hashtable--------
//---------------------------------------------
void add_movie(int id, double rating) {
    struct movie* res;
    HASH_FIND_INT(movie_hashtable, &id, res);
    if (res == NULL) {
        struct movie* new_movie = (struct movie*)malloc(sizeof(struct movie));
        new_movie -> id = id;
        new_movie -> rating_num = 1;
        new_movie -> num = 1;
        new_movie -> total_rating = rating;
        HASH_ADD_INT(movie_hashtable, id, new_movie);
    } else {
        res -> rating_num ++;
        res -> num ++;
        res -> total_rating += rating;
    }
}

void delete_movie(int id) {
    struct movie* res;
    HASH_FIND_INT(movie_hashtable, &id, res);
    res -> rating_num --;
}

int find_movie(int id) {
    struct movie *result;
    HASH_FIND_INT(movie_hashtable, &id, result);
    if (result == NULL) {
        return 0;
    } else {
        return result -> rating_num;
    }
}

double get_movie_average(int id) {
    struct movie *result;
    HASH_FIND_INT(movie_hashtable, &id, result);
    if (result == NULL) {
        return 0;
    } else {
        return result -> total_rating / ((double)result -> num);
    }
}

void add_user(int id) {
    struct user* res;
    HASH_FIND_INT(user_hashtable, &id, res);
    if (res == NULL) {
        struct user* new_user = (struct user*)malloc(sizeof(struct user));
        new_user -> id = id;
        new_user -> rating_num = 1;
        HASH_ADD_INT(user_hashtable, id, new_user);
    } else {
        res -> rating_num ++;
    }
}

void delete_user(int id) {
    struct user* res;
    HASH_FIND_INT(user_hashtable, &id, res);
    res -> rating_num --;
}

int find_user(int id) {
    struct user *result;
    HASH_FIND_INT(user_hashtable, &id, result);
    if (result == NULL) {
        return 0;
    } else {
        return result -> rating_num;
    }
}

//------------------------------------------------
//-----------process input data-------------------
//------------------------------------------------

/*
 * get the input relevant stat
 */
void getInputStat(char* inputFilename) {
    userNum = 0;
    movieNum = 0;
    ratingNum = 0;

    FILE* fp = fopen(inputFilename, "r");
    if (fp == NULL) {
        fprintf(stderr, "Failed to read input file %s\n", inputFilename);
        exit(-1);
    }

    ssize_t read;
    size_t len = 0;
    char* line = NULL;
    int uid = 0;
    int mid = 0;
    double rating = 0;
    if ((read = getline(&line, &len, fp)) != -1) {
        printf("Start reading file, get input stat\n");
    }

    while ((read = getline(&line, &len, fp)) != -1) {
        ratingNum++;

        // split the data line with comma
        char* tmp = line;
        char* last_pos = line;
        int word_len = 0;
        int comma_cnt = 0;

        while (comma_cnt < 3) {
            if (*tmp == ',') {
                word_len = tmp - last_pos;
                char tmpbuf[MAX_NAME_LEN] = "";
                strncpy(tmpbuf, last_pos, word_len);
                if (comma_cnt == 0) {
                    uid = atoi(tmpbuf);
                } else if (comma_cnt == 1) {
                    mid = atoi(tmpbuf);
                } else {
                    rating = atof(tmpbuf);
                }
                comma_cnt++;
                last_pos = tmp + 1;
            }
            tmp++;
        }
        add_movie(mid, rating);  // keep record of rating number for each movie
        add_user(uid);           // keep record of rating number of each user
        userNum = max(userNum, uid);
        movieNum = max(movieNum, mid);
    }
    printf("There are %d ratings, %d users, %d movies\n", ratingNum, userNum, movieNum);
}

void initUserStartIndex() {
    int i;
    int sum = 0;
    for (i = 1; i <= userNum; ++i) {
        int user_rating = find_user(i);
        userStartIdx[i - 1] = sum;
        sum += user_rating;
    }
    userStartIdx[userNum] = ratingNum;
}

void initMovieStartIndex() {
    int i;
    int sum = 0;
    for (i = 1; i <= movieNum; ++i) {
        int movie_rating = find_movie(i);
        movieStartIdx[i - 1] = sum;
        sum += movie_rating;
    }
    movieStartIdx[movieNum] = ratingNum;
}

void initRatingMatrix(int uid, int mid, double rating) {
    int movieIdx = movieStartIdx[mid] - find_movie(mid);
    int userIdx = userStartIdx[uid] - find_user(uid);
    userId[movieIdx] = uid;
    userRating[movieIdx] = rating;
    movieId[userIdx] = mid;
    movieRating[userIdx] = rating;
    delete_movie(mid);
    delete_user(uid);
}

void initMatrix() {
    printf("Init movie matrix\n");
    int i, j;
    for (i = 0; i < movieNum; ++i) {
        int start_idx = i * info.numFeatures;
        movieMatrix[start_idx] = get_movie_average(i + 1);
        // initialize the random small value
        for (j = 1; j < info.numFeatures; ++j) {
            movieMatrix[start_idx + j] = (rand() % 1000) / 1000.0;
        }
    }
}

// read input file
// initializa the rating matrix & movie matrix
void readInput(char* inputFilename) {
    FILE* fp = fopen(inputFilename, "r");
    if (fp == NULL) {
        fprintf(stderr, "Failed to read input file %s\n", inputFilename);
        exit(-1);
    }

    ssize_t read;
    size_t len = 0;
    char* line = NULL;
    int uid = 0;
    int mid = 0;
    double rating = 0;
    if ((read = getline(&line, &len, fp)) != -1) {
        printf("Start reading file, get rating details\n");
    }

    while ((read = getline(&line, &len, fp)) != -1) {
        // split the data line with comma
        char* tmp = line;
        char* last_pos = line;
        int word_len = 0;
        int comma_cnt = 0;

        while (comma_cnt < 3) {
            if (*tmp == ',') {
                word_len = tmp - last_pos;
                char tmpbuf[MAX_NAME_LEN] = "";
                strncpy(tmpbuf, last_pos, word_len);
                if (comma_cnt == 0) {
                    uid = atoi(tmpbuf);
                } else if (comma_cnt == 1) {
                    mid = atoi(tmpbuf);
                } else {
                    rating = atof(tmpbuf);
                }
                comma_cnt++;
                last_pos = tmp + 1;
            }
            tmp++;
        }
        // initialize the rating matrix (use compression)
        initRatingMatrix(uid, mid, rating);
    }
}

// init basic arg
void init(int numFeatures, int numIter, double lambda) {
    info.numIter = numIter;
    info.numFeatures = numFeatures;
    info.lambda = lambda;
}

/*
 * allocate data structure memory
 */
void init_data(int userNum, int movieNum, int ratingNum, int nproc) {
    printf("allocating data memory\n");

    int movieSpan = (movieNum + nproc - 1) / nproc;
    int userSpan = (userNum + nproc - 1) / nproc;

    userStartIdx = (int*)calloc(userNum + 1, sizeof(int));
    movieId =      (int*)calloc(ratingNum + 1, sizeof(int));
    movieRating =  (double*)calloc(ratingNum + 1, sizeof(double));

    movieStartIdx = (int*)calloc(movieNum + 1, sizeof(int));
    userId =        (int*)calloc(ratingNum + 1, sizeof(int));
    userRating =    (double*)calloc(ratingNum + 1, sizeof(double));

    movieMatrix = (feature_t *)calloc(info.numFeatures * movieSpan * nproc, sizeof(feature_t));
    userMatrix =  (feature_t *)calloc(info.numFeatures * userSpan * nproc, sizeof(feature_t));
}

//-----------------------------------------------------------------
//---------------------real computation----------------------------
//-----------------------------------------------------------------
void compute(int procID, int nproc, char* inputFilename,
             int numFeatures, int numIterations, double lambda)
{
    const int root = 0; // set the rank 0 processor as the root
    int span;

    // initialize
    init(numFeatures, numIterations, lambda);

    /* Read the input file and initialization */
    if (procID == root) {
        getInputStat(inputFilename);
    }

    /* broadcast the information of basic stat*/
    MPI_Bcast(&userNum, 1, MPI_INT, root, MPI_COMM_WORLD);
    MPI_Bcast(&movieNum, 1, MPI_INT, root, MPI_COMM_WORLD);
    MPI_Bcast(&ratingNum, 1, MPI_INT, root, MPI_COMM_WORLD);

    // initialize the data structure
    init_data(userNum, movieNum, ratingNum, nproc);

    if (procID == root) {
        initUserStartIndex();
        initMovieStartIndex();
    }

    /* broadcast the information about start index */
    MPI_Bcast(userStartIdx, userNum + 1, MPI_INT, root, MPI_COMM_WORLD);
    MPI_Bcast(movieStartIdx, movieNum + 1, MPI_INT, root, MPI_COMM_WORLD);

    /* Read the input file and get detailed rating */
    if (procID == root) {
        readInput(inputFilename);
    }

    MPI_Bcast(movieId, ratingNum + 1, MPI_INT, root, MPI_COMM_WORLD);
    MPI_Bcast(movieRating, ratingNum + 1, MPI_DOUBLE, root, MPI_COMM_WORLD);
    MPI_Bcast(userId, ratingNum + 1, MPI_INT, root, MPI_COMM_WORLD);
    MPI_Bcast(userRating, ratingNum + 1, MPI_DOUBLE, root, MPI_COMM_WORLD);

    // initialize movie/ user matrix
    if (procID == root) {
        initMatrix();
    }

    MPI_Bcast(movieMatrix, numFeatures * movieNum, MPI_DOUBLE, root, MPI_COMM_WORLD);

    // initialize n_user and n_movie
    int* n_user =  (int *)malloc(sizeof(int) * userNum);
    int* n_movie = (int *)malloc(sizeof(int) * movieNum);
    int i, j;
    for (i = 0; i < userNum; ++i) {
        n_user[i] = userStartIdx[i+1] - userStartIdx[i];
    }
    for (i = 0; i < movieNum; ++i) {
        n_movie[i] = movieStartIdx[i+1] - movieStartIdx[i];
    }

    // create gsl objects
    gsl_matrix * M = gsl_matrix_alloc (numFeatures, movieNum);
    gsl_matrix * U = gsl_matrix_alloc (numFeatures, userNum);
    gsl_matrix * R = gsl_matrix_alloc (userNum, movieNum);

    gsl_matrix * A = gsl_matrix_alloc (numFeatures, numFeatures);
    gsl_matrix * Ainv = gsl_matrix_alloc (numFeatures, numFeatures);

    gsl_matrix * E = gsl_matrix_alloc (numFeatures, numFeatures);
    for (i = 0; i < numFeatures; ++i)
        gsl_matrix_set(E, i, i, 1);
    gsl_matrix * F = gsl_matrix_alloc (numFeatures, numFeatures);

    gsl_vector * V = gsl_vector_alloc (numFeatures);

    gsl_vector * Rm = gsl_vector_alloc (movieNum);
    gsl_vector * Ru = gsl_vector_alloc (userNum);

    gsl_permutation *p = gsl_permutation_alloc(numFeatures);
    gsl_vector * O = gsl_vector_alloc (numFeatures);

    /* Start */
    // start iteration
    if (procID == root) {
        printf("Start iteration\n");
        computePredictionRMSE(M, U, R);
    }

    int iter = 0;
    for (iter = 0; iter < numIterations; ++iter) {

        // Each processor solve user feature
        if (procID == root) {
            printf("Updating user (iter = %d)\n", iter);
        }

        int user_idx;
        int user_num = 0;

        // check span
        // NOTE THE NUMBER SHOULD BE DIVISIBLE
        span = (userNum + nproc - 1) / nproc;
        int user_start_idx = min(procID * span, userNum);
        int user_end_idx = min(user_start_idx + span, userNum);

        for (user_idx = user_start_idx; user_idx < user_end_idx; ++user_idx) {
            user_num++;

            gsl_matrix_set_zero (M);
            for (i = userStartIdx[user_idx]; i < userStartIdx[user_idx+1]; ++i)
                for (j = 0; j < numFeatures; ++j)
                    gsl_matrix_set(M, j, i - userStartIdx[user_idx], movieMatrix[(movieId[i]-1) * numFeatures + j]);

            gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, M, M, 0, A);
            gsl_matrix_memcpy(F, E);
            gsl_matrix_scale(F, lambda * n_user[user_idx]);
            gsl_matrix_add(A, F); // A = M M^T + lambda * n * E

            gsl_vector_set_zero (Rm);
            for (i = userStartIdx[user_idx]; i < userStartIdx[user_idx+1]; ++i)
                gsl_vector_set(Rm, i - userStartIdx[user_idx], movieRating[i]);

            gsl_blas_dgemv(CblasNoTrans, 1, M, Rm, 0, V); // V = M R^T

            int s;
            gsl_linalg_LU_decomp(A, p, &s);
            gsl_linalg_LU_invert(A, p, Ainv);

            gsl_blas_dgemv(CblasNoTrans, 1, Ainv, V, 0, O); // O = A^-1 V

            for (j = 0; j < numFeatures; ++j)
                userMatrix[user_idx * numFeatures + j] = gsl_vector_get(O, j);
        }

        // allgather (everyone gets a local copy of U)
        MPI_Allgather(&userMatrix[user_start_idx * numFeatures],
                      span * numFeatures, MPI_DOUBLE,
                      userMatrix, span * numFeatures,
                      MPI_DOUBLE, MPI_COMM_WORLD);

        // solver movie feature matrix
        if (procID == root) {
            printf("Updating movie (iter = %d)\n", iter);
        }

        int movie_idx;
        int movie_num = 0;

        span = (movieNum + nproc - 1) / nproc;
        int movie_start_idx = min(procID * span, movieNum);
        int movie_end_idx = min(movie_start_idx + span, movieNum);

        for (movie_idx = movie_start_idx; movie_idx < movie_end_idx; ++movie_idx) {
            if (movieStartIdx[movie_idx] == movieStartIdx[movie_idx+1])
                continue;

            movie_num++;
            gsl_matrix_set_zero (U);
            for (i = movieStartIdx[movie_idx]; i < movieStartIdx[movie_idx+1]; ++i)
                for (j = 0; j < numFeatures; ++j)
                    gsl_matrix_set(U, j, i - movieStartIdx[movie_idx], userMatrix[(userId[i]-1) * numFeatures + j]);

            gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1, U, U, 0, A);
            gsl_matrix_memcpy(F, E);
            gsl_matrix_scale(F, lambda * n_movie[movie_idx]);
            gsl_matrix_add(A, F); // A = U U^T + lambda * n * E

            gsl_vector_set_zero (Ru);
            for (i = movieStartIdx[movie_idx]; i < movieStartIdx[movie_idx+1]; ++i)
                gsl_vector_set(Ru, i - movieStartIdx[movie_idx], userRating[i]);

            gsl_blas_dgemv(CblasNoTrans, 1, U, Ru, 0, V); // V = U R^T

            int s;
            gsl_linalg_LU_decomp(A, p, &s);
            gsl_linalg_LU_invert(A, p, Ainv);

            gsl_blas_dgemv(CblasNoTrans, 1, Ainv, V, 0, O); // O = A^-1 V

            for (j = 0; j < numFeatures; ++j)
                movieMatrix[movie_idx * numFeatures + j] = gsl_vector_get(O, j);
        }

        // allgather (everyone gets a local copy of M)
        MPI_Allgather(&movieMatrix[movie_start_idx * numFeatures],
                      span * numFeatures, MPI_DOUBLE,
                      movieMatrix, span * numFeatures,
                      MPI_DOUBLE, MPI_COMM_WORLD);

        // compute prediction and rmse
        if (procID == root) {
            computePredictionRMSE(M, U, R);
        }

    }

    gsl_matrix_free (M);
    gsl_matrix_free (U);
    gsl_matrix_free (R);
    gsl_matrix_free (A);
    gsl_matrix_free (Ainv);
    gsl_matrix_free (E);
    gsl_matrix_free (F);

    gsl_vector_free (V);
    gsl_vector_free (Rm);
    gsl_vector_free (Ru);
    gsl_vector_free (O);

    gsl_permutation_free (p);

    free(n_user);
    free(n_movie);

    free(userStartIdx);
    free(movieId);
    free(movieRating);
    free(movieStartIdx);
    free(userId);
    free(userRating);
    free(movieMatrix);
    free(userMatrix);
}

void computePredictionRMSE(gsl_matrix * M, gsl_matrix * U, gsl_matrix * R) {
    int i, j, user_idx;

    for (i = 0; i < movieNum; ++i)
        for (j = 0; j < info.numFeatures; ++j)
            gsl_matrix_set(M, j, i, movieMatrix[i * info.numFeatures + j]);

    for (i = 0; i < userNum; ++i)
        for (j = 0; j < info.numFeatures; ++j)
            gsl_matrix_set(U, j, i, userMatrix[i * info.numFeatures + j]);

    gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, U, M, 0, R);

    double rmse = 0;

    for (user_idx = 0; user_idx < userNum; ++user_idx)
        for (i = userStartIdx[user_idx]; i < userStartIdx[user_idx+1]; ++i) {
            double diff = movieRating[i] - gsl_matrix_get(R, user_idx, movieId[i] - 1);
            rmse += diff * diff;
        }

    printf("RMSE = %f\n", rmse / ratingNum);
}

ext: c


text:
#pragma once

#include <memory_resource>
#include <vector>

#include <gsl-lite/gsl-lite.hpp>

#include <thrust/sort.h>
#include <cuda/define_specifiers.hpp>

#include <thrustshift/memory-resource.h>

namespace thrustshift {

//! Col indices must be ordered
template <typename DataType, typename IndexType>
class CSR {
   public:
	using value_type = DataType;
	using index_type = IndexType;

   private:
	bool cols_are_sorted() {
		for (size_t row_id = 0; row_id < this->num_rows(); ++row_id) {
			if (!std::is_sorted(col_indices_.begin() + row_ptrs_[row_id],
			                    col_indices_.begin() + row_ptrs_[row_id + 1])) {
				return false;
			}
		}
		return true;
	}

   public:
	CSR() : row_ptrs_(1, 0, &pmr::default_resource), num_cols_(0) {
	}

	template <class DataRange,
	          class ColIndRange,
	          class RowPtrsRange,
	          class MemoryResource>
	CSR(DataRange&& values,
	    ColIndRange&& col_indices,
	    RowPtrsRange&& row_ptrs,
	    size_t num_cols,
	    MemoryResource& memory_resource)
	    : values_(values.begin(), values.end(), &memory_resource),
	      col_indices_(col_indices.begin(),
	                   col_indices.end(),
	                   &memory_resource),
	      row_ptrs_(row_ptrs.begin(), row_ptrs.end(), &memory_resource),
	      num_cols_(num_cols) {
		gsl_Expects(values.size() == col_indices.size());
		gsl_Expects(row_ptrs.size() > 0);
		gsl_Expects(row_ptrs[0] == 0);
		gsl_ExpectsAudit(cols_are_sorted());
	}

	template <class DataRange, class ColIndRange, class RowPtrsRange>
	CSR(DataRange&& values,
	    ColIndRange&& col_indices,
	    RowPtrsRange&& row_ptrs,
	    size_t num_cols)
	    : CSR(std::forward<DataRange>(values),
	          std::forward<ColIndRange>(col_indices),
	          std::forward<RowPtrsRange>(row_ptrs),
	          num_cols,
	          pmr::default_resource) {
	}

	// The copy constructor is declared explicitly to ensure
	// managed memory is used per default.
	CSR(const CSR& other)
	    : CSR(other.values(),
	          other.col_indices(),
	          other.row_ptrs(),
	          other.num_cols()) {
	}

	CSR(CSR&& other) = default;

	CSR& operator=(const CSR& other) = default;

	gsl_lite::span<DataType> values() {
		return gsl_lite::make_span(values_);
	}

	gsl_lite::span<const DataType> values() const {
		return gsl_lite::make_span(values_);
	}

	gsl_lite::span<IndexType> col_indices() {
		return gsl_lite::make_span(col_indices_);
	}

	gsl_lite::span<const IndexType> col_indices() const {
		return gsl_lite::make_span(col_indices_);
	}

	gsl_lite::span<IndexType> row_ptrs() {
		return gsl_lite::make_span(row_ptrs_);
	}

	gsl_lite::span<const IndexType> row_ptrs() const {
		return gsl_lite::make_span(row_ptrs_);
	}

	size_t num_rows() const {
		return row_ptrs_.size() - 1;
	}

	size_t num_cols() const {
		return num_cols_;
	}

	void sort_column_indices() {
		for (size_t row_id = 0; row_id < num_rows(); ++row_id) {
			const auto nns_start = row_ptrs_[row_id];
			const auto nns_end = row_ptrs_[row_id + 1];
			thrust::sort_by_key(thrust::host,
			                    col_indices_.begin() + nns_start,
			                    col_indices_.begin() + nns_end,
			                    values_.begin() + nns_start);
		}
	}

	/*! \brief Add additional elements to each row of the matrix.
	 *
	 *  A row might already have the maximum number of element. Therefore
	 *  not every row might be extended by exactly `num_max_additional_elements_per_row`.
	 *
	 */
	void extend_rows(int num_max_additional_elements_per_row,
	                 DataType value = 0) {
		if (num_max_additional_elements_per_row == 0) {
			return;
		}
		gsl_Expects(num_max_additional_elements_per_row >= 0);
		std::vector<DataType> tmp_values(values_.begin(), values_.end());
		std::vector<IndexType> tmp_col_indices(col_indices_.begin(),
		                                       col_indices_.end());
		size_t num_additional_elements = 0;
		const int nc = gsl_lite::narrow<int>(num_cols());
		auto get_num_additional_elements_per_row = [&](size_t row_id) {
			const int num_elements_curr_row =
			    row_ptrs_[row_id + 1] - row_ptrs_[row_id];
			return std::max(std::min(num_max_additional_elements_per_row,
			                         nc - num_elements_curr_row),
			                0);
		};
		for (size_t row_id = 0; row_id < num_rows(); ++row_id) {
			num_additional_elements +=
			    get_num_additional_elements_per_row(row_id);
		}

		const size_t new_nnz = values_.size() + num_additional_elements;
		values_.resize(new_nnz);
		col_indices_.resize(new_nnz);
		int nns_offset = 0;
		for (size_t row_id = 0; row_id < num_rows(); ++row_id) {
			const int num_additional_elements_curr_row =
			    get_num_additional_elements_per_row(row_id);
			int nns_id = row_ptrs_[row_id];
			const int nns_end = row_ptrs_[row_id + 1];
			const int nnz_curr_row = nns_end - nns_id;
			int curr_col_id = 0;
			int num_added_elements = 0;
			for (int new_nns_id = nns_offset;
			     new_nns_id < nns_offset + num_additional_elements_curr_row +
			                      row_ptrs_[row_id + 1] - row_ptrs_[row_id];
			     ++new_nns_id) {
				const int other_col_id = nns_id == nns_end
				                             ? std::numeric_limits<int>::max()
				                             : tmp_col_indices[nns_id];
				if (curr_col_id < other_col_id &&
				    num_added_elements < num_additional_elements_curr_row) {
					values_[new_nns_id] = value;
					col_indices_[new_nns_id] = curr_col_id;
					++curr_col_id;
					++num_added_elements;
				}
				else {
					values_[new_nns_id] = tmp_values[nns_id];
					col_indices_[new_nns_id] = other_col_id;
					curr_col_id = other_col_id + 1;
					++nns_id;
				}
			}
			row_ptrs_[row_id] = nns_offset;
			nns_offset += nnz_curr_row + num_additional_elements_curr_row;
		}
		row_ptrs_.back() = values_.size();
		gsl_ExpectsAudit(std::is_sorted(row_ptrs_.begin(), row_ptrs_.end()));
		gsl_ExpectsAudit(cols_are_sorted());
	}

   private:
	std::pmr::vector<DataType> values_;
	std::pmr::vector<IndexType> col_indices_;
	std::pmr::vector<IndexType> row_ptrs_;
	size_t num_cols_;
};

template <typename DataType, typename IndexType>
bool operator==(const CSR<DataType, IndexType>& a,
                const CSR<DataType, IndexType>& b) {
	return equal(a.values(), b.values()) &&
	       equal(a.col_indices(), b.col_indices()) &&
	       equal(a.row_ptrs(), b.row_ptrs()) &&
	       a.num_rows() == b.num_rows() && a.num_cols() == b.num_cols();
}

template <typename DataType, typename IndexType>
class CSR_view {

   public:
	using value_type = DataType;
	using index_type = IndexType;

	template <typename OtherDataType, typename OtherIndexType>
	CSR_view(CSR<OtherDataType, OtherIndexType>& owner)
	    : values_(owner.values()),
	      col_indices_(owner.col_indices()),
	      row_ptrs_(owner.row_ptrs()),
	      num_cols_(owner.num_cols()) {
	}

	template <typename OtherDataType, typename OtherIndexType>
	CSR_view(const CSR<OtherDataType, OtherIndexType>& owner)
	    : values_(owner.values()),
	      col_indices_(owner.col_indices()),
	      row_ptrs_(owner.row_ptrs()),
	      num_cols_(owner.num_cols()) {
	}

	template <class DataRange,
	          class ColIndRange,
	          class RowPtrsRange>
	CSR_view(DataRange&& values,
	    ColIndRange&& col_indices,
	    RowPtrsRange&& row_ptrs,
	    size_t num_cols)
	    : values_(values),
	      col_indices_(col_indices),
	      row_ptrs_(row_ptrs),
	      num_cols_(num_cols) {
		gsl_Expects(values.size() == col_indices.size());
		gsl_Expects(row_ptrs.size() > 0);
		gsl_ExpectsAudit(row_ptrs[0] == 0);
		gsl_ExpectsAudit(row_ptrs.back() == values.size());
		gsl_ExpectsAudit(std::is_sorted(row_ptrs.begin(), row_ptrs.end()));
	}

	CSR_view(const CSR_view& other) = default;
	CSR_view(CSR_view&& other) = default;

	CUDA_FHD gsl_lite::span<DataType> values() {
		return values_;
	}

	CUDA_FHD gsl_lite::span<const DataType> values() const {
		return values_;
	}

	CUDA_FHD gsl_lite::span<IndexType> col_indices() {
		return col_indices_;
	}

	CUDA_FHD gsl_lite::span<const IndexType> col_indices() const {
		return col_indices_;
	}

	CUDA_FHD gsl_lite::span<IndexType> row_ptrs() {
		return row_ptrs_;
	}

	CUDA_FHD gsl_lite::span<const IndexType> row_ptrs() const {
		return row_ptrs_;
	}

	CUDA_FHD size_t num_rows() const {
		return row_ptrs_.size() - 1;
	}

	CUDA_FHD size_t num_cols() const {
		return num_cols_;
	}

	CUDA_FHD auto max_row_nnz() const {
		using I = typename std::remove_const<IndexType>::type;
		I mnnz = 0;
		for (size_t row_id = 1; row_id < row_ptrs_.size(); ++row_id) {
			mnnz = std::max(row_ptrs_[row_id] - row_ptrs_[row_id - 1], mnnz);
		}
		return mnnz;
	}

   private:
	gsl_lite::span<DataType> values_;
	gsl_lite::span<IndexType> col_indices_;
	gsl_lite::span<IndexType> row_ptrs_;
	size_t num_cols_;
};

} // namespace thrustshift

ext: h

max_stars_repo_name: codecircuit/thrustshift

max_issues_repo_name: codecircuit/thrustshift

max_forks_repo_name: codecircuit/thrustshift


text:
#ifndef _H_LIMBER
#define _H_LIMBER

#include <gsl/gsl_spline.h>
#include <stdbool.h>
#include "interp2d.h"
//#include <gsl/gsl_interp2d.h>
//#include <gsl/gsl_spline2d.h>

// The integrated C_ell are in general allowed to be zero or negative if
// they describe cross-correlations. We use these statuses to describe
// errors where a log was requested also.
// LIMBER_STATUS_NEGATIVE is probably always an error
// but LIMBER_STATUS_ZERO is not necessarily
#define LIMBER_STATUS_OK 0
#define LIMBER_STATUS_ZERO 1
#define LIMBER_STATUS_NEGATIVE 2
#define LIMBER_STATUS_ERROR 3

// These are the options you can set for
// the Limber integrator.
typedef struct limber_config{
	bool xlog;  // The output spline will be in terms of log(ell) not ell
	bool ylog;
	int n_ell;  // Number of ell values you want in the spline
	double * ell;  // The chosen ell values you want
	double prefactor; //Scaling prefactor
    int status; // did everything go okay?
    double absolute_tolerance;
    double relative_tolerance;
} limber_config;

gsl_spline * limber_integral(limber_config * config,
	gsl_spline * WbX, gsl_spline * WfX, gsl_spline * WmX, gsl_spline * WbY, gsl_spline * WfY, gsl_spline * WmY, Interpolator2D * P, Interpolator2D * f, Interpolator2D * D, Interpolator2D * BB);

#endif


text:
//! \file reaction.h
//! Data for an incident neutron reaction

#ifndef OPENMC_REACTION_H
#define OPENMC_REACTION_H

#include <string>

#include <gsl/gsl>
#include "hdf5.h"

#include "openmc/reaction_product.h"
#include "openmc/vector.h"

namespace openmc {

//==============================================================================
//! Data for a single reaction including cross sections (possibly at multiple
//! temperatures) and reaction products (with secondary angle-energy
//! distributions)
//==============================================================================

class Reaction {
public:
  //! Construct reaction from HDF5 data
  //! \param[in] group HDF5 group containing reaction data
  //! \param[in] temperatures Desired temperatures for cross sections
  explicit Reaction(hid_t group, const vector<int>& temperatures);

  //! \brief Calculate reaction rate based on group-wise flux distribution
  //
  //! \param[in] i_temp Temperature index
  //! \param[in] energy Energy group boundaries in [eV]
  //! \param[in] flux Flux in each energy group (not normalized per eV)
  //! \param[in] grid Nuclide energy grid
  //! \return Reaction rate
  double collapse_rate(gsl::index i_temp, gsl::span<const double> energy,
    gsl::span<const double> flux, const vector<double>& grid) const;

  //! Cross section at a single temperature
  struct TemperatureXS {
    int threshold;
    vector<double> value;
  };

  int mt_;             //!< ENDF MT value
  double q_value_;     //!< Reaction Q value in [eV]
  bool scatter_in_cm_; //!< scattering system in center-of-mass?
  bool redundant_;     //!< redundant reaction?
  vector<TemperatureXS> xs_;         //!< Cross section at each temperature
  vector<ReactionProduct> products_; //!< Reaction products
};

//==============================================================================
// Non-member functions
//==============================================================================

//! Return reaction name given an ENDF MT value
//
//! \param[in] mt  ENDF MT value
//! \return Name of the corresponding reaction
std::string reaction_name(int mt);

//! Return reaction type (MT value) given a reaction name
//
//! \param[in] name  Reaction name
//! \return Corresponding reaction type (MT value)
int reaction_type(std::string name);

} // namespace openmc

#endif // OPENMC_REACTION_H

max_stars_repo_name: cjwyett/openmc

max_issues_repo_name: cjwyett/openmc

max_forks_repo_name: cjwyett/openmc


text:
/* testutilities.c
 *
 * 24-May-2008 Initial write: Ting Zhao
 */

#include <stdlib.h>
#include <string.h>
#include <utilities.h>
#include <math.h>
#include "tz_utilities.h"
#include "tz_error.h"
#include "tz_file_list.h"
#include "tz_image_lib_defs.h"
#include "tz_stack_utils.h"

/*
#include "tz_matlabio.h"
#include "tz_rastergeom.h"
#include <gsl/gsl_math.h>
*/
int main(int argc, char *argv[])
{
  static char *Spec[] = {"[-t]", NULL};

  Process_Arguments(argc, argv, Spec, 1);

  if (Is_Arg_Matched("-t")) {
    /* example test */

    /* extract file name without extension from a path */
    const char *path = "../data/test.tif";
    char *name = fname(path, NULL);
    if (strcmp(name, "test") != 0) {
      PRINT_EXCEPTION("Bug?", "Unexpected file name.");
      return 1;
    }

    /* The user is responsible for freeing the returned value */
    free(name);

    /* The routine returns NULL when no name can be extracted */
    path = ".";
    if (fname(path, NULL) != NULL) {
      PRINT_EXCEPTION("Bug?", "Unexpected file name.");
      return 1;
    }

    path = "../data/";
    if (fname(path, NULL) != NULL) {
      PRINT_EXCEPTION("Bug?", "Unexpected file name.");
      return 1;
    }

    /* Different parts can be put together to form a path */
    char new_path[500];
    fullpath_e("../data/", "test", "tif", new_path);
    if (strcmp(new_path, "../data/test.tif") != 0) {
      printf("%s\n", new_path);
      PRINT_EXCEPTION("Bug?", "Unexpected file path.");
      return 1;
    }

    fullpath_e("", "test", "tif", new_path);
    if (strcmp(new_path, "test.tif") != 0) {
      printf("%s\n", new_path);
      PRINT_EXCEPTION("Bug?", "Unexpected file path.");
      return 1;
    }

    { /* unit test */
      char buffer[500];

      if (fname(NULL, buffer) != NULL) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (fname("\0", buffer) != NULL) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (fname(".", buffer) != NULL) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (strcmp(fname("t", buffer), "t") != 0) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (fname("/", buffer) != NULL) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (fname("test/", buffer) != NULL) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (strcmp(fname("test", buffer), "test") != 0) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (strcmp(fname("./test", buffer), "test") != 0) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (strcmp(fname("/test", buffer), "test") != 0) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (strcmp(fname("/.test", buffer), ".test") != 0) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (strcmp(fname("t.est", buffer), "t") != 0) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (strcmp(fname("/t.est", buffer), "t") != 0) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (fname("/.", buffer) != NULL) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (fname("data/.", buffer) != NULL) {
	PRINT_EXCEPTION("Bug?", "Unexpected file name.");
	return 1;
      }

      if (strcmp(fullpath_e("dir", "file", "ext", buffer),
		 "dir/file.ext") != 0) {
	PRINT_EXCEPTION("Bug?", "Unexpected path name.");
	return 1;
      }

      if (strcmp(fullpath_e("dir/", "file", "ext", buffer),
		 "dir/file.ext") != 0) {
	PRINT_EXCEPTION("Bug?", "Unexpected path name.");
	return 1;
      }

      if (strcmp(fullpath_e("dir", "file", NULL, buffer),
		 "dir/file") != 0) {
	PRINT_EXCEPTION("Bug?", "Unexpected path name.");
	return 1;
      }

      if (strcmp(fullpath_e("dir", "file", ".", buffer),
		 "dir/file.") != 0) {
	printf("%s\n", buffer);
	PRINT_EXCEPTION("Bug?", "Unexpected path name.");
	return 1;
      }

      if (strcmp(fullpath_e("dir", "file", ".ext", buffer),
		 "dir/file.ext") != 0) {
	printf("%s\n", buffer);
	PRINT_EXCEPTION("Bug?", "Unexpected path name.");
	return 1;
      }

    }

    printf(":) Testing passed.\n");
    return 0;
  }

#if 0
  FILE *fp = fopen("../data/test.fbd", "r");

  size_t len;
  char *line;
  while (!feof(fp)) {
    if ((line = fgetln(fp, &len)) != NULL) {
      printf("read %lu: ", len);
      line[len-1] = '\0';
      printf("%s\n", line);
    }
  }

  fclose(fp);
#endif

#if 0
  char path[100];
  Find_Matlab_Binpath(path);
  printf("%s : %lu\n", path, strlen(path));
#endif

#if 0
  int a = 10000;
  int b = 11;
  XOR_SWAP(a, b);

  printf("%d %d\n", a, b);
#endif

#if 0
  printf("%d\n", 2 << 2);
  printf("%d\n", Compare_Float(1.01, 1.002, 0.0005));
  printf("%d\n", gsl_fcmp(1.01, 1.002, 0.0005));
#endif

#if 0
  char *s1 = "hello";
  printf("%p\n", s1);
  s1[0] = 't';
  printf("%s\n", s1);
#endif

#if 0
  double x = -100.4350;
  printf("%g\n", Cube_Root(x) * Cube_Root(x) * Cube_Root(x));
#endif

#if 0
  printf("%d\n", Raster_Line_Map(852, 451, 296));
  return 1;

  int m = 33;
  int n = 10;
  int i;
  for (i = 0; i < m; i++) {
    printf("%d ", Raster_Line_Map(m, n, i));
  }
  printf("\n");
#endif

#if 0
  int w1 = 302;
  int h1 = 586;
  int nw2;
  int nh2;
  Raster_Ratio_Scale(w1, h1, 425, 236, &nw2, &nh2);
  printf("%d x %d\n", nw2, nh2);
  printf("%g => %g\n", (double) w1 / h1, (double) nw2 / nh2);
#endif

#if 0
  //printf("%u\n", Hexstr_To_Uint("BA9B"));
  char str[9];
  printf("%s\n", Uint_To_Hexstr(262, str));
#endif

#if 0
  FILE *fp = fopen("/Users/zhaot/Work/neurolabi/data/movie/a0.BMP", "r");
  Fprint_File_Binary(fp, 20, stdout);
  fclose(fp);
#endif

#if 0
  printf("%d\n", Round_Div(110, -19));
#endif

#if 0
  printf("%d\n", Raster_Linear_Map(31, 10, 100, 0, 231));
#endif

#if 0
  printf("%d\n", Raster_Point_Zoom_Offset(0, 2, 100, 512, 500, 1));
#endif

#if 0
  fcopy("../data/test.tif", "../data/test2.tif");
#endif

#if 0
  printf("%d bytes\n", fsize("../data/test.tif"));
#endif

#if 0
  printf("%g %g %g\n", Infinity, -Infinity, NaN);
#endif

#if 0
  size_t x = 1000;
  printf("%zd\n", (size_t) 1000 * 1000 * 1000 * 4);
#endif

#if 0
  size_t index = Stack_Util_Offset(1000, 1000, 1000, 2000, 2000, 2000);
  printf("%zd\n", index);
  int x, y, z;
  Stack_Util_Coord(index, 2000, 2000, &x, &y, &z);
  printf("%d %d %d\n", x, y, z);
#endif

#if 0
  printf("%s\n", fname("../data/test.dot", NULL));
  printf("%s\n", dname("../data/test.dot", NULL));
  printf("%s\n", dname("/test.dot", NULL));
  printf("%s\n", dname("test.dot", NULL));
#endif

#if 0
  char *line = malloc(2048 * 2048 * 4);
  char *data = malloc(1024 * 1024 * 4);
  int i;
  int j;
  tic();
  for (i = 0; i < 2048; i++) {
    for (j = 0; j < 2048; j++) {
      if (*line != *data++) {
	*line++ = *data;
	*line++ = *data;
	*line = *data;
	line += 2;
      }
      /*
      memset(line, *data++, 3);
      line += 4;
      */
    }
  }
  printf("%lld\n", toc());
#endif

#if 0
  char buf[11];
  Memset_Pattern4( buf, "1234", 10 );
  buf[10] = '\0';
  printf("%s\n", buf);

  Memset_Pattern4( buf, "4321", 8 );
  buf[8] = '\0';
  printf("%s\n", buf);

  Memset_Pattern4( buf, "4321", 2);
  buf[2] = '\0';
  printf("%s\n", buf);
#endif

#if 0
  char *block[1000000];
  int i;

  tic();
  for (i = 0; i < 1000000; i++) {
    block[i] = (char*) malloc(i % 991 + 10);
  }
  printf("%lld\n", toc());

  int index[1000000];
  int j;
  int offset = 0;
  for (i = 0; i < 100; i++) {
    for (j = 0; j < 10000; j++) {
      index[offset++] = j * 100 + i;
    }
  }

  tic();
  for (i = 999999; i >= 0; i--) {
    free(block[index[i]]);
  }
  printf("%lld\n", toc());
#endif

#if 0
  int (*test)[10];
  int test2[10];
  test = &test2;
  printf("%p\n", test);
  printf("%p\n", *test);
  printf("%p\n", test2);
#endif

#if 0
  int a[100];
  int i = 0;
  i[a] = 0;
  printf("%d\n", i[a]);
#endif

#if 0
  free(NULL);
#endif

#if 0
  int i;
  for (i = 0; i < 10; i++) {
    int n = 0;
    n++;
    printf("%d\n", n);
  }
#endif

#if 0
  printf("%g\n", round(0.5));
  printf("%g\n", floor(0.5));
#endif

#if 0
  if (fcmp("../data/test2.tb", "../data/test3.tb")) {
    printf("The files are different.\n");
  } else {
    printf("The files are the same.\n");
  }
#endif

#if 0
  char str[65];
  printf("%s\n", double_binstr(0.0, str));
#endif

#if 0
  printf("%d\n", MEDIAN3(0, 1, 2));
  printf("%d\n", MEDIAN3(0, 2, 1));
  printf("%d\n", MEDIAN3(1, 0, 2));
  printf("%d\n", MEDIAN3(2, 1, 0));
  printf("%d\n", MEDIAN3(1, 2, 0));
  printf("%d\n", MEDIAN3(2, 0, 1));
  printf("%d\n", MEDIAN3(0, 1, 1));
  printf("%d\n", MEDIAN3(1, 0, 1));
  printf("%d\n", MEDIAN3(1, 1, 0));
  printf("%d\n", MEDIAN3(2, 1, 1));
  printf("%d\n", MEDIAN3(1, 2, 1));
  printf("%d\n", MEDIAN3(1, 1, 2));
#endif

#if 0
  unsigned char ch1 = 255;
  char ch = ch1;
  printf("%x\n", ch);

  ch1 = -1;
  printf("%x\n", ch1);
#endif

#if 0
  //size_t fs = fsize("../data/spikes/pm002.nev");
  FILE *fp = fopen("../data/spikes/pm002.nev", "r");
  /*
  fseek(fp, fs - 104*5, SEEK_SET);
  uint32_t tmp1;
  uint16_t tmp2;
  uint8_t tmp3;
  fread(&tmp1, 4, 1, fp);
  fread(&tmp2, 2, 1, fp);
  fread(&tmp3, 1, 1, fp);
  printf("%u, %u, %u\n", tmp1, tmp2, tmp3);
  */

  uint32_t head_size;
  fseek(fp, 12, SEEK_SET);
  fread(&head_size, 4, 1, fp);
  fseek(fp, head_size, SEEK_SET);
  int packetcnt = 4325366;
  char *buffer = (char*) malloc(packetcnt * 104);
  fread(buffer, 1, packetcnt * 104, fp);
  int offset = 0;
  uint32_t tmp1;
  uint16_t tmp2;
  uint8_t tmp3;

  int i;
  for (i = 0; i < packetcnt; i++) {
    tmp1 = *((uint32_t*) (buffer + offset));
    offset += 4;
    tmp2 = *((uint16_t*) (buffer + offset));
    offset += 2;
    tmp3 = *((uint8_t*) (buffer + offset));
    if (i > 4325360) {
      printf("%u %u %u\n", tmp1, tmp2, tmp3);
    }
    offset += 1;
    offset += 97;
  }

  fclose(fp);
#endif

#if 0
  File_List *list = New_File_List();

  //File_List_Load_Dir("/home/zhaot/Work/neutube/neurolabi/data/test", "tif", list);
  File_List_Load_Dir("../data", "tif", list);

  File_List_Sort_By_Number(list);

  Print_File_List(list);
#endif

#if 0
  int start, end;
  dir_fnum_pair("../data/ting_example_stack/superpixel_maps", ".*\\.png",
      &start, &end);
  printf("%d, %d\n", start, end);
#endif

#if 0
  color_t color;
  uint8_t overflow = Value_To_Color(0, color);

  printf("%u %u %u %u\n", overflow, color[2], color[1], color[0]);
#endif

#if 0
  printf("%*sD\n", 3, "t");
#endif

#if 1
#endif

  return 0;
}


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

#include <gsl/gsl_errno.h>

#include "ccl.h"

/* ------- ROUTINE: ccl_linear spacing ------
INPUTS: [xmin,xmax] of the interval to be divided in N bins
OUTPUT: bin edges in range [xmin,xmax]
*/

double * ccl_linear_spacing(double xmin, double xmax, int N)
{
  double dx = (xmax-xmin)/(N -1.);

  double * x = malloc(sizeof(double)*N);
  if (x==NULL) {
    ccl_raise_warning(
      CCL_ERROR_MEMORY,
      "ERROR: Could not allocate memory for linear-spaced array (N=%d)\n", N);
    return x;
  }

  for (int i=0; i<N; i++) {
    x[i] = xmin + dx*i;
  }
  x[0]=xmin; //Make sure roundoff errors don't spoil edges
  x[N-1]=xmax; //Make sure roundoff errors don't spoil edges

  return x;
}

/* ------- ROUTINE: ccl_linlog spacing ------
 * INPUTS: [xminlog,xmax] of the interval to be divided in bins
 *         xmin when linear spacing starts
 *         Nlog number of logarithmically spaced bins
 *         Nlin number of linearly spaced bins
 * OUTPUT: bin edges in range [xminlog,xmax]
 * */

double * ccl_linlog_spacing(double xminlog, double xmin, double xmax, int Nlog, int Nlin)
{
  if (Nlog<2) {
    ccl_raise_warning(
      CCL_ERROR_LINLOGSPACE,
      "ERROR: Cannot make log-spaced array with %d points - need at least 2\n", Nlog);
    return NULL;
  }

  if (!(xminlog>0 && xmin>0)) {
    ccl_raise_warning(
      CCL_ERROR_LINLOGSPACE,
      "ERROR: Cannot make log-spaced array xminlog or xmin non-positive (had %le, %le)\n", xminlog, xmin);
    return NULL;
  }

  if (xminlog>xmin){
    ccl_raise_warning(CCL_ERROR_LINLOGSPACE, "ERROR: xminlog must be smaller as xmin");
    return NULL;
  }

  if (xmin>xmax){
    ccl_raise_warning(CCL_ERROR_LINLOGSPACE, "ERROR: xmin must be smaller as xmax");
    return NULL;
  }

  double * x = malloc(sizeof(double)*(Nlin+Nlog-1));
  if (x==NULL) {
    ccl_raise_warning(
      CCL_ERROR_MEMORY,
      "ERROR: Could not allocate memory for array of size (Nlin+Nlog-1)=%d)\n", (Nlin+Nlog-1));
    return x;
  }

  double dx = (xmax-xmin)/(Nlin -1.);
  double log_xchange = log(xmin);
  double log_xmin = log(xminlog);
  double dlog_x = (log_xchange - log_xmin) /  (Nlog-1.);

  for (int i=0; i<Nlin+Nlog-1; i++) {
    if (i<Nlog)
        x[i] = exp(log_xmin + dlog_x*i);
    if (i>=Nlog)
        x[i] = xmin + dx*(i-Nlog+1);
  }

  x[0]=xminlog; //Make sure roundoff errors don't spoil edges
  x[Nlog-1]=xmin; //Make sure roundoff errors don't spoil edges
  x[Nlin+Nlog-2]=xmax; //Make sure roundoff errors don't spoil edges

  return x;
}

/* ------- ROUTINE: ccl_log spacing ------
INPUTS: [xmin,xmax] of the interval to be divided logarithmically in N bins
TASK: divide an interval in N logarithmic bins
OUTPUT: bin edges in range [xmin,xmax]
*/

double * ccl_log_spacing(double xmin, double xmax, int N)
{
  if (N<2) {
    ccl_raise_warning(
      CCL_ERROR_LOGSPACE,
      "ERROR: Cannot make log-spaced array with %d points - need at least 2\n", N);
    return NULL;
  }

  if (!(xmin>0 && xmax>0)) {
    ccl_raise_warning(
      CCL_ERROR_LOGSPACE,
      "ERROR: Cannot make log-spaced array xmax or xmax non-positive (had %le, %le)\n", xmin, xmax);
    return NULL;
  }

  double log_xmax = log(xmax);
  double log_xmin = log(xmin);
  double dlog_x = (log_xmax - log_xmin) /  (N-1.);

  double * x = malloc(sizeof(double)*N);
  if (x==NULL) {
    ccl_raise_warning(
      CCL_ERROR_MEMORY,
      "ERROR: Could not allocate memory for log-spaced array (N=%d)\n", N);
    return x;
  }

  double xratio = exp(dlog_x);
  x[0] = xmin; //Make sure roundoff errors don't spoil edges
  for (int i=1; i<N-1; i++) {
    x[i] = x[i-1] * xratio;
  }
  x[N-1]=xmax; //Make sure roundoff errors don't spoil edges

  return x;
}

//Spline creator
//n     -> number of points
//x     -> x-axis
//y     -> f(x)-axis
//y0,yf -> values of f(x) to use beyond the interpolation range
SplPar *ccl_spline_init(int n,double *x,double *y,double y0,double yf)
{
  SplPar *spl=malloc(sizeof(SplPar));
  if(spl==NULL)
    return NULL;

  spl->intacc=gsl_interp_accel_alloc();
  spl->spline=gsl_spline_alloc(gsl_interp_cspline,n);
  int parstatus=gsl_spline_init(spl->spline,x,y,n);
  if(parstatus) {
    gsl_interp_accel_free(spl->intacc);
    gsl_spline_free(spl->spline);
    return NULL;
  }

  spl->x0=x[0];
  spl->xf=x[n-1];
  spl->y0=y0;
  spl->yf=yf;

  return spl;
}

//Evaluates spline at x checking for bound errors
double ccl_spline_eval(double x,SplPar *spl)
{
  if(x<=spl->x0)
    return spl->y0;
  else if(x>=spl->xf)
    return spl->yf;
  else {
    double y;
    int stat=gsl_spline_eval_e(spl->spline,x,spl->intacc,&y);
    if (stat!=GSL_SUCCESS) {
      ccl_raise_gsl_warning(stat, "ccl_utils.c: ccl_splin_eval():");
      return NAN;
    }
    return y;
  }
}

#define CCL_GAMMA1 2.6789385347077476336556 //Gamma(1/3)
#define CCL_GAMMA2 1.3541179394264004169452 //Gamma(2/3)
#define CCL_ROOTPI12 21.269446210866192327578 //12*sqrt(pi)
double ccl_j_bessel(int l,double x)
{
  double jl;
  double ax=fabs(x);
  double ax2=x*x;
  if(l<0) {
    fprintf(stderr,"CosmoMas: l>0 for Bessel function");
    exit(1);
  }

  if(l<7) {
    if(l==0) {
      if(ax<0.1) jl=1-ax2*(1-ax2/20.)/6.;
      else jl=sin(x)/x;
    }
    else if(l==1) {
      if(ax<0.2) jl=ax*(1-ax2*(1-ax2/28)/10)/3;
      else jl=(sin(x)/ax-cos(x))/ax;
    }
    else if(l==2) {
      if(ax<0.3) jl=ax2*(1-ax2*(1-ax2/36)/14)/15;
      else jl=(-3*cos(x)/ax-sin(x)*(1-3/ax2))/ax;
    }
    else if(l==3) {
      if(ax<0.4)
	jl=ax*ax2*(1-ax2*(1-ax2/44)/18)/105;
      else
	jl=(cos(x)*(1-15/ax2)-sin(x)*(6-15/ax2)/ax)/ax;
    }
    else if(l==4) {
      if(ax<0.6)
	jl=ax2*ax2*(1-ax2*(1-ax2/52)/22)/945;
      else
	jl=(sin(x)*(1-(45-105/ax2)/ax2)+cos(x)*(10-105/ax2)/ax)/ax;
    }
    else if(l==5) {
      if(ax<1.0)
	jl=ax2*ax2*ax*(1-ax2*(1-ax2/60)/26)/10395;
      else {
	jl=(sin(x)*(15-(420-945/ax2)/ax2)/ax-
	    cos(x)*(1-(105-945/ax2)/ax2))/ax;
      }
    }
    else {
      if(ax<1.0)
	jl=ax2*ax2*ax2*(1-ax2*(1-ax2/68)/30)/135135;
      else {
	jl=(sin(x)*(-1+(210-(4725-10395/ax2)/ax2)/ax2)+
	    cos(x)*(-21+(1260-10395/ax2)/ax2)/ax)/ax;
      }
    }
  }
  else {
    double nu=l+0.5;
    double nu2=nu*nu;

    if(ax<1.0E-40) jl=0;
    else if((ax2/l)<0.5) {
      jl=(exp(l*log(ax/nu)-M_LN2+nu*(1-M_LN2)-(1-(1-3.5/nu2)/(30*nu2))/(12*nu))/nu)*
	(1-ax2/(4*nu+4)*(1-ax2/(8*nu+16)*(1-ax2/(12*nu+36))));
    }
    else if((l*l/ax)<0.5) {
      double beta=ax-0.5*M_PI*(l+1);
      jl=(cos(beta)*(1-(nu2-0.25)*(nu2-2.25)/(8*ax2)*(1-(nu2-6.25)*(nu2-12.25)/(48*ax2)))-
	  sin(beta)*(nu2-0.25)/(2*ax)*(1-(nu2-2.25)*(nu2-6.25)/(24*ax2)*
				       (1-(nu2-12.25)*(nu2-20.25)/(80*ax2))))/ax;
    }
    else {
      double l3=pow(nu,0.325);
      if(ax<nu-1.31*l3) {
	double cosb=nu/ax;
	double sx=sqrt(nu2-ax2);
	double cotb=nu/sx;
	double secb=ax/nu;
	double beta=log(cosb+sx/ax);
	double cot3b=cotb*cotb*cotb;
	double cot6b=cot3b*cot3b;
	double sec2b=secb*secb;
	double expterm=((2+3*sec2b)*cot3b/24
			-((4+sec2b)*sec2b*cot6b/16
			  +((16-(1512+(3654+375*sec2b)*sec2b)*sec2b)*cot3b/5760
			    +(32+(288+(232+13*sec2b)*sec2b)*sec2b)*sec2b*cot6b/(128*nu))*
			  cot6b/nu)/nu)/nu;
	jl=sqrt(cotb*cosb)/(2*nu)*exp(-nu*beta+nu/cotb-expterm);
      }
      else if(ax>nu+1.48*l3) {
	double cosb=nu/ax;
	double sx=sqrt(ax2-nu2);
	double cotb=nu/sx;
	double secb=ax/nu;
	double beta=acos(cosb);
	double cot3b=cotb*cotb*cotb;
	double cot6b=cot3b*cot3b;
	double sec2b=secb*secb;
	double trigarg=nu/cotb-nu*beta-0.25*M_PI-
	  ((2+3*sec2b)*cot3b/24+(16-(1512+(3654+375*sec2b)*sec2b)*sec2b)*
	   cot3b*cot6b/(5760*nu2))/nu;
	double expterm=((4+sec2b)*sec2b*cot6b/16-
			(32+(288+(232+13*sec2b)*sec2b)*sec2b)*
			sec2b*cot6b*cot6b/(128*nu2))/nu2;
	jl=sqrt(cotb*cosb)/nu*exp(-expterm)*cos(trigarg);
      }
      else {
	double beta=ax-nu;
	double beta2=beta*beta;
	double sx=6/ax;
	double sx2=sx*sx;
	double secb=pow(sx,0.3333333333333333333333);
	double sec2b=secb*secb;

	jl=(CCL_GAMMA1*secb+beta*CCL_GAMMA2*sec2b
	    -(beta2/18-1.0/45.0)*beta*sx*secb*CCL_GAMMA1
	    -((beta2-1)*beta2/36+1.0/420.0)*sx*sec2b*CCL_GAMMA2
	    +(((beta2/1620-7.0/3240.0)*beta2+1.0/648.0)*beta2-1.0/8100.0)*sx2*secb*CCL_GAMMA1
	    +(((beta2/4536-1.0/810.0)*beta2+19.0/11340.0)*beta2-13.0/28350.0)*beta*sx2*sec2b*CCL_GAMMA2
	    -((((beta2/349920-1.0/29160.0)*beta2+71.0/583200.0)*beta2-121.0/874800.0)*
	      beta2+7939.0/224532000.0)*beta*sx2*sx*secb*CCL_GAMMA1)*sqrt(sx)/CCL_ROOTPI12;
      }
    }
  }
  if((x<0)&&(l%2!=0)) jl=-jl;

  return jl;
}

//Spline destructor
void ccl_spline_free(SplPar *spl)
{
  gsl_spline_free(spl->spline);
  gsl_interp_accel_free(spl->intacc);
  free(spl);
}

max_stars_repo_name: Russell-Jones-OxPhys/CCL

max_issues_repo_name: Russell-Jones-OxPhys/CCL

max_forks_repo_name: Russell-Jones-OxPhys/CCL


text:
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <err.h>

#include <gsl/gsl_rng.h>
#include <tskit/tables.h>

#define check_tsk_error(val) if (val < 0) {\
    errx(EXIT_FAILURE, "line %d: %s", __LINE__, tsk_strerror(val));\
}

void
simulate(tsk_table_collection_t *tables, int N, int T, int simplify_interval, gsl_rng *rng)
{
    tsk_id_t *buffer, *parents, *children, child, left_parent, right_parent;
    double breakpoint;
    int ret, j, t, b;

    assert(simplify_interval != 0); // leads to division by zero
    buffer = malloc(2 * N * sizeof(tsk_id_t));
    if (buffer == NULL) {
        errx(EXIT_FAILURE, "Out of memory");
    }
    tables->sequence_length = 1.0;
    parents = buffer;
    for (j = 0; j < N; j++) {
        parents[j] = tsk_node_table_add_row(&tables->nodes, 0, T,
                TSK_NULL, TSK_NULL, NULL, 0);
        check_tsk_error(parents[j]);
    }
    b = 0;
    for (t = T - 1; t >= 0; t--) {
        /* Alternate between using the first and last N values in the buffer */
        parents = buffer + (b * N);
        b = (b + 1) % 2;
        children = buffer + (b * N);
        for (j = 0; j < N; j++) {
            child = tsk_node_table_add_row(&tables->nodes, 0, t,
                    TSK_NULL, TSK_NULL, NULL, 0);
            check_tsk_error(child);
            left_parent = parents[gsl_rng_uniform_int(rng, N)];
            right_parent = parents[gsl_rng_uniform_int(rng, N)];
            do {
                breakpoint = gsl_rng_uniform(rng);
            } while (breakpoint == 0); /* tiny proba of breakpoint being 0 */
            ret = tsk_edge_table_add_row(&tables->edges, 0, breakpoint, left_parent, child);
            check_tsk_error(ret);
            ret = tsk_edge_table_add_row(&tables->edges, breakpoint, 1, right_parent, child);
            check_tsk_error(ret);
            children[j] = child;
        }
        if (t % simplify_interval == 0) {
            printf("Simplify at generation %d: (%d nodes %d edges)", t,
                    tables->nodes.num_rows, tables->edges.num_rows);
            /* Note: Edges must be sorted for simplify to work, and we use a brute force
             * approach of sorting each time here for simplicity. This is inefficient. */
            ret = tsk_table_collection_sort(tables, NULL, 0);
            check_tsk_error(ret);
            ret = tsk_table_collection_simplify(tables, children, N, 0, NULL);
            check_tsk_error(ret);
            printf(" -> (%d nodes %d edges)\n", tables->nodes.num_rows,
                    tables->edges.num_rows);
            for (j = 0; j < N; j++) {
                children[j] = j;
            }
        }
    }
    free(buffer);
}

int
main(int argc, char **argv)
{
    int ret;
    tsk_table_collection_t tables;
    gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

    if (argc != 5) {
        errx(EXIT_FAILURE, "usage: N T simplify-interval output-file");
    }
    ret = tsk_table_collection_init(&tables, 0);
    check_tsk_error(ret);
    simulate(&tables, atoi(argv[1]), atoi(argv[2]), atoi(argv[3]), rng);
    ret = tsk_table_collection_dump(&tables, argv[4], 0);
    check_tsk_error(ret);

    tsk_table_collection_free(&tables);
    gsl_rng_free(rng);
    return 0;
}


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

#pragma once

#include <gsl.h>
#include <nlohmann/json_fwd.hpp>
#include <cstdint>
#include <string>

/**
 * The Event class represents the information needed for a single
 * audit event entry.
 */
class Event {
public:
    Event() = delete;

    /**
     * Construct and initialize a new Event structure based off the
     * provided JSON. See ../README.md for information about the
     * layout of the JSON element.
     *
     * @param entry
     * @throws std::runtime_error for errors accessing the expected
     *                            elements
     */
    explicit Event(const nlohmann::json& json);

    /// The identifier for this entry
    uint32_t id;
    /// The name of the entry
    std::string name;
    /// The full description of the entry
    std::string description;
    /// Set to true if this entry should be handled synchronously
    bool sync;
    /// Set to true if this entry is enabled (or should be dropped)
    bool enabled;
    /// Set to true if the user may enable filtering for the enry
    bool filtering_permitted;
    /// The textual representation of the JSON describing mandatory
    /// fields in the event (NOTE: this is currently not enforced
    /// by the audit daemon)
    std::string mandatory_fields;
    /// The textual representation of the JSON describing the optional
    /// fields in the event (NOTE: this is currently not enforced
    /// by the audit daemon)
    std::string optional_fields;
};


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

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_blas.h>

int
gsl_multifit_gradient (const gsl_matrix * J, const gsl_vector * f,
                       gsl_vector * g)
{
  int status = gsl_blas_dgemv (CblasTrans, 1.0, J, f, 0.0, g);
  return status;
}


text:
/* examples/C/ssmfe/hermitian.c - Example code for SPRAL_SSMFE package */
/* Hermitian operator example */
#include "spral.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <cblas.h>

/* central differences for i d/dx */
void apply_idx(int n, int m, const double complex *x_ptr, double complex *y_ptr) {
   const double complex (*x)[n] = (const double complex (*)[n]) x_ptr;
   double complex (*y)[n] = (double complex (*)[n]) y_ptr;
   for(int j=0; j<m; j++) {
      for(int i=0; i<n; i++) {
         int il = (i==0)   ? n-1 : i-1;
         int ir = (i==n-1) ? 0   : i+1;
         y[j][i] = _Complex_I * (x[j][ir] - x[j][il]);
      }
   }
}

/* main routine */
int main(void) {
   const int n   = 80;                 /* problem size */
   const int nep = 5;                  /* eigenpairs wanted */

   double lambda[nep];                 /* eigenvalues */
   double complex X[nep][n];           /* eigenvectors */
   struct spral_ssmfe_rciz rci;        /* reverse communication data */
   struct spral_ssmfe_options options; /* options */
   void *keep;                         /* private data */
   struct spral_ssmfe_inform inform;   /* information */

   /* Initialize options to default values */
   spral_ssmfe_default_options(&options);

   rci.job = 0; keep = NULL;
   while(true) { /* reverse communication loop */
      spral_ssmfe_standard_double_complex(&rci, nep, nep, lambda, n,
         &X[0][0], n, &keep, &options, &inform);
      switch ( rci.job ) {
      case 1:
         apply_idx(n, rci.nx, rci.x, rci.y);
         break;
      case 2:
         // No preconditioning
         break;
      default:
         goto finished;
      }
   }
finished:
   printf("%d eigenpairs converged in %d iterations\n", inform.left, inform.iteration);
   for(int i=0; i<inform.left; i++)
      printf(" lambda[%1d] = %13.7e\n", i, lambda[i]);
   spral_ssmfe_free_double_complex(&keep, &inform);

   /* Success */
   return 0;
}

max_stars_repo_name: mjacobse/spral

max_issues_repo_name: mjacobse/spral

max_forks_repo_name: mjacobse/spral


text:
#pragma once

#include <vector>

#include <gsl-lite/gsl-lite.hpp>

#include <Eigen/Sparse>

#include <thrustshift/CSR.h>
#include <thrustshift/container-conversion.h>
#include <thrustshift/managed-vector.h>
#include <thrustshift/memory-resource.h>

namespace thrustshift {

namespace eigen {

//! Return a CSR matrix with the default memory resource of the CSR class
template <class EigenSparseMatrix>
auto sparse_mtx2csr(EigenSparseMatrix&& m_) {

	using DataType =
	    typename std::remove_reference<EigenSparseMatrix>::type::value_type;
	using StorageIndex =
	    typename std::remove_reference<EigenSparseMatrix>::type::StorageIndex;
	// Create copy because we might modify the container with `makeCompressed`
	Eigen::SparseMatrix<DataType, Eigen::RowMajor, StorageIndex> m = m_;
	m.makeCompressed();
	auto m_data = m.data();
	auto nnz = m.nonZeros();
	auto rows = m.rows();

	const DataType* A = &m_data.value(0);
	const StorageIndex* IA = m.outerIndexPtr();
	const StorageIndex* JA = &m_data.index(0);

	managed_vector<DataType> seq_A(A, A + nnz);
	managed_vector<StorageIndex> seq_JA(JA, JA + nnz);
	managed_vector<StorageIndex> seq_IA(IA, IA + rows + 1);

	return CSR<DataType, StorageIndex>(
	    seq_A, seq_JA, seq_IA, gsl_lite::narrow<size_t>(m.cols()));
}

// Forward declaration to avoid co-dependent headers
template <class EigenSparseMatrix, class COO_C>
EigenSparseMatrix coo2sparse_mtx(COO_C&& coo);

template <class EigenSparseMatrix, class CSR_C>
EigenSparseMatrix csr2sparse_mtx(CSR_C&& csr) {

	using DataType =
	    typename std::remove_reference<EigenSparseMatrix>::type::value_type;
	using StorageIndex =
	    typename std::remove_reference<EigenSparseMatrix>::type::StorageIndex;

	return coo2sparse_mtx<EigenSparseMatrix>(
	    csr2coo<thrustshift::COO<DataType, StorageIndex>>(
	        std::forward<CSR_C>(csr), pmr::default_resource));
}

} // namespace eigen

} // namespace thrustshift

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
#ifdef USE_MKL
#include <mkl.h>
#else
#include <cblas.h>
#endif

// contains all calls to C-interface BLAS functions
// e.g. http://www.netlib.org/lapack/explore-html/d1/dff/cblas__example1_8c_source.html

// returns the dot product of x and y
// http://www.netlib.org/lapack/explore-html/d5/df6/ddot_8f_source.html
double specex_dot(int n, const double *x, const double *y){
  return cblas_ddot(n, x, 1, y, 1);
}

// y += alpha*x
// http://www.netlib.org/lapack/explore-html/d9/dcd/daxpy_8f_source.html
void specex_axpy(int n, const double *alpha, const double *x, const double *y){
  cblas_daxpy(n, *alpha, x, 1, y, 1);
}

// A += alpha*x*x**T, where A is a symmetric matrx (only lower half is filled)
// http://www.netlib.org/lapack/explore-html/d3/d60/dsyr_8f_source.html
void specex_syr(int n, const double *alpha, const double *x, const double *A){
  cblas_dsyr(CblasColMajor, CblasLower, n, *alpha, x, 1, A, n);
}

// C = alpha*A**T + beta*C, where A is a symmetric matrx
// http://www.netlib.org/lapack/explore-html/dc/d05/dsyrk_8f_source.html
void specex_syrk(int n, int k, const double *alpha, const double *A, const double *beta,
		 const double *C){
  cblas_dsyrk(CblasColMajor, CblasLower, CblasNoTrans, n, k, *alpha, A, n, *beta, C, n);
}

// y = alpha*A*x + beta*y
// http://www.netlib.org/lapack/explore-html/dc/da8/dgemv_8f_source.html
void specex_gemv(int m, int n, const double *alpha, const double *A, const double *x,
		 const double *beta, const double *y){
  cblas_dgemv(CblasColMajor, CblasNoTrans, m, n, *alpha, A, m, x, 1, *beta, y, 1);

}

// C = alpha*A*B + beta*C
// http://www.netlib.org/lapack/explore-html/d7/d2b/dgemm_8f_source.html
void specex_gemm(int m, int n, int k, const double *alpha, const double *A, const double *B,
		 const double *beta, const double *C){
  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, *alpha, A, m, B, k, *beta, C, m);
}

max_stars_repo_name: tskisner/specex

max_issues_repo_name: tskisner/specex

max_forks_repo_name: tskisner/specex


text:
/* randist/logistic.c
 *
 * 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.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* The logistic distribution has the form,

   p(x) dx = (1/a) exp(-x/a) / (1 + exp(-x/a))^2 dx

   for -infty < x < infty */

double
gsl_ran_logistic (const gsl_rng * r, const double a)
{
  double x, z;

  do
    {
      x = gsl_rng_uniform_pos (r);
    }
  while (x == 1);

  z = log (x / (1 - x));

  return a * z;
}

double
gsl_ran_logistic_pdf (const double x, const double a)
{
  double u = exp (-fabs(x)/a);
  double p = u / (fabs(a) * (1 + u) * (1 + u));
  return p;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>

#include "tests.h"

void
test_hpr2 (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   int N = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float Ap[] = { 0.159f, -0.13f };
   float X[] = { 0.854f, 0.851f };
   int incX = -1;
   float Y[] = { 0.526f, -0.267f };
   int incY = -1;
   float Ap_expected[] = { -0.284974f, 0.0f };
   cblas_chpr2(order, uplo, N, alpha, X, incX, Y, incY, Ap);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], flteps, "chpr2(case 1458) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], flteps, "chpr2(case 1458) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float Ap[] = { 0.159f, -0.13f };
   float X[] = { 0.854f, 0.851f };
   int incX = -1;
   float Y[] = { 0.526f, -0.267f };
   int incY = -1;
   float Ap_expected[] = { -0.284974f, 0.0f };
   cblas_chpr2(order, uplo, N, alpha, X, incX, Y, incY, Ap);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], flteps, "chpr2(case 1459) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], flteps, "chpr2(case 1459) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float Ap[] = { 0.159f, -0.13f };
   float X[] = { 0.854f, 0.851f };
   int incX = -1;
   float Y[] = { 0.526f, -0.267f };
   int incY = -1;
   float Ap_expected[] = { -0.284974f, 0.0f };
   cblas_chpr2(order, uplo, N, alpha, X, incX, Y, incY, Ap);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], flteps, "chpr2(case 1460) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], flteps, "chpr2(case 1460) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float Ap[] = { 0.159f, -0.13f };
   float X[] = { 0.854f, 0.851f };
   int incX = -1;
   float Y[] = { 0.526f, -0.267f };
   int incY = -1;
   float Ap_expected[] = { -0.284974f, 0.0f };
   cblas_chpr2(order, uplo, N, alpha, X, incX, Y, incY, Ap);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], flteps, "chpr2(case 1461) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], flteps, "chpr2(case 1461) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int N = 1;
   double alpha[2] = {-0.3, 0.1};
   double Ap[] = { 0.772, 0.997 };
   double X[] = { -0.173, -0.839 };
   int incX = -1;
   double Y[] = { 0.941, -0.422 };
   int incY = -1;
   double Ap_expected[] = { 0.829742, 0.0 };
   cblas_zhpr2(order, uplo, N, alpha, X, incX, Y, incY, Ap);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], dbleps, "zhpr2(case 1462) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], dbleps, "zhpr2(case 1462) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 1;
   double alpha[2] = {-0.3, 0.1};
   double Ap[] = { 0.772, 0.997 };
   double X[] = { -0.173, -0.839 };
   int incX = -1;
   double Y[] = { 0.941, -0.422 };
   int incY = -1;
   double Ap_expected[] = { 0.829742, 0.0 };
   cblas_zhpr2(order, uplo, N, alpha, X, incX, Y, incY, Ap);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], dbleps, "zhpr2(case 1463) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], dbleps, "zhpr2(case 1463) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 1;
   double alpha[2] = {-0.3, 0.1};
   double Ap[] = { 0.772, 0.997 };
   double X[] = { -0.173, -0.839 };
   int incX = -1;
   double Y[] = { 0.941, -0.422 };
   int incY = -1;
   double Ap_expected[] = { 0.829742, 0.0 };
   cblas_zhpr2(order, uplo, N, alpha, X, incX, Y, incY, Ap);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], dbleps, "zhpr2(case 1464) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], dbleps, "zhpr2(case 1464) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 1;
   double alpha[2] = {-0.3, 0.1};
   double Ap[] = { 0.772, 0.997 };
   double X[] = { -0.173, -0.839 };
   int incX = -1;
   double Y[] = { 0.941, -0.422 };
   int incY = -1;
   double Ap_expected[] = { 0.829742, 0.0 };
   cblas_zhpr2(order, uplo, N, alpha, X, incX, Y, incY, Ap);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(Ap[2*i], Ap_expected[2*i], dbleps, "zhpr2(case 1465) real");
       gsl_test_rel(Ap[2*i+1], Ap_expected[2*i+1], dbleps, "zhpr2(case 1465) imag");
     };
   };
  };

}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc