text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdint.h>
#include <limits.h> // for int limits etc.
#include <assert.h>
#include <time.h>
#include <gsl/gsl_histogram.h>
#include <omp.h>

#include "main.h"
#include "dSFMT-src-2.1/dSFMT.h"
//#include "ziggurat/ziggurat.h"

//#define STATS
#define MULTIPLIER 1000000  // seperates the coordiantes (poor man's hash)

int main (int argc, char * argv[]) {
	int num_particles = DSFMT_N64*4;
	int num_defects = 1024*4;
	int box_size = 128;
	long long nsteps = 1024*1024*2;

	//omp_init_lock(&omplock);
	//omp_set_num_threads(1);
	if (argc == 5) {
		num_particles = atoi(argv[1]);
		num_defects = atoi(argv[2]);
		box_size = atoi(argv[3]);
		nsteps = atoll(argv[4]);
	}
	else {
		printf("\n# ***** Using default values! *****\n\n");
		printf("# usage: rw #particles #defects #box #steps\n");

	}

	// Information:
	printf("#-------------------- Paramters --------------------\n");
	printf("# Particles: %i\n", num_particles);
	printf("# Box size: %i\n", box_size);
	printf("# Defects: %i (Density: %.2e)\n", num_defects, (float) num_defects/pow(box_size,3) );
	printf("# Steps: %lli\n", nsteps);
	printf("#---------------------------------------------------\n");

	assert( (box_size & (box_size - 1)) == 0 ); // check if box_size is power of two
	//zigset(1);
	dsfmt_t dsfmt;

	int seed = 1;

	// int num_random  = (num_particles < 1024) ? 1024: num_particles;
	// int  *correlation_times = malloc(num_particles * sizeof(int)); // random distribution of correlation times
	// char *directions = malloc(nsteps*sizeof(char)); // directions
	int **particles = malloc2D_i(num_particles, 4); // 2d array for particles: x,y,z,mag
	int **defects = malloc2D_i(num_defects, 3); // 2d array for defect coordinates

	//int *particle;
	//int direction;

	// init random number generator
	dsfmt_init_gen_rand(&dsfmt, seed);

	// check if we can create the hashes
	//assert(box_size < MULTIPLIER);
	assert(num_defects < pow(box_size,3));

#ifdef STATS
	// statistics
	// histogram of diretions
	gsl_histogram * h = gsl_histogram_alloc (6);
	gsl_histogram_set_ranges_uniform (h, 0, 6);

	// histogram of visits
	gsl_histogram * hvisits = gsl_histogram_alloc (10000);
	gsl_histogram_set_ranges_uniform (hvisits, 0, 10000);

	// gsl histogram seems not to be thread safe, allow only 1 thread
	omp_set_num_threads(1);
#endif

	// Start simulation

	// distribute particles from 0 to +box_size
	for (int j = 0; j < 3; j++) {
		for (int i=0 ;i < num_particles; i++) {
			particles[i][j] = (int) (dsfmt_genrand_close_open(&dsfmt)*box_size);
		}
	}

	// distribute defects from 0 to +box_size
	for (int j = 0; j < 3; j++) {
		for (int i=0 ; i < num_defects; i++) {
			int val = (int) (dsfmt_genrand_close_open(&dsfmt)*box_size);
			defects[i][j] = val;
		}
	}

	// METHOD 1: now create a hashed list to find them later
	// This will be a fallback to METHOD 2, in case memory is not enough
	int64_t *hash_list = malloc( num_defects*sizeof(int64_t) );
	for (int i = 0; i < num_defects; i++) {
		hash_list[i] = hash(defects[i][0],defects[i][1],defects[i][2]);
	}

	qsort(hash_list, num_defects, sizeof(int64_t), int64_cmp);

	// METHOD 2: create lookup table for x -> lookup table for y -> lookup table for z
	// The smart thing is that the pointers are NULL if there is no defect in the corresponding slab,
	// so only the z coordinates are really arrays
	//
	// If I need more space one could use the chars as bit fields,
	// One needs to calculate the offset (or index of char) to get to the proper group though.
	// get offset: offset = coord/sizeof(char) oder coord >> log2(sizeof(char))
	// set bit: array[offset] |= 1 << coord%sizeof(char)

	// in check_defekt_3d:
	// check bit: array[offset] & 1 << coord%sizeof(char)  > 1
	//
/*
	int mem_size_table = 0;
	// first coordinate (x) will be an array of pointers to an array of pointers
	char ***lookup_table = malloc(box_size * sizeof(char**));
	mem_size_table += box_size * sizeof(char**);
	// initialize the arrays to NULL pointer
	for (int i = 0; i < box_size; i++) {
		lookup_table[i]=NULL;
	}

	for (int i = 0; i < num_defects; i++) {
		int x_index = defects[i][0];
		int y_index = defects[i][1];
		int z_index = defects[i][2] / sizeof(char);
		// check if there  is already an array at x_index ...
		if (lookup_table[x_index] == NULL) {
			// ... it's not! Create an array of pointers for the second coordinate
			lookup_table[x_index] = malloc(box_size * sizeof(char*)); // malloc second coordinate pointers
			mem_size_table += box_size * sizeof(char*);
			for (int i = 0; i < box_size; i++)
				lookup_table[x_index][i]=NULL; // initialize the second coordiante pointers to NULL
		}

		// check if there is already an array at [x_index][y_index]
		if (lookup_table[x_index][y_index] == NULL) { // check if third coordinate array exists
			lookup_table[x_index][y_index] = malloc(box_size * sizeof(char)); // malloc third coordinate array
			mem_size_table += box_size * sizeof(char);
			for (int i = 0; i < box_size; i++)
				lookup_table[x_index][y_index][i]=0; // initialize the third array to zero

		}
		// set the defect coordinate
		lookup_table[x_index][y_index][z_index] = 1;

		 }

	int test_particle[4];
	double start = omp_get_wtime();
	for (int i = 0; i < 1024*1024*128; i++) {
		for (int j = 0; j < 3; j++) {
			test_particle[j] = (int) (dsfmt_genrand_close_open(&dsfmt)*box_size);
		}
		check_defect_3d(test_particle, lookup_table, 10);
	}
	double stop = omp_get_wtime();
	printf("# Time Method 2: %.2fs\n", stop-start);
	printf("# Lookup table size M1: %10.1fkB %i %i\n",  mem_size_table/((float) (1024)), num_defects, box_size);

*/

	// Method 3: Using the scheme above but for x,y,z seperately
	int ltb_N = (int) ceil( ( (double)box_size ) / sizeof(int));
	int xltb[ltb_N];
	int yltb[ltb_N];
	int zltb[ltb_N];
	for (int i = 0; i < ltb_N; i++) {
		xltb[i]=0;
		yltb[i]=0;
		zltb[i]=0;
	}
	for (int i = 0; i < num_defects; i++) {
		int xi =  defects[i][0] / sizeof(int);
		int xbit = defects[i][0] % sizeof(int);
		//printf("x %i %i %i %i\n",i, defects[i][0], xi, xbit);

		xltb[xi] |= 1 << xbit;

		int yi = defects[i][1] / sizeof(int);
		int ybit = defects[i][1] % sizeof(int);
		yltb[yi] |= 1 << ybit;
		//printf("y %i %i %i %i\n",i, defects[i][1], yi, ybit);

		int zi = defects[i][2] / sizeof(int);
		int zbit = defects[i][2] % sizeof(int);
		zltb[zi] |= 1 << zbit;
		//printf("z %i %i %i %i\n",i, defects[i][2], zi, zbit);
	}

/*
	start = omp_get_wtime();
	for (int i = 0; i < 1024*1024*128; i++) {
		for (int j = 0; j < 3; j++) {
			test_particle[j] = (int) (dsfmt_genrand_close_open(&dsfmt)*box_size);
		}
		check_defect_ltb(test_particle, xltb, yltb, zltb , 0);
	}
	stop = omp_get_wtime();
	printf("Time Method 3: %.2fs\n", stop-start);
*/

	printf("# Lookup table size M2: %10.1fkB %i %i\n",  3*ltb_N*sizeof(int) / ((float) (1024)), num_defects, box_size);

	// check if the lookup table is correct
	for (int i = 0; i < num_defects; i++) {

		int x = defects[i][0];
		int y = defects[i][1];
		int z = defects[i][2];
		//printf("Test: %i\n",lookup_table[x][y][z]);
/*
		assert(lookup_table[x][y][z] == 1); // Method 2
*/
		// Method 3
		int xi = x/ sizeof(int);
		int xbit = x % sizeof(int);
		assert( (xltb[xi] & (1<<xbit)) != 0 );

		int yi = y/ sizeof(int);
		int ybit = y % sizeof(int);
		assert( (yltb[yi] & (1<<ybit)) != 0 );

		int zi = z/ sizeof(int);
		int zbit  = z % sizeof(int);
		assert( (zltb[zi] & (1<<zbit)) != 0 );

	}

	/*
	 for (int i = 0; i < num_particles; i++) {
	 if (particles[i][3] == 1) printf("%i\n", i);
	 }
	 */

	/******************************* loop *********************************/
	// exchange outer with inner loop
	printf("\n# Starting ...\n");

	int  *mags = malloc(nsteps * sizeof(int)); // magnetization per step
	for (int i = 0; i < nsteps; i++) {
		mags[i] = 0;
	}

	// loop over particles
	double calc_time=0;
	printf("MinSize: %i\n", DSFMT_N64);
	double *dir_pool = malloc(DSFMT_N64 * sizeof(double));
#pragma omp parallel for reduction(+:calc_time) firstprivate(dir_pool)
	for (int i = 0; i < num_particles; i+=DSFMT_N64) {
		// every thread gets its own RNG
		dsfmt_t dsfmt;
		dsfmt_init_gen_rand(&dsfmt, i);

		/*
		 double *random_numbers_steps = malloc(nsteps * sizeof(double));
		 // create random numbers for the movements (directions 1..6)
		 dsfmt_fill_array_open_close(&dsfmt, random_numbers_steps, nsteps);

		 // scale the to 0,1,2,3,4,5 (the 6 directions)
		 for (int i=0 ;i < nsteps; i++) {
		 directions[i] = (short) (random_numbers_steps[i]*6);
		 }
		 free(random_numbers_steps);
		 */

		// distribution of correlation times, rexp,rnor are NOT thread safe!
		//for (int i=0 ;i < num_particles; i++) {
		//	correlation_times[i] = (int) rexp()*30;
		//}

		// loop over steps
		for (int step = 0; step < nsteps; step++) {
			dsfmt_fill_array_open_close(&dsfmt, dir_pool, DSFMT_N64);
			double start = omp_get_wtime();
			// doing batches of particles
			for (int j = 0; j < DSFMT_N64; j++) {
				int* particle = particles[i+j];
				int direction = dir_pool[j];
				//int direction = (int) (dsfmt_genrand_close_open(&dsfmt)*6);
				// only move particles which have not met defect yet == 0
				// or see how often they met a defefct >= 0
				if (particle[3] == 0) {

					// random step
					move_particle(particle, direction);

					// obey periodic boundary conditions, i.e. fold back
					check_pbc(particle, box_size);

					int tc = 10;

					// check_defect(particle, tc , hash_list, num_defects);
					// check_defect_tlb(particle, tc, hash_min, span, defekt_ltb);
					// check_defect_3d(particle, lookup_table, tc);
					check_defect_ltb(particle, xltb, yltb, zltb , tc);

					// ref_check_defect(particle, defects, num_defects);
				}
				else { // particle is trapped, decrease the residual waiting time
					particle[3] -= 1;
				}
				#pragma omp atomic
				mags[step] += particle[3];
				//gsl_histogram_increment (hvisits, particle[3]);
				//if (magnetization == num_particles) main_loop_break = 1;
				//int tid = omp_get_thread_num();
				//printf("Thread %i:  %i %i\n", id, i, direction);
				/*if (step%2000 == 0) {
				 printf("# Step: %8i (MAG: %5i)\r", step, magnetization);
				 fflush(stdout);
				 }*/
				//printf("%8i %8i %8i %8i\n",particle[0],particle[1],particle[2],particle[3]);

			} // end sub particle loop
			double stop = omp_get_wtime();
			calc_time += (stop-start);

		} // end steps loop
		/*
		 if (i%32 == 0) {
		 double stop = omp_get_wtime();
		 printf("# Particle: %8i (%8.3f s) Magnetization: %8i\r",i , (stop-start)/32 , magnetization);
		 fflush(stdout);
		 }
		 */
		// open the file we are writing to

		//#pragma omp critical
		//printf("# Particle: %8i (%8.3f s) \n",i , stop-start);
	} // end particle loop
	printf("Speed: %.2e s/particle \n", calc_time/num_particles);

	FILE *outFile;
	char fname[] = "binout.omp";
	sprintf(fname, "binout.om%i",0);
	outFile = fopen(fname, "w");
	// use fwrite to write binary data to the file
	fwrite(mags, sizeof(mags[0]), nsteps, outFile);
	fclose(outFile);
	print_array(particles, 10, 3);

	//print_array(particles, num_particles, 4);
	free(mags);
	free2D_i(particles);
	free2D_i(defects);
#ifdef STATS
	printf("Directions drawn:\n");
	gsl_histogram_fprintf (stdout, h, "%g", "%g");
	printf("\n");
	gsl_histogram_free (h);
#endif
	//omp_destroy_lock(&omplock);
	return 0;
}

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

/* qsort C-string comparison function */
int cstring_cmp(const void *a, const void *b)
{
	//    const char **ia = (const char **)a;
	//    const char **ib = (const char **)b;
	//    return strcmp(*ia, *ib);
	return strcmp ( (const char*)a, (const char*)b);

	/* strcmp functions works exactly as expected from
	 comparison function */
}

/* asm long long comparison function */
/*
 int asm64_comp(const void *a, const void *b) {
 int i=0;
 __asm__(
 "mov (%%rdi), %%rdx\n\t"    // Subtract low word
 "sub (%%rsi), %%rdx\n\t"
 "mov 8(%%rdi), %%rdi\n\t"    // Subtract high word
 "sbb 8(%%rsi), %%rdi\n\t"
 "sbb %%eax, %%eax\n\t"    // %eax = -1 if below, zero otherwise
 "or %%rdx, %%rdi\n\t"    // %rdi is non-zero if comparison is non-zero
 "neg %%rdi\n\t"    // carry flag is 1 if comparison is non-zero
 "adc %%eax, %%eax\n\t" // Result in %eax
 "movl %%eax, %0\n\t"
 : "=a" (i)
 :"r" (a), "r" (b)
 );
 return i;
 }
 */

