text:
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

/*
 Contains collection of low-level routines for handling lists,
 one-hot matrices, evaluation metrics (f1-scores), and other.

 Dimitris Berberidis
 University of Minnesota 2017-2018
*/

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <inttypes.h>
#include <cblas.h>
#include <sys/stat.h>
#include <stdbool.h>

#include "my_utils.h"
#include "my_defs.h"

void assert_all_nodes_present(csr_graph graph, const sz_long* seed_indices, sz_med num_seeds){

	for(sz_med i=0;i<num_seeds;i++){
		if(seed_indices[i]>graph.num_nodes){
			printf("ERROR: Seed node index does not appear in edgelist (probably an isolated node)\n");
		        exit(EXIT_FAILURE);
		}
	}

}

sz_long rand_lim(sz_long limit) {
	/* return a random number between 0 and limit inclusive.
	 */
	int divisor = RAND_MAX/((int)limit+1);
	int retval;
	do {
		retval = rand() / divisor;
	} while (retval > limit);
	return retval;
}

//Draw random samples with replacement from 0 to N-1
void random_sample( sz_long* seeds, abstract_labels labels, abstract_labels all_labels, sz_med num_seeds, sz_long N){

	sz_long temp;
	sz_short flag;

	//Draw seed indexes
	seeds[0]=rand_lim(N-1);
	for(sz_med i=1;i<num_seeds;i++){
		do{
			temp=rand_lim(N-1);
			flag=0;
			for(sz_med j=0;j<i;j++){
				if(temp==seeds[j])
					flag=1;
			}
		}while(flag==1);
		seeds[i]=temp;
	}

	//Draw corresponding labels
	if(all_labels.is_multilabel){
		for(sz_med i=0;i<num_seeds;i++){
			for(sz_short j=0; j<all_labels.mlabel.num_class;j++)
				labels.mlabel.bin[j][i] = all_labels.mlabel.bin[j][seeds[i]];
		}
	}else{

		for(sz_med i=0;i<num_seeds;i++){
			labels.mclass[i]=all_labels.mclass[seeds[i]];
			//printf("%d\n",seeds[i]);
		}
	}

	// +1 seed indexes !
	for(sz_med i=0;i<num_seeds;i++){seeds[i]+=1;}

}

//My comparator for two collumn array. Sorts second col according to first
int compare ( const void *pa, const void *pb )
{
	const sz_long *a = *(const sz_long **)pa;
	const sz_long *b = *(const sz_long **)pb;
	if(a[0] == b[0])
		return a[1] - b[1];
	else
		return a[0] - b[0];
}

//My comparator for a struct with double and index
int compare2 ( const void *pa, const void *pb )
{
  val_and_ind *a1 = (val_and_ind*)pa;
  val_and_ind *a2 = (val_and_ind*)pb;
  if((*a1).val>(*a2).val)return -1;
  else if((*a1).val<(*a2).val)return 1;
  else return 0;
}

//frobenious norm of double-valued square matrix
double frob_norm(double* A, sz_med dim){
	double norm=0.0f;

	for(sz_med i=0;i<dim*dim;i++){
		norm+=pow(A[i],2.0f);
	}

	return sqrt(norm);
}

//mean of double array
double mean(double* a, int len){
	double sum=0;
	for(int i=0;i<len;i++) sum+=a[i];
	return sum/(double)len;
}

//Simple linear system solver using the LU decomposition
/* INPUT: A,P filled in LUPDecompose; b - rhs vector; N - dimension
 * OUTPUT: x - solution vector of A*x=b
 */
void LUPSolve(double **A, int *P, double *b, int N, double *x){

    for (int i = 0; i < N; i++) {
        x[i] = b[P[i]];

        for (int k = 0; k < i; k++)
            x[i] -= A[i][k] * x[k];
    }

    for (int i = N - 1; i >= 0; i--) {
        for (int k = i + 1; k < N; k++)
            x[i] -= A[i][k] * x[k];

        x[i] = x[i] / A[i][i];
    }
}

//Simple LU decomposition routine
/* INPUT: A - array of pointers to rows of a square matrix having dimension N
 *        Tol - small tolerance number to detect failure when the matrix is near degenerate
 * OUTPUT: Matrix A is changed, it contains both matrices L-E and U as A=(L-E)+U such that P*A=L*U.
 *        The permutation matrix is not stored as a matrix, but in an integer vector P of size N+1
 *        containing column indexes where the permutation matrix has "1". The last element P[N]=S+N,
 *        where S is the number of row exchanges needed for determinant computation, det(P)=(-1)^S
 */
int LUPDecompose(double **A, int N, double Tol, int *P) {

    int i, j, k, imax;
    double maxA, *ptr, absA;

    for (i = 0; i <= N; i++)
        P[i] = i; //Unit permutation matrix, P[N] initialized with N

    for (i = 0; i < N; i++) {
        maxA = 0.0;
        imax = i;

        for (k = i; k < N; k++)
            if ((absA = fabs(A[k][i])) > maxA) {
                maxA = absA;
                imax = k;
            }

        if (maxA < Tol) return 0; //failure, matrix is degenerate

        if (imax != i) {
            //pivoting P
            j = P[i];
            P[i] = P[imax];
            P[imax] = j;

            //pivoting rows of A
            ptr = A[i];
            A[i] = A[imax];
            A[imax] = ptr;

            //counting pivots starting from N (for determinant)
            P[N]++;
        }

        for (j = i + 1; j < N; j++) {
            A[j][i] /= A[i][i];

            for (k = i + 1; k < N; k++)
                A[j][k] -= A[j][i] * A[i][k];
        }
    }

    return 1;  //decomposition done
}

//Interface for CBLAS matrix vector product
void matvec(double*y, double* A, double* x, sz_med M, sz_med N ){

	for(int i=0;i<M;i++){y[i]=0.0f;}

	cblas_dgemv( CblasRowMajor , CblasNoTrans , (int)M , (int)N, 1.0f, A, (int)M, x, 1, 0.0f, y, 1);

}

//Interface for CBLAS matrix vector product
void matvec_trans(double*y, double* A, double* x, sz_med M, sz_med N ){

	for(sz_med i=0;i<M;i++){y[i]=0.0f;}

	cblas_dgemv( CblasRowMajor , CblasTrans , (int)M , (int)N, 1.0f, A, (int)M, x, 1, 0.0f, y, 1);

}

//Interface for CBLAS trnaspose-matrix vector product
void matvec_trans_long( double* y , double* A, double* x, sz_long N, sz_med p ){

	cblas_dgemv (CblasRowMajor, CblasTrans, (int) p, (int)N , 1.0f, A,  (int) N , x, 1 , 0.0f, y, 1);

}

//Interface for CBLAS mutrix matrix product
void matrix_matrix_product(double*C, double* A, double* B, sz_long m, sz_med k , sz_short n){

	cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, (int)m,
		     (int)n, (int)k, 1.0f, A, (int)k, B, (int)n, 0.0f, C, (int)n);

}

//Project vector onto simplex by alternating projections onto line and positive quadrant
//Operation happens in place
void project_to_simplex( double* x, sz_med N ){
	double sum,a;
	sz_short flag;

	do{
		flag=0;
		sum=0.0f;

		for(sz_med i=0; i<N; i++) sum+=x[i];

		a=(sum - 1.0f)/(double)N;

		for(sz_med i=0; i<N; i++){
			x[i]-=a;
			if(x[i]<= - PROJ_TOL){
				x[i]=0.0f;
				flag=1;}
		}

	}while(flag==1);
}

//find l_max norm between vecors of known length
double max_diff(double* x, double* y , sz_med L){
	double dif;
	double max_dif=0.0f;
	for(sz_med i=0;i<L;i++){
		dif=fabs( x[i] - y[i] );
		max_dif = ( dif > max_dif ) ? dif : max_dif ;
	}
	return max_dif;
}

//Print array (double)
void print_array_1D(double* arr, sz_long N, sz_long M){
	sz_long i,j;
	printf("Array: \n");
	for(i=0;i<N;i++){
		printf("\n");
		for(j=0;j<M;j++) printf("%lf  ",arr[i*M + j]);
	}
	printf("\n");
}

//Pradict labels from largest soft label
void predict_labels( class_t* label_out, double* soft_labels, class_t* class,
	             sz_long graph_size, sz_short num_class ){
	sz_short max_ind;
	double max_val;

	for(sz_long i=0;i<graph_size;i++){
		max_val=-100.0f;
		for(sz_short j=0;j<num_class;j++){
			if(max_val<soft_labels[i*num_class + j]){
				max_val=soft_labels[i*num_class + j];
				max_ind=j;}
		}
		label_out[i]=class[max_ind];
	}
}

//Pradict labels from largest soft label (Soft labels are transposed array)
void predict_labels_type2( class_t* label_out, double* soft_labels, class_t* class,
			   sz_long graph_size, sz_short num_class ){

	sz_short max_ind;
	double max_val;

	for(sz_long j=0;j<graph_size;j++){
		max_val=-100.0f;
		for(sz_short i=0;i<num_class;i++){
			if(max_val<soft_labels[i*graph_size + j]){
				max_val=soft_labels[i*graph_size + j];
				max_ind=i;}
		}
		label_out[j]=class[max_ind];
	}
}

//Print the two collumn array of edges
void print_edge_list(const sz_long** edge_list, sz_long len){
	printf("EDGE LIST: \n");
	for(sz_long i=0;i<len;i++)
		printf("%"PRIu64"  %"PRIu64"\n", (uint64_t) edge_list[i][0], (uint64_t) edge_list[i][1]);
}

//Print array (double)
void print_predictions(class_t* arr, sz_long N){
	printf("Predicted labels: \n");
	for(sz_long i=0;i<N;i++) printf("%"PRIi16"\n", (int16_t) arr[i]);
}

//Check if file is valid
int file_isreg(char *path) {
    struct stat st;

    if (stat(path, &st) < 0)
        return -1;

    return S_ISREG(st.st_mode);
}

//Elementwise subtract b from a and store in c
void my_array_sub(double* c, double* a, double* b, sz_long N ){
	for(sz_long i=0;i<N;i++) c[i] = a[i] -b[i];
}

//Max of array

sz_short max_u8( sz_short* a, sz_long N ){
	sz_short max= 0;
	for(sz_long i=0; i<N ; i++) max = ( a[i]>max ) ? a[i] : max ;
	return max;
}

sz_long max_u64( sz_long* a, sz_long N ){
	sz_long max= 0;
	for(sz_long i=0; i<N ; i++) max = ( a[i]>max ) ? a[i] : max ;
	return max;
}

//Return edge list and count to main

sz_long** give_edge_list( char* file_name, sz_long* count ){

	sz_long** buffer= (sz_long **)malloc(EDGE_BUFF_SIZE * sizeof(sz_long *));

	for(sz_long i=0;i<EDGE_BUFF_SIZE;i++)
		buffer[i]= (sz_long*) malloc(2*sizeof(sz_long));

	FILE* file= fopen(file_name, "r");

	if(!file) printf("ERROR: Cannot open graph file");

	// Read adjacency into buffer into buffer and return length count=edges
	*count= read_adjacency_to_buffer(buffer,file);
	printf("Number of edges: %"PRIu64"\n", (uint64_t) *count);

	//print_edge_list( buffer, *count);

	//Free excess memory
	for(sz_long i=*count+1;i<EDGE_BUFF_SIZE;i++)
	{free(buffer[i]);}
	buffer=realloc(buffer,(*count)*sizeof(sz_long*));

	return buffer;
}

//Read .txt file into buffer
sz_long read_adjacency_to_buffer(sz_long** buffer, FILE* file){
        sz_long count = 0;
	for (; count < EDGE_BUFF_SIZE; ++count)
	{
		int got = fscanf(file, "%"SCNu64"\t%"SCNu64"\n", &buffer[count][0] , &buffer[count][1]);
		if ((got != 2)||( (buffer[count][0]==0) && (buffer[count][1]==0))) break;
		// Stop scanning if wrong number of tokens (maybe end of file) or zero input
	}
	fclose(file);
	return count;
}

//Read class.txt file into buffer (ignore first collumn here)
class_t* read_labels(char* filename, sz_long* label_count){
	sz_long count = 0;
	sz_long* indexes = (sz_long*) malloc(sizeof(sz_long)*CLASS_BUFF_SIZE);
	class_t* buffer = (class_t*) malloc(sizeof(class_t)*CLASS_BUFF_SIZE);
	FILE* file= fopen(filename, "r");

	if(!file) printf("ERROR: Cannot open label file");

	for (; count < CLASS_BUFF_SIZE; ++count)
	{
		int got = fscanf(file, "%"SCNu64"%"SCNu8"", &indexes[count] , &buffer[count]);
		if (got != 2) break; // wrong number of tokens - maybe end of file
	}
	fclose(file);

	*label_count = count;
	buffer = realloc(buffer,sizeof(class_t)*count);
	indexes = realloc(indexes,sizeof(sz_long)*count);

	my_relative_sorting( indexes, buffer, count );

	free(indexes);
	return buffer;
}

//Read seed and label file when in operation mode
sz_long* read_seed_file( char* filename, sz_med* num_seeds, sz_short* num_class, abstract_labels* label_in ){

        //First read file into buffers
	sz_long count = 0;
	sz_long* index_buffer = (sz_long*) malloc(sizeof(sz_long)*CLASS_BUFF_SIZE);
	class_t* label_buffer = (class_t*) malloc(sizeof(class_t)*CLASS_BUFF_SIZE);
	FILE* file= fopen(filename, "r");

	if(!file) printf("ERROR: Cannot open seed file");

	for (; count < CLASS_BUFF_SIZE; ++count)
	{
		int got = fscanf(file, "%"SCNu64"%"SCNd8"", &index_buffer[count] , &label_buffer[count]);
		if (got != 2) break; // wrong number of tokens - maybe end of file
	}
	fclose(file);
	label_buffer = realloc(label_buffer,sizeof(class_t)*count);
	index_buffer = realloc(index_buffer,sizeof(sz_long)*count);

	//prepare input labels and seeds for is_multilabel or multi_class

	sz_long* seed_indices;

	if(label_in->is_multilabel){
		*num_class = max_u8( (sz_short*) label_buffer, count);
		my_relative_sorting( index_buffer, label_buffer, count );
		seed_indices = find_unique_from_sorted( index_buffer, count , num_seeds );

		label_in->mlabel = init_one_hot(*num_class , *num_seeds);
		sz_long j=0;
		for(sz_long i=0;i<count;i++){
			if(i>0) j = (index_buffer[i] == index_buffer[i-1] ) ? j : j+1;
			label_in->mlabel.bin[label_buffer[i]-1][j]=1;
		}

	        free(index_buffer);
		free(label_buffer);
	}else{
		seed_indices = index_buffer;
		label_in->mclass = label_buffer;
		*num_seeds = (sz_med) count;
	}

	return seed_indices;
}