int int64_cmp(const void *a, const void *b)
{
	const int64_t *x = a, *y = b;
	if(*x > *y)
		return 1;
	else
		return (*x < *y) ? -1 : 0;
}

/* qsort int comparison function */
int int_cmp(const void *a, const void *b)
{
	//    const int *ia = (const int *)a; // casting pointer types
	//    const int *ib = (const int *)b;
    //return *ia  - *ib;
	/* integer comparison: returns negative if b > a
	 and positive if a > b */
	return ( *(int*)a - *(int*)b );
}

static inline int64_t hash(int x, int y, int z) {
	return  (int64_t) MULTIPLIER* (int64_t) MULTIPLIER * (int64_t) x + (int64_t) MULTIPLIER * (int64_t) y + (int64_t) z;
}

int** malloc2D_i(long nrows, long ncolumns){
	int **array = malloc(nrows * sizeof(int *));
	array[0] = malloc(nrows * ncolumns * sizeof(int));
	if (array[0] == NULL) printf("Could not allocate memory");
	for(int i = 1; i < nrows; i++)
		array[i] = array[0] + i * ncolumns;

	// set all elements to 0
	for (int i = 0; i < nrows; i++) {
		for (int j = 0; j < ncolumns; j++) {
			array[i][j] = 0;
		}
	}
	return array;
}

/*
 char** malloc2D_char(long nrows, long ncolumns){
 char **array = malloc(nrows * sizeof(char *));
 array[0] = malloc(nrows * ncolumns * sizeof(char));
 if (array[0] == NULL) printf("Could not allocate memory");
 for(int i = 1; i < nrows; i++)
 array[i] = array[0] + i * ncolumns;
 return array;
 }
 */
char** malloc2D_char(long nrows, long ncolumns){
	char **array = malloc(nrows * sizeof(char *));
	for(int i = 0; i < nrows; i++)
		array[i] = malloc(ncolumns * sizeof(char));
	return array;
}

void free2D_i(int** array) {
	//free(&array[0]);
	free(array);
}

void print_array(int **array, int nrows, int ncolumns){
	for (int i = 0; i < nrows; i++) {
		for (int j = 0; j < ncolumns-1; j++) {
			printf("%i ",array[i][j]);
		}
		printf("%i\n",array[i][ncolumns-1]);
	}
}

void move_particle(int *particle, int direction){
	switch (direction) {
		case 0:
			particle[0] += 1;
			break;
		case 1:
			particle[0] -= 1;
			break;
		case 2:
			particle[1] += 1;
			break;
		case 3:
			particle[1] -= 1;
			break;
		case 4:
			particle[2] += 1;
			break;
		case 5:
			particle[2] -= 1;
			break;
	} // end switch statement
}

static inline void check_pbc(int* particle, int box_size) {
	for (int i = 0; i < 3; i++) {
		// % is NOT the mod operator, but the REMAINDER, it is not working for negative numbers (of course in C only)
		//particle[i] = particle[i] % box_size + (particle[i]<0?box_size:0);
		particle[i] &= (box_size - 1);
	}
}
/* binary search */
void check_defect(int* particle, int correlation_time, int64_t* hash_list, int num_defects ){
	int64_t hash_val;
	hash_val = hash(particle[0],particle[1],particle[2]);
	int * ptr;
	ptr = bsearch( &hash_val, hash_list, num_defects , sizeof(int64_t), int64_cmp);
	if (ptr != NULL)
		particle[3] = correlation_time;
}

/* lookup table 1 */

void check_defect_hash(int* particle, int correlation_time, int64_t hash_list_min, int64_t span, char* defekt_ltb){
	int64_t hash_val, offset;
	hash_val = hash(particle[0],particle[1],particle[2]);
	offset = hash_val -  hash_list_min;
	if ((offset >= 0) && (offset < span)) {
		if (defekt_ltb[ offset ] == 1)
			particle[3] = correlation_time;
	}
}

/* check Method 3 */

void check_defect_3d(int* particle, char*** lookup, int correlation_time){
	int x = particle[0];
	int y = particle[1];
	int z = particle[2];
	if (lookup[x] != NULL) {
		if (lookup[x][y] != NULL) {
			if (lookup[x][y][z] == 1) {
				particle[3] = correlation_time;
			}
		}
	}
}

void check_defect_ltb(int* particle, int* x, int* y, int* z, int correlation_time){
	int i = particle[0] / sizeof(int);
	int bit = 1<< particle[0] % sizeof(int);
	if (  (x[i] & 1 << bit) > 0) {
		i = particle[1] / sizeof(int);
		bit = 1<< particle[1] % sizeof(int);
		if (  (y[i] & 1 << bit) > 0) {
			i = particle[2] / sizeof(int);
			bit = 1<< particle[2] % sizeof(int);
			if (  (z[i] & 1 << bit) > 0) {
				particle[3] = correlation_time;
			}
		}
	}
		/*
	 switch ( (x[particle[0] / sizeof(int)]) & (1<< particle[0] % sizeof(int))  ) {
		 case 0:
			 break;
		 default:
			 switch ( (y[particle[1] / sizeof(int)]) & (1<< particle[1] % sizeof(int))  ) {
				 case 0:
					 break;
				 default:
					 switch ( (z[particle[2] / sizeof(int)]) & (1<< particle[2] % sizeof(int))  ) {
						 case 0:
							 break;
						 default:
							 particle[3] = correlation_time;
					 }
			 }
	 }*/
}

/*
 void check_defect(int* particle, int correlation_time, int64_t* hash_list, int num_defects ){
 int64_t hash_val;
 hash_val = hash(particle[0],particle[1],particle[2]);

 int bsearch = 0;
 int left = 0;
 int right = num_defects-1;

 while (bsearch == 0 && left <= right) {
 // int middle = (left + right) / 2;
 // better: avoid integer overflow
 int middle = left + (right-left) / 2;

 if (hash_val == hash_list[middle]) {
 bsearch = 1;
 particle[3] = correlation_time;
 }
 else {
 if (hash_val < hash_list[middle]) right = middle - 1;
 if (hash_val > hash_list[middle]) left = middle + 1;
 }
 }
 }

 */

// REFERENCE METHOD
void ref_check_defect(int* particle, int** defect_coords, int num_defects){
	int isDefect = 0;
	for (int i = 0; (i < num_defects) && (isDefect == 0); i++) {
		for (int j=0; j<3; j++ ){
			if (particle[j] != defect_coords[i][j]) {
				break; // coordinate mismatch, go to next particle (break loop over coordinates)
			}
			else {
				if (j==2) { // x,y and z ccordinate match
					particle[3] = (int) (rexp()*3) ; // set scalar value
					isDefect = 1; // break outer loop
				}
			}
		}
	}
}

/*

 void HT_check_defect(int* particle, Fnv64_t* hash_list, int num_defects ){
 int* pItem;
 char hash_string[18];
 Fnv64_t hash_val;
 snprintf(hash_string, sizeof(hash_string), "%5i %5i %5i", particle[0], particle[1], particle[2]);
 hash_val = fnv_64_str(hash_string, FNV0_64_INIT);
 pItem = (int*) bsearch (&hash_val, hash_list, num_defects, sizeof (Fnv64_t), fnv64_cmp);
 if (pItem != 0) particle[3] = 1; // set scalar value

 }

 */

/*
 // now create a hashed list to find them easier (hopefully)
 Fnv64_t hash_val;

 Fnv64_t *hash_list = malloc( num_defects*sizeof(Fnv64_t) );
 for (int i = 0; i < num_defects; i++) {
 snprintf(hash_string, sizeof(hash_string), "%5i %5i %5i", defects[i][0], defects[i][1], defects[i][2]);
 hash_val = fnv_64_str(hash_string, FNV1_64_INIT);
 hash_list[i] = hash_val;
 }
 qsort(hash_list, num_defects,sizeof(Fnv64_t), fnv64_cmp);
 */

/* mit chars
 char **hash_list = malloc2D_char(num_defects, 18);
 //char hash_list[20][18];
 for (int i = 0; i < num_defects; i++) {
 snprintf(hash_string, sizeof(hash_string), "%5i %5i %5i", defects[i][0], defects[i][1], defects[i][2]);
 hash_list[i] = hash_string;
 printf("%03i %s\n",i, hash_list[i]);
 }

 printf("Sorted\n");

 qsort(hash_list, num_defects, sizeof(hash_list[0]), cstring_cmp);//cmpstring_up);
 */

max_stars_repo_name: mrosenstihl/projects

max_issues_repo_name: mrosenstihl/projects

max_forks_repo_name: mrosenstihl/projects


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */
/*
 *
 *
 *       Filename:  hamil.c
 *
 *    Description:  Create an hamiltonian generator in Bloch form, from a given
 *    			vector omega
 *
 *        Version:  1.0
 *        Created:  03/05/2014 14:16:00
 *       Revision:  none
 *        License:  BSD
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 *
 */

#include "funcs.h"
#include <gsl/gsl_matrix.h>

/*
 *      FUNCTION
 *         Name:  ham_gen
 *  Description:
 *
 */
int ham_gen ( gsl_matrix* h, const double* o )
{
	gsl_matrix_set (h, 1, 2, o[3] ) ;
	gsl_matrix_set (h, 1, 3, -o[2] ) ;
	gsl_matrix_set (h, 2, 1, -o[3] ) ;
	gsl_matrix_set (h, 2, 3, o[1] ) ;
	gsl_matrix_set (h, 3, 1, o[2] ) ;
	gsl_matrix_set (h, 3, 2, -o[1] ) ;

	/*  Multiplies times -4 to obtain the Bloch matrix */
	gsl_matrix_scale ( h, -4 ) ;

	return 0;
}		/* -----  end of function ham_gen  ----- */

ext: c


text:
#ifndef __GSLEXTRA_H__
#define __GSLEXTRA_H__

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>

typedef struct gsl_quad_tensor
{
    int i;
    int j;
    int k;
    int l;

    gsl_matrix *** element;
}gsl_quad_tensor;

gsl_quad_tensor * gsl_quad_tensor_alloc(int i, int j, int k, int l);
gsl_quad_tensor * gsl_quad_tensor_calloc(int i,int j, int k, int l);
void gsl_quad_tensor_free(gsl_quad_tensor * q);

double gsl_quad_tensor_get(gsl_quad_tensor * q, int i, int j, int k, int l);
void gsl_quad_tensor_get_matrix(gsl_matrix * m, gsl_quad_tensor * q, int i, int j);

void gsl_quad_tensor_set(gsl_quad_tensor * q, int i, int j, int k, int l, double x);
void gsl_quad_tensor_set_matrix(gsl_quad_tensor * q, int i, int j, const gsl_matrix * m);

void gsl_quad_tensor_add(gsl_quad_tensor * A, gsl_quad_tensor * b);
void gsl_quad_tensor_sub(gsl_quad_tensor * A, gsl_quad_tensor * b);
void gsl_quad_tensor_scale(gsl_quad_tensor * q, double x);

void gsl_vector_complex_convert(gsl_vector * source, gsl_vector_complex * target, int length);
void gsl_matrix_complex_convert(gsl_matrix * source, gsl_matrix_complex * target, int rows, int columns);

void gsl_vector_complex_extract(gsl_vector_complex * source, gsl_vector * real, gsl_vector * imag, int length);
void gsl_matrix_complex_extract(gsl_vector_complex * source, gsl_matrix * real,gsl_matrix * imag, int rows, int columns);

void gsl_vector_complex_combine(gsl_vector * real, gsl_vector * imag, gsl_vector_complex * target);
void gsl_matrix_complex_combine(gsl_matrix * real,gsl_matrix * imag, gsl_matrix_complex * target);

void gsl_matrix_diag(gsl_matrix * target, gsl_vector * diag, int length);
void gsl_matrix_complex_diag(gsl_matrix_complex * target, gsl_vector_complex * diag, int length);

void gsl_matrix_mul(gsl_matrix * A, gsl_matrix *B, gsl_matrix * Result,int Acolumn,int Arow,int Bcolumn);
void gsl_matrix_complex_mul(gsl_matrix_complex * A, gsl_matrix_complex *B, gsl_matrix_complex * Result,int Acolumn,int Arow,int Bcolumn);

double gsl_vector_inner_product(gsl_vector * A, gsl_vector * B,int length);
gsl_complex gsl_vector_complex_product(gsl_vector_complex * A, gsl_vector_complex * B, int length);
gsl_complex gsl_vector_complex_inner_product(gsl_vector_complex * A, gsl_vector_complex * B, int length);
void gsl_vector_transform(gsl_vector * vec,gsl_matrix * trf,int length);
void gsl_vector_complex_transform(gsl_vector_complex * vec, gsl_matrix_complex * trf,int length);
void gsl_matrix_unitmatrix(gsl_matrix * m,int length);
void gsl_matrix_complex_unitmatrix(gsl_matrix_complex * m,int length);
void gsl_vector_complex_conjugate(gsl_vector_complex * v, int length);
void gsl_matrix_complex_conjugate(gsl_matrix_complex * m, int rows, int columns);

void gsl_matrix_square_root(gsl_matrix * dest, gsl_matrix * src, int length);
void gsl_matrix_inverse_square_root(gsl_matrix* dest, gsl_matrix * src, int length);
void gsl_eigen_Lowdin_diag(gsl_matrix * m, gsl_matrix * S, gsl_vector * eigen, gsl_matrix * eigenvec, int length);
void gsl_matrix_normalize(gsl_matrix * coef,int length, int columns);
void gsl_matrix_flip(gsl_matrix * dest,int height,int width);

#endif

ext: h

max_stars_repo_name: Walter-Feng/Overlap

max_issues_repo_name: Walter-Feng/Overlap

max_forks_repo_name: Walter-Feng/Overlap