//write predicted labels (or ranking in multilabel case) to output file
void save_predictions(char* filename, abstract_label_output label_out, sz_long len, sz_short num_class){

	FILE* file = fopen(filename, "w");

	if(!file) printf("ERROR: Cannot open outfile");

	if(label_out.is_multilabel){
		val_and_ind line_of_out[num_class];
		for(sz_long i=0; i<len; i++){
			for(sz_short j=0; j<num_class; j++)
				line_of_out[j] = (val_and_ind) {.val = label_out.mlabel[j*len + i], .ind=(int)j};

			qsort( line_of_out, num_class, sizeof(line_of_out[0]), compare2);

			fprintf(file, "%"SCNu64":\t", (uint64_t) i+1 );

			for(sz_short j=0; j<num_class; j++)
				fprintf(file, "%"PRIu16" ", (uint16_t) line_of_out[j].ind +1 );

			fprintf(file, "\n");
		}
	}else{
		for(sz_long i=0; i<len; i++)
			fprintf(file, "%"SCNu64"\t%"SCNd16"\n", (uint64_t) i+1 , (int16_t) label_out.mclass[i]);
	}

	fclose(file);
}

//Read class.txt file into one_hot_matrix (ignore first collumn here)
one_hot_mat read_one_hot_mat(char* filename, sz_long* label_count){
	sz_long count = 0;
	sz_long* indexes = (sz_long*) malloc(sizeof(sz_long)*CLASS_BUFF_SIZE);
	sz_short* buffer = (sz_short*) malloc(sizeof(sz_short)*CLASS_BUFF_SIZE);
	FILE* file = fopen(filename, "r");

	if(!file) printf("ERROR: Cannot open label file");

	for (; count < CLASS_BUFF_SIZE; ++count)
	{
		int got = fscanf(file, "%"SCNu64"\t%"SCNu8"", &indexes[count] , &buffer[count]);
		if (got != 2) break; // wrong number of tokens - maybe end of file
	}
	fclose(file);

	buffer = realloc(buffer,sizeof(sz_short)*count);
	indexes = realloc(indexes,sizeof(sz_long)*count);

	sz_short num_class = max_u8( buffer, count);

	sz_long length = max_u64( indexes, count);

//	printf("LENGTH %"PRIu64"\n",length);

	one_hot_mat all_labels = init_one_hot(num_class , length);

	for(sz_long i=0;i<count;i++) all_labels.bin[buffer[i]-1][indexes[i]-1]=1;

	*label_count = length;

	free(indexes);
	free(buffer);
	return all_labels;
}

// Sort A and ind with respect to indexes in ind
void my_relative_sorting( sz_long* ind, class_t* A, sz_long len ){
	sz_long* temp = (sz_long*)malloc(len*sizeof(sz_long));
	sz_long* sorted_inds = (sz_long*)malloc(len*sizeof(sz_long));

	for(sz_long i=0;i<len;i++){ temp[i] = i ;}

	sz_long** zipped_arrays = zip_arrays(ind,temp,len);

        qsort( zipped_arrays , len, sizeof(zipped_arrays[0]), compare);

        unzip_array(zipped_arrays, temp, sorted_inds, len );

        rearange(sorted_inds, (void*) A, "class_t" , len );

	free(temp);
	free(sorted_inds);

}

//Zip two arrays into a list of length-2 lists
sz_long** zip_arrays(sz_long* A, sz_long* B, sz_long len){
	sz_long** zipped = (sz_long**)malloc(len*sizeof(sz_long*));

	for(sz_long i=0;i<len;i++){
		zipped[i] = (sz_long*)malloc(2*sizeof(sz_long));
                ** (zipped +i) = A[i];
                *(* (zipped +i)+1) = B[i];
	}
	return zipped;
}

//Unzip two arrays into input pointers and destroy zipped array
void unzip_array(sz_long** zipped, sz_long* unzip_1, sz_long* unzip_2, sz_long len ){

	for(sz_long i=0;i<len;i++){
		unzip_1[i] = zipped[i][0];
		unzip_2[i] = zipped[i][1];
		free(*(zipped+i));
	}

	free(zipped);
}

//Rearange elements of array A accortding to given indexes
//Works for any type of array as long as the type is provided
void rearange(sz_long* ind, void* A, char* type , sz_long len ){

	if(strcmp(type,"double")==0){
		double A_temp[len];
		memcpy(A_temp, A, len*sizeof(double));
		for(sz_long i=0;i<len;i++){ *((double*)A + i) = A_temp[ind[i]];}
	}else if(strcmp(type,"sz_long")==0){
		sz_long A_temp[len];
		memcpy(A_temp, A, len*sizeof(sz_long));
		for(sz_long i=0;i<len;i++){ *((sz_long*)A + i) = A_temp[ind[i]];}
	}else if(strcmp(type,"uint32_t")==0){
		uint32_t A_temp[len];
		memcpy(A_temp, A, len*sizeof(uint32_t));
		for(sz_long i=0;i<len;i++){ *((uint32_t*)A + i) = A_temp[ind[i]];}
	}else if(strcmp(type,"sz_med")==0){
		sz_med A_temp[len];
		memcpy(A_temp, A, len*sizeof(sz_med));
		for(sz_long i=0;i<len;i++){ *((sz_med*)A + i) = A_temp[ind[i]];}
	}else if(strcmp(type,"sz_short")==0){
		sz_short A_temp[len];
		memcpy(A_temp, A, len*sizeof(sz_short));
		for(sz_long i=0;i<len;i++){ *((sz_short*)A + i) = A_temp[ind[i]];}
	}else if(strcmp(type,"int32_t")==0){
		int32_t A_temp[len];
		memcpy(A_temp, A, len*sizeof(int32_t));
		for(sz_long i=0;i<len;i++){ *((int32_t*)A + i) = A_temp[ind[i]];}
	}else if(strcmp(type,"int16_t")==0){
		int16_t A_temp[len];
		memcpy(A_temp, A, len*sizeof(int16_t));
		for(sz_long i=0;i<len;i++){ *((int16_t*)A + i) = A_temp[ind[i]];}
	}else if(strcmp(type,"class_t")==0){
		class_t A_temp[len];
		memcpy(A_temp, A, len*sizeof(class_t));
		for(sz_long i=0;i<len;i++){ *((class_t*)A + i) = A_temp[ind[i]];}
	}else if(strcmp(type,"long double")==0){
		long double A_temp[len];
		memcpy(A_temp, A, len*sizeof(long double));
		for(sz_long i=0;i<len;i++){ *((long double*)A + i) = A_temp[ind[i]];}
	}else if(strcmp(type,"int")==0){
		int A_temp[len];
		memcpy(A_temp, A, len*sizeof(int));
		for(sz_long i=0;i<len;i++){ *((int*)A + i) = A_temp[ind[i]];}
	}else{
		printf("Unknown rearange type\n");
		exit(EXIT_FAILURE);
	}

}

//Remove items of given indexes from array (list)
// Return result in NEW  list
sz_long* remove_from_list(const sz_long* list, const sz_long* indexes_to_be_removed,
			   sz_long len, sz_long num_removed ){

	sz_long* new_list = (sz_long*) malloc((len-num_removed)*sizeof(sz_long));

	int mask[len];

	memset(mask, 0, len*sizeof(int));

	for(sz_long i =0; i<num_removed; i++){ mask[indexes_to_be_removed[i]] =1 ;}

	sz_long k=0;
	for(sz_long i =0; i<len; i++){
		if(mask[i]==0){
		    new_list[k++] = list[i];
		}
	}

	return new_list;
}

//Find unique elements among known number of entries in buffer
//Also eturn number of unique entries
sz_short find_unique(class_t* unique_elements, const class_t* buffer,sz_med N){

	sz_short new,end=1;

	unique_elements[0]=buffer[0];
	for(sz_med i=0;i<N;i++){
		new=1;
		for(sz_short j=0;j<end;j++)
		{if(unique_elements[j]==buffer[i])
			new=0;}

		if(new)
			unique_elements[end++]=buffer[i];

	}
	return end;
}

//Find unique entries from list of sorted (unisgned indexes)
sz_long* find_unique_from_sorted( sz_long* sorted_buffer, sz_long len , sz_med* num_unique ){
	sz_long* unique = (sz_long*) malloc(len*sizeof(sz_long));
	*num_unique = 1;

	unique[0] = sorted_buffer[0];
	for(sz_long i=1;i<len;i++){
		if(!(sorted_buffer[i] == sorted_buffer[i-1] )){
			*num_unique += 1;
			unique[*num_unique-1] = sorted_buffer[i];
		}
	}
	unique = realloc(unique,*num_unique*sizeof(sz_long));
	return unique;
}

// Converts list to one hot binary matrix of one_hot_mat type
// Label count is number of indexes with known labels
// label_count<=length
// If label_count<=length, then rows of the one_hot matrix
// without a coresponding labeled index will be [0...0]
one_hot_mat list_to_one_hot( sz_long* ind , class_t* labels, sz_short num_class ,
			     class_t* class ,  sz_long label_count ,sz_long length){

	one_hot_mat one_hot = init_one_hot( num_class , length);

	for(sz_long i=0;i<label_count;i++){
		for(sz_short j=0;j<num_class;j++){
			if( labels[i] == class[j] ){
				one_hot.bin[j][ind[i]]=1;
				break;
			}
		}
	}

	return one_hot;
}

//Allocate one_hot_mat
one_hot_mat init_one_hot(sz_short num_class ,sz_long length){
	one_hot_mat one_hot;
	one_hot.bin = (sz_short**) malloc(num_class*sizeof(sz_short*));

	for(sz_short i=0;i<num_class;i++){
		one_hot.bin[i] = (sz_short*) malloc(length*sizeof(sz_short));
		for(sz_long j=0;j<length;j++){
			one_hot.bin[i][j] =0 ;
		}
	}
	one_hot.num_class = num_class;
	one_hot.length = length;
	return one_hot;
}

// Destroy (free) one_hot matrix
void destroy_one_hot(one_hot_mat one_hot){
	for(sz_short i=0;i<one_hot.num_class;i++) free(one_hot.bin[i]);
	free(one_hot.bin);
}

// Return array with number of non-zero entries in each row of one-hot-matrix
sz_short* return_num_labels_per_node( one_hot_mat one_hot ){
	sz_short* num_labels = (sz_short*)malloc(one_hot.length*sizeof(sz_short));
	for(sz_long i=0;i<one_hot.length;i++){
		num_labels[i]=0;
		for(sz_short j=0;j<one_hot.num_class;j++) num_labels[i]+= one_hot.bin[j][i];
	}
	return num_labels;
}

// Return a one_hot_mat with non-zeros for each row given to
// the k-largest (specified by num_lpn) corresponding soft labels
// Does not do sorting and has O(kC) complexity per node instead
// of O(C logC). Will be slow if C and k large
one_hot_mat top_k_mlabel( double* soft_labels , sz_short* num_lpn, sz_long length, sz_short num_class ){
	double max_max,max,val;
	sz_short max_ind;
	one_hot_mat label_out = init_one_hot(num_class ,length);

	for(sz_long i=0; i<length;i++ ){
		max_max =1.1f;
		for(sz_short k = 0; k<num_lpn[i];k++){
			max = 0.0f;
			for(sz_short j=0; j<num_class ;j++ ){
				val = soft_labels[j*length +i];
				if( ( val < max_max) && (val > max) ){
					max = val;
					max_ind = j;
				}
			}
			max_max = max;
			label_out.bin[max_ind][i]=1;
		}
	}

	return label_out;
}

// Computes rate of true positive, False positive and false negative
// A binary mask determines which values I am interested in (usually unlabeled entries..)
detector_stats get_detector_stats( sz_short* true_val, sz_short* pred_val, sz_long len, int* mask ){

	detector_stats stats = {.true_pos = 0.0,
			        .true_neg =0.0,
			        .false_pos = 0.0,
			        .false_neg = 0.0};

	sz_long k=0;
	for(sz_long i=0;i<len;i++){
	    if(mask[i]!=0){
	    	k++;
		if(true_val[i]==1){
			if(pred_val[i]==1){
				stats.true_pos+=1.0;
			}else{
				stats.false_neg+=1.0;
			}
		}else{
			if(pred_val[i]==1){
				stats.false_pos+=1.0;
			}else{
				stats.true_neg+=1.0;
			}
		}
	    }
	}

	stats.true_pos/=(double) k;
	stats.true_neg/=(double) k;
	stats.false_pos/=(double) k;
	stats.false_neg/=(double) k;

	return stats;
}

// Harmonic mean of 2 real numbers
double harmonic_mean( double x, double y ){
	if(x+y == 0.0f){
		printf("ERROR: Harmonic mean evaluated at 0\n ");
		exit(EXIT_FAILURE);
	}
	return (2.0f*(x*y))/(x+y);
}

// Produce micro and macro averaged precision and recall
// Input is the more simple detector_stats for different classes
classifier_stats  get_class_stats(detector_stats* all_stats , sz_short num_class ){

	classifier_stats stats = {.micro_precision=0.0f,
				  .macro_precision=0.0f,
				  .micro_recall=0.0f,
				  .macro_recall=0.0f };

	double sum_of_true_pos = 0.0f, micro_prec_denom = 0.0f, micro_recall_denom = 0.0f ;

	for(sz_short i=0; i<num_class;i++){
		sum_of_true_pos += all_stats[i].true_pos;
		micro_prec_denom += all_stats[i].true_pos + all_stats[i].false_pos;
		micro_recall_denom += all_stats[i].true_pos + all_stats[i].false_neg;

		stats.macro_precision += all_stats[i].true_pos/( all_stats[i].true_pos + all_stats[i].false_pos );
		stats.macro_recall += all_stats[i].true_pos/( all_stats[i].true_pos + all_stats[i].false_neg );
	}

	stats.micro_precision = sum_of_true_pos/micro_prec_denom;
	stats.micro_recall = sum_of_true_pos/micro_recall_denom;

	stats.macro_precision /= (double) num_class;
	stats.macro_recall /= (double) num_class;

	return stats;
}

// Takes array of detector statistics as input and returns array with f-1 scores for every class

void get_per_class_f1_scores(double* f1_per_class, detector_stats* all_stats, sz_short num_class ){

	for(sz_short i=0;i<num_class;i++){
		double TP = all_stats[i].true_pos;
		double FN = all_stats[i].false_neg;
		double FP = all_stats[i].false_pos;
		f1_per_class[i] = 2.0*TP / ( 2.0*TP + FN + FP );
	}
}

// Computes Micro and Macro f1-score
// Inputs must be one_hot structs
// f1_scores are only evaluated on the indexes denoted by unlabeled

f1_scores get_averaged_f1_scores( one_hot_mat true_labels, one_hot_mat pred_labels,
	                          sz_long* unlabeled , sz_long num_unlabeled ){

	if( (true_labels.num_class != pred_labels.num_class) || (true_labels.length != pred_labels.length) ){
		printf("Classifier stats ERROR: One-hot matrixes dimensions dont match\n");
		exit(EXIT_FAILURE);
	}

	int mask[true_labels.length];
	memset(mask, 0 , true_labels.length*sizeof(int));
	for(sz_long i=0;i< num_unlabeled; i++ ){ mask[unlabeled[i]]=1; }

        // Start by gathering detection statistics per each class

	detector_stats all_stats[true_labels.num_class];

	for(sz_short i=0;i<true_labels.num_class;i++){
		all_stats[i] = get_detector_stats( true_labels.bin[i] , pred_labels.bin[i] , true_labels.length , mask );
	}

	// Use detection stats to get Micro and Macro averaged F1 scores

	f1_scores scores = {.micro = 0.0f,
			    .macro = 0.0f };

	classifier_stats class_stats = get_class_stats(all_stats , true_labels.num_class );

	double f1_per_class[true_labels.num_class];

	get_per_class_f1_scores(f1_per_class, all_stats, true_labels.num_class);

	scores.micro = harmonic_mean(class_stats.micro_precision, class_stats.micro_recall );
        scores.macro = mean( f1_per_class, (int)true_labels.num_class );
/*
        This definition of Macro F1 (using macro averaged precision and recall) does not seem to work well
	scores.macro = harmonic_mean(class_stats.macro_precision, class_stats.macro_recall ); //
*/
	return scores;
}

//Compute error rate

double accuracy(class_t* true_label, class_t* pred,  sz_long* unlabeled, sz_long num_unlabeled){
	double sum=0.0;
	for(sz_long i=0;i<num_unlabeled;i++){
			if(true_label[unlabeled[i]]==pred[unlabeled[i]])
				sum+=1.0;
		}

	return sum/(double)num_unlabeled;
}


text:
//#include "readArray.h"
//#include "tensor_train.h"

#include "matrix.h"
#include "tensor.h"
#include "sketch.h"
#include "tt.h"
#include "VTime.h"
#include "readArray_MPI.h"
#include "PSTT.h"
#include "my_io.h"
#include "SSTT.h"

#include <omp.h>

#include <stdio.h>
#include <cblas.h>
#include <lapacke.h>
#include <math.h>
#include <string.h>
#include <gperftools/heap-profiler.h>

#ifndef HEAD
#define HEAD ((int) 0)
#endif

#ifndef FILE_BUF
#define FILE_BUF ((int) 256)
#endif

VTime_gathered* get_TT(tensor_train* tt, MPI_tensor* ten, int mid, int solvetype, int *heap_profile, char* path, long* total_alloc_ptr);

void MyHeapProfilerStart(int heap_profile, char* path, char* prefix, int rank)
{
    if (heap_profile){
        char* profile_path = (char*) calloc(FILE_BUF, sizeof(char));
        snprintf(profile_path, FILE_BUF, "%s%s%d", path, prefix, rank);
        HeapProfilerStart(profile_path);
    }
}

void MyHeapProfilerStop(int* heap_profile, char* path, char* prefix, int rank, long* total_alloc_ptr)
{
    if (*heap_profile){
        HeapProfilerDump("");
        HeapProfilerStop();

        char file_name[256];
        sprintf(file_name, "%s%s%d.000%d.heap", path, prefix, rank, *heap_profile);
        FILE *f = fopen(file_name, "r");
        while(f){
            *heap_profile = *heap_profile + 1;
            sprintf(file_name, "%s%s%d.000%d.heap", path, prefix, rank, *heap_profile);
            f = fopen(file_name, "r");
        }
        sprintf(file_name, "%s%s%d.000%d.heap", path, prefix, rank, *heap_profile - 1);
        f = fopen(file_name, "r");

        char total_alloc_str[256];
        char c = fgetc(f);
        while( c != '['){
            c = fgetc(f);
        }
        while( c != ':'){
            c = fgetc(f);
        }
        c = fgetc(f);
        int ii = 0;
        while( c != ']'){
            if (c != ' '){
                total_alloc_str[ii] = c;
                ii = ii+1;
            }
            c = fgetc(f);
        }
        total_alloc_str[ii] = '\0';

        int info = sscanf(total_alloc_str, "%ld", total_alloc_ptr);
//        printf("total_alloc_str = %s, total_alloc = %ld\n", total_alloc_str, *total_alloc_ptr);

//        long MB = 1;
//        for (int ii = 0; ii < 20; ++ii){
//            MB = MB * 2;
//        }
//        printf("size = %f MB\n ", (double) *total_alloc_ptr / MB);
    }
}

long* gather_heap_data(int heap_profile, long* total_allocs, int size, int rank, MPI_Comm comm)
{

    int head = 0;
    long* total_allocs_matrix = NULL;

    if (heap_profile){
        if (rank == head){
            total_allocs_matrix = (long*) calloc((size + 1)*3, sizeof(long));
        }
        MPI_Gather(total_allocs, 3, MPI_LONG, total_allocs_matrix + 3, 3, MPI_LONG, head, comm);

        for (int ii = 0; ii < size; ++ii){
            for (int jj = 0; jj < 3; ++jj){
                total_allocs[jj + ii*3];
            }
        }

        if (rank == head){
            for (int jj = 1; jj < size + 1; ++jj){
                for (int ii = 0; ii < 3; ++ii){
                    total_allocs_matrix[ii] += total_allocs_matrix[ii + jj*3];
                }
            }
            long MB = 1;
            for (int ii = 0; ii < 20; ++ii){
                MB = MB*2;
            }
            printf("TOTAL ALGORITHM SIZE = %f\n", (double) total_allocs_matrix[2]/MB);
        }
    }

    return total_allocs_matrix;
}

//Input:
// argv[1] = solvetype
// argv[2] = tensortype
// -d      = dimension
// -n      = sizes
// -r      = ranks
// -p      = path
// -e      = nps
// -m      = mid (optional)
// -h      = heap_profile
// -M      = M for gaussian bumps (only for tensortype == 2)
// -g      = gamma for gaussian bumps (only for tensortype == 2)
int main (int argc, char *argv[])
{
    MPI_Init(&argc, &argv);

    int world_rank;
    int world_size;
    MPI_Comm comm_world = MPI_COMM_WORLD;
    MPI_Comm_rank(comm_world, &world_rank);
    MPI_Comm_size(comm_world, &world_size);
    if (world_rank == HEAD){ printf("Starting MPI_timing\n");}

    int solvetype = atoi(argv[1]);
    int tensortype = atoi(argv[2]);

    int* dvec = NULL;
    int d = 0;

    int* midvec = NULL;
    int mid = -1;

    int* heap_profile_vec = NULL;
    int heap_profile = 0;

    int* nps = NULL;

    int* n = NULL;
    int* r = NULL;

    void (*f_ten)(double* restrict, int*, int*, const void*);
    void* parameters = NULL;

    long* total_allocs = NULL;

    double acc = 1e-10;

    // Taken chars: 'd', 'm', 'n', 'r', 'p', 'e', 'g'

    dvec = get_int_star_input('d', 1, argc, argv);
    if (dvec == NULL){
        printf("No input d for Hilbert Tensor\n");
        MPI_Finalize();
        return 0;
    }
    else{
        d = dvec[0];
    }

    midvec = get_int_star_input('m', 1, argc, argv);
    if (midvec != NULL){
        mid = midvec[0];
    }

    // 'e' is for 'eeeeeeeeh-nps'
    nps = get_int_star_input('e', d, argc, argv);
    if (nps == NULL){
        printf("No input nps\n");
        MPI_Finalize();
        return 0;
    }

    n = get_int_star_input('n', d, argc, argv);
    if (n == NULL){
        printf("No input n\n");
        MPI_Finalize();
        return 0;
    }

    r = get_int_star_input('r', d+1, argc, argv);
    if (r == NULL){
        printf("No input r\n");
        MPI_Finalize();
        return 0;
    }

    char* path = get_string_input ('p', argc, argv);
    if (path == NULL){
        printf("Must input path with -p\n");
        MPI_Finalize();
        return 0;
    }

    heap_profile_vec = get_int_star_input('h', 1, argc, argv);
    if (heap_profile_vec != NULL){
        heap_profile = *heap_profile_vec;
        total_allocs = (long*) calloc(3, sizeof(long));
    }

    // Get the tensor pointer (I have _no_ idea how to make a function for this anymore)
    if (tensortype == 1){
        f_ten = &f_hilbert;
        parameters = p_hilbert_init(d, n);
    }
    else if (tensortype == 2){
        int* Mvec = NULL;
        int M;
        Mvec = get_int_star_input('M', 1, argc, argv);
        if (Mvec != NULL){
            M = Mvec[0];
        }
        else{
            printf("Must input M with -M when tensortype == 2\n");
            MPI_Finalize();
            return 0;
        }

        double* gammavec = NULL;
        double gamma;
        gammavec = get_double_star_input('g', 1, argc, argv);
        if (gammavec != NULL){
            gamma = gammavec[0];
        }
        else{
            printf("Must input gamma with -g when tensortype == 2\n");
            MPI_Finalize();
            return 0;
        }

        int seed = 168234591; // Arbitrary number
        parameters = unit_random_p_gaussian_bumps_init(d, n, M, gamma, seed);
        f_ten = &f_gaussian_bumps;

        free(gammavec); gammavec = NULL;
        free(Mvec);     Mvec = NULL;
    }

    MyHeapProfilerStart(heap_profile, path, ".tt", world_rank);
    tensor_train* tt = TT_init_rank(d, n, r);
    MyHeapProfilerStop(&heap_profile, path, ".tt", world_rank, total_allocs);

    MyHeapProfilerStart(heap_profile, path, ".ten", world_rank);
    MPI_tensor* ten = MPI_tensor_init(d, n, nps, MPI_COMM_WORLD, f_ten, parameters);
    MyHeapProfilerStop(&heap_profile, path, ".ten", world_rank, total_allocs + 1);

    VTime_gathered* tms = get_TT(tt, ten, mid, solvetype, &heap_profile, path, total_allocs + 2);
    if (world_rank == HEAD){ printf("Got the tensor train\n");}

    if (solvetype == 1){
        ten = MPI_tensor_init(d, n, nps, MPI_COMM_WORLD, f_ten, parameters);
    }
    double err = tt_error(tt, ten);
    if (world_rank == HEAD){
        printf("ERROR = %e\n", err);
        printf("TIME = %f\n", tms->t[0]);
    }
    MPI_tensor_free(ten);

    long* total_allocs_matrix = gather_heap_data(heap_profile, total_allocs, world_size, world_rank, comm_world);
    if (world_rank == HEAD){
        my_tt_write(tt, tms, err, nps, tensortype, solvetype, path, total_allocs_matrix);
    }

    if (tms != NULL){;
        VTime_gathered_free(tms); tms = NULL;
    }
    TT_free(tt); tt = NULL;

    if (tensortype == 1){
        p_hilbert_free(parameters); parameters = NULL;
    }
    else if (tensortype == 2){
        p_gaussian_bumps_free(parameters); parameters = NULL;
    }

    free(r);      r = NULL;
    free(n);      n = NULL;
    free(nps);    nps = NULL;
    free(midvec); midvec = NULL;
    free(dvec);   dvec = NULL;
    free(total_allocs); total_allocs = NULL;
    free(heap_profile_vec); heap_profile_vec = NULL;

    MPI_Finalize();

    return 0;
}

/* Key for argv[1]
 *  1 : streaming_SSTT
 *  2 : streaming_PSTT2
 *  3 : streaming_PSTT2_onepass
*/
VTime_gathered* get_TT(tensor_train* tt, MPI_tensor* ten, int mid, int solvetype, int *heap_profile, char* path, long* total_alloc_ptr)
{
    VTime* tm;
    int world_rank = ten->rank;
    MPI_Comm comm = ten->comm;

    if (world_rank==HEAD){ printf("Getting the tensor train ("); }
    if (solvetype == 1){
        if (world_rank==HEAD){ printf("SSTT)\n"); }
        MyHeapProfilerStart(*heap_profile, path, ".SSTT", world_rank);
        tm = SSTT(tt, ten);
        MyHeapProfilerStop(heap_profile, path, ".SSTT", world_rank, total_alloc_ptr);
    }
    if (solvetype == 2){
        if (world_rank==HEAD){ printf("PSTT2)\n"); }
        MyHeapProfilerStart(*heap_profile, path, ".PSTT", world_rank);
        tm = PSTT2(tt, ten, mid);
        MyHeapProfilerStop(heap_profile, path, ".PSTT", world_rank, total_alloc_ptr);
    }
    if (solvetype == 3){
        if (world_rank==HEAD){ printf("PSTT2_onepass)\n"); }
        MyHeapProfilerStart(*heap_profile, path, ".PSTT_onepass", world_rank);
        tm = PSTT2_onepass(tt, ten, mid);
        MyHeapProfilerStop(heap_profile, path, ".PSTT_onepass", world_rank, total_alloc_ptr);
    }
    if (world_rank==HEAD){ printf("Got tensor train\n");}

    VTime_gathered* tms = VTime_gather(tm, comm);

    if (tm != NULL){
        VTime_free(tm); tm = NULL;
    }

    return tms;
}

ext: c


text:
#include <assert.h>
#include <gsl/gsl_sf_gamma.h>
#include <inttypes.h>
#include <math.h>
#include <stdlib.h>

#include "index.h"
#include "parse_mmaps.h"
#include "probability.h"