text:
#include <stdio.h>
     #include <math.h>
     #include <gsl/gsl_sort.h>
     #include <gsl/gsl_wavelet.h>
     #include <gsl/gsl_math.h>

     int
     main (int argc, char **argv)
     {
       int i, j, k, n = 256, nc = 20;
       double *data = malloc (n * sizeof (double));
       double *abscoeff = malloc (n * sizeof (double));
       size_t *p = malloc (n * sizeof (size_t));

       FILE * f;
       gsl_wavelet *w;
       gsl_wavelet_workspace *work;

       w = gsl_wavelet_alloc (gsl_wavelet_daubechies, 4);
       work = gsl_wavelet_workspace_alloc (n);

       f = fopen (argv[1], "r");
       for (i = 0; i < n; i++)
         {
           fscanf (f, "%lg", &data[i]);
         }
       fclose (f);

       gsl_wavelet_transform_forward (w, data, 1, n, work);

       for (i = 0; i < n; i++)
         {
           abscoeff[i] = fabs (data[i]);
	   printf ("abscoeff[%d] = %g\n", i, abscoeff[i]);
         }

       printf ("(-1,0) = %g\n\n", abscoeff[0]);
       i = 1;
       for (j = 0; j < 8; j++)
	 {
	   for (k = 0; k < gsl_pow_int(2, j); k++)
	     {
	       printf ("(%d,%d) = %g\n", j, k, abscoeff[i]);
	       i++;
	     }
	   printf("\n");
	 }

       /* ************** */
       j = 0;
       for (i = 0; i < n; ++i)
	 j += abscoeff[i];
       for (i = 0; i < n; ++i)
	 abscoeff[i] /= j;
       /* ************** */

       gsl_sort_index (p, abscoeff, 1, n);

       for (i = 0; (i + nc) < n; i++)
	 {
	   printf ("p[%d] = %ld\n", i, p[i]);
	   data[p[i]] = 0;
	 }

       gsl_wavelet_transform_inverse (w, data, 1, n, work);

       for (i = 0; i < n; i++)
         {
           printf ("%g\n", data[i]);
         }

       gsl_wavelet_free (w);
       gsl_wavelet_workspace_free (work);

       free (data);
       free (abscoeff);
       free (p);
       return 0;
     }

max_stars_repo_name: klaricmn/snippets

max_issues_repo_name: klaricmn/snippets

max_forks_repo_name: klaricmn/snippets


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"

namespace Library
{
	class Mesh;
}

namespace Rendering
{
	class TexturedModelDemo final : public Library::DrawableGameComponent
	{
	public:
		TexturedModelDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);

		bool AnimationEnabled() const;
		void SetAnimationEnabled(bool enabled);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		struct CBufferPerObject
		{
			DirectX::XMFLOAT4X4 WorldViewProjection{ Library::MatrixHelper::Identity };
		};

		void CreateVertexBuffer(const Library::Mesh& mesh, gsl::not_null<ID3D11Buffer**> vertexBuffer) const;

		inline static const float RotationRate{ DirectX::XM_PI };

		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		CBufferPerObject mCBufferPerObject;
		winrt::com_ptr<ID3D11VertexShader> mVertexShader;
		winrt::com_ptr<ID3D11PixelShader> mPixelShader;
		winrt::com_ptr<ID3D11InputLayout> mInputLayout;
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		winrt::com_ptr<ID3D11Buffer> mConstantBuffer;
		winrt::com_ptr<ID3D11ShaderResourceView> mColorTexture;
		std::uint32_t mIndexCount{ 0 };
		float mRotationAngle{ 0.0f };
		bool mAnimationEnabled{ true };
		bool mUpdateConstantBuffer{ true };
	};
}


text:
#ifndef MSNHINFERENCECFG_H
#define MSNHINFERENCECFG_H
#include <stdint.h>
#include <float.h>
#include <string>
#include <vector>
#include <chrono>
#include "Msnhnet/utils/MsnhException.h"
#include "Msnhnet/config/MsnhnetMacro.h"
#include <math.h>
#include <string.h>

#ifdef USE_OMP
#include <omp.h>
#endif

#ifdef USE_NEON
#include <arm_neon.h>
#endif

#ifdef USE_OPEN_BLAS
#include <cblas.h>
#endif

#ifdef USE_OPENGL
#include <Msnhnet/config/MsnhnetOpenGL.h>
#endif

#ifndef OMP_THREAD
#define OMP_THREAD omp_get_max_threads()
#endif

#define CUDA_THREADS 512

#define MIN_OMP_DATA 10000

#define MSNHNET_VERSION 1200

#define EFFCIENT_ALIGN 16

enum ActivationType
{
    LOGISTIC,
    RELU,
    RELU6,
    RELIE,
    RAMP,
    TANH,
    PRELU,
    PLSE,
    LEAKY,
    ELU,
    LOGGY,
    STAIR,
    HARDTAN,
    LHTAN,
    SOFT_PLUS,
    SELU,
    SWISH,
    HARD_SWISH,
    MISH,
    NORM_CHAN,
    NORM_CHAN_SOFTMAX,
    NORM_CHAN_SOFTMAX_MAXVAL,
    NONE
};

enum LayerType
{
    CONVOLUTIONAL,
    DECONVOLUTIONAL,
    CONNECTED,
    MAXPOOL,
    LOCAL_AVGPOOL,
    GLOBAL_AVGPOOL,
    SOFTMAX,
    CROP,
    ROUTE,
    VARIABLE_OP,
    NORMALIZATION,
    AVGPOOL,
    ACTIVE,
    BATCHNORM,
    NETWORK,
    YOLO,
    YOLO_OUT,
    GAUSSIAN_YOLO,
    UPSAMPLE,
    L2NORM,
    EMPTY,
    VIEW,
    PERMUTE,
    PIXEL_SHUFFLE,
    SLICE,
    REDUCTION,
    CONFIG,
    RES_BLOCK,
    RES_2_BLOCK,
    CONCAT_BLOCK,
    ADD_BLOCK,
    PADDING
};

enum Arithmetic
{
    ARITH_ADD = 0,
    ARITH_SUB,
    ARITH_SUB_INV,
    ARITH_MUL,
    ARITH_DIV,
    ARITH_DIV_INV
};

enum Scientific
{
    SCI_ABS=0,
    SCI_ACOS,
    SCI_ASIN,
    SCI_ATAN,
    SCI_COS,
    SCI_COSH,
    SCI_SIN,
    SCI_SINH,
    SCI_TAN,
    SCI_TANH,
    SCI_EXP,
    SCI_POW,
    SCI_LOG,
    SCI_LOG10,
    SCI_SQRT
};

enum ReductionType
{
    REDUCTION_SUM,
    REDUCTION_MEAN
};

enum WeightsType
{
    NO_WEIGHTS,
    PER_FEATURE,
    PER_CHANNEL
};

enum WeightsNorm
{
    NO_NORM,
    RELU_NORM,
    SOFTMAX_NORM
};

#endif


text:
#define PyGSL_NUMERIC
#define PyGSL_NUMARRAY
#include "numeric_numarray.h"
#include <gsl/gsl_errno.h>
#include <pygsl/general_helpers.h>

static int
cp_numeric_ptrs(void)
{
     return GSL_SUCCESS;
}

static int
cp_numarray_ptrs(void)
{
     return GSL_SUCCESS;
}