double reference_probability(const struct mmap_info* mmap_info,
			     struct revision_assignment* parent,
			     struct revision_assignment* child,
			     int disagrees,
			     const struct index_patch* index_patch) {
  double action_prob = 1.0;
  if (parent != NULL && parent->topic >= 0 && parent->pov >= 0
      && child->topic >= 0 && child->pov >= 0) {
    struct revision_assignment_header* revision_assignment_header
      = (struct revision_assignment_header*)mmap_info->revision_assignment_mmap;
    struct topic_summary* topic_summary;
    struct pov_summary* pov_dist;
    get_topic_summary(mmap_info, child->topic, &topic_summary, &pov_dist, NULL);
    if (parent->topic == child->topic) {
      if (parent->pov != child->pov) {
	// POV reference
	struct pov_summary* pov_summary = get_ant_pov(mmap_info, pov_dist,
						      child->pov, parent->pov);
	double denom = (double)(pov_summary->norevert_count
				+ pov_summary->revert_count)
	  + revision_assignment_header->psi_alpha
	  + revision_assignment_header->psi_beta;
	int correction = 0;
	if (index_patch != NULL
	    && index_patch->topic == child->topic
	    && index_patch->pov == child->pov
	    && index_patch->parent_pov == parent->pov
	    && index_patch->parent_topic == index_patch->topic) {
	  assert(denom -= 1.0 >= 0.0);
	  if (!disagrees == !(index_patch->disagrees)) {
	    correction += 1;
	  }
	}
	if (index_patch != NULL
	    && index_patch->child_topic == child->topic
	    && index_patch->child_pov == child->pov
	    && index_patch->pov == parent->pov
	    && index_patch->topic == index_patch->child_topic) {
	  assert(denom -= 1.0 >= 0.0);
	  if (!disagrees == !(index_patch->child_disagrees)) {
	    correction += 1;
	  }
	}
	if (disagrees) {
	  action_prob *= ((double)(pov_summary->revert_count
				   - correction)
			  + revision_assignment_header->psi_alpha) / denom;
	} else {
	  action_prob *= ((double)(pov_summary->norevert_count
				   - correction)
			  + revision_assignment_header->psi_beta) / denom;
	}
      } else {
	// topic reference
	double denom = (double)(topic_summary->revert_topic_count
				+ topic_summary->norevert_topic_count
				+ revision_assignment_header->gamma_alpha
				+ revision_assignment_header->gamma_beta);
	int correction = 0;
	if (index_patch != NULL
	    && index_patch->topic == child->topic
	    && index_patch->parent_pov == index_patch->pov
	    && index_patch->parent_topic == index_patch->topic) {
	  assert(denom -= 1.0 >= 0.0);
	  if (!disagrees == !(index_patch->disagrees)) {
	    correction += 1;
	  }
	}
	if (index_patch != NULL
	    && index_patch->child_topic == child->topic
	    && index_patch->pov == index_patch->child_pov
	    && index_patch->topic == index_patch->child_topic) {
	  assert(denom -= 1.0 >= 0.0);
	  if (!disagrees == !(index_patch->child_disagrees)) {
	    correction += 1;
	  }
	}
	if (disagrees) {
	  action_prob *= ((double)(topic_summary->revert_topic_count
				   - correction)
			  + revision_assignment_header->gamma_alpha) / denom;
	} else {
	  action_prob *= ((double)(topic_summary->norevert_topic_count
				   - correction)
			  + revision_assignment_header->gamma_beta) / denom;
	}
      }
    } else {
      // general reference
      double denom = (double)(topic_summary->revert_general_count
			      + topic_summary->norevert_general_count)
	+ revision_assignment_header->gamma_alpha
	+ revision_assignment_header->gamma_beta;
      int correction = 0;
      if (index_patch != NULL
	  && index_patch->parent_topic >= 0
	  && index_patch->topic == child->topic
	  && index_patch->parent_topic != index_patch->topic) {
	assert(denom -= 1.0 >= 0.0);
	if (!disagrees == !(index_patch->disagrees)) {
	  correction += 1;
	}
      }
      if (index_patch != NULL
	  && index_patch->topic >= 0
	  && index_patch->child_topic == child->topic
	  && index_patch->topic != index_patch->child_topic) {
	assert(denom -= 1.0 >= 0.0);
	if (!disagrees == !(index_patch->child_disagrees)) {
	  correction = 1;
	}
      }
      if (disagrees) {
	action_prob *= ((double)(topic_summary->revert_general_count
				 - correction)
			+ revision_assignment_header->gamma_alpha) / denom;
      } else {
	action_prob *= ((double)(topic_summary->norevert_general_count
				 - correction)
			+ revision_assignment_header->gamma_beta) / denom;
      }
    }
  } // Else, a constant not dependant on POV or topic
  assert(action_prob > 0.0);
  return action_prob;
}

double revision_probability(const struct mmap_info* mmap_info, int64_t revision_id,
			    int32_t topic, int32_t pov,
			    int include_child,
			    int include_user_topic_pov,
			    const struct index_patch* index_patch) {
  const struct revision* revision = get_revision(mmap_info, revision_id);
  struct revision_assignment current_assignment;
  current_assignment.topic = topic;
  current_assignment.pov = pov;
  double ret = 1.0;
  // Revert probability
  if (revision->parent >= 0) {
    ret *= reference_probability(mmap_info,
				 get_revision_assignment(mmap_info, revision->parent),
				 &current_assignment, revision->disagrees,
				 index_patch);
  }
  if (revision->child >= 0 && include_child) {
    ret *= reference_probability(mmap_info, &current_assignment,
				 get_revision_assignment(mmap_info, revision->child),
				 get_revision(mmap_info, revision->child)->disagrees,
				 index_patch);
  }

  // Topic, pov probability
  struct revision_assignment_header* revision_assignment_header
    = (struct revision_assignment_header*)mmap_info->revision_assignment_mmap;
  if (include_user_topic_pov) {
    int pov_correction = 0;
    if (index_patch
	&& index_patch->user == revision->user
	&& pov == index_patch->pov
	&& topic == index_patch->topic) {
      pov_correction = 1;
    }

    double* user_pov_dist;
    get_user_topics(mmap_info, revision->user, &user_pov_dist);
    // Alpha is already added to this array. The normalizing constant is
    // independant of topic/pov, so we exclude it here.
    ret *= (user_pov_dist[topic * revision_assignment_header->pov_per_topic + pov]
	    - pov_correction);
  }

  // Page probability
  struct topic_summary_header* topic_summary_header
    = (struct topic_summary_header*)mmap_info->topic_index_mmap;
  int64_t* page_dist;
  struct topic_summary* topic_summary;
  get_topic_summary(mmap_info, topic, &topic_summary, NULL, &page_dist);
  int64_t topic_page_revisions = page_dist[revision->article];
  int64_t topic_revisions = topic_summary->total_revisions;
  if (index_patch && index_patch->topic == topic) {
    topic_revisions -= 1;
    if (index_patch->page == revision->article) {
      topic_page_revisions -= 1;
    }
  }
  ret *= ((double)(topic_page_revisions) + revision_assignment_header->beta)
    / ((double)(topic_revisions)
       + revision_assignment_header->beta * topic_summary_header->num_pages);
  assert(ret > 0.0);
  return ret;
}

void fill_index_patch(const struct mmap_info* mmap_info,
		      int64_t revision_id,
		      struct index_patch* index_patch) {
  const struct revision* revision = get_revision(mmap_info, revision_id);
  index_patch->user = revision->user;
  index_patch->page = revision->article;
  const struct revision_assignment* revision_assignment
    = get_revision_assignment(mmap_info, revision_id);
  index_patch->topic = revision_assignment->topic;
  index_patch->pov = revision_assignment->pov;
  index_patch->disagrees = revision->disagrees;
  if (revision->parent >= 0) {
    const struct revision_assignment* parent_assignment
      = get_revision_assignment(mmap_info, revision->parent);
    index_patch->parent_topic = parent_assignment->topic;
    index_patch->parent_pov = parent_assignment->pov;
  } else {
    index_patch->parent_topic = -1;
    index_patch->parent_pov = -1;
  }
  if (revision->child >= 0) {
    const struct revision_assignment* child_assignment
      = get_revision_assignment(mmap_info, revision->child);
    index_patch->child_topic = child_assignment->topic;
    index_patch->child_pov = child_assignment->pov;
    index_patch->child_disagrees = get_revision(mmap_info,
						revision->child)->disagrees;
  } else {
    index_patch->child_topic = -1;
    index_patch->child_pov = -1;
    index_patch->child_disagrees = -1;
  }
}

double users_pages_probability_modn(const struct mmap_info* mmap_info, int sample, int modn) {
  struct revision_assignment_header* revision_assignment_header
    = (struct revision_assignment_header*)mmap_info->revision_assignment_mmap;
  struct user_topic_header* user_topic_header = (struct user_topic_header*)mmap_info->user_topic_mmap;
  double log_likelihood = 0.0;
  int64_t user_edits;
  double* user_topic_pov_dist;
  int topic_pov_count
    = revision_assignment_header->num_topics * revision_assignment_header->pov_per_topic;
  for (int64_t user_num = sample; user_num < user_topic_header->num_users; user_num += modn) {
    get_user(mmap_info, user_num, &user_edits, NULL);
    get_user_topics(mmap_info, user_num, &user_topic_pov_dist);
    if (user_edits <= 0) {
      // This user does not actually exist
      continue;
    }
    log_likelihood -= gsl_sf_lngamma(user_edits + revision_assignment_header->alpha * topic_pov_count);
    for (int topic = 0; topic < revision_assignment_header->num_topics; ++topic) {
      for (int pov = 0; pov < revision_assignment_header->pov_per_topic; ++pov) {
	log_likelihood += gsl_sf_lngamma(user_topic_pov_dist[topic * revision_assignment_header->pov_per_topic
							     + pov]);
      }
    }
  }
  struct topic_summary_header* topic_summary_header
    = (struct topic_summary_header*)mmap_info->topic_index_mmap;
  int64_t* page_dist;
  for (int topic = 0; topic < revision_assignment_header->num_topics; ++topic) {
    get_topic_summary(mmap_info, topic, NULL, NULL, &page_dist);
    for (int64_t page = sample; page < topic_summary_header->num_pages; page += modn) {
      log_likelihood += gsl_sf_lngamma(page_dist[page] + revision_assignment_header->beta);
    }
  }
  return log_likelihood;
}

double log_likelihood(const struct mmap_info* mmap_info) {
  return log_likelihood_gamma(mmap_info, 1);
}

double log_likelihood_gamma(const struct mmap_info* mmap_info, int include_users_pages) {
  struct user_topic_header* user_topic_header = (struct user_topic_header*)mmap_info->user_topic_mmap;
  struct revision_assignment_header* revision_assignment_header
    = (struct revision_assignment_header*)mmap_info->revision_assignment_mmap;
  int topic_pov_count
    = revision_assignment_header->num_topics * revision_assignment_header->pov_per_topic;
  double log_likelihood = 0.0;
  log_likelihood += user_topic_header->num_users
    * gsl_sf_lngamma(revision_assignment_header->alpha * topic_pov_count);
  log_likelihood -= user_topic_header->num_users * topic_pov_count
    * gsl_sf_lngamma(revision_assignment_header->alpha);
  if (include_users_pages) {
    log_likelihood += users_pages_probability_modn(mmap_info, 0, 1);
  }

  struct topic_summary_header* topic_summary_header
    = (struct topic_summary_header*)mmap_info->topic_index_mmap;
  struct topic_summary* topic_summary;
  struct pov_summary* pov_dist;
  struct pov_summary* pov_summary;
  log_likelihood += revision_assignment_header->num_topics
    * gsl_sf_lngamma(revision_assignment_header->beta * topic_summary_header->num_pages);
  log_likelihood -= revision_assignment_header->num_topics * topic_summary_header->num_pages
    * gsl_sf_lngamma(revision_assignment_header->beta);
  for (int topic = 0; topic < revision_assignment_header->num_topics; ++topic) {
    get_topic_summary(mmap_info, topic, &topic_summary, &pov_dist, NULL);
    log_likelihood -= gsl_sf_lngamma(topic_summary->total_revisions + revision_assignment_header->beta
				     * topic_summary_header->num_pages);
    // General reverts
    log_likelihood += gsl_sf_lngamma(topic_summary->revert_general_count
				     + revision_assignment_header->gamma_alpha);
    log_likelihood += gsl_sf_lngamma(topic_summary->norevert_general_count
				     + revision_assignment_header->gamma_beta);
    log_likelihood -= gsl_sf_lngamma(topic_summary->revert_general_count
				     + topic_summary->norevert_general_count
				     + revision_assignment_header->gamma_alpha
				     + revision_assignment_header->gamma_beta);
    // Topic reverts
    log_likelihood += gsl_sf_lngamma(topic_summary->revert_topic_count
				     + revision_assignment_header->gamma_alpha);
    log_likelihood += gsl_sf_lngamma(topic_summary->norevert_topic_count
				     + revision_assignment_header->gamma_beta);
    log_likelihood -= gsl_sf_lngamma(topic_summary->revert_topic_count
				     + topic_summary->norevert_topic_count
				     + revision_assignment_header->gamma_alpha
				     + revision_assignment_header->gamma_beta);

    // POV reverts
    for (int pov = 0; pov < revision_assignment_header->pov_per_topic; ++pov) {
      for (int ant_pov = 0; ant_pov < revision_assignment_header->pov_per_topic - 1; ++ant_pov) {
	pov_summary = pov_dist + pov * (revision_assignment_header->pov_per_topic - 1) + ant_pov;
	log_likelihood += gsl_sf_lngamma(pov_summary->revert_count
					 + revision_assignment_header->psi_alpha);
	log_likelihood += gsl_sf_lngamma(pov_summary->norevert_count
					 + revision_assignment_header->psi_beta);
	log_likelihood -= gsl_sf_lngamma(pov_summary->revert_count
					 + pov_summary->norevert_count
					 + revision_assignment_header->psi_alpha
					 + revision_assignment_header->psi_beta);
      }
    }
  }
  log_likelihood += 2 * revision_assignment_header->num_topics
    * gsl_sf_lngamma(revision_assignment_header->gamma_alpha + revision_assignment_header->gamma_beta);
  log_likelihood -= 2 * revision_assignment_header->num_topics
    * gsl_sf_lngamma(revision_assignment_header->gamma_alpha);
  log_likelihood -= 2 * revision_assignment_header->num_topics
    * gsl_sf_lngamma(revision_assignment_header->gamma_beta);

  log_likelihood += revision_assignment_header->pov_per_topic
    * (revision_assignment_header->pov_per_topic - 1)
    * revision_assignment_header->num_topics
    * gsl_sf_lngamma(revision_assignment_header->psi_alpha + revision_assignment_header->psi_beta);
  log_likelihood -= revision_assignment_header->pov_per_topic
    * (revision_assignment_header->pov_per_topic - 1)
    * revision_assignment_header->num_topics
    * gsl_sf_lngamma(revision_assignment_header->psi_alpha);
  log_likelihood -= revision_assignment_header->pov_per_topic
    * (revision_assignment_header->pov_per_topic - 1)
    * revision_assignment_header->num_topics
    * gsl_sf_lngamma(revision_assignment_header->psi_beta);

  return log_likelihood;
}

max_stars_repo_name: allenlavoie/topic-pov

max_issues_repo_name: allenlavoie/topic-pov

max_forks_repo_name: allenlavoie/topic-pov