void
initnumx(void)
{
     init_pygsl();
     cp_numeric_ptrs();
     cp_numarray_ptrs();
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#ifndef LIBTENSOR_CBLAS_H_H
#define LIBTENSOR_CBLAS_H_H

#ifdef USE_CBLAS
extern "C" { // Fixes older cblas.h versions without extern "C"
#include <cblas.h>
}
#endif // USE_CBLAS

#ifdef USE_GSL
#include <gsl/gsl_cblas.h>
#endif // USE_GSL

#endif // LIBTENSOR_CBLAS_H_H

max_stars_repo_name: pjknowles/libtensor

max_issues_repo_name: pjknowles/libtensor

max_forks_repo_name: pjknowles/libtensor


text:
#ifndef _B_SYSTEM_H_
#define _B_SYSTEM_H_

#include <iostream>
#include <vector>

using namespace std;

#include "Ode.h"
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

class B_system{
	public:
		B_system(double ini_vx_pass, double ini_vy_pass,
				 double ini_xx_pass, double ini_xy_pass,
				 double t_start_pass, double t_end_pass);
		~B_system();
		void solve();                  //solve the ODE and do the interpolation
		double vx(double t_now);       //return velocity in x direction at t_now
		double vy(double t_now);       //return velocity in y direction at t_now
		double xx(double t_now);       //return location in x direction at t_now
		double xy(double t_now);       //return location in y direction at t_now

	private:
		double ini_vx;                 //initial velocity in x direction
		double ini_vy;                 //initial velocity in x direction
		double ini_xx;                 //initial location in x direction
		double ini_xy;                 //initial location in x direction
		double t_start;                //the start time
		double t_end;                  //the end time
		int N = 1000;      			   //number of dots used in solving ODE
		ODE ode;                       //the initinization of the ODE class see the B_system constructor
		vector<double> ti_ode;         //this vector stores all
		vector<double> vx_ode;         //this vector stores all velocity in x dirction solved from ODE
		vector<double> vy_ode;         //this vector stores all velocity in y dirction solved from ODE
		gsl_interp_accel *vx_acc;      //workspace for interpolation loop up
		gsl_interp_accel *vy_acc;      //workspace for interpolation loop up
		gsl_spline *vx_spline;         //workspace for interpolation interface
		gsl_spline *vy_spline;         //workspace for interpolation interface
};

#endif


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"
#include "PointLight.h"

namespace Rendering
{
	class EnvironmentMappingMaterial;

	class EnvironmentMappingDemo final : public Library::DrawableGameComponent
	{
	public:
		EnvironmentMappingDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		EnvironmentMappingDemo(const EnvironmentMappingDemo&) = delete;
		EnvironmentMappingDemo(EnvironmentMappingDemo&&) = default;
		EnvironmentMappingDemo& operator=(const EnvironmentMappingDemo&) = default;
		EnvironmentMappingDemo& operator=(EnvironmentMappingDemo&&) = default;
		~EnvironmentMappingDemo();

		float AmbientLightIntensity() const;
		void SetAmbientLightIntensity(float intensity);

		float EnvironmentIntensity() const;
		void SetEnvironmentIntensity(float intensity);

		float ReflectionAmount() const;
		void SetReflectionAmount(float amount);

		virtual void Initialize() override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		inline static const float RotationRate{ DirectX::XM_PI };

		std::shared_ptr<EnvironmentMappingMaterial> mMaterial;
		DirectX::XMFLOAT4X4 mWorldMatrix{ Library::MatrixHelper::Identity };
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::uint32_t mIndexCount{ 0 };
		bool mUpdateMaterial{ true };
	};
}


text:
/*
 * accelerate_includes.h
 *
 *  Created on: Sep 29, 2014
 *      Author: Tama
 */

#ifndef ACCELERATE_INCLUDES_H_
#define ACCELERATE_INCLUDES_H_

#ifdef __linux__
/*LINUX IMPLEMENTATION*/
/*the libraries:

*liblapacke.a
*liblapack.a
*libcblas.a
*librefblas.a

are supposed to be in some standard place*/
//!!! The constant MKL_ENABLED must be defined at compile time
#ifdef CUDA_ENABLED
#include <cuda.h>
#include <cuda_runtime_api.h>
#include <cublas.h>
#include <cula_lapack.h>
#include <cula_lapack_device.h>
#include <curand.h>
#endif

#ifdef MKL_ENABLED
//#include <mkl_lapacke.h>
#include <mkl_cblas.h>
#include <mkl.h>

//#include <lapacke_mangling.h>
#else
#include "lapacke.h"
#include "cblas.h"
//#include <lapacke.h>
//#include <cblas.h>
//#include <lapacke_mangling.h>
#endif

typedef lapack_int integer;

typedef float __lpkreal;
typedef double __lpkdoublereal;
typedef lapack_complex_float __lpkcomplex;
typedef lapack_complex_double __lpkdoublecomplex;

template <typename T>
struct constants{
//	const __lpkcomplex complex_one ={1.f,0.f};
//		const __lpkcomplex complex_zero = {0.f,0.f};
//		const __lpkreal real_one = 1.f;
//		const __lpkreal real_zero = 0.f;
//		const __lpkreal typed_one = 1.f;
//		const __lpkreal typed_zero = 0.f;
//		const __lpkreal typed_mone = -1.f;
};

template <>
struct constants<__lpkreal>{
	const __lpkcomplex complex_one ={1.f,0.f};
	const __lpkcomplex complex_zero = {0.f,0.f};
	const __lpkreal real_one = 1.f;
	const __lpkreal real_zero = 0.f;
	const __lpkreal typed_one = 1.f;
	const __lpkreal typed_zero = 0.f;
	const __lpkreal typed_mone = -1.f;
	const char typeCheck = 's';
};

template <>
struct constants<__lpkdoublereal>{
	const __lpkdoublecomplex complex_one ={1.,0.};
	const __lpkdoublecomplex complex_zero = {0.,0.};
	const __lpkdoublereal real_one = 1.;
	const __lpkdoublereal real_zero = 0.;
	const __lpkdoublereal typed_one = 1.;
	const __lpkdoublereal typed_zero = 0.;
	const __lpkdoublereal typed_mone = -1.;
	const char typeCheck = 'd';
};

template <>
struct constants<__lpkcomplex>{
	const __lpkcomplex complex_one ={1.f,0.f};
	const __lpkcomplex complex_zero = {0.f,0.f};
	const __lpkreal real_one = 1.f;
	const __lpkreal real_zero = 0.f;
	const __lpkcomplex typed_one = {1.f,0.f};
	const __lpkcomplex typed_zero = {0.f,0.f};
	const __lpkcomplex typed_mone = {-1.f,0.f};
	const char typeCheck = 'c';
};

template <>
struct constants<__lpkdoublecomplex>{
	const __lpkdoublecomplex complex_one ={1.,0.};
	const __lpkdoublecomplex complex_zero = {0.,0.};
	const __lpkdoublereal real_one = 1.;
	const __lpkdoublereal real_zero = 0.;
	const __lpkdoublecomplex typed_one = {1.,0.};
	const __lpkdoublecomplex typed_zero = {0.,0.};
	const __lpkdoublecomplex typed_mone = {-1.,0.};
const char typeCheck = 'z';
};

//#ifdef SINGLEPREC
//typedef float lpkreal;
//typedef lapack_complex_float lpkcomplex;
//const lpkcomplex complex_one = {1.f,0.f};
//const lpkcomplex complex_zero = {0.f,0.f};
//const lpkreal real_one = 1.0f;
//const lpkreal real_zero = 0.0f;
////Real Norm 2 single prec
//#define NORM2REAL cblas_snrm2
//
//#ifdef REAL
////SINGLE PREC REAL
//const lpkreal typed_one = 1.0f;
//const lpkreal typed_zero = 0.0f;
//const lpkreal typed_mone = -1.0f;
//#define GEQRF_ sgeqrf_
//#define UNGQR_ sungqr_
//#define GEMM cblas_sgemm
//#define AXPY cblas_saxpy
//#define DOT cblas_sdotc
//#define COPY cblas_scopy
//#define GESVD_ sgesvd_
//#define DAG CblasTrans
//
//#else
//
//const lpkcomplex typed_one = {1.0f,0.f};
//const lpkcomplex typed_mone = {-1.0f,0.f};
//const lpkcomplex typed_zero = {0.0f,0.0f};
////SINGLE PREC COMPLEX
//#define GEQRF_ cgeqrf_
//#define UNGQR_ cungqr_
//#define GEMM cblas_cgemm
//#define AXPY cblas_caxpy
//#define DOT cblas_cdotc_sub
//#define COPY cblas_ccopy
//#define GESVD_ cgesvd_
//#define DAG CblasConjTrans
//#endif //REAL
//
//
//#else //DOUBLEPREC
//typedef double  lpkreal;
//typedef lapack_complex_double lpkcomplex;
//const lpkcomplex complex_one = {1.,0.};
//const lpkcomplex complex_zero = {0.,0.};
//const lpkreal real_one = 1.0;
//const lpkreal real_zero = 0.0;
//
////Real Norm 2 double prec
//#define NORM2REAL cblas_dnrm2
//#ifdef REAL
////DOULBE PREC REAL
//const lpkreal typed_one = 1.0;
//const lpkreal typed_mone = -1.0;
//const lpkreal typed_zero = 0.0;
//#define GEQRF_ dgeqrf_
//#define UNGQR_ dungqr_
//#define GEMM cblas_dgemm
//#define AXPY cblas_daxpy
//#define DOT cblas_ddotc
//#define COPY cblas_dcopy
//#define GESVD_ dgesvd_
//#define DAG CblasTrans
//#else
////DOUBLE PREC COMPLEX
//const lpkcomplex typed_one = {1.0,0.};
//const lpkcomplex typed_mone = {-1.0,0.};
//const lpkcomplex typed_zero = {0.0,0.0};
//
//
//#define GEQRF_ zgeqrf_
//#define UNGQR_ zungqr_
//#define GEMM cblas_zgemm
//#define AXPY cblas_zaxpy
//#define DOT cblas_zdotc_sub
//#define COPY cblas_zcopy
//#define GESVD_ zgesvd_
//#define DAG CblasConjTrans
//#endif
//
//#endif //END DOUBLE PREC

#endif //LINUX

#ifdef __APPLE__
#include "/System/Library/Frameworks/Accelerate.framework/Headers/Accelerate.h"
//#include <Accelerate.framework/Headers/Accelerate.h>
#include "/System/Library/Frameworks/Accelerate.framework/Frameworks/vecLib.framework/Headers/clapack.h"
#include "/System/Library/Frameworks/Accelerate.framework/Frameworks/vecLib.framework/Headers/cblas.h"

/******APPLE IMPLEMENTATION******/
typedef __CLPK_integer integer;

typedef __CLPK_real __lpkreal;
typedef __CLPK_doublereal __lpkdoublereal;
typedef __CLPK_complex __lpkcomplex;
typedef __CLPK_doublecomplex __lpkdoublecomplex;

template <typename T>
struct constants{
//	const __lpkcomplex complex_one ={1.f,0.f};
//		const __lpkcomplex complex_zero = {0.f,0.f};
//		const __lpkreal real_one = 1.f;
//		const __lpkreal real_zero = 0.f;
//		const __lpkreal typed_one = 1.f;
//		const __lpkreal typed_zero = 0.f;
//		const __lpkreal typed_mone = -1.f;
};

template <>
struct constants<__lpkreal>{
	const __lpkcomplex complex_one ={1.f,0.f};
	const __lpkcomplex complex_zero = {0.f,0.f};
	const __lpkreal real_one = 1.f;
	const __lpkreal real_zero = 0.f;
	const __lpkreal typed_one = 1.f;
	const __lpkreal typed_zero = 0.f;
	const __lpkreal typed_mone = -1.f;
	const char typeCheck = 's';
};

template <>
struct constants<__lpkdoublereal>{
	const __lpkdoublecomplex complex_one ={1.,0.};
	const __lpkdoublecomplex complex_zero = {0.,0.};
	const __lpkdoublereal real_one = 1.;
	const __lpkdoublereal real_zero = 0.;
	const __lpkdoublereal typed_one = 1.;
	const __lpkdoublereal typed_zero = 0.;
	const __lpkdoublereal typed_mone = -1.;
	const char typeCheck = 'd';
};

template <>
struct constants<__lpkcomplex>{
	const __lpkcomplex complex_one ={1.f,0.f};
	const __lpkcomplex complex_zero = {0.f,0.f};
	const __lpkreal real_one = 1.f;
	const __lpkreal real_zero = 0.f;
	const __lpkcomplex typed_one = {1.f,0.f};
	const __lpkcomplex typed_zero = {0.f,0.f};
	const __lpkcomplex typed_mone = {-1.f,0.f};
	const char typeCheck = 'c';
};

template <>
struct constants<__lpkdoublecomplex>{
	const __lpkdoublecomplex complex_one ={1.,0.};
	const __lpkdoublecomplex complex_zero = {0.,0.};
	const __lpkdoublereal real_one = 1.;
	const __lpkdoublereal real_zero = 0.;
	const __lpkdoublecomplex typed_one = {1.,0.};
	const __lpkdoublecomplex typed_zero = {0.,0.};
	const __lpkdoublecomplex typed_mone = {-1.,0.};
const char typeCheck = 'z';
};

//#ifdef SINGLEPREC
//typedef __CLPK_real lpkreal;
//typedef __CLPK_complex lpkcomplex;
//const lpkcomplex complex_one = {1.f,0.f};
//const lpkcomplex complex_zero = {0.f,0.f};
//const lpkreal real_one = 1.0f;
//const lpkreal real_zero = 0.0f;
////Real Norm 2 single prec
//#define NORM2REAL cblas_snrm2
//
//#ifdef REAL
////SINGLE PREC REAL
//const lpkreal typed_one = 1.0f;
//const lpkreal typed_zero = 0.0f;
//const lpkreal typed_mone = -1.0f;
//
//#define GEQRF_ sgeqrf_
//#define UNGQR_ sungqr_
//#define GEMM cblas_sgemm
//#define AXPY cblas_saxpy
//#define DOT cblas_sdotc
//#define COPY cblas_scopy
//#define GESVD_ sgesvd_
//
//#else
//
//const lpkcomplex typed_one = {1.0f,0.f};
//const lpkcomplex typed_mone = {-1.0f,0.f};
//const lpkcomplex typed_zero = {0.0f,0.0f};
////SINGLE PREC COMPLEX
//#define GEQRF_ cgeqrf_
//#define UNGQR_ cungqr_
//#define GEMM cblas_cgemm
//#define AXPY cblas_caxpy
//#define DOT cblas_cdotc_sub
//#define COPY cblas_ccopy
//#define GESVD_ cgesvd_
//#endif //REAL

//#else //DOUBLEPREC
//typedef __CLPK_doublereal lpkreal;
//typedef __CLPK_doublecomplex lpkcomplex;
//const lpkcomplex complex_one = {1.,0.};
//const lpkcomplex complex_zero = {0.,0.};
//const lpkreal real_one = 1.0;
//const lpkreal real_zero = 0.0;
//
////Real Norm 2 double prec
////#define NORM2REAL cblas_dnrm2
////#ifdef REAL
////DOULBE PREC REAL
//template <>
//struct constants<lpkreal>{
//	const lpkcomplex complex_one ={1.,0.};
//	const lpkcomplex complex_zero = {0.,0.};
//	const lpkreal real_one = 1.;
//	const lpkreal real_zero = 0.;
//	const lpkreal typed_one = 1.;
//	const lpkreal typed_zero = 0.;
//	const lpkreal typed_mone = -1.;
//	const char typeCheck = 'd';
//};
//const lpkreal typed_one = 1.0;
//const lpkreal typed_mone = -1.0;
//const lpkreal typed_zero = 0.0;
//#define GEQRF_ dgeqrf_
//#define UNGQR_ dorgqr_
//#define GEMM cblas_dgemm
//#define AXPY cblas_daxpy
//#define DOT cblas_ddotc
//#define COPY cblas_dcopy
//#define GESVD_ dgesvd_

//#else
//DOUBLE PREC COMPLEX
//Questa struttura dovra` essere modificata alla fine
//template<>
//struct constants<lpkcomplex>{
//	const lpkcomplex complex_one ={1.,0.};
//	const lpkcomplex complex_zero = {0.,0.};
//	const lpkreal real_one = 1.;
//	const lpkreal real_zero = 0.;
//	const lpkcomplex typed_one = {1.,0.};
//	const lpkcomplex typed_zero = {0.,0.};
//	const lpkcomplex typed_mone = {-1.,0.};
//	const char typeCheck = 'z';
//};

//const lpkcomplex typed_one = {1.0,0.};
//const lpkcomplex typed_mone = {-1.0,0.};
//const lpkcomplex typed_zero = {0.0,0.0};

//#define GEQRF_ zgeqrf_
//#define UNGQR_ zungqr_
//#define GEMM cblas_zgemm
//#define AXPY cblas_zaxpy
//#define DOT cblas_zdotc_sub
//#define COPY cblas_zcopy
//#define GESVD_ zgesvd_

//#endif

//#endif //END DOUBLE PREC
#endif //APPLE
/******Tama's constants******/

#endif /* ACCELERATE_INCLUDES_H_ */

max_stars_repo_name: kindaguy/RRSVD

max_issues_repo_name: kindaguy/RRSVD

max_forks_repo_name: kindaguy/RRSVD


text:
#include <config.h>
#include <math.h>
#include <gsl/gsl_ntuple.h>
#include <gsl/gsl_histogram.h>

struct data
{
  double x;
  double y;
  double z;
};

int sel_func (void *ntuple_data, void *params);
double val_func (void *ntuple_data, void *params);

int
main (void)
{
  struct data ntuple_row;
  int i;

  gsl_ntuple *ntuple = gsl_ntuple_open ("test.dat", &ntuple_row,
                                        sizeof (ntuple_row));

  gsl_histogram *h = gsl_histogram_calloc_uniform (100, 0., 10.);

  gsl_ntuple_select_fn S;
  gsl_ntuple_value_fn V;

  double scale = 1.5;

  S.function = &sel_func;
  S.params = &scale;

  V.function = &val_func;
  V.params = 0;

  gsl_ntuple_project (h, ntuple, &V, &S);

  gsl_histogram_fprintf (stdout, h, "%f", "%f");

  gsl_histogram_free (h);

  gsl_ntuple_close (ntuple);
}

int
sel_func (void *ntuple_data, void *params)
{
  double x, y, z, E, scale;
  scale = *(double *) params;

  x = ((struct data *) ntuple_data)->x;
  y = ((struct data *) ntuple_data)->y;
  z = ((struct data *) ntuple_data)->z;

  E = x * x + y * y + z * z;

  return E / scale > 1;
}

double
val_func (void *ntuple_data, void *params)
{
  double x, y, z;

  x = ((struct data *) ntuple_data)->x;
  y = ((struct data *) ntuple_data)->y;
  z = ((struct data *) ntuple_data)->z;

  return x * x + y * y + z * z;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#ifndef ALG_LOMV
#define ALG_LOMV

#include "string.h"
#include "assert.h"
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>

double linf(int d, gsl_vector * a, gsl_vector * b);

void psi(int p, int q, gsl_vector * theta, gsl_matrix * B, gsl_matrix * V, gsl_vector * Delta, gsl_vector * x);

gsl_vector * ffp(int p, int q, gsl_vector * theta, gsl_matrix * B, gsl_matrix * V, gsl_vector * Delta);

gsl_vector * lo_minvar(int p, int q, gsl_matrix * B, gsl_matrix * V, gsl_vector * Delta);

double * ffp_C_interface(int p, int q, double* theta, double** B, double** V, double* Delta);

double* lo_minvar_C_interface(int p, int q, double**B, double** V, double* Delta);

double* psi_C_interface(int p, int q, double**B, double** V, double* Delta);

#endif


text:
/* Copyright (c) 2020 Stijn Hinterding, Utrecht University
 * This sofware is licensed under the MIT license (see the LICENSE file)
*/

#ifndef FAKECOMMUNICATOR_H
#define FAKECOMMUNICATOR_H

#include <phast_gui/interfaces/itimetaggercommunicator.h>

#include <gsl/gsl_rng.h>
#include <chrono>
#include <mutex>

class FakeCommunicator : public ITimeTaggerCommunicator
{
private:
    gsl_rng* r;
    double timeunit_seconds;
    double decayrate_inv_seconds;
    double pulse_period_seconds;
    double detection_probability;

    std::chrono::time_point<std::chrono::high_resolution_clock> prev_time;
    bool first_run;

    int64_t last_largest_time;
    int64_t last_pulse_time;
    int64_t pulse_period_timeunits;
    int64_t last_pulse_published;
    uint64_t sync_counter;
    uint64_t sync_divider;

    int64_t chan1_induced_delay;
    int64_t chan2_induced_delay;

    std::string descriptor;

    bool chan0_enabled;
    bool chan1_enabled;
    bool chan2_enabled;

    bool is_running;

    double sim_speedup_factor;
    std::mutex m;

    std::map<chan_id,timestamp> first_times;
    std::map<chan_id,timestamp> last_times;
    std::map<chan_id,uint64_t> n_events;

    int64_t n_emitters;
private:
    void gen_pulses(std::vector<int64_t> *timestamps,
                    std::vector<uint8_t> *channel_IDs, int64_t duration);

public:
    FakeCommunicator(const FakeCommunicator&) = delete;
    FakeCommunicator& operator=(const FakeCommunicator&) = delete;
    ~FakeCommunicator() override;

    FakeCommunicator(double timeunit_seconds,
                     double decayrate_inv_seconds,
                     double pulse_period_seconds,
                     double detection_probability,
                     uint64_t sync_divider=128,
                     int64_t chan1_induced_delay=0,
                     int64_t chan2_induced_delay=0,
                     int64_t n_emitters=1,
                     double sim_speedup_factor=1);

    void SetSimSpeedupFactor(double value)
    {
        std::lock_guard<std::mutex> l(this->m);
        if (value <= 0.001)
            value = 0.001;

        this->sim_speedup_factor = value;
    }

    double GetSimSpeedupFactor() const
    {
        return this->sim_speedup_factor;
    }

    void SetDecayRate(double value)
    {
        std::lock_guard<std::mutex> l(this->m);
        this->decayrate_inv_seconds = value;
    }

    int64_t GetNumEmitters() const
    {
        return this->n_emitters;
    }

    double GetDecayRate() const
    {
        return this->decayrate_inv_seconds;
    }

    void SetNumEmitters(int64_t value)
    {
        std::lock_guard<std::mutex> l(this->m);

        if (value == 0)
            value = 1;

        this->n_emitters = value;
    }
    void SetPulsePeriod(double value)
    {
        std::lock_guard<std::mutex> l(this->m);
        this->pulse_period_seconds = value;
        this->pulse_period_timeunits = this->pulse_period_seconds/this->timeunit_seconds;
    }

    double GetPulsePeriod() const
    {
        return this->pulse_period_seconds;
    }

    void SetDetectionProbability(double value)
    {
        std::lock_guard<std::mutex> l(this->m);
        if (value > 1.0)
            value = 1.0;

        if (value < 0.0)
            value = 0.0;

        this->detection_probability = value;
    }

    double GetDetectionProbability() const
    {
        return this->detection_probability;
    }

    void SetSyncDivider(uint64_t value)
    {
        std::lock_guard<std::mutex> l(this->m);
        this->sync_divider = value;
    }

    virtual uint64_t ReceiveData(std::vector<int64_t>* timestamps,
                             std::vector<uint8_t>* channel_IDs) override;

    virtual uint64_t GetSyncDivider(chan_id channel_id) override
    {
        if (channel_id == 0)
            return this->sync_divider;

        return 1;
    }

    virtual void DisableChan(uint64_t ID);

    virtual bool ConnectedToDevice() override
    {
        return true;
    }

    virtual const std::string& DeviceDescriptor() const override
    {
        return this->descriptor;
    }

    virtual bool TryToConnect(int64_t) override
    {
        return true;
    }

    virtual bool DataLossSinceLastCall() override
    {
        return false;
    }

    virtual double TimeUnit() const override
    {
        return this->timeunit_seconds;
    }

    virtual bool AnyDeviceAvailable() const override
    {
        return true;
    }

    virtual bool Disconnect() override
    {
        return true;
    }

    virtual uint64_t GetNumDevicesConnected() override
    {
        return 1;
    }

    virtual bool IsRealDevice() const
    {
        return false;
    }
};

#endif // FAKECOMMUNICATOR_H

max_stars_repo_name: stijnhinterding/phast

max_issues_repo_name: stijnhinterding/phast

max_forks_repo_name: stijnhinterding/phast


text:
//////////////////////////////////////////// Spin glass phase transition ////////////////////////////////////////

/*

Phase space of the SK model.

Please cite the following paper when you use this code.

[Ezaki T, Fonseca dos Reis E, Watanabe T, Sakaki M, Masuda N. Closer to critical resting-state neural dynamics in individuals with higher fluid intelligence. Commun Biol 3:1 (2020).](https://www.nature.com/articles/s42003-020-0774-y)

If you find a bug in this code, please contact the authors.

Author: Elohim Fonseca dos Reis

Date: July 8, 2020

Parameters:
----------
   This program maps the SK model phase space using the following parameters that have to be set by the user before compiling:
       * Total number of spins (N);
       * Number of time steps to calculate the averages given a fixed configuration of $J_{ij}$ (tdim);
       * Number of $J_{ij}$ configurations for the configurational average (conf_num);
       * Number of sweeps in the transient period (thermal);
       * Array with the mean values of $J_{ij}$ (mu). The user has to set the minimum value (mu_min), maximum value (mu_max) and the step (mu_step). The size of the array with the mean values of $J_{ij}$ is given by 1+(mu_max-mu_min)/mu_step.
       * Array with the standard deviation values of $J_{ij}$ (sd). The user has to set the minimum value (sd_min), the maximum value (sd_max) and the step (sd_step). The size of the array with the standard deviation values of $J_{ij}$ is given by 1+(max_sd-min_sd)/sd_step.

Outputs:
-------
The program outputs the following measures:
 * spin glass susceptibility (Xsg);
 * uniform susceptibility (Xuni);
 * spin glass order parameter (q);
 * magnetization (m);
 * specific heat (c).

The outputs of the program are multiple .txt files, one for each pair of (mu, sd), and each .txt file contains one line with seven values in this
order: mu, sd, Xsg, Xuni, q, m, c. For example, if the user sets *n* values for the mu array and *m* values for the sd array, the program outputs *nm* .txt files. The files are named as: file_0_0.txt, file_0_1.txt, ..., file_1_0.txt, file_1_1.txt... The reason for generating multiple .txt files is for parallelizing the program.

The program has a built-in naive parallelization for the mu and sd arrays, i.e., after compiling the code, if the user runs the same program on 10 different instances,
for example, the total time will decrease 10 times. So the user might be interested in running an array of jobs, each job running the same program for a given set of parameters.

*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <unistd.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include "mt19937ar.h"

// Parameters
#define N 264          // Total number of spins
#define tdim 10000     // Thermal average dimension size
#define conf_num 1000  // Number of interaction configurations
#define thermal 10000  // Equilibration sweeps

#define mu_min -0.002   // Minimum Average Interaction
#define mu_max 0.01     // Maximum Average Interaction
#define mu_step 0.0005  // Average interaction step
#define sd_min 0        // Minimum interaction standard deviation
#define sd_max 0.15     // Maximum interaction standard deviation
#define sd_step 0.0075  // Interaction sd step

// Matrix Serializations
#define SSArr(SS,i,j) SS[ tdim*i + j ] // Spin series array
#define JArr(J,i,j) J[ N*i + j ] // Interactions array
#define CovArr(var,i,j) var[N*i + j] // Covariance array

// Global Variables
int *s; // Spins
int *SS; // Spin series array (N x T)
double *J; // Interactions (Jij) drawn from a gaussian distribution

// Functions
double randnum(); // Real random number [0,1)
void init_randnum(); // Seed of the random number generator
void sweep(); // Sweep of the system with N metropolis steps
double cov(int *A, int k, int l); // Covariance cov[A(k),A(l)]
void elapsed_time();

int main()
{

///////////////////////////// Summary /////////////////////////////
	printf("\n\nSummary:\n\n");
	printf("Number of spins: %d\n", N);
	printf("Thermal average dimension: %d\n",tdim);
	printf("Number of configurations: %d\n",conf_num);
	printf("Equilibration sweeps: %d\n",thermal);
	printf("Number of parameter values: %d mean(J),  %d std(J)", (int)round(((mu_max - mu_min) / mu_step + 1)), (int)round(((sd_max - sd_min) / sd_step + 1)));
	printf("\n\n");

///////////////////////////// Model Parameters /////////////////////////////

	/* Average interaction array (main function variable) */
	int mu_size = (int)round( (mu_max - mu_min) / mu_step + 1); //printf("\n\nmu size: %d\n\n", mu_size);
	double *mu=malloc(mu_size*sizeof(double));
	if(!mu){
		printf("Out of memory.1\n");
		exit(1);
	}
	for(int i=0; i<mu_size; i++){
		mu[i] = (mu_min + i*mu_step); //printf("m[%d]=%f\n\n",i,mu[i]);
	}

	/* Interaction standard deviation array (main function variable) */
	int sd_size = (int)round( (sd_max - sd_min) / sd_step + 1); //printf("\n\nsd size: %d\n\n", sd_size);
	double *sd=malloc(sd_size*sizeof(double));
	if(!sd){
		printf("Out of memory.2\n");
		exit(1);
	}
	for(int i=0; i<sd_size; i++){
		sd[i] = (sd_min + i*sd_step); //printf("sd[%d]=%f\n\n",i,sd[i]);
	}

//////////////////////////////// Variables ///////////////////////////////

	/* Spin array (global) */
	s=malloc(N*sizeof(int));
	if(!s){
		printf("Out of memory.\n");
		exit(1);
	}

	/* Spin series array (global) */
	SS=malloc(N*tdim*sizeof(int));
	if(!SS){
		printf("Out of memory.\n");
		exit(1);
	}

	/* Interaction matrix (global) */
	J=malloc(N*N*sizeof(double));
	if(!J){
		printf("Out of memory.\n");
		exit(1);
	}

	/* Covariance matrix configuration average sum */
	double *SumC=calloc(N*N,sizeof(double));
	if(!SumC){
		printf("Out of memory.\n");
		exit(1);
	}

	/* Covariance matrix configuration average product sum */
	double *ProdC=calloc(N*N,sizeof(double));
	if(!ProdC){
		printf("Out of memory.\n");
		exit(1);
	}

	/* Magnetization */
	double *m=malloc(mu_size*sd_size*sizeof(double));
	if(!m){
		printf("Out of memory.\n");
		exit(1);
	}

	/* SG order parameter */
	double *q=malloc(mu_size*sd_size*sizeof(double));
	if(!q){
		printf("Out of memory.\n");
		exit(1);
	}

	/* Uniform Susceptibiliy */
	double *Xuni=calloc(mu_size*sd_size,sizeof(double));
	if(!Xuni){
		printf("Out of memory.\n");
		exit(1);
	}

	/* Spin Glass Susceptibiliy */
	double *Xsg=calloc(mu_size*sd_size,sizeof(double));
	if(!Xsg){
		printf("Out of memory.\n");
		exit(1);
	}

	/* Specific Heat */
	double *c=calloc(mu_size*sd_size,sizeof(double));
	if(!c){
		printf("Out of memory.\n");
		exit(1);
	}

/////////////////////////////// Program ///////////////////////////////

	/* Initialize mt rng */
	init_randnum();

	/* Interaction array variables */
	const gsl_rng_type * Q;
	gsl_rng * r;
	Q = gsl_rng_mt19937;
	r = gsl_rng_alloc (Q);

	//////////////////////////////////////////////////////////////////////////////////
	/////////////////////////////// AVERAGE LOOP (A) /////////////////////////////////
	//////////////////////////////////////////////////////////////////////////////////
	for(int A=0; A<mu_size; A++){

	//////////////////////////////////////////////////////////////////////////////////
	///////////////////////// STANDARD DEVIATION LOOP (B) ////////////////////////////
	//////////////////////////////////////////////////////////////////////////////////
	for(int B=0; B<sd_size; B++){

		// Output file
		FILE *file;
		char *fname=malloc(100*sizeof(char));
		snprintf(fname,100,"file_%d_%d.txt",A,B);

		// Check if file exists
		if( access( fname, F_OK ) != -1 ) {
		    // file exists
		    continue;
		}

		// file doesn't exist: perform calculations

		// Save blank file
		if((file = fopen(fname,"w"))==NULL){
			printf("Cannot open file.\n");
		        exit(1);
		}
		fclose(file);
		printf("\nfile_%d_%d.txt  mu=%f  sd=%f",A,B, mu[A], sd[B]);

		//~~~~~~~~~ Beginning of calculations

		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		// Initial condition
		for(int i=0; i<N; i++) s[i]=1;

		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		// Reset convariance array
		memset(SumC, 0, N*N*sizeof(*SumC));
		memset(ProdC, 0, N*N*sizeof(*ProdC));

		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		// Reset specific heat variables
		double ProdSumE = 0;
		double SumProdE = 0;

		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		// Configuration Average
		for(int D=0; D<conf_num; D++){

			//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~;
			// Interaction array
			for(int i=0; i<N; i++){
			for(int j=i; j<N; j++){
				JArr(J,i,j) = mu[A] + gsl_ran_gaussian_ziggurat(r,sd[B]);
				JArr(J,j,i) = JArr(J,i,j);
				JArr(J,i,i) = 0;
			}}

			//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
			// Thermalization
			for(int i=0; i<thermal; i++){
				sweep();
			}

			//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
			// Sampling
			for (int j=0; j<tdim; j++){
				sweep();
				for (int i=0; i<N; i++){
					SSArr(SS,i,j)=s[i];
				}
			}

			//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
			// Magnetization and SG order parameter
			double mag=0;
			double sgop=0;
			double aux=0;
			for(int i=0; i<N; i++){
				aux=0;
				for(int t=0; t<tdim; t++){
					aux += SSArr(SS,i,t);
				}
				mag += aux;
				sgop += aux*aux;
			}
			m[sd_size*A + B] += mag;
			q[sd_size*A + B] += sgop;

			//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
			// Covariance Matrix: configuration average sum (SumC) for Xuni and product (ProdC) for Xsg
			double cov_ij;
			for(int i=0; i<N; i++){
			for(int j=i; j<N; j++){
				cov_ij = cov(SS,i,j);
				CovArr(SumC,i,j) += cov_ij;
				CovArr(SumC,j,i) = CovArr(SumC,i,j);
				CovArr(ProdC,i,j) += cov_ij*cov_ij;
				CovArr(ProdC,j,i) = CovArr(ProdC,i,j);
			}}

			//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
			// Specific heat:
			double sumE=0, sumE2=0, E=0;
			for(int t=1; t<tdim; t++){
				E=0;
				for(int i=0; i<N; i++){
				for(int j=i; j<N; j++){
					E -= JArr(J,i,j) * SSArr(SS,i,t) *SSArr(SS,j,t);
				}}
				sumE += E;
				sumE2 += E*E;
			}
			ProdSumE += sumE*sumE;
			SumProdE += sumE2;

		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		} // End of Configuration Average

		m[sd_size*A+B] = fabs(m[sd_size*A+B]/(1.0*N*tdim*conf_num));
		q[sd_size*A+B] = q[sd_size*A+B]/(1.0*N*tdim*tdim*conf_num);

		int NN=N*N;
		for(int i=0; i<NN; i++){
			Xuni[sd_size*A+B] += SumC[i];
			Xsg[sd_size*A+B] += ProdC[i];
		}
		Xuni[sd_size*A+B] = Xuni[sd_size*A+B]/(1.0*N*conf_num);
		Xsg[sd_size*A+B] = Xsg[sd_size*A+B]/(1.0*N*conf_num);

		c[sd_size*A+B] = (SumProdE/(1.0*tdim) - ProdSumE/(1.0*tdim*tdim))/(1.0*N*conf_num);

		// Save file with results
		if((file = fopen(fname,"w"))==NULL){
			printf("Cannot open file.\n");
        		exit(1);
        	}

		fprintf(file,"%f\t%f\t%f\t%f\t%f\t%f\t%f", mu[A],
							   sd[B],
							   Xsg[A*sd_size+B],
							   Xuni[A*sd_size+B],
							   q[A*sd_size+B],
							   m[A*sd_size+B],
							   c[A*sd_size+B]);

        	fclose(file);
		free(fname); fname = NULL;

	}/////////////////////////////////////////////////////////////////////////////////
	///////////////////////// END OF STANDARD DEVIATION LOOP (B) /////////////////////
	//////////////////////////////////////////////////////////////////////////////////

	}/////////////////////////////////////////////////////////////////////////////////
	////////////////////////////// END OF  AVERAGE LOOP (A) //////////////////////////
	//////////////////////////////////////////////////////////////////////////////////

// !!!! FREE POINTERS !!!!
	free(mu); mu=NULL;
	free(sd); sd=NULL;
	free(s); s=NULL;
	free(SS); SS=NULL;
	free(J); J=NULL;
	free(m); m=NULL;
	free(q); q=NULL;
	free(SumC); SumC=NULL;
	free(ProdC); ProdC=NULL;
	free(Xuni); Xuni=NULL;
	free(Xsg); Xsg=NULL;
	free(c); c=NULL;
	gsl_rng_free (r);

//-------- END OF THE PROGRAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	elapsed_time();
	return 0;

}

// 	FUNCTIONS	//

/*********************************************************************************************
		                      Metropolis algorithm
**********************************************************************************************/

void sweep()
{
	int k;
	double Heff, delta;

	for(int i=0; i<N; i++)
	{
		/* Choose a site at random */
		k = N*randnum();

		/* Calculate the effective field */
		Heff=0;
		for(int j=0; j<N; j++) Heff += JArr(J,k,j)*s[j];

		/* Calculate the change in energy */
		delta = s[k]*(Heff);

		/* Decide whether to flip the spin */
		if (delta <= 0){
			s[k] = -s[k];
		}
		else if ( randnum() < exp(-2*delta) ){
			s[k] = -s[k];
		}
	}
}

/*********************************************************************************************
                      Covariance: function to calculate correlations
**********************************************************************************************/

double cov(int *A, int k, int l)
{
	long long int P=0, Qk=0, Ql=0;
	double R=0;

	for(int t=0; t<tdim; t++){
		P += SSArr(A,k,t) * SSArr(A,l,t);
		Qk += SSArr(A,k,t);
		Ql += SSArr(A,l,t);
	}

	R = P/(1.0*tdim) - Qk * Ql/(1.0*tdim*tdim);

	return R;
}

/*********************************************************************************************
                 Mersenne Twister random number generator: real number in [0,1)
**********************************************************************************************/

/*Seeding the MT random number generator*/
void init_randnum()
{
	//init_genrand(time(NULL));
	init_genrand(0);
}
/*Generates a real number "rand()" between [0,1)*/
double randnum()
{
	double num;
	num=genrand_real2();
	return num;
}

/*********************************************************************************************
                                       Execution time
**********************************************************************************************/

void elapsed_time(void)
{
	printf("\n.\n.\n.\nElapsed time: %f min\n\n", clock()/(60.0*CLOCKS_PER_SEC));
}


text:
/* ode-initval/gear1.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.
 */

/* Gear 1 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

typedef struct
{
  double *k;
  double *y0;
}
gear1_state_t;

static void *
gear1_alloc (size_t dim)
{
  gear1_state_t *state = (gear1_state_t *) malloc (sizeof (gear1_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for gear1_state", GSL_ENOMEM);
    }

  state->k = (double *) malloc (dim * sizeof (double));

  if (state->k == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k", GSL_ENOMEM);
    }

  state->y0 = (double *) malloc (dim * sizeof (double));

  if (state->y0 == 0)
    {
      free (state->k);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
    }

  return state;
}

static int
gear1_apply(void * vstate,
            size_t dim,
            double t,
            double h,
            double y[],
            double yerr[],
            const double dydt_in[],
            double dydt_out[],
            const gsl_odeiv_system * sys)
{
  gear1_state_t *state = (gear1_state_t *) vstate;

  const int iter_steps = 3;
  int status = 0;
  int nu;
  size_t i;

  double * const k = state->k;
  double * const y0 = state->y0;

  DISCARD_POINTER(dydt_in); /* prevent warning about unused parameter */

  DBL_MEMCPY(y0, y, dim);

  /* iterative solution */
  for(nu=0; nu<iter_steps; nu++) {
    int s = GSL_ODEIV_FN_EVAL(sys, t + h, y, k);
    GSL_STATUS_UPDATE(&status, s);
    for(i=0; i<dim; i++) {
      y[i] = y0[i] + h * k[i];
    }
  }

  /* fudge the error estimate */
  for(i=0; i<dim; i++) {
    yerr[i] = h * h * k[i];
  }

  if(dydt_out != NULL) {
    DBL_MEMCPY(dydt_out, k, dim);
  }

  return status;
}