text:
/* Copyright 2019-2020 Canaan 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/gsl-lite.hpp>
#include <type_traits>

#if defined(_MSC_VER)
#ifdef NNCASE_DLL
#define NNCASE_API __declspec(dllexport)
#elif defined(NNCASE_SHARED_LIBS)
#define NNCASE_API __declspec(dllimport)
#else
#define NNCASE_API
#endif
#else
#define NNCASE_API
#endif

#if defined(_MSC_VER)
#define NNCASE_UNREACHABLE() __assume(0)
#else
#define NNCASE_UNREACHABLE() __builtin_unreachable()
#endif

#if gsl_CPP17_OR_GREATER
#define NNCASE_INLINE_VAR inline
#define NNCASE_UNUSED [[maybe_unused]]
namespace nncase
{
template <class Callable, class... Args>
using invoke_result_t = std::invoke_result_t<Callable, Args...>;
}
#else
#define NNCASE_INLINE_VAR
#if defined(_MSC_VER)
#define NNCASE_UNUSED
#else
#define NNCASE_UNUSED __attribute__((unused))
#endif
namespace nncase
{
template <class Callable, class... Args>
using invoke_result_t = std::result_of_t<Callable(Args...)>;
}
#endif

#define NNCASE_LITTLE_ENDIAN 1

#define NNCASE_HAVE_STD_BYTE gsl_CPP17_OR_GREATER
#define NNCASE_NODISCARD gsl_NODISCARD
#define NNCASE_NORETURN gsl_NORETURN

#define BEGIN_NS_NNCASE_RUNTIME \
    namespace nncase            \
    {                           \
        namespace runtime       \
        {
#define END_NS_NNCASE_RUNTIME \
    }                         \
    }

#define BEGIN_NS_NNCASE_RT_STACKVM \
    namespace nncase               \
    {                              \
        namespace runtime          \
        {                          \
            namespace stackvm      \
            {
#define END_NS_NNCASE_RT_STACKVM \
    }                            \
    }                            \
    }

#define BEGIN_NS_NNCASE_KERNELS \
    namespace nncase            \
    {                           \
        namespace kernels       \
        {

#define END_NS_NNCASE_KERNELS \
    }                         \
    }

#ifndef DEFINE_ENUM_BITMASK_OPERATORS
#define DEFINE_ENUM_BITMASK_OPERATORS(ENUMTYPE) gsl_DEFINE_ENUM_BITMASK_OPERATORS(ENUMTYPE)
#endif

namespace nncase
{
struct default_init_t
{
};

NNCASE_INLINE_VAR constexpr default_init_t default_init {};
}

ext: h


text:
#include <math.h>
#include <stdlib.h>
#include <drfftw_mpi.h>
#include <mpi.h>
#include <gsl/gsl_rng.h>

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

int main(int argc, char **argv)
{
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
  MPI_Comm_size(MPI_COMM_WORLD, &NTask);

  if(argc < 2)
    {
      if(ThisTask == 0)
	{
	  fprintf(stdout, "\nParameters are missing.\n");
	  fprintf(stdout, "Call with <ParameterFile>\n\n");
	}
      MPI_Finalize();
      exit(0);
    }

  read_parameterfile(argv[1]);

  set_units();

  initialize_powerspectrum();

  initialize_ffts();

  read_glass(GlassFile);

  displacement_fields();

  write_particle_data();

  if(NumPart)
    free(P);

  free_ffts();

  if(ThisTask == 0)
    {
      printf("\nIC's generated.\n\n");
      printf("Initial scale factor = %g\n", InitTime);
      printf("\n");
    }

  MPI_Barrier(MPI_COMM_WORLD);
  print_spec();

  MPI_Finalize();		/* clean up & finalize MPI */
  exit(0);
}

void displacement_fields(void)
{
  MPI_Request request;
  MPI_Status status;
  gsl_rng *random_generator;
  int i, j, k, ii, jj, kk, axes;
  int n;
  int sendTask, recvTask;
  double fac, vel_prefac;
  double kvec[3], kmag, kmag2, p_of_k;
  double delta, phase, ampl, hubble_a;
  double u, v, w;
  double f1, f2, f3, f4, f5, f6, f7, f8;
  double dis, maxdisp, max_disp_glob;
  unsigned int *seedtable;

#ifdef CORRECT_CIC
  double fx, fy, fz, ff, smth;
#endif

  if(ThisTask == 0)
    {
      printf("\nstart computing displacement fields...\n");
      fflush(stdout);
    }

  hubble_a =
    Hubble * sqrt(Omega / pow(InitTime, 3) + (1 - Omega - OmegaLambda) / pow(InitTime, 2) + OmegaLambda);

  vel_prefac = InitTime * hubble_a * F_Omega(InitTime);

  vel_prefac /= sqrt(InitTime);	/* converts to Gadget velocity */

  if(ThisTask == 0)
    printf("vel_prefac= %g  hubble_a=%g fom=%g \n", vel_prefac, hubble_a, F_Omega(InitTime));

  fac = pow(2 * PI / Box, 1.5);

  maxdisp = 0;

  random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);

  gsl_rng_set(random_generator, Seed);

  if(!(seedtable = malloc(Nmesh * Nmesh * sizeof(unsigned int))))
    FatalError(4);

  for(i = 0; i < Nmesh / 2; i++)
    {
      for(j = 0; j < i; j++)
	seedtable[i * Nmesh + j] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i + 1; j++)
	seedtable[j * Nmesh + i] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i; j++)
	seedtable[(Nmesh - 1 - i) * Nmesh + j] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i + 1; j++)
	seedtable[(Nmesh - 1 - j) * Nmesh + i] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i; j++)
	seedtable[i * Nmesh + (Nmesh - 1 - j)] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i + 1; j++)
	seedtable[j * Nmesh + (Nmesh - 1 - i)] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i; j++)
	seedtable[(Nmesh - 1 - i) * Nmesh + (Nmesh - 1 - j)] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i + 1; j++)
	seedtable[(Nmesh - 1 - j) * Nmesh + (Nmesh - 1 - i)] = 0x7fffffff * gsl_rng_uniform(random_generator);
    }

  for(axes = 0; axes < 3; axes++)
    {
      if(ThisTask == 0)
	{
	  printf("\nstarting axes=%d...\n", axes);
	  fflush(stdout);
	}

      /* first, clean the array */
      for(i = 0; i < Local_nx; i++)
	for(j = 0; j < Nmesh; j++)
	  for(k = 0; k <= Nmesh / 2; k++)
	    {
	      Cdata[(i * Nmesh + j) * (Nmesh / 2 + 1) + k].re = 0;
	      Cdata[(i * Nmesh + j) * (Nmesh / 2 + 1) + k].im = 0;
	    }

      for(i = 0; i < Nmesh; i++)
	{
	  ii = Nmesh - i;
	  if(ii == Nmesh)
	    ii = 0;
	  if((i >= Local_x_start && i < (Local_x_start + Local_nx)) ||
	     (ii >= Local_x_start && ii < (Local_x_start + Local_nx)))
	    {
	      for(j = 0; j < Nmesh; j++)
		{
		  gsl_rng_set(random_generator, seedtable[i * Nmesh + j]);

		  for(k = 0; k < Nmesh / 2; k++)
		    {
		      phase = gsl_rng_uniform(random_generator) * 2 * PI;
		      do
			ampl = gsl_rng_uniform(random_generator);
		      while(ampl == 0);

		      if(i == Nmesh / 2 || j == Nmesh / 2 || k == Nmesh / 2)
			continue;
		      if(i == 0 && j == 0 && k == 0)
			continue;

		      if(i < Nmesh / 2)
			kvec[0] = i * 2 * PI / Box;
		      else
			kvec[0] = -(Nmesh - i) * 2 * PI / Box;

		      if(j < Nmesh / 2)
			kvec[1] = j * 2 * PI / Box;
		      else
			kvec[1] = -(Nmesh - j) * 2 * PI / Box;

		      if(k < Nmesh / 2)
			kvec[2] = k * 2 * PI / Box;
		      else
			kvec[2] = -(Nmesh - k) * 2 * PI / Box;

		      kmag2 = kvec[0] * kvec[0] + kvec[1] * kvec[1] + kvec[2] * kvec[2];
		      kmag = sqrt(kmag2);

		      if(SphereMode == 1)
			{
			  if(kmag * Box / (2 * PI) > Nsample / 2)	/* select a sphere in k-space */
			    continue;
			}
		      else
			{
			  if(fabs(kvec[0]) * Box / (2 * PI) > Nsample / 2)
			    continue;
			  if(fabs(kvec[1]) * Box / (2 * PI) > Nsample / 2)
			    continue;
			  if(fabs(kvec[2]) * Box / (2 * PI) > Nsample / 2)
			    continue;
			}

		      p_of_k = PowerSpec(kmag);

		      p_of_k *= -log(ampl);

		      delta = fac * sqrt(p_of_k) / Dplus;	/* scale back to starting redshift */

#ifdef CORRECT_CIC
		      /* do deconvolution of CIC interpolation */
		      fx = fy = fz = 1;
		      if(kvec[0] != 0)
			{
			  fx = (kvec[0] * Box / 2) / Nmesh;
			  fx = sin(fx) / fx;
			}
		      if(kvec[1] != 0)
			{
			  fy = (kvec[1] * Box / 2) / Nmesh;
			  fy = sin(fy) / fy;
			}
		      if(kvec[2] != 0)
			{
			  fz = (kvec[2] * Box / 2) / Nmesh;
			  fz = sin(fz) / fz;
			}
		      ff = 1 / (fx * fy * fz);
		      smth = ff * ff;

		      delta *= smth;
		      /* end deconvolution */
#endif
		      if(k > 0)
			{
			  if(i >= Local_x_start && i < (Local_x_start + Local_nx))
			    {
			      Cdata[((i - Local_x_start) * Nmesh + j) * (Nmesh / 2 + 1) + k].re =
				-kvec[axes] / kmag2 * delta * sin(phase);
			      Cdata[((i - Local_x_start) * Nmesh + j) * (Nmesh / 2 + 1) + k].im =
				kvec[axes] / kmag2 * delta * cos(phase);
			    }
			}
		      else	/* k=0 plane needs special treatment */
			{
			  if(i == 0)
			    {
			      if(j >= Nmesh / 2)
				continue;
			      else
				{
				  if(i >= Local_x_start && i < (Local_x_start + Local_nx))
				    {
				      jj = Nmesh - j;	/* note: j!=0 surely holds at this point */

				      Cdata[((i - Local_x_start) * Nmesh + j) * (Nmesh / 2 + 1) + k].re =
					-kvec[axes] / kmag2 * delta * sin(phase);
				      Cdata[((i - Local_x_start) * Nmesh + j) * (Nmesh / 2 + 1) + k].im =
					kvec[axes] / kmag2 * delta * cos(phase);

				      Cdata[((i - Local_x_start) * Nmesh + jj) * (Nmesh / 2 + 1) + k].re =
					-kvec[axes] / kmag2 * delta * sin(phase);
				      Cdata[((i - Local_x_start) * Nmesh + jj) * (Nmesh / 2 + 1) + k].im =
					-kvec[axes] / kmag2 * delta * cos(phase);
				    }
				}
			    }
			  else	/* here comes i!=0 : conjugate can be on other processor! */
			    {
			      if(i >= Nmesh / 2)
				continue;
			      else
				{
				  ii = Nmesh - i;
				  if(ii == Nmesh)
				    ii = 0;
				  jj = Nmesh - j;
				  if(jj == Nmesh)
				    jj = 0;

				  if(i >= Local_x_start && i < (Local_x_start + Local_nx))
				    {
				      Cdata[((i - Local_x_start) * Nmesh + j) * (Nmesh / 2 + 1) + k].re =
					-kvec[axes] / kmag2 * delta * sin(phase);
				      Cdata[((i - Local_x_start) * Nmesh + j) * (Nmesh / 2 + 1) + k].im =
					kvec[axes] / kmag2 * delta * cos(phase);
				    }

				  if(ii >= Local_x_start && ii < (Local_x_start + Local_nx))
				    {
				      Cdata[((ii - Local_x_start) * Nmesh + jj) * (Nmesh / 2 + 1) +
					    k].re = -kvec[axes] / kmag2 * delta * sin(phase);
				      Cdata[((ii - Local_x_start) * Nmesh + jj) * (Nmesh / 2 + 1) +
					    k].im = -kvec[axes] / kmag2 * delta * cos(phase);
				    }
				}
			    }
			}
		    }
		}
	    }
	}

      rfftwnd_mpi(Inverse_plan, 1, Disp, Workspace, FFTW_NORMAL_ORDER);		/** FFT **/

      /* now get the plane on the right side from neighbour on the right,
         and send the left plane */

      recvTask = ThisTask;
      do
	{
	  recvTask--;
	  if(recvTask < 0)
	    recvTask = NTask - 1;
	}
      while(Local_nx_table[recvTask] == 0);

      sendTask = ThisTask;
      do
	{
	  sendTask++;
	  if(sendTask >= NTask)
	    sendTask = 0;
	}
      while(Local_nx_table[sendTask] == 0);

      /* use non-blocking send */

      if(Local_nx > 0)
	{
	  MPI_Isend(&Disp[0],
		    sizeof(fftw_real) * Nmesh * (2 * (Nmesh / 2 + 1)),
		    MPI_BYTE, recvTask, 10, MPI_COMM_WORLD, &request);

	  MPI_Recv(&Disp[(Local_nx * Nmesh) * (2 * (Nmesh / 2 + 1))],
		   sizeof(fftw_real) * Nmesh * (2 * (Nmesh / 2 + 1)),
		   MPI_BYTE, sendTask, 10, MPI_COMM_WORLD, &status);

	  MPI_Wait(&request, &status);
	}

      /* read-out displacements */

      for(n = 0; n < NumPart; n++)
	{
	  {
	    u = P[n].Pos[0] / Box * Nmesh;
	    v = P[n].Pos[1] / Box * Nmesh;
	    w = P[n].Pos[2] / Box * Nmesh;

	    i = (int) u;
	    j = (int) v;
	    k = (int) w;

	    if(i == (Local_x_start + Local_nx))
	      i = (Local_x_start + Local_nx) - 1;
	    if(i < Local_x_start)
	      i = Local_x_start;
	    if(j == Nmesh)
	      j = Nmesh - 1;
	    if(k == Nmesh)
	      k = Nmesh - 1;

	    u -= i;
	    v -= j;
	    w -= k;

	    i -= Local_x_start;
	    ii = i + 1;
	    jj = j + 1;
	    kk = k + 1;

	    if(jj >= Nmesh)
	      jj -= Nmesh;
	    if(kk >= Nmesh)
	      kk -= Nmesh;

	    f1 = (1 - u) * (1 - v) * (1 - w);
	    f2 = (1 - u) * (1 - v) * (w);
	    f3 = (1 - u) * (v) * (1 - w);
	    f4 = (1 - u) * (v) * (w);
	    f5 = (u) * (1 - v) * (1 - w);
	    f6 = (u) * (1 - v) * (w);
	    f7 = (u) * (v) * (1 - w);
	    f8 = (u) * (v) * (w);

	    dis = Disp[(i * Nmesh + j) * (2 * (Nmesh / 2 + 1)) + k] * f1 +
	      Disp[(i * Nmesh + j) * (2 * (Nmesh / 2 + 1)) + kk] * f2 +
	      Disp[(i * Nmesh + jj) * (2 * (Nmesh / 2 + 1)) + k] * f3 +
	      Disp[(i * Nmesh + jj) * (2 * (Nmesh / 2 + 1)) + kk] * f4 +
	      Disp[(ii * Nmesh + j) * (2 * (Nmesh / 2 + 1)) + k] * f5 +
	      Disp[(ii * Nmesh + j) * (2 * (Nmesh / 2 + 1)) + kk] * f6 +
	      Disp[(ii * Nmesh + jj) * (2 * (Nmesh / 2 + 1)) + k] * f7 +
	      Disp[(ii * Nmesh + jj) * (2 * (Nmesh / 2 + 1)) + kk] * f8;

	    P[n].Vel[axes] = dis;

	    if(dis > maxdisp)
	      maxdisp = dis;
	  }
	}
    }

  /* now add displacement to Lagrangian coordinates, and multiply velocities by correct factor */
  for(n = 0; n < NumPart; n++)
    {
      for(axes = 0; axes < 3; axes++)
	{
	  P[n].Pos[axes] += P[n].Vel[axes];
	  P[n].Vel[axes] *= vel_prefac;
	  P[n].Pos[axes] = periodic_wrap(P[n].Pos[axes]);
	}
    }

  gsl_rng_free(random_generator);

  MPI_Reduce(&maxdisp, &max_disp_glob, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);

  if(ThisTask == 0)
    {
      printf("\nMaximum displacement: %g kpc/h, in units of the part-spacing= %g\n",
	     max_disp_glob, max_disp_glob / (Box / Nmesh));
    }
}

double periodic_wrap(double x)
{
  while(x >= Box)
    x -= Box;

  while(x < 0)
    x += Box;

  return x;
}

void set_units(void)		/* ... set some units */
{
  UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s;

  G = GRAVITY / pow(UnitLength_in_cm, 3) * UnitMass_in_g * pow(UnitTime_in_s, 2);
  Hubble = HUBBLE * UnitTime_in_s;
}

void initialize_ffts(void)
{
  int total_size, i, additional;
  int local_ny_after_transpose, local_y_start_after_transpose;
  int *slab_to_task_local;
  size_t bytes;

  Inverse_plan = rfftw3d_mpi_create_plan(MPI_COMM_WORLD,
					 Nmesh, Nmesh, Nmesh, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);

  rfftwnd_mpi_local_sizes(Inverse_plan, &Local_nx, &Local_x_start,
			  &local_ny_after_transpose, &local_y_start_after_transpose, &total_size);

  Local_nx_table = malloc(sizeof(int) * NTask);
  MPI_Allgather(&Local_nx, 1, MPI_INT, Local_nx_table, 1, MPI_INT, MPI_COMM_WORLD);

  if(ThisTask == 0)
    {
      for(i = 0; i < NTask; i++)
	printf("Task=%d Local_nx=%d\n", i, Local_nx_table[i]);
      fflush(stdout);
    }

  Slab_to_task = malloc(sizeof(int) * Nmesh);
  slab_to_task_local = malloc(sizeof(int) * Nmesh);

  for(i = 0; i < Nmesh; i++)
    slab_to_task_local[i] = 0;

  for(i = 0; i < Local_nx; i++)
    slab_to_task_local[Local_x_start + i] = ThisTask;

  MPI_Allreduce(slab_to_task_local, Slab_to_task, Nmesh, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

  free(slab_to_task_local);

  additional = (Nmesh) * (2 * (Nmesh / 2 + 1));	/* additional plane on the right side */

  Disp = (fftw_real *) malloc(bytes = sizeof(fftw_real) * (total_size + additional));
  Workspace = (fftw_real *) malloc(bytes += sizeof(fftw_real) * total_size);

  if(Disp && Workspace)
    {
      if(ThisTask == 0)
	printf("\nallocated %g Mbyte on Task %d for FFT's\n", bytes / (1024.0 * 1024.0), ThisTask);
    }
  else
    {
      printf("failed to allocate %g Mbyte on Task %d\n", bytes / (1024.0 * 1024.0), ThisTask);
      printf("bailing out.\n");
      FatalError(1);
    }

  Cdata = (fftw_complex *) Disp;	/* transformed array */
}

void free_ffts(void)
{
  free(Workspace);
  free(Disp);
  free(Slab_to_task);
  rfftwnd_mpi_destroy_plan(Inverse_plan);
}

int FatalError(int errnum)
{
  printf("FatalError called with number=%d\n", errnum);
  fflush(stdout);
  MPI_Abort(MPI_COMM_WORLD, errnum);
  exit(0);
}

static double A, B, alpha, beta, V, gf;

double fnl(double x)		/* Peacock & Dodds formula */
{
  return x * pow((1 + B * beta * x + pow(A * x, alpha * beta)) /
		 (1 + pow(pow(A * x, alpha) * gf * gf * gf / (V * sqrt(x)), beta)), 1 / beta);
}

void print_spec(void)
{
  double k, knl, po, dl, dnl, neff, kf, kstart, kend, po2, po1, DDD;
  char buf[1000];
  FILE *fd;

  if(ThisTask == 0)
    {
      sprintf(buf, "%s/inputspec_%s.txt", OutputDir, FileBase);

      fd = fopen(buf, "w");

      gf = GrowthFactor(0.001, 1.0) / (1.0 / 0.001);

      DDD = GrowthFactor(1.0 / (Redshift + 1), 1.0);

      fprintf(fd, "%12g %12g\n", Redshift, DDD);	/* print actual starting redshift and
							   linear growth factor for this cosmology */

      kstart = 2 * PI / (1000.0 * (3.085678e24 / UnitLength_in_cm));	/* 1000 Mpc/h */
      kend = 2 * PI / (0.001 * (3.085678e24 / UnitLength_in_cm));	/* 0.001 Mpc/h */

      for(k = kstart; k < kend; k *= 1.025)
	{
	  po = PowerSpec(k);
	  dl = 4.0 * PI * k * k * k * po;

	  kf = 0.5;

	  po2 = PowerSpec(1.001 * k * kf);
	  po1 = PowerSpec(k * kf);

	  if(po != 0 && po1 != 0 && po2 != 0)
	    {
	      neff = (log(po2) - log(po1)) / (log(1.001 * k * kf) - log(k * kf));

	      if(1 + neff / 3 > 0)
		{
		  A = 0.482 * pow(1 + neff / 3, -0.947);
		  B = 0.226 * pow(1 + neff / 3, -1.778);
		  alpha = 3.310 * pow(1 + neff / 3, -0.244);
		  beta = 0.862 * pow(1 + neff / 3, -0.287);
		  V = 11.55 * pow(1 + neff / 3, -0.423) * 1.2;

		  dnl = fnl(dl);
		  knl = k * pow(1 + dnl, 1.0 / 3);
		}
	      else
		{
		  dnl = 0;
		  knl = 0;
		}
	    }
	  else
	    {
	      dnl = 0;
	      knl = 0;
	    }

	  fprintf(fd, "%12g %12g    %12g %12g\n", k, dl, knl, dnl);
	}
      fclose(fd);
    }
}

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


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

/* Author:  G. Jungman */
/* Miscellaneous support functions for Bessel function evaluations.
 */
#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_airy.h"
#include "gsl_sf_elementary.h"
#include "gsl_sf_exp.h"
#include "gsl_sf_gamma.h"
#include "gsl_sf_trig.h"

#include "error.h"

#include "bessel_amp_phase.h"
#include "bessel_temme.h"
#include "bessel.h"

#define CubeRoot2_  1.25992104989487316476721060728

/* Debye functions [Abramowitz+Stegun, 9.3.9-10] */

inline static double
debye_u1(const double * tpow)
{
  return (3.0*tpow[1] - 5.0*tpow[3])/24.0;
}

inline static double
debye_u2(const double * tpow)
{
  return (81.0*tpow[2] - 462.0*tpow[4] + 385.0*tpow[6])/1152.0;
}

inline
static double debye_u3(const double * tpow)
{
  return (30375.0*tpow[3] - 369603.0*tpow[5] + 765765.0*tpow[7] - 425425.0*tpow[9])/414720.0;
}

inline
static double debye_u4(const double * tpow)
{
  return (4465125.0*tpow[4] - 94121676.0*tpow[6] + 349922430.0*tpow[8] -
          446185740.0*tpow[10] + 185910725.0*tpow[12])/39813120.0;
}

inline
static double debye_u5(const double * tpow)
{
  return (1519035525.0*tpow[5]     - 49286948607.0*tpow[7] +
          284499769554.0*tpow[9]   - 614135872350.0*tpow[11] +
          566098157625.0*tpow[13]  - 188699385875.0*tpow[15])/6688604160.0;
}

#if 0
inline
static double debye_u6(const double * tpow)
{
  return (2757049477875.0*tpow[6] - 127577298354750.0*tpow[8] +
          1050760774457901.0*tpow[10] - 3369032068261860.0*tpow[12] +
          5104696716244125.0*tpow[14] - 3685299006138750.0*tpow[16] +
          1023694168371875.0*tpow[18])/4815794995200.0;
}
#endif

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_bessel_IJ_taylor_e(const double nu, const double x,
                             const int sign,
                             const int kmax,
                             const double threshold,
                             gsl_sf_result * result
                             )
{
  /* CHECK_POINTER(result) */

  if(nu < 0.0 || x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    if(nu == 0.0) {
      result->val = 1.0;
      result->err = 0.0;
    }
    else {
      result->val = 0.0;
      result->err = 0.0;
    }
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result prefactor;   /* (x/2)^nu / Gamma(nu+1) */
    gsl_sf_result sum;

    int stat_pre;
    int stat_sum;
    int stat_mul;

    if(nu == 0.0) {
      prefactor.val = 1.0;
      prefactor.err = 0.0;
      stat_pre = GSL_SUCCESS;
    }
    else if(nu < INT_MAX-1) {
      /* Separate the integer part and use
       * y^nu / Gamma(nu+1) = y^N /N! y^f / (N+1)_f,
       * to control the error.
       */
      const int    N = (int)floor(nu + 0.5);
      const double f = nu - N;
      gsl_sf_result poch_factor;
      gsl_sf_result tc_factor;
      const int stat_poch = gsl_sf_poch_e(N+1.0, f, &poch_factor);
      const int stat_tc   = gsl_sf_taylorcoeff_e(N, 0.5*x, &tc_factor);
      const double p = pow(0.5*x,f);
      prefactor.val  = tc_factor.val * p / poch_factor.val;
      prefactor.err  = tc_factor.err * p / poch_factor.val;
      prefactor.err += fabs(prefactor.val) / poch_factor.val * poch_factor.err;
      prefactor.err += 2.0 * GSL_DBL_EPSILON * fabs(prefactor.val);
      stat_pre = GSL_ERROR_SELECT_2(stat_tc, stat_poch);
    }
    else {
      gsl_sf_result lg;
      const int stat_lg = gsl_sf_lngamma_e(nu+1.0, &lg);
      const double term1  = nu*log(0.5*x);
      const double term2  = lg.val;
      const double ln_pre = term1 - term2;
      const double ln_pre_err = GSL_DBL_EPSILON * (fabs(term1)+fabs(term2)) + lg.err;
      const int stat_ex = gsl_sf_exp_err_e(ln_pre, ln_pre_err, &prefactor);
      stat_pre = GSL_ERROR_SELECT_2(stat_ex, stat_lg);
    }

    /* Evaluate the sum.
     * [Abramowitz+Stegun, 9.1.10]
     * [Abramowitz+Stegun, 9.6.7]
     */
    {
      const double y = sign * 0.25 * x*x;
      double sumk = 1.0;
      double term = 1.0;
      int k;

      for(k=1; k<=kmax; k++) {
        term *= y/((nu+k)*k);
        sumk += term;
        if(fabs(term/sumk) < threshold) break;
      }

      sum.val = sumk;
      sum.err = threshold * fabs(sumk);

      stat_sum = ( k >= kmax ? GSL_EMAXITER : GSL_SUCCESS );
    }

    stat_mul = gsl_sf_multiply_err_e(prefactor.val, prefactor.err,
                                        sum.val, sum.err,
                                        result);

    return GSL_ERROR_SELECT_3(stat_mul, stat_pre, stat_sum);
  }
}

/* x >> nu*nu+1
 * error ~ O( ((nu*nu+1)/x)^3 )
 *
 * empirical error analysis:
 *   choose  GSL_ROOT3_MACH_EPS * x > (nu*nu + 1)
 *
 * This is not especially useful. When the argument gets
 * large enough for this to apply, the cos() and sin()
 * start loosing digits. However, this seems inevitable
 * for this particular method.
 */
int
gsl_sf_bessel_Jnu_asympx_e(const double nu, const double x, gsl_sf_result * result)
{
  double mu   = 4.0*nu*nu;
  double mum1 = mu-1.0;
  double mum9 = mu-9.0;
  double chi = x - (0.5*nu + 0.25)*M_PI;
  double P   = 1.0 - mum1*mum9/(128.0*x*x);
  double Q   = mum1/(8.0*x);
  double pre = sqrt(2.0/(M_PI*x));
  double c   = cos(chi);
  double s   = sin(chi);
  double r   = mu/x;
  result->val  = pre * (c*P - s*Q);
  result->err  = pre * GSL_DBL_EPSILON * (fabs(c*P) + fabs(s*Q));
  result->err += pre * fabs(0.1*r*r*r);
  return GSL_SUCCESS;
}

/* x >> nu*nu+1
 */
int
gsl_sf_bessel_Ynu_asympx_e(const double nu, const double x, gsl_sf_result * result)
{
  double ampl;
  double theta;
  double alpha = x;
  double beta  = -0.5*nu*M_PI;
  int stat_a = gsl_sf_bessel_asymp_Mnu_e(nu, x, &ampl);
  int stat_t = gsl_sf_bessel_asymp_thetanu_corr_e(nu, x, &theta);
  double sin_alpha = sin(alpha);
  double cos_alpha = cos(alpha);
  double sin_chi   = sin(beta + theta);
  double cos_chi   = cos(beta + theta);
  double sin_term     = sin_alpha * cos_chi + sin_chi * cos_alpha;
  double sin_term_mag = fabs(sin_alpha * cos_chi) + fabs(sin_chi * cos_alpha);
  result->val  = ampl * sin_term;
  result->err  = fabs(ampl) * GSL_DBL_EPSILON * sin_term_mag;
  result->err += fabs(result->val) * 2.0 * GSL_DBL_EPSILON;

  if(fabs(alpha) > 1.0/GSL_DBL_EPSILON) {
    result->err *= 0.5 * fabs(alpha);
  }
  else if(fabs(alpha) > 1.0/GSL_SQRT_DBL_EPSILON) {
    result->err *= 256.0 * fabs(alpha) * GSL_SQRT_DBL_EPSILON;
  }

  return GSL_ERROR_SELECT_2(stat_t, stat_a);
}

/* x >> nu*nu+1
 */
int
gsl_sf_bessel_Inu_scaled_asympx_e(const double nu, const double x, gsl_sf_result * result)
{
  double mu   = 4.0*nu*nu;
  double mum1 = mu-1.0;
  double mum9 = mu-9.0;
  double pre  = 1.0/sqrt(2.0*M_PI*x);
  double r    = mu/x;
  result->val = pre * (1.0 - mum1/(8.0*x) + mum1*mum9/(128.0*x*x));
  result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + pre * fabs(0.1*r*r*r);
  return GSL_SUCCESS;
}

/* x >> nu*nu+1
 */
int
gsl_sf_bessel_Knu_scaled_asympx_e(const double nu, const double x, gsl_sf_result * result)
{
  double mu   = 4.0*nu*nu;
  double mum1 = mu-1.0;
  double mum9 = mu-9.0;
  double pre  = sqrt(M_PI/(2.0*x));
  double r    = nu/x;
  result->val = pre * (1.0 + mum1/(8.0*x) + mum1*mum9/(128.0*x*x));
  result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + pre * fabs(0.1*r*r*r);
  return GSL_SUCCESS;
}

/* nu -> Inf; uniform in x > 0  [Abramowitz+Stegun, 9.7.7]
 *
 * error:
 *   The error has the form u_N(t)/nu^N  where  0 <= t <= 1.
 *   It is not hard to show that |u_N(t)| is small for such t.
 *   We have N=6 here, and |u_6(t)| < 0.025, so the error is clearly
 *   bounded by 0.025/nu^6. This gives the asymptotic bound on nu
 *   seen below as nu ~ 100. For general MACH_EPS it will be
 *                     nu > 0.5 / MACH_EPS^(1/6)
 *   When t is small, the bound is even better because |u_N(t)| vanishes
 *   as t->0. In fact u_N(t) ~ C t^N as t->0, with C ~= 0.1.
 *   We write
 *                     err_N <= min(0.025, C(1/(1+(x/nu)^2))^3) / nu^6
 *   therefore
 *                     min(0.29/nu^2, 0.5/(nu^2+x^2)) < MACH_EPS^{1/3}
 *   and this is the general form.
 *
 * empirical error analysis, assuming 14 digit requirement:
 *   choose   x > 50.000 nu   ==>  nu >   3
 *   choose   x > 10.000 nu   ==>  nu >  15
 *   choose   x >  2.000 nu   ==>  nu >  50
 *   choose   x >  1.000 nu   ==>  nu >  75
 *   choose   x >  0.500 nu   ==>  nu >  80
 *   choose   x >  0.100 nu   ==>  nu >  83
 *
 * This makes sense. For x << nu, the error will be of the form u_N(1)/nu^N,
 * since the polynomial term will be evaluated near t=1, so the bound
 * on nu will become constant for small x. Furthermore, increasing x with
 * nu fixed will decrease the error.
 */
int
gsl_sf_bessel_Inu_scaled_asymp_unif_e(const double nu, const double x, gsl_sf_result * result)
{
  int i;
  double z = x/nu;
  double root_term = sqrt(1.0 + z*z);
  double pre = 1.0/sqrt(2.0*M_PI*nu * root_term);
  double eta = root_term + log(z/(1.0+root_term));
  double ex_arg = ( z < 1.0/GSL_ROOT3_DBL_EPSILON ? nu*(-z + eta) : -0.5*nu/z*(1.0 - 1.0/(12.0*z*z)) );
  gsl_sf_result ex_result;
  int stat_ex = gsl_sf_exp_e(ex_arg, &ex_result);
  if(stat_ex == GSL_SUCCESS) {
    double t = 1.0/root_term;
    double sum;
    double tpow[16];
    tpow[0] = 1.0;
    for(i=1; i<16; i++) tpow[i] = t * tpow[i-1];
    sum = 1.0 + debye_u1(tpow)/nu + debye_u2(tpow)/(nu*nu) + debye_u3(tpow)/(nu*nu*nu)
          + debye_u4(tpow)/(nu*nu*nu*nu) + debye_u5(tpow)/(nu*nu*nu*nu*nu);
    result->val  = pre * ex_result.val * sum;
    result->err  = pre * ex_result.val / (nu*nu*nu*nu*nu*nu);
    result->err += pre * ex_result.err * fabs(sum);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    result->val = 0.0;
    result->err = 0.0;
    return stat_ex;
  }
}

/* nu -> Inf; uniform in x > 0  [Abramowitz+Stegun, 9.7.8]
 *
 * error:
 *   identical to that above for Inu_scaled
 */
int
gsl_sf_bessel_Knu_scaled_asymp_unif_e(const double nu, const double x, gsl_sf_result * result)
{
  int i;
  double z = x/nu;
  double root_term = sqrt(1.0 + z*z);
  double pre = sqrt(M_PI/(2.0*nu*root_term));
  double eta = root_term + log(z/(1.0+root_term));
  double ex_arg = ( z < 1.0/GSL_ROOT3_DBL_EPSILON ? nu*(z - eta) : 0.5*nu/z*(1.0 + 1.0/(12.0*z*z)) );
  gsl_sf_result ex_result;
  int stat_ex = gsl_sf_exp_e(ex_arg, &ex_result);
  if(stat_ex == GSL_SUCCESS) {
    double t = 1.0/root_term;
    double sum;
    double tpow[16];
    tpow[0] = 1.0;
    for(i=1; i<16; i++) tpow[i] = t * tpow[i-1];
    sum = 1.0 - debye_u1(tpow)/nu + debye_u2(tpow)/(nu*nu) - debye_u3(tpow)/(nu*nu*nu)
          + debye_u4(tpow)/(nu*nu*nu*nu) - debye_u5(tpow)/(nu*nu*nu*nu*nu);
    result->val  = pre * ex_result.val * sum;
    result->err  = pre * ex_result.err * fabs(sum);
    result->err += pre * ex_result.val / (nu*nu*nu*nu*nu*nu);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    result->val = 0.0;
    result->err = 0.0;
    return stat_ex;
  }
}

/* Evaluate J_mu(x),J_{mu+1}(x) and Y_mu(x),Y_{mu+1}(x)  for |mu| < 1/2
 */
int
gsl_sf_bessel_JY_mu_restricted(const double mu, const double x,
                               gsl_sf_result * Jmu, gsl_sf_result * Jmup1,
                               gsl_sf_result * Ymu, gsl_sf_result * Ymup1)
{
  /* CHECK_POINTER(Jmu) */
  /* CHECK_POINTER(Jmup1) */
  /* CHECK_POINTER(Ymu) */
  /* CHECK_POINTER(Ymup1) */

  if(x < 0.0 || fabs(mu) > 0.5) {
    Jmu->val   = 0.0;
    Jmu->err   = 0.0;
    Jmup1->val = 0.0;
    Jmup1->err = 0.0;
    Ymu->val   = 0.0;
    Ymu->err   = 0.0;
    Ymup1->val = 0.0;
    Ymup1->err = 0.0;
    GSL_ERROR ("error", GSL_EDOM);
  }
  else if(x == 0.0) {
    if(mu == 0.0) {
      Jmu->val   = 1.0;
      Jmu->err   = 0.0;
    }
    else {
      Jmu->val   = 0.0;
      Jmu->err   = 0.0;
    }
    Jmup1->val = 0.0;
    Jmup1->err = 0.0;
    Ymu->val   = 0.0;
    Ymu->err   = 0.0;
    Ymup1->val = 0.0;
    Ymup1->err = 0.0;
    GSL_ERROR ("error", GSL_EDOM);
  }
  else {
    int stat_Y;
    int stat_J;

    if(x < 2.0) {
      /* Use Taylor series for J and the Temme series for Y.
       * The Taylor series for J requires nu > 0, so we shift
       * up one and use the recursion relation to get Jmu, in
       * case mu < 0.
       */
      gsl_sf_result Jmup2;
      int stat_J1 = gsl_sf_bessel_IJ_taylor_e(mu+1.0, x, -1, 100, GSL_DBL_EPSILON,  Jmup1);
      int stat_J2 = gsl_sf_bessel_IJ_taylor_e(mu+2.0, x, -1, 100, GSL_DBL_EPSILON, &Jmup2);
      double c = 2.0*(mu+1.0)/x;
      Jmu->val  = c * Jmup1->val - Jmup2.val;
      Jmu->err  = c * Jmup1->err + Jmup2.err;
      Jmu->err += 2.0 * GSL_DBL_EPSILON * fabs(Jmu->val);
      stat_J = GSL_ERROR_SELECT_2(stat_J1, stat_J2);
      stat_Y = gsl_sf_bessel_Y_temme(mu, x, Ymu, Ymup1);
      return GSL_ERROR_SELECT_2(stat_J, stat_Y);
    }
    else if(x < 1000.0) {
      double P, Q;
      double J_ratio;
      double J_sgn;
      const int stat_CF1 = gsl_sf_bessel_J_CF1(mu, x, &J_ratio, &J_sgn);
      const int stat_CF2 = gsl_sf_bessel_JY_steed_CF2(mu, x, &P, &Q);
      double Jprime_J_ratio = mu/x - J_ratio;
      double gamma = (P - Jprime_J_ratio)/Q;
      Jmu->val = J_sgn * sqrt(2.0/(M_PI*x) / (Q + gamma*(P-Jprime_J_ratio)));
      Jmu->err = 4.0 * GSL_DBL_EPSILON * fabs(Jmu->val);
      Jmup1->val = J_ratio * Jmu->val;
      Jmup1->err = fabs(J_ratio) * Jmu->err;
      Ymu->val = gamma * Jmu->val;
      Ymu->err = fabs(gamma) * Jmu->err;
      Ymup1->val = Ymu->val * (mu/x - P - Q/gamma);
      Ymup1->err = Ymu->err * fabs(mu/x - P - Q/gamma) + 4.0*GSL_DBL_EPSILON*fabs(Ymup1->val);
      return GSL_ERROR_SELECT_2(stat_CF1, stat_CF2);
    }
    else {
      /* Use asymptotics for large argument.
       */
      const int stat_J0 = gsl_sf_bessel_Jnu_asympx_e(mu,     x, Jmu);
      const int stat_J1 = gsl_sf_bessel_Jnu_asympx_e(mu+1.0, x, Jmup1);
      const int stat_Y0 = gsl_sf_bessel_Ynu_asympx_e(mu,     x, Ymu);
      const int stat_Y1 = gsl_sf_bessel_Ynu_asympx_e(mu+1.0, x, Ymup1);
      stat_J = GSL_ERROR_SELECT_2(stat_J0, stat_J1);
      stat_Y = GSL_ERROR_SELECT_2(stat_Y0, stat_Y1);
      return GSL_ERROR_SELECT_2(stat_J, stat_Y);
    }
  }
}

int
gsl_sf_bessel_J_CF1(const double nu, const double x,
                    double * ratio, double * sgn)
{
  const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  const int maxiter = 10000;
  int n = 1;
  double Anm2 = 1.0;
  double Bnm2 = 0.0;
  double Anm1 = 0.0;
  double Bnm1 = 1.0;
  double a1 = x/(2.0*(nu+1.0));
  double An = Anm1 + a1*Anm2;
  double Bn = Bnm1 + a1*Bnm2;
  double an;
  double fn = An/Bn;
  double dn = a1;
  double s  = 1.0;

  while(n < maxiter) {
    double old_fn;
    double del;
    n++;
    Anm2 = Anm1;
    Bnm2 = Bnm1;
    Anm1 = An;
    Bnm1 = Bn;
    an = -x*x/(4.0*(nu+n-1.0)*(nu+n));
    An = Anm1 + an*Anm2;
    Bn = Bnm1 + an*Bnm2;

    if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
      An /= RECUR_BIG;
      Bn /= RECUR_BIG;
      Anm1 /= RECUR_BIG;
      Bnm1 /= RECUR_BIG;
      Anm2 /= RECUR_BIG;
      Bnm2 /= RECUR_BIG;
    }

    old_fn = fn;
    fn = An/Bn;
    del = old_fn/fn;

    dn = 1.0 / (2.0*(nu+n)/x - dn);
    if(dn < 0.0) s = -s;

    if(fabs(del - 1.0) < 2.0*GSL_DBL_EPSILON) break;
  }

  *ratio = fn;
  *sgn   = s;

  if(n >= maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

/* Evaluate the continued fraction CF1 for J_{nu+1}/J_nu
 * using Gautschi (Euler) equivalent series.
 * This exhibits an annoying problem because the
 * a_k are not positive definite (in fact they are all negative).
 * There are cases when rho_k blows up. Example: nu=1,x=4.
 */
#if 0
int
gsl_sf_bessel_J_CF1_ser(const double nu, const double x,
                        double * ratio, double * sgn)
{
  const int maxk = 20000;
  double tk   = 1.0;
  double sum  = 1.0;
  double rhok = 0.0;
  double dk = 0.0;
  double s  = 1.0;
  int k;

  for(k=1; k<maxk; k++) {
    double ak = -0.25 * (x/(nu+k)) * x/(nu+k+1.0);
    rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
    tk  *= rhok;
    sum += tk;

    dk = 1.0 / (2.0/x - (nu+k-1.0)/(nu+k) * dk);
    if(dk < 0.0) s = -s;

    if(fabs(tk/sum) < GSL_DBL_EPSILON) break;
  }

  *ratio = x/(2.0*(nu+1.0)) * sum;
  *sgn   = s;

  if(k == maxk)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}
#endif

/* Evaluate the continued fraction CF1 for I_{nu+1}/I_nu
 * using Gautschi (Euler) equivalent series.
 */
int
gsl_sf_bessel_I_CF1_ser(const double nu, const double x, double * ratio)
{
  const int maxk = 20000;
  double tk   = 1.0;
  double sum  = 1.0;
  double rhok = 0.0;
  int k;

  for(k=1; k<maxk; k++) {
    double ak = 0.25 * (x/(nu+k)) * x/(nu+k+1.0);
    rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
    tk  *= rhok;
    sum += tk;
    if(fabs(tk/sum) < GSL_DBL_EPSILON) break;
  }

  *ratio = x/(2.0*(nu+1.0)) * sum;

  if(k == maxk)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

int
gsl_sf_bessel_JY_steed_CF2(const double nu, const double x,
                           double * P, double * Q)
{
  const int max_iter = 10000;
  const double SMALL = 1.0e-100;

  int i = 1;

  double x_inv = 1.0/x;
  double a = 0.25 - nu*nu;
  double p = -0.5*x_inv;
  double q = 1.0;
  double br = 2.0*x;
  double bi = 2.0;
  double fact = a*x_inv/(p*p + q*q);
  double cr = br + q*fact;
  double ci = bi + p*fact;
  double den = br*br + bi*bi;
  double dr = br/den;
  double di = -bi/den;
  double dlr = cr*dr - ci*di;
  double dli = cr*di + ci*dr;
  double temp = p*dlr - q*dli;
  q = p*dli + q*dlr;
  p = temp;
  for (i=2; i<=max_iter; i++) {
    a  += 2*(i-1);
    bi += 2.0;
    dr = a*dr + br;
    di = a*di + bi;
    if(fabs(dr)+fabs(di) < SMALL) dr = SMALL;
    fact = a/(cr*cr+ci*ci);
    cr = br + cr*fact;
    ci = bi - ci*fact;
    if(fabs(cr)+fabs(ci) < SMALL) cr = SMALL;
    den = dr*dr + di*di;
    dr /= den;
    di /= -den;
    dlr = cr*dr - ci*di;
    dli = cr*di + ci*dr;
    temp = p*dlr - q*dli;
    q = p*dli + q*dlr;
    p = temp;
    if(fabs(dlr-1.0)+fabs(dli) < GSL_DBL_EPSILON) break;
  }

  *P = p;
  *Q = q;

  if(i == max_iter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

/* Evaluate continued fraction CF2, using Thompson-Barnett-Temme method,
 * to obtain values of exp(x)*K_nu and exp(x)*K_{nu+1}.
 *
 * This is unstable for small x; x > 2 is a good cutoff.
 * Also requires |nu| < 1/2.
 */
int
gsl_sf_bessel_K_scaled_steed_temme_CF2(const double nu, const double x,
                                       double * K_nu, double * K_nup1,
                                       double * Kp_nu)
{
  const int maxiter = 10000;

  int i = 1;
  double bi = 2.0*(1.0 + x);
  double di = 1.0/bi;
  double delhi = di;
  double hi    = di;

  double qi   = 0.0;
  double qip1 = 1.0;

  double ai = -(0.25 - nu*nu);
  double a1 = ai;
  double ci = -ai;
  double Qi = -ai;

  double s = 1.0 + Qi*delhi;

  for(i=2; i<=maxiter; i++) {
    double dels;
    double tmp;
    ai -= 2.0*(i-1);
    ci  = -ai*ci/i;
    tmp  = (qi - bi*qip1)/ai;
    qi   = qip1;
    qip1 = tmp;
    Qi += ci*qip1;
    bi += 2.0;
    di  = 1.0/(bi + ai*di);
    delhi = (bi*di - 1.0) * delhi;
    hi += delhi;
    dels = Qi*delhi;
    s += dels;
    if(fabs(dels/s) < GSL_DBL_EPSILON) break;
  }

  hi *= -a1;

  *K_nu   = sqrt(M_PI/(2.0*x)) / s;
  *K_nup1 = *K_nu * (nu + x + 0.5 - hi)/x;
  *Kp_nu  = - *K_nup1 + nu/x * *K_nu;
  if(i == maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

int gsl_sf_bessel_cos_pi4_e(double y, double eps, gsl_sf_result * result)
{
  const double sy = sin(y);
  const double cy = cos(y);
  const double s = sy + cy;
  const double d = sy - cy;
  const double abs_sum = fabs(cy) + fabs(sy);
  double seps;
  double ceps;
  if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
    const double e2 = eps*eps;
    seps = eps * (1.0 - e2/6.0 * (1.0 - e2/20.0));
    ceps = 1.0 - e2/2.0 * (1.0 - e2/12.0);
  }
  else {
    seps = sin(eps);
    ceps = cos(eps);
  }
  result->val = (ceps * s - seps * d)/ M_SQRT2;
  result->err = 2.0 * GSL_DBL_EPSILON * (fabs(ceps) + fabs(seps)) * abs_sum / M_SQRT2;

  /* Try to account for error in evaluation of sin(y), cos(y).
   * This is a little sticky because we don't really know
   * how the library routines are doing their argument reduction.
   * However, we will make a reasonable guess.
   * FIXME ?
   */
  if(y > 1.0/GSL_DBL_EPSILON) {
    result->err *= 0.5 * y;
  }
  else if(y > 1.0/GSL_SQRT_DBL_EPSILON) {
    result->err *= 256.0 * y * GSL_SQRT_DBL_EPSILON;
  }

  return GSL_SUCCESS;
}

int gsl_sf_bessel_sin_pi4_e(double y, double eps, gsl_sf_result * result)
{
  const double sy = sin(y);
  const double cy = cos(y);
  const double s = sy + cy;
  const double d = sy - cy;
  const double abs_sum = fabs(cy) + fabs(sy);
  double seps;
  double ceps;
  if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
    const double e2 = eps*eps;
    seps = eps * (1.0 - e2/6.0 * (1.0 - e2/20.0));
    ceps = 1.0 - e2/2.0 * (1.0 - e2/12.0);
  }
  else {
    seps = sin(eps);
    ceps = cos(eps);
  }
  result->val = (ceps * d + seps * s)/ M_SQRT2;
  result->err = 2.0 * GSL_DBL_EPSILON * (fabs(ceps) + fabs(seps)) * abs_sum / M_SQRT2;

  /* Try to account for error in evaluation of sin(y), cos(y).
   * See above.
   * FIXME ?
   */
  if(y > 1.0/GSL_DBL_EPSILON) {
    result->err *= 0.5 * y;
  }
  else if(y > 1.0/GSL_SQRT_DBL_EPSILON) {
    result->err *= 256.0 * y * GSL_SQRT_DBL_EPSILON;
  }

  return GSL_SUCCESS;
}

/************************************************************************
 *                                                                      *
  Asymptotic approximations 8.11.5, 8.12.5, and 8.42.7 from
  G.N.Watson, A Treatise on the Theory of Bessel Functions,
  2nd Edition (Cambridge University Press, 1944).
  Higher terms in expansion for x near l given by
  Airey in Phil. Mag. 31, 520 (1916).

  This approximation is accurate to near 0.1% at the boundaries
  between the asymptotic regions; well away from the boundaries
  the accuracy is better than 10^{-5}.
 *                                                                      *
 ************************************************************************/
#if 0
double besselJ_meissel(double nu, double x)
{
  double beta = pow(nu, 0.325);
  double result;

  /* Fitted matching points.   */
  double llimit = 1.1 * beta;
  double ulimit = 1.3 * beta;

  double nu2 = nu * nu;

  if (nu < 5. && x < 1.)
    {
      /* Small argument and order. Use a Taylor expansion. */
      int k;
      double xo2 = 0.5 * x;
      double gamfactor = pow(nu,nu) * exp(-nu) * sqrt(nu * 2. * M_PI)
	* (1. + 1./(12.*nu) + 1./(288.*nu*nu));
      double prefactor = pow(xo2, nu) / gamfactor;
      double C[5];

      C[0] = 1.;
      C[1] = -C[0] / (nu+1.);
      C[2] = -C[1] / (2.*(nu+2.));
      C[3] = -C[2] / (3.*(nu+3.));
      C[4] = -C[3] / (4.*(nu+4.));

      result = 0.;
      for(k=0; k<5; k++)
	result += C[k] * pow(xo2, 2.*k);

      result *= prefactor;
    }
  else if(x < nu - llimit)
    {
      /* Small x region: x << l.    */
      double z = x / nu;
      double z2 = z*z;
      double rtomz2 = sqrt(1.-z2);
      double omz2_2 = (1.-z2)*(1.-z2);

      /* Calculate Meissel exponent. */
      double term1 = 1./(24.*nu) * ((2.+3.*z2)/((1.-z2)*rtomz2) -2.);
      double term2 = - z2*(4. + z2)/(16.*nu2*(1.-z2)*omz2_2);
      double V_nu = term1 + term2;

      /* Calculate the harmless prefactor. */
      double sterlingsum = 1. + 1./(12.*nu) + 1./(288*nu2);
      double harmless = 1. / (sqrt(rtomz2*2.*M_PI*nu) * sterlingsum);

      /* Calculate the logarithm of the nu dependent prefactor. */
      double ln_nupre = rtomz2 + log(z) - log(1. + rtomz2);

      result = harmless * exp(nu*ln_nupre - V_nu);
    }
  else if(x < nu + ulimit)
    {
      /* Intermediate region 1: x near nu. */
      double eps = 1.-nu/x;
      double eps_x = eps * x;
      double eps_x_2 = eps_x * eps_x;
      double xo6 = x/6.;
      double B[6];
      static double gam[6] = {2.67894, 1.35412, 1., 0.89298, 0.902745, 1.};
      static double sf[6] = {0.866025, 0.866025, 0., -0.866025, -0.866025, 0.};

      /* Some terms are identically zero, because sf[] can be zero.
       * Some terms do not appear in the result.
       */
      B[0] = 1.;
      B[1] = eps_x;
      /* B[2] = 0.5 * eps_x_2 - 1./20.; */
      B[3] = eps_x * (eps_x_2/6. - 1./15.);
      B[4] = eps_x_2 * (eps_x_2 - 1.)/24. + 1./280.;
      /* B[5] = eps_x * (eps_x_2*(0.5*eps_x_2 - 1.)/60. + 43./8400.); */

      result  = B[0] * gam[0] * sf[0] / pow(xo6, 1./3.);
      result += B[1] * gam[1] * sf[1] / pow(xo6, 2./3.);
      result += B[3] * gam[3] * sf[3] / pow(xo6, 4./3.);
      result += B[4] * gam[4] * sf[4] / pow(xo6, 5./3.);

      result /= (3.*M_PI);
    }
  else
    {
      /* Region of very large argument. Use expansion
       * for x>>l, and we need not be very exacting.
       */
      double secb = x/nu;
      double sec2b= secb*secb;

      double cotb = 1./sqrt(sec2b-1.);      /* cotb=cot(beta) */

      double beta = acos(nu/x);
      double trigarg = nu/cotb - nu*beta - 0.25 * M_PI;

      double cot3b = cotb * cotb * cotb;
      double cot6b = cot3b * cot3b;

      double sum1, sum2, expterm, prefactor, trigcos;

      sum1  = 2.0 + 3.0 * sec2b;
      trigarg -= sum1 * cot3b / (24.0 * nu);

      trigcos = cos(trigarg);

      sum2 = 4.0 + sec2b;
      expterm = sum2 * sec2b * cot6b / (16.0 * nu2);

      expterm = exp(-expterm);
      prefactor = sqrt(2. * cotb / (nu * M_PI));

      result = prefactor * expterm * trigcos;
    }

  return  result;
}
#endif


text:
/*
 * Gsl_blas_dgemm: test gsl_blas_dgemm (matrix . matrix)
 *
 * Copyright (c) 2012 Jérémie Decock
 *
 * Required: GSL library (libgsl0-dev)
 * Usage: gcc gsl_blas_dgemm.c -lgsl -lgslcblas -lm
 *    or: gcc gsl_blas_dgemm.c $(pkg-config --libs gsl)
 */

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

int main (void)
{
    gsl_matrix * m1 = gsl_matrix_alloc(3, 3);
    gsl_matrix * m2 = gsl_matrix_alloc(3, 3);
    gsl_matrix * m3 = gsl_matrix_alloc(3, 3);

    gsl_matrix_set_all(m1, 3.0);
    gsl_matrix_set_all(m2, 2.0);

    /* Compute m3 = m1.m2 */

    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,
                   1.0, m1, m2,
                   0.0, m3);

    int i, j;
    for(i=0 ; i<3 ; i++) {
        printf("[");
        for(j=0 ; j<3 ; j++) {
            printf(" %f ", gsl_matrix_get(m3, i, j));
        }
        printf("]\n");
    }

    gsl_matrix_free(m1);
    gsl_matrix_free(m2);
    gsl_matrix_free(m3);

    return 0;
}

max_stars_repo_name: jeremiedecock/snippets

max_issues_repo_name: jeremiedecock/snippets

max_forks_repo_name: jeremiedecock/snippets


text:
/*
* -----------------------------------------------------------------
*  ODE Solver Library --- ode_lib.h
*  Version: 1.6180
*  Date: Feb 19, 2010
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  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.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the header file for a library
*  with the ODE solution tools.
* -----------------------------------------------------------------
*/

#ifndef __ODESOLVER_LIB_H__
#define __ODESOLVER_LIB_H__

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <cvode/cvode.h>

#undef  RTOL
#define RTOL 1.0e-6

#undef  ATOL
#define ATOL 1.0e-15

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

double uround(void);

double wnorm(int n,
             double *v1,
             double *v2);

void ewtset(int n,
            double *v,
            double atol,
            double rtol,
            double *ewt);

int jacobian(CVRhsFn f,
             void *f_data,
             gsl_vector *Fx,
             gsl_vector *x,
             double t,
             double atol,
             double rtol,
             gsl_matrix *J);

int gradient(CVRhsFn f,
             void *f_data,
             void *cvode_mem,
             double t0,
             double delta_t,
             double atol,
             double rtol,
             gsl_vector *phi,
             gsl_vector *Rphi,
             gsl_matrix *A);

void linear_approx(gsl_vector *x,
                   gsl_vector *x0,
                   gsl_vector *Fx0,
                   gsl_matrix *DFx0,
                   gsl_vector *Fx);

int odesolver_init(CVRhsFn f,
                   void *data,
                   double t0,
                   int mxsteps,
		           double atol,
                   double rtol,
                   gsl_vector *x,
                   void *cvode_mem);

int odesolver_reinit(CVRhsFn f,
                     void *f_data,
                     double t0,
		             double atol,
                     double rtol,
                     gsl_vector *x,
                     void *cvode_mem);

int odesolver(void *cvode_mem,
              double tf,
              gsl_vector *Fx);

#endif /* __ODESOLVER_LIB_H__ */

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
#if !defined(PETIGAPART_H)
#define PETIGAPART_H

#include <petsc.h>

PETSC_EXTERN PetscErrorCode IGA_Partition(PetscInt,PetscInt,
                                          PetscInt,const PetscInt[],
                                          PetscInt[],PetscInt[]);
PETSC_EXTERN PetscErrorCode IGA_Distribute(PetscInt,
                                           const PetscInt[],const PetscInt[],
                                           const PetscInt[],PetscInt[],PetscInt[]);

#endif/*PETIGAPART_H*/

max_stars_repo_name: otherlab/petiga

max_issues_repo_name: otherlab/petiga

max_forks_repo_name: otherlab/petiga