static int
gear1_reset (void *vstate, size_t dim)
{
  gear1_state_t *state = (gear1_state_t *) vstate;

  DBL_ZERO_MEMSET (state->k, dim);
  DBL_ZERO_MEMSET (state->y0, dim);

  return GSL_SUCCESS;
}

static unsigned int
gear1_order (void *vstate)
{
  gear1_state_t *state = (gear1_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 2;
}

static void
gear1_free (void *vstate)
{
  gear1_state_t *state = (gear1_state_t *) vstate;
  free (state->k);
  free (state->y0);
  free (state);
}

static const gsl_odeiv_step_type gear1_type = { "gear1",        /* name */
  1,                            /* can use dydt_in */
  0,                            /* gives exact dydt_out */
  &gear1_alloc,
  &gear1_apply,
  &gear1_reset,
  &gear1_order,
  &gear1_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_gear1 = &gear1_type;


text:
#ifndef PRK_PETSC_H_
#define PRK_PETSC_H_

#define PRK_PETSC_USE_MPI 1

#ifdef PRK_PETSC_USE_MPI
#include <mpi.h>
#endif

#include <petsc.h>
#include <petscsys.h>
#include <petscvec.h>
#include <petscmat.h>

#endif // PRK_PETSC_H_

max_stars_repo_name: hattom/Kernels

max_issues_repo_name: hattom/Kernels

max_forks_repo_name: hattom/Kernels


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

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_linalg.h>

#define L 1 //1 2 3 4
#define LL 2 //2^(1 2 3 4)
#define LLL 4 //2^(2 3 4 5)

#define N 8

#define M 1 //# of trials
#define R 1 //the number of different initial conditions
#define S 100 //the number of different rules
#define T1 100 //transient
#define T2 1000 //time window used
int **mat;
int component[N],cnum;

void findcomponent(int node)
{
	int i;

	for (i=0;i<N;i++){
		if (mat[node][i]==1 && component[i]==0){
			component[i]=cnum;
			findcomponent(i);
		}
	}
}

double dabs(double x)
{
	if (x<0.0)return -x;
	else return x;
}

int main()
{
	FILE *fp,*fq,*fe,*fg,*fh,*fc,*fa;

	int i,j,k,l,t,m,r,s;
	int flag,flag1,flag2,tmp,tmp1,tmp2,tmp3,tmp4,ttmp;
	double dtmp;

	int x[N],xx[N][L+1],y[N];
	int in[N];
	int **nb,**rule;

	int check[N],sum,check_sum;

	double p1[N][LL],p2[N][LLL];
	double ****p3,****p4;
	int num1,num2,num3;

	double indeg[N],outdeg[N];

	double p=0.5;

	double q;
	int kav=2;
	int **adj,**adjflag;

	double *llap;
	double **lap,**edge,**grad,**circ;
	double div[N],sol[N];
	double diag[N];
	double **uu,**vv,**lap_inv;

	double el2_av,gl2_av,cl2_av,cratio_av;
	double el2,gl2,cl2,cratio;

	double *curl,*vp;
	int trinum;
	int **tri;

	double *ccurcur;
	double **curcur;
	double *diag2;
	double **uu2,**vv2,**curcur_inv;

	double **lcirc,**gcirc;

	double hl2_av,lcl2_av,gcratio_av,lcratio_av;
	double hl2,lcl2,gcratio,lcratio;

	double et,gt,ht,ct;
	double et_av,gt_av,ht_av,ct_av;

	char filename[100];

	gsl_matrix *v_gsl=gsl_matrix_alloc(N,N);
	gsl_vector *s_gsl=gsl_vector_alloc(N);
	gsl_vector *work=gsl_vector_alloc(N);
	gsl_matrix_view lap_gsl;

	gsl_matrix_view cur_gsl;

	const gsl_rng_type * TYPE;
	gsl_rng * ran;

	gsl_rng_env_setup();
	TYPE=gsl_rng_default;
	ran=gsl_rng_alloc(TYPE);
	gsl_rng_set(ran,2);

	adj=malloc(sizeof(int*)*N);
	adjflag=malloc(sizeof(int*)*N);

	mat=malloc(sizeof(int*)*N);
	lap=malloc(sizeof(double*)*N);
	edge=malloc(sizeof(double*)*N);
	grad=malloc(sizeof(double*)*N);
	circ=malloc(sizeof(double*)*N);
	lcirc=malloc(sizeof(double*)*N);
	gcirc=malloc(sizeof(double*)*N);

	uu=malloc(sizeof(double*)*N);
	vv=malloc(sizeof(double*)*N);
	lap_inv=malloc(sizeof(double*)*N);

	llap=malloc(sizeof(double)*(N*N));

	p3=malloc(sizeof(double***)*N);
	p4=malloc(sizeof(double***)*N);

	for (i=0;i<N;i++){
		adj[i]=malloc(sizeof(int)*N);
		adjflag[i]=malloc(sizeof(int)*N);

		mat[i]=malloc(sizeof(int)*N);

		lap[i]=malloc(sizeof(double)*N);
		edge[i]=malloc(sizeof(double)*N);
		grad[i]=malloc(sizeof(double)*N);
		circ[i]=malloc(sizeof(double)*N);
		lcirc[i]=malloc(sizeof(double)*N);
		gcirc[i]=malloc(sizeof(double)*N);

		uu[i]=malloc(sizeof(double)*N);
		vv[i]=malloc(sizeof(double)*N);
		lap_inv[i]=malloc(sizeof(double)*N);

		p3[i]=malloc(sizeof(double**)*N);
		p4[i]=malloc(sizeof(double**)*N);
		for (j=0;j<N;j++){
			p3[i][j]=malloc(sizeof(double*)*LL);
			p4[i][j]=malloc(sizeof(double*)*LLL);
			for (k=0;k<LL;k++){
				p3[i][j][k]=malloc(sizeof(double)*LL);
			}
			for (k=0;k<LLL;k++){
				p4[i][j][k]=malloc(sizeof(double)*LL);
			}
		}
	}

	sprintf(filename,"rtn_hodge_sw_av_n%d_k%d.txt",N,kav);
	fp=fopen(filename,"w");

	sprintf(filename,"rtn_hodge_sw_trial_n%d_k%d.txt",N,kav);
	fq=fopen(filename,"w");

	for (q=0.1;q<0.11;q=q+0.1){

		el2_av=0.0; gl2_av=0.0; cl2_av=0.0;
		cratio_av=0.0;

		hl2_av=0.0; lcl2_av=0.0;
		gcratio_av=0.0; lcratio_av=0.0;

		et_av=0.0; gt_av=0.0; ht_av=0.0; ct_av=0.0;

		for (m=0;m<M;m++){

			for (i=0;i<N;i++){
				for (j=0;j<N;j++){
					adj[i][j]=0;
					adjflag[i][j]=0;
				}
			}
			for (i=0;i<N;i++){
				for (k=-kav;k<kav+1;k++){
					if (k!=0){
						adj[(i+k+N)%N][i]=1;
						adj[i][(i+k+N)%N]=1;
					}
				}
			}

			for (i=0;i<N;i++){
				for (k=-kav;k<kav+1;k++){
					if (k!=0 && gsl_rng_uniform(ran)<q){
						adjflag[(i+k+N)%N][i]=1;
						adjflag[i][(i+k+N)%N]=1;
					}
				}
			}

			for (i=0;i<N;i++){
				for (j=i+1;j<N;j++){
					if (adjflag[i][j]==1){
						if (gsl_rng_uniform(ran)<0.5){
							while (1){
								tmp=(int)((double)N*gsl_rng_uniform(ran));

								if (j!=tmp && adj[tmp][j]==0){
									adj[i][j]=0;
									adj[j][i]=0;
									adj[tmp][j]=1;
									adj[j][tmp]=1;
									break;
								}
							}
						} else{
							while (1){
								tmp=(int)((double)N*gsl_rng_uniform(ran));

								if (i!=tmp && adj[i][tmp]==0){
									adj[i][j]=0;
									adj[j][i]=0;
									adj[tmp][i]=1;
									adj[i][tmp]=1;
									break;
								}
							}
						}
					}
				}
			}

			for (i=0;i<N;i++){
				for (j=i+1;j<N;j++){
					if (adj[i][j]==1){
						if (gsl_rng_uniform(ran)<0.5)adj[i][j]=0;
						else adj[j][i]=0;
					}
				}
			}

			for (i=0;i<N;i++){
				in[i]=0;
				for (j=0;j<N;j++){
					in[i]+=adj[j][i];
				}
			}

			nb=malloc(sizeof(int*)*N);
			rule=malloc(sizeof(int*)*N);
			for (i=0;i<N;i++){
				nb[i]=malloc(sizeof(int)*in[i]);
				rule[i]=malloc(sizeof(int)*in[i]);
			}

			//in-neighbors
			for (i=0;i<N;i++){
				k=0;
				for (j=0;j<N;j++){
					if (adj[j][i]==1){
						nb[i][k]=j;
						k++;
					}
				}
			}
			//

			for (i=0;i<N;i++){
				for (j=0;j<N;j++){
					edge[i][j]=0.0;
				}
			}

			for (s=0;s<S;s++){
				//random generation of rules
				for (i=0;i<N;i++){
					for (j=0;j<in[i];j++){
					//	rule[i][j]=gsl_ran_gaussian(ran,sigma);
						if (gsl_rng_uniform(ran)<0.5)rule[i][j]=1;
						else rule[i][j]=-1;
					}
				}

				for (k=0;k<N;k++){
					for (i=0;i<LL;i++){
						p1[k][i]=0.0;
					}

					for (i=0;i<LLL;i++){
						p2[k][i]=0.0;
					}

					for (l=0;l<in[k];l++){
						for (i=0;i<LL;i++){
							for (j=0;j<LL;j++){
								p3[k][nb[k][l]][i][j]=0.0;
							}
						}
						for (i=0;i<LLL;i++){
							for (j=0;j<LL;j++){
								p4[k][nb[k][l]][i][j]=0.0;
							}
						}
					}
				}

				for (r=0;r<R;r++){

					//initial state
					for (i=0;i<N;i++){
						if (gsl_rng_uniform(ran)<p)x[i]=1;
						else x[i]=-1;
					}

					//transient
					for (t=0;t<T1;t++){

						for (i=0;i<N;i++){
							tmp=0;
							for (j=0;j<in[i];j++){
								tmp+=rule[i][j]*x[nb[i][j]];
							}
							if (tmp<0)y[i]=-1;
							else y[i]=1;
						}

						for (i=0;i<N;i++){
							x[i]=y[i];
						}
					}

					for (i=0;i<N;i++){
						for (j=0;j<L;j++){
							xx[i][j]=0;
						}
						xx[i][L]=x[i];
					}

					//time window used to calculation
					for (t=0;t<T2+L-1;t++){
						for (i=0;i<N;i++){
							tmp=0;
							for (j=0;j<in[i];j++){
								tmp+=rule[i][j]*x[nb[i][j]];
							}
							if (tmp<0)y[i]=-1;
							else y[i]=1;

							for (j=0;j<L;j++){
								xx[i][j]=xx[i][j+1];
							}
							xx[i][L]=y[i];
						}

						if (t>=L-1){
							for (i=0;i<N;i++){
								num1=0; num2=0;
								for (k=0;k<L;k++){
									num1+=((xx[i][k]+1)/2)*(int)pow(2.0,k);
									num2+=((xx[i][k]+1)/2)*(int)pow(2.0,k);
								}
								num2+=((xx[i][L]+1)/2)*(int)pow(2.0,L);

								p1[i][num1]+=1.0;
								p2[i][num2]+=1.0;

								for (j=0;j<in[i];j++){
									num3=0;
									for (k=0;k<L;k++){
										num3+=((xx[nb[i][j]][k]+1)/2)*(int)pow(2.0,k);
									}

									p3[i][nb[i][j]][num1][num3]+=1.0;
									p4[i][nb[i][j]][num2][num3]+=1.0;
								}
							}
						}

						for (i=0;i<N;i++){
							x[i]=y[i];
						}
					}
				}//end of r-loop

				for (k=0;k<N;k++){
					for (i=0;i<LL;i++){
						p1[k][i]=p1[k][i]/(double)(T2*R);
					}

					for (i=0;i<LLL;i++){
						p2[k][i]=p2[k][i]/(double)(T2*R);
					}

					for (l=0;l<in[k];l++){
						for (i=0;i<LL;i++){
							for (j=0;j<LL;j++){
								p3[k][nb[k][l]][i][j]=p3[k][nb[k][l]][i][j]/(double)(T2*R);
							}
						}
						for (i=0;i<LLL;i++){
							for (j=0;j<LL;j++){
								p4[k][nb[k][l]][i][j]=p4[k][nb[k][l]][i][j]/(double)(T2*R);
							}
						}
					}
				}

				//construction of edge flow
				for (i=0;i<N;i++){
					for (j=0;j<in[i];j++){
						for (num2=0;num2<LLL;num2++){
							for (num3=0;num3<LL;num3++){
								if (p4[i][nb[i][j]][num2][num3]>0.0)edge[nb[i][j]][i]+=p4[i][nb[i][j]][num2][num3]*(log(p4[i][nb[i][j]][num2][num3])/log(2.0) + log(p1[i][num2%LL])/log(2.0) - log(p2[i][num2])/log(2.0) - log(p3[i][nb[i][j]][num2%LL][num3])/log(2.0));
							}
						}
					}
				}
			}//end of s-loop

			for (i=0;i<N;i++){
				for (j=0;j<N;j++){
					edge[i][j]=edge[i][j]/(double)S;
				}
			}

			for (i=0;i<N;i++){
				for (j=i+1;j<N;j++){
					dtmp=edge[i][j]-edge[j][i];
					edge[i][j]=dtmp;
					edge[j][i]=-dtmp;
				}
			}

			//finding the number of components
			for (i=0;i<N;i++){
				for (j=0;j<N;j++){
					mat[i][j]=0;
				}
			}

			for (i=0;i<N;i++){
				for (j=0;j<in[i];j++){
					if (i!=nb[i][j]){
						mat[nb[i][j]][i]=1;
						mat[i][nb[i][j]]=1;
					}
				}
			}

			tmp=0;
			for (i=0;i<N;i++){
				for (j=i+1;j<N;j++){
					if (mat[i][j]==1)tmp++;
				}
			}
			et=(double)tmp;

			for (i=0;i<N;i++){
				component[i]=0;
			}

			cnum=1; component[0]=cnum;
			findcomponent(0);
			for (i=1;i<N;i++){
				if (component[i]==0){
					cnum++;
					component[i]=cnum;
					findcomponent(i);
				}
			}

			gt=(double)(N-cnum);

			//graph Laplacian
			for (i=0;i<N;i++){
				for (j=0;j<N;j++){
					lap[i][j]=0.0;
				}
			}

			for (i=0;i<N;i++){
				lap[i][i]=0.0;
				for (j=0;j<N;j++){
					lap[i][i]+=(double)mat[i][j];
				}
			}

			for (i=0;i<N;i++){
				for (j=0;j<in[i];j++){
					if (i!=nb[i][j]){
						lap[nb[i][j]][i]=-1.0;
						lap[i][nb[i][j]]=-1.0;
					}
				}
			}

			//generalized inverse of lap
			for (i=0;i<N;i++){
				for (j=0;j<N;j++){
					llap[i*N+j]=lap[i][j];
				}
			}

			lap_gsl=gsl_matrix_view_array(llap,N,N);
			gsl_linalg_SV_decomp(&lap_gsl.matrix,v_gsl,s_gsl,work);

			for (i=0;i<N-cnum;i++){
				diag[i]=1.0/gsl_vector_get(s_gsl,i);
			}
			for (i=N-cnum;i<N;i++){
				diag[i]=0.0;
			}

			for (i=0;i<N;i++){
				for (j=0;j<N;j++){
					uu[i][j]=gsl_matrix_get(&lap_gsl.matrix,i,j);
					vv[i][j]=gsl_matrix_get(v_gsl,i,j);
				}
			}

			for (i=0;i<N;i++){
				for (j=0;j<N;j++){
					lap_inv[i][j]=0.0;
					for (k=0;k<N;k++){
						lap_inv[i][j]+=vv[i][k]*diag[k]*uu[j][k];
					}
				}
			}
			//

			//divergence of edge flow
			for (i=0;i<N;i++){
				div[i]=0.0;
				for (j=0;j<N;j++){
					div[i]+=edge[i][j];
				}
			}

			//solution (negative potential)
			for (i=0;i<N;i++){
				sol[i]=0.0;
				for (j=0;j<N;j++){
					sol[i]-=lap_inv[i][j]*div[j];
				}
			}

			//gradient flow
			for (i=0;i<N;i++){
				for (j=0;j<N;j++){
					grad[i][j]=sol[j]-sol[i];
				}
			}

			//circular flow
			for (i=0;i<N;i++){
				for (j=0;j<N;j++){
					circ[i][j]=edge[i][j]-grad[i][j];
				}
			}

			el2=0.0; gl2=0.0; cl2=0.0;
			for (i=0;i<N;i++){
				for (j=i+1;j<N;j++){
					if (lap[i][j]<0.0){
						el2+=edge[i][j]*edge[i][j];
						gl2+=grad[i][j]*grad[i][j];
						cl2+=circ[i][j]*circ[i][j];
					}
				}
			}

			if (el2>0.0)cratio=cl2/el2;
			else cratio=0.0;

			el2_av+=el2/(double)M;
			gl2_av+=gl2/(double)M;
			cl2_av+=cl2/(double)M;
			cratio_av+=cratio/(double)M;

			//decomposing circular flow into harmonic and curl flows
			trinum=0;
			for (i=0;i<N;i++){
				for (j=i+1;j<N;j++){
					for (k=j+1;k<N;k++){
						if (mat[i][j]*mat[j][k]*mat[k][i]!=0)trinum++;
					}
				}
			}

			if (trinum>0){
				curl=malloc(sizeof(double)*trinum);
				tri=malloc(sizeof(int*)*trinum);
				for (i=0;i<trinum;i++){
					tri[i]=malloc(sizeof(int)*3);
				}

				l=0;
				for (i=0;i<N;i++){
					for (j=i+1;j<N;j++){
						for (k=j+1;k<N;k++){
							if (mat[i][j]*mat[j][k]*mat[k][i]!=0){
								tri[l][0]=i; tri[l][1]=j; tri[l][2]=k;
								curl[l]=circ[i][j]+circ[j][k]-circ[i][k];
								l++;
							}
						}
					}
				}

				ccurcur=malloc(sizeof(double)*trinum*trinum);
				diag2=malloc(sizeof(double)*trinum);
				vp=malloc(sizeof(double)*trinum);

				curcur=malloc(sizeof(double*)*trinum);
				uu2=malloc(sizeof(double*)*trinum);
				vv2=malloc(sizeof(double*)*trinum);
				curcur_inv=malloc(sizeof(double*)*trinum);
				for (i=0;i<trinum;i++){
					curcur[i]=malloc(sizeof(double)*trinum);
					uu2[i]=malloc(sizeof(double)*trinum);
					vv2[i]=malloc(sizeof(double)*trinum);
					curcur_inv[i]=malloc(sizeof(double)*trinum);
				}

				for (i=0;i<trinum;i++){
					for (j=0;j<trinum;j++){
						curcur[i][j]=0.0;
						if (tri[i][0]==tri[j][0] && tri[i][1]==tri[j][1])curcur[i][j]+=1.0;
						if (tri[i][0]==tri[j][0] && tri[i][1]==tri[j][2])curcur[i][j]-=1.0;
						if (tri[i][0]==tri[j][1] && tri[i][1]==tri[j][2])curcur[i][j]+=1.0;
						if (tri[i][1]==tri[j][0] && tri[i][2]==tri[j][1])curcur[i][j]+=1.0;
						if (tri[i][1]==tri[j][0] && tri[i][2]==tri[j][2])curcur[i][j]-=1.0;
						if (tri[i][1]==tri[j][1] && tri[i][2]==tri[j][2])curcur[i][j]+=1.0;
						if (tri[i][0]==tri[j][0] && tri[i][2]==tri[j][1])curcur[i][j]-=1.0;
						if (tri[i][0]==tri[j][0] && tri[i][2]==tri[j][2])curcur[i][j]+=1.0;
						if (tri[i][0]==tri[j][1] && tri[i][2]==tri[j][2])curcur[i][j]-=1.0;
					}
				}

				gsl_matrix *v_gsl2=gsl_matrix_alloc(trinum,trinum);
				gsl_vector *s_gsl2=gsl_vector_alloc(trinum);
				gsl_vector *work2=gsl_vector_alloc(trinum);

				//generalized inverse of curcur
				for (i=0;i<trinum;i++){
					for (j=0;j<trinum;j++){
						ccurcur[i*trinum+j]=curcur[i][j];
					}
				}

				cur_gsl=gsl_matrix_view_array(ccurcur,trinum,trinum);
				gsl_linalg_SV_decomp(&cur_gsl.matrix,v_gsl2,s_gsl2,work2);

				tmp=0;
				for (i=0;i<trinum;i++){
					dtmp=gsl_vector_get(s_gsl2,i);
					if (dtmp>1e-8)diag2[i]=1.0/dtmp;
					else {
						diag2[i]=0.0;
						tmp++;
					}
				}

				ct=(double)(trinum-tmp);
				ht=et-gt-ct;

				for (i=0;i<trinum;i++){
					for (j=0;j<trinum;j++){
						uu2[i][j]=gsl_matrix_get(&cur_gsl.matrix,i,j);
						vv2[i][j]=gsl_matrix_get(v_gsl2,i,j);
					}
				}

				for (i=0;i<trinum;i++){
					for (j=0;j<trinum;j++){
						curcur_inv[i][j]=0.0;
						for (k=0;k<trinum;k++){
							curcur_inv[i][j]+=vv2[i][k]*diag2[k]*uu2[j][k];
						}
					}
				}

				for (i=0;i<trinum;i++){
					vp[i]=0.0;
					for (j=0;j<trinum;j++){
						vp[i]+=curcur_inv[i][j]*curl[j];
					}
				}

				for (i=0;i<N;i++){
					lcirc[i][i]=0.0;
					for (j=i+1;j<N;j++){
						lcirc[i][j]=0.0;
						if (mat[i][j]!=0){
							for (k=0;k<trinum;k++){
								if (tri[k][0]==i && tri[k][1]==j)lcirc[i][j]+=vp[k];
								if (tri[k][0]==i && tri[k][2]==j)lcirc[i][j]-=vp[k];
								if (tri[k][1]==i && tri[k][2]==j)lcirc[i][j]+=vp[k];
							}
						}
					}
				}

				for (i=0;i<N;i++){
					for (j=0;j<i;j++){
						lcirc[i][j]=-lcirc[j][i];
					}
				}

				for (i=0;i<N;i++){
					for (j=0;j<N;j++){
						gcirc[i][j]=circ[i][j]-lcirc[i][j];
					}
				}

				hl2=0.0; lcl2=0.0;
				for (i=0;i<N;i++){
					for (j=i+1;j<N;j++){
						if (lap[i][j]<0.0){
							hl2+=gcirc[i][j]*gcirc[i][j];
							lcl2+=lcirc[i][j]*lcirc[i][j];
						}
					}
				}

				fe=fopen("edgeflow.txt","w");
				fg=fopen("gradientflow.txt","w");
				fh=fopen("harmonicflow.txt","w");
				fc=fopen("curlflow.txt","w");
				fa=fopen("adjmatrix.txt","w");
				for (i=0;i<N;i++){
					for (j=0;j<N;j++){
						if (adj[i][j]==1)fprintf(fa,"1 ");
						else fprintf(fa,"0 ");

						if (lap[i][j]<0.0){
							fprintf(fe,"%lf ",edge[i][j]);
							fprintf(fg,"%lf ",grad[i][j]);
							fprintf(fh,"%lf ",gcirc[i][j]);
							fprintf(fc,"%lf ",lcirc[i][j]);
						} else{
							fprintf(fe,"%lf ",0.0);
							fprintf(fg,"%lf ",0.0);
							fprintf(fh,"%lf ",0.0);
							fprintf(fc,"%lf ",0.0);
						}
					}
					fprintf(fe,"\n");
					fprintf(fg,"\n");
					fprintf(fh,"\n");
					fprintf(fc,"\n");
					fprintf(fa,"\n");
				}
				fclose(fe);
				fclose(fg);
				fclose(fh);
				fclose(fc);
				fclose(fa);

				gsl_matrix_free(v_gsl2);
				gsl_vector_free(s_gsl2);
				gsl_vector_free(work2);

				free(ccurcur);
				free(diag2);
				free(vp);

				if (curcur){
					for (i=0;i<trinum;i++){
						free(curcur[i]);
					}
				}
				if (uu2){
					for (i=0;i<trinum;i++){
						free(uu2[i]);
					}
				}
				if (vv2){
					for (i=0;i<trinum;i++){
						free(vv2[i]);
					}
				}
				if (curcur_inv){
					for (i=0;i<trinum;i++){
						free(curcur_inv[i]);
					}
				}

			} else{
				ct=0.0;
				ht=et-gt-ct;

				for (i=0;i<N;i++){
					for (j=0;j<N;j++){
						gcirc[i][j]=circ[i][j];
					}
				}

				hl2=0.0; lcl2=0.0;
				for (i=0;i<N;i++){
					for (j=i+1;j<N;j++){
						if (lap[i][j]<0.0){
							hl2+=gcirc[i][j]*gcirc[i][j];
						}
					}
				}
			}

			if (cl2>0.0)gcratio=hl2/el2;
			else gcratio=0.0;

			hl2_av+=hl2/(double)M;
			lcl2_av+=lcl2/(double)M;
			gcratio_av+=gcratio/(double)M;

			lcratio=cratio-gcratio;
			lcratio_av+=lcratio/(double)M;

			et_av+=et/(double)M;
			gt_av+=gt/(double)M;
			ht_av+=ht/(double)M;
			ct_av+=ct/(double)M;

			if (m%1==0)printf("%lf %d %d %lf %lf %lf %lf %lf %lf\n",q,m,cnum,el2,gl2,hl2,lcl2,gcratio,lcratio);
			fprintf(fq,"%lf %d %d %.15lf %.15lf %.15lf %.15lf %.15lf %.15lf %.15lf %.15lf %.15lf %.15lf\n",q,m,cnum,el2,gl2,hl2,lcl2,gcratio,lcratio,et,gt,ht,ct);

			if (nb){
				for (i=0;i<N;i++){
					if (nb[i])free(nb[i]);
				}
				free(nb);
			}

			if (rule){
				for (i=0;i<N;i++){
					if (rule[i])free(rule[i]);
				}
				free(rule);
			}
		}//end of m-loop

		fprintf(fp,"%lf %.15lf %.15lf %.15lf %.15lf %.15lf %.15lf %.15lf %.15lf %.15lf %.15lf\n",q,el2_av,gl2_av,hl2_av,lcl2_av,gcratio_av,lcratio_av,et_av,gt_av,ht_av,ct_av);
	}//end of q-loop

	fclose(fp);
	fclose(fq);

	return 0;
}


text:
/****************************************************************
  wigner_gsl.h / wigner_gsl_twice.h

  Defines Wigner coupling and recoupling symbols as wrappers for GSL
  angular momentum functions:

    - wigner_gsl.h -- takes HalfInt angular momentum arguments (RECOMMENDED)
    - wigner_gsl_twice.h -- takes integer "twice value" angular momentum arguments

  Naming convention:
    - Function names *not* ending in '2' accept HalfInt arguments J.
    - Function names ending in '2' accept integer arguments 2*J.

  See, e.g., appendix to de Shalit and Talmi for underlying formulas.

  DOCUMENTATION: See wigner_gsl.h (rather than wigner_gsl_twice.h) for
  more complete function comments.

  Language: C++

  Mark A. Caprio
  University of Notre Dame

  + 02/16/10 (mac): Created.
  + 11/13/15 (mac): Add unitary 6J for (12)3-(13)2 recoupling
    and Racah reduction factor.
  + 02/27/16 (mac): Update includes for restructured header files.
  + 03/8/16 (mac): Enclose in namespace.
  + 06/8/16 (mac): Update #define guard directive.
  + 06/21/16 (mac): Remove Racah reduction factor. Update comments.
  + 10/18/16 (mac): Update Unitary6J comment. Rename wigner2_gsl.h to
    wigner_gsl_twice.h.
  + 04/28/18 (mac): Restore missing Hat2 and ParitySign2 to
    wigner_gsl_twice.h.
  + 04/10/20 (pjf): Replace assertions with exceptions.

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

#ifndef WIGNER2_GSL_H_
#define WIGNER2_GSL_H_

#include <stdexcept>

#include <gsl/gsl_sf_coupling.h>

#include "am.h"

namespace am {

  // angular momentum hat symbol

  inline
    double Hat2(int two_j)
  {
    return sqrt(static_cast<double>(two_j+1));
  }

  // phase sign (-)^sum

  inline
    int ParitySign2(int two_sum)
  {
    if ((two_sum%2)!=0) throw std::domain_error("two_sum not even");
    int remainder = abs(two_sum/2) % 2;
    int sign = 1 - 2*remainder;
    return sign;
  }

  // Wigner3J(ja,jb,jc,ma,mb,mc)
  //   returns Wigner 3-J symbol
  //   wrapper for gsl_sf_coupling_3j

  inline
    double Wigner3J2(
		     int two_ja, int two_jb, int two_jc,
		     int two_ma, int two_mb, int two_mc
		     )
  {
    return gsl_sf_coupling_3j(
			      two_ja, two_jb, two_jc,
			      two_ma, two_mb, two_mc
			      );
  }

  // ClebschGordan(ja,ma,jb,mb,jc,mc)
  //   returns Clebsch-Gordan coefficient
  //   wrapper for gsl_sf_coupling_3j

  inline
    double ClebschGordan2(
			  int two_ja, int two_ma,
			  int two_jb, int two_mb,
			  int two_jc, int two_mc
			  )
  {
    return Hat2(two_jc)*ParitySign2(two_ja-two_jb+two_mc)
      *gsl_sf_coupling_3j(
			  two_ja, two_jb, two_jc,
			  two_ma, two_mb, -two_mc
			  );
  }

  // Wigner6J(ja,jb,jc,jd,je,jf)
  //   returns Wigner 6-J symbol
  //   wrapper for gsl_sf_coupling_6j

  inline
    double Wigner6J2(
		     int two_ja, int two_jb, int two_jc,
		     int two_jd, int two_je, int two_jf
		     )
  {
    return gsl_sf_coupling_6j(
			      two_ja, two_jb, two_jc,
			      two_jd, two_je, two_jf
			      );
  }

  // Unitary6J(ja,jb,jc,jd,je,jf)
  //   wrapper for gsl_sf_coupling_6j
  //   returns unitary recoupling symbol for (12)3-1(23) recoupling

  inline
    double Unitary6J2(
		      int two_ja, int two_jb, int two_jc,
		      int two_jd, int two_je, int two_jf
		      )
  {
    return ParitySign2(two_ja+two_jb+two_jd+two_je)*Hat2(two_jc)*Hat2(two_jf)*gsl_sf_coupling_6j(
												 two_ja, two_jb, two_jc,
												 two_jd, two_je, two_jf
												 );
  }

  // Wigner9J(ja,jb,jc,jd,je,jf,jg,jh,ji)
  //   returns Wigner 9-J symbol
  //   wrapper for gsl_sf_coupling_9j

  inline
    double Wigner9J2(
		     int two_ja, int two_jb, int two_jc,
		     int two_jd, int two_je, int two_jf,
		     int two_jg, int two_jh, int two_ji
		     )
  {
    return gsl_sf_coupling_9j(
			      two_ja, two_jb, two_jc,
			      two_jd, two_je, two_jf,
			      two_jg, two_jh, two_ji
			      );
  }

  // Unitary9J(ja,jb,jc,jd,je,jf,jg,jh,ji)
  //   returns unitary 9-J symbol
  //   wrapper for gsl_sf_coupling_9j

  inline
    double Unitary9J2(
		      int two_ja, int two_jb, int two_jc,
		      int two_jd, int two_je, int two_jf,
		      int two_jg, int two_jh, int two_ji
		      )
  {
    return Hat2(two_jc)*Hat2(two_jf)*Hat2(two_jg)*Hat2(two_jh)
      *gsl_sf_coupling_9j(
			  two_ja, two_jb, two_jc,
			  two_jd, two_je, two_jf,
			  two_jg, two_jh, two_ji
			  );
  }

} // namespace
#endif

max_stars_repo_name: mark-caprio/am

max_issues_repo_name: mark-caprio/am

max_forks_repo_name: mark-caprio/am


text:
/* dwt_undec.c
 * Paul Demorest, 2007/10
 *
 * Undecimated wavelet transform and inverse.  These
 * do the usual wavelet filter tree, but don't decimate
 * by 2 after each step.  This makes the transform
 * translation-invariant (to translations which are
 * integer numbers of samples anyways).
 *
 * Ends up with n*log2(n) sized array, from size n input.
 *
 * Wavelet coeffs g and h need to be externally supplied.
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_wavelet.h>

/* Do one UDWT step on input data, out-of-place (original
 * data array is unaltered).
 * Vars:
 *   in  pointer to input data (n doubles)
 *   outl pointer to low-pass output data (n doubles)
 *   outh pointer to high-pass output data (n doubles)
 *   n   size of in/out
 *   h   Wavelet coeffs, low-pass (nc doubles)
 *   g   Wavelet coeffs, high-pass (nc doubles)
 *   nc  size of h/g
 *   s   scale (filters expand by 2^s)
 */
void udwt_forward_step(double *in, double *outl, double *outh, size_t n,
        const double h[], const double g[], size_t nc, int s) {

    int i,j,idx;

    /* Zero output array */
    for (i=0; i<n; i++) { outl[i]=0.0; outh[i]=0.0; }

    /* Filter */
    for (i=0; i<n; i++) {
        for (j=0; j<nc; j++) {
            idx = (i + (1<<s)*j) % n;
            outl[i] += h[j] * in[idx];
            outh[i] += g[j] * in[idx];
        }
    }
}

/* Inverse UDWT step similar to above.  Remember to stay
 * consistent about which half contains h and which
 * contains g chunks!
 */
void udwt_inverse_step(double *inl, double *inh, double *out, size_t n,
        const double h[], const double g[], size_t nc, int s) {

    int i,j,idx;

    /* Zero output */
    for (i=0; i<n; i++) { out[i]=0.0; }

    /* filter, add */
    for (i=0; i<n; i++) {
        for (j=0; j<nc; j++) {
            idx = (i - (1<<s)*j) % n;
            out[i] += h[j] * inl[idx];
            out[i] += g[j] * inh[idx];
        }
    }

    /* Inverse needs a rescaling */
    for (i=0; i<n; i++) { out[i] *= 0.5; }
}

/* Return integer log 2 of input.
 * Error if not a power of 2?
 */
static int log2i(int n) {
    int rv=0;
    while ((n>>rv)>1) { rv++; }
    return(rv);
}

/* Do a undec wavelet transform on input data n.
 * Vars:
 *   in  pointer to input data (n doubles)
 *   out pointer to output data (n*log2(n) doubles)
 *   n   size of in
 *   h   Wavelet coeffs, high-pass (nc doubles)
 *   g   Wavelet coeffs, low-pass (nc doubles)
 *   nc  size of h/g
 *
 * This can be called using a GSL wavelet definition struct like:
 * gsl_wavelet *w;
 * (etc, init wavelet...)
 * dwt_undec_transform(in, &out, n, w->h1, w->g1, w->nc);
 */
int _dwt_undec_transform(double *in, double *out, size_t n,
        const double h[], const double g[], size_t nc) {

    /* Determine number of levels of decomp */
    int nl=0;
    nl = log2i(n);

    /* maybe do some value checks on n, nl .. */

    /* Temporary low-pass outputs */
    double *tmp1 = (double *)malloc(sizeof(double)*n);
    double *tmp2 = (double *)malloc(sizeof(double)*n);

    /* Loop to do dwts */
    int i;
    double *curi, *curol, *curoh, *swp;
    curi = in;
    curol = tmp1;
    curoh = out;
    for (i=0; i<nl; i++) {
        udwt_forward_step(curi, curol, curoh, n, h, g, nc, i);
        if (i==0) { curi=tmp2; }
        swp = curi;
        curi = curol;
        curol = swp;
        curoh = &(out[(i+1)*n]);
    }

    free(tmp1);
    free(tmp2);
    return(nl);
}

/* Inverse, similar to above, but switch in/out
 * This can be called using a GSL wavelet definition struct like:
 * gsl_wavelet *w;
 * (etc, init wavelet...)
 * dwt_undec_inverse(in, out, n, w->h2, w->g2, w->nc);
 */
int _dwt_undec_inverse(double *in, double *out, size_t n,
        const double h[], const double g[], size_t nc) {

    /* Determine number of levels of decomp */
    int nl=0;
    nl = log2i(n);

    /* maybe do some value checks on n, nl .. */

    /* Temporary low-pass outputs */
    int i;
    double *tmp1 = (double *)malloc(sizeof(double)*n);
    double *tmp2 = (double *)malloc(sizeof(double)*n);
    for (i=0; i<n; i++) { tmp1[i]=0.0; }

    /* Loop to do dwts */
    double *curil, *curih, *curo, *swp;
    curil = tmp1;
    curih = &in[(nl-1)*n];
    curo = tmp2;
    for (i=nl-1; i>=0; i--) {
        udwt_inverse_step(curil, curih, curo, n, h, g, nc, i);
        swp = curil;
        curil = curo;
        curo = swp;
        curih = &in[(i-1)*n];
    }

    /* We end up with final output in curil */
    for (i=0; i<n; i++) { out[i] = curil[i]; }

    free(tmp1);
    free(tmp2);
    return(nl);
}

/* Wrappers using GSL wavelet definitions */
int dwt_undec_transform(double *in, double *out, size_t n, const gsl_wavelet *w) {
    return _dwt_undec_transform(in, out, n, w->h1, w->g1, w->nc);
}
int dwt_undec_inverse(double *in, double *out, size_t n, const gsl_wavelet *w) {
    return _dwt_undec_inverse(in, out, n, w->h2, w->g2, w->nc);
}