text:
// Sammlung verwendeter Funktionen
#include "global.h"
#include "targets.h"
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>
#include <math.h>
#include <stdio.h>
#include "utility.h"
#include <string.h>
#include <time.h>
#include <omp.h>

void VVerlet_Step(const int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t, //mach einzelnen update Schritt nach Velocity-Verlet ohne Hinderniss
					void (*derivmethod) (double *y, double *ans, double t,int N)) {
	double temp_vec[N], abl[N];
	int i,j;
	for (j = 0; j < DIM; j++ )		// heavy particle outside parallel region
	{

		v[j] = v[j] + 0.5 * a[j] * MAXSTEPSIZE;
		x[j] = x[j] + v[j] * MAXSTEPSIZE;
		temp_vec[j] = x[j];					//temp_vec zur abl- berechnung
		temp_vec[DIM + j] = v[j] * mass;
	}
	#pragma omp parallel for 			//parallelisiere über Osszillatoren
	for (i= 0 ; i < OSSZI; i++)
	{
		for (j = 0; j < DIM; j++ )
		{
			v[DIM + j + i*DIM] = v[DIM + j + i*DIM] + 0.5 * a[DIM + j + i*DIM] * MAXSTEPSIZE;
			x[DIM + j + i*DIM] = x[DIM + j + i*DIM] + v[DIM + j + i*DIM] * MAXSTEPSIZE;
			temp_vec[2*DIM + i * 2*DIM + j] = x[DIM + j + i*DIM];	//temp_vec zur abl- berechnung
			temp_vec[3*DIM + i * 2*DIM + j] = v[DIM + j + i*DIM] ;
		}
	}

	derivmethod(temp_vec, abl, *t, N); //bilde Ableitung zu t
										// gehe wieder in x-v Form für a
	for (j = 0; j < DIM; j++ )
	{
		a[j] = abl[DIM + j] / mass;
		v[j] = v[j] + 0.5 * a[j] * MAXSTEPSIZE;
	}

	#pragma omp parallel for
	for (i= 0 ; i < OSSZI; i++)
	{
		for (j = 0; j < DIM; j++ )
		{
			a[DIM +j + i*DIM] = abl[3*DIM + j + i * 2*DIM]/massq[i];
			v[DIM +j + i*DIM] = v[DIM +j + i*DIM]\
									+ 0.5 * a[DIM +j + i*DIM] * MAXSTEPSIZE;
		}
	}
	//#pragma omp parallel for
	//for (i= 0 ; i < N/2; i++)
	//{
	//	v[i] = v[i] + 0.5 * a[i] * MAXSTEPSIZE;
	//}
	Update_Lattice_Position(lattice_position, x, LATTICE_SPACING);
	*t = *t + MAXSTEPSIZE;
}

void VVerlet_Step_deriv(const int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t, //mach einzelnen update Schritt nach Velocity-Verlet ohne Hinderniss
					void (*derivmethod) (double *y, double *ans, double t,int N)) {

	int i,j;
	SETNUMTHREADS
	double sum[DIM];
	#pragma omp parallel 			// Values for Bath particles parralized///////////////////
	{
		int ID = omp_get_thread_num();
		int MAX_THREADS = omp_get_num_threads();

		if (MAX_THREADS < DIM)  // Catch case : not enough threads
		{	for (j= 0; j<DIM; j++)
			{
				sum[j] = 0.0f;				// prepare sum for heavy particle
			}
		}else if (ID < DIM)  // heavy particle dimensions go to threads
		{
			sum[ID] = 0.0f;				// prepare sum for heavy particle
		}
		#pragma omp for
		for (i= 0 ; i < N/2; i++)	// Erstes Halbupdate ////
		{
			double v_temp = v[i] = v[i] + 0.5 * a[i] * MAXSTEPSIZE;		//vnext = v + a/2 *dt
			x[i] = x[i] + v_temp * MAXSTEPSIZE;			//xnext = x + vnext *dt
		}
		#pragma omp for 		// Values for Bath particles parralized
		for (i = 0 ; i < OSSZI; i++)
		{
			double coup = coupling[i];
			double om = ommega[i];
			for(j = 0; j < DIM; j++)
			{		// Ableitung /
				double a_temp = a[DIM + j + i*DIM] = - om * om * x[DIM + j + i*DIM]\
											  + coup * x[j]; 		// dp_q/dt = -w^2 * q+ gamma * x
				v[DIM +j + i*DIM] = v[DIM +j + i*DIM]\
										+ 0.5 * a_temp * MAXSTEPSIZE; // Zweites Halbupdate ////

			}
		}
		double private_sum[DIM];
		double private_c[DIM] = {0.0,0.0};			// Kahan correction
		double private_t[DIM] = {0.0,0.0};
		double private_y[DIM] = {0.0,0.0};
		for(j = 0; j < DIM; j++)		// Values for heavy particle
		{
			private_sum[j] = 0.0; 			// JEder thread berechnet Partialsumme
		}
		#pragma omp for 			// berechne die Summe gamma_i *q_i - gamma_i^2/ommega_i^2 * x
		for (i = 0 ; i < OSSZI; i++)
		{
			double coup = coupling[i];
			double om = ommega[i];
			for (j=0; j<DIM; j++)
			{	// Kahan summation scheme
				private_y[j] = coup * x[DIM + j + i*DIM]\
						 - pow(coup,2.0)/pow(om,2.0) * x[j]  - private_c[j];
				private_t[j] = private_sum[j] + private_y[j];
				private_c[j] = (private_t[j] - private_sum[j]) - private_y[j];
				private_sum[j] = private_t[j];
			}
		}
        for(j=0; j<DIM; j++)
        {
        	#pragma omp atomic		//addiere threadsummen auf
            sum[j] =sum[j] + private_sum[j];
        }
	    #pragma omp barrier
	    if ( MAX_THREADS < DIM)  // Catch case : not enough threads
		{	for (j= 0; j<DIM; j++)
			{
				a[j] = sum[j] /mass ;   	//p-Ableitung
			}
		}else if (ID < DIM)  // heavy particle dimensions go to threads
		{
			a[ID] = sum[ID] /mass ;   	//p-Ableitung
		}
	}

	for(j = 0; j< DIM; j++)
	{
		v[j] = v[j] + 0.5 * a[j] * MAXSTEPSIZE;
	}
	*t = *t + MAXSTEPSIZE;
}

void VVerlet_Step_hardsphere_reflect(const int N,  double  x[] , double v[], double a[], double *t, //mach einzelnen update Schritt nach Velocity-Verlet ohne Hinderniss
					void (*derivmethod) (double *y, double *ans, double t,int N)) {

	int i,j;
	SETNUMTHREADS
	double sum[DIM];						// variable stepsize maximum MAXSTEPSIZE
	double target_vec[20];
	int reflec_flag = 0; 					// set to 1of target if root found
	int target_nr_hit = 0;					// number of actual target hit
	double contact_stepsize;
	// upper search bound
	int number_targets;
	double lower_t_global[2*DIM]; 			// array containing contact times for all targets in cell
	int reflec_flag_global[2*DIM];			// global array tro indicate reflection on target
	for (i=0; i < 2*DIM; i++)
	{
		lower_t_global[i] = MAXSTEPSIZE;	// all times a equal or smaller stepsize, so search for minimum later
		reflec_flag_global[i] = 0;

	}
	// determine time to collsion --------------------------------------------------------------------------------------------
	contact_stepsize = MAXSTEPSIZE;
	Update_Lattice_Position(lattice_position, x, LATTICE_SPACING);
	if 	(	(number_targets =  TARGET_METHOD(lattice_position, target_vec, 	// Find Targets
		LATTICE_SPACING, ESCAPE_CELLS)\
		)>0
		)
	{
		#pragma omp parallel for
		for (int target_i = 0; target_i < number_targets; target_i++)
		{
			double lower_t; 						// lower search bound
			double upper_t;
			int signold; int signnew;
			lower_t = 0.0; 							// lower search bound
			upper_t = MAXSTEPSIZE;					// upper search bound
			double steps = 500.0;
			int root_flag = 1;						// set to one if root found and also one at fitst loop
			int root_loop = 0;						// number of current loop
			int root_loop_max = 3;
			while ((root_flag == 1)	&&	(root_loop < root_loop_max))
			{
				root_flag = 0;
				root_loop += 1;
			// get sign lower bound before loop
				double eval_value = 0.0;
				for(i = 0; i< DIM; i++)
				{
					double v_temp = v[i] + 0.5 * a[i] * lower_t;
					double x_temp = x[i] + v_temp * lower_t;
					eval_value += pow( (x_temp - target_vec[i + DIM * target_i]) , 2.0);		// f = sum (x-R)^2
				}
				eval_value = eval_value - pow(TARGET_LENGTH , 2.0);
				if (eval_value > 0)
				{
					signold = 1;
				}else
				{
					signold = -1;
				}
				double dt = (upper_t - lower_t)/ steps;
			// check all signs in intervall till first change
				for (double t_check = lower_t; t_check < upper_t; t_check += dt)
				{
					eval_value = 0.0;
					for(i = 0; i< DIM; i++)
					{
						double v_temp = v[i] + 0.5 * a[i] * t_check;
						double x_temp = x[i] + v_temp * t_check;
						eval_value += pow( (x_temp - target_vec[i + DIM * target_i]) , 2.0);		// f = sum (x-R)^2
					}
					eval_value = eval_value - pow(TARGET_LENGTH , 2.0);
					//printf("%2.2e\n", eval_value);
					if (eval_value > 0)
					{
						signnew = 1;
					}else
					{
						signnew = -1;
					}
					if (signnew != signold)				// at first sign change adjust intervall
					{
						lower_t = t_check - dt;
						steps = 10000.0;				// refine only 4 orders of magnitude per loop
						upper_t = t_check;
						root_flag = 1;					// redo loop to refine root
						lower_t_global[target_i] = lower_t;
						break;
					}
					signold = signnew;
				}
			}
		}
	}// collision detection end--------------------------------------------------------------------------------------------------
	double minimum_t = MAXSTEPSIZE;
	for (i=0; i < number_targets; i++)				//find smallest reflection time
	{
		if(lower_t_global[i] < minimum_t)
		{
			reflec_flag = 1;
			target_nr_hit = i;
			minimum_t = lower_t_global[i];
		}
	}
	if (reflec_flag ==1 )
	{
		if (TARGET_FLAG)				// check for first contact
		{
			time_first_contact += *t + minimum_t;
			TARGET_FLAG = 0;// no need to test again after 1st contact
		}
		contact_stepsize =  minimum_t;				// pick lower boud so particle is JUST outside the boundary - Else zero search converges badly
	}else
	{
		contact_stepsize = MAXSTEPSIZE;
	}

	#pragma omp parallel 							// Values for Bath particles parralized//
	{
		int ID = omp_get_thread_num();
		int MAX_THREADS = omp_get_num_threads();

		if (MAX_THREADS < DIM) 						// Catch case : not enough threads
		{	for (j= 0; j<DIM; j++)
			{
				sum[j] = 0.0f;						// prepare sum for heavy particle
			}
		}else if (ID < DIM)  						// heavy particle dimensions go to threads
		{
			sum[ID] = 0.0f;							// prepare sum for heavy particle
		}
		#pragma omp for
		for (i= 0 ; i < N/2; i++)					// Erstes Halbupdate ////
		{
			double v_temp = v[i] = v[i] + 0.5 * a[i] * contact_stepsize;		//vnext = v + a/2 *dt
			x[i] = x[i] + v_temp * contact_stepsize;			//xnext = x + vnext *dt
		}
		#pragma omp for 		// Values for Bath particles parralized
		for (i = 0 ; i < OSSZI; i++)
		{
			double coup = coupling[i];
			double om = ommega[i];
			for(j = 0; j < DIM; j++)
			{		// Ableitung /
				double a_temp = a[DIM + j + i*DIM] = - om * om * x[DIM + j + i*DIM]\
											  + coup * x[j]; 		// dp_q/dt = -w^2 * q+ gamma * x
				v[DIM +j + i*DIM] = v[DIM +j + i*DIM]\
										+ 0.5 * a_temp * contact_stepsize; // Zweites Halbupdate ////

			}
		}
		double private_sum[DIM];
		double private_c[DIM] = {0.0,0.0};			// Kahan correction
		double private_t[DIM] = {0.0,0.0};
		double private_y[DIM] = {0.0,0.0};
		for(j = 0; j < DIM; j++)					// Values for heavy particle
		{
			private_sum[j] = 0.0; 					// ech thread calcs partial sum
		}
		#pragma omp for 							// calc sum gamma_i *q_i - gamma_i^2/ommega_i^2 * x
		for (i = 0 ; i < OSSZI; i++)
		{
			double coup = coupling[i];
			double om = ommega[i];
			for (j=0; j<DIM; j++)
			{										// Kahan summation scheme
				private_y[j] = coup * x[DIM + j + i*DIM]\
						 - pow(coup,2.0)/pow(om,2.0) * x[j]  - private_c[j];
				private_t[j] = private_sum[j] + private_y[j];
				private_c[j] = (private_t[j] - private_sum[j]) - private_y[j];
				private_sum[j] = private_t[j];
			}
		}
        for(j=0; j<DIM; j++)
        {
        	#pragma omp atomic		//addiere threadsummen auf
            sum[j] =sum[j] + private_sum[j];
        }
	    #pragma omp barrier
	    if ( MAX_THREADS < DIM)  // Catch case : not enough threads
		{	for (j= 0; j<DIM; j++)
			{
				a[j] = sum[j] /mass ;   	//p-Ableitung
			}
		}else if (ID < DIM)  // heavy particle dimensions go to threads
		{
			a[ID] = sum[ID] /mass ;   	//p-Ableitung
		}
	}

	for(j = 0; j< DIM; j++)
	{
		v[j] = v[j] + 0.5 * a[j] * contact_stepsize;
	}
	*t = *t + contact_stepsize;
	double norm = 0.0;
	for(j = 0; j< DIM; j++)
	{
			norm  += pow(target_vec[j + DIM * target_nr_hit] - x[j],2.0);
	}
	norm = sqrt(norm);
	if (number_targets)
	{
		if(norm < TARGET_LENGTH * 0.98)									// catch case particle entered sphere, only nonempty cells
		{
			FLAG_DUMP_VALUE  = 1;
			printf("\nsphere entered at t=%2.2e, lattice point %d %d", *t, lattice_position[0], lattice_position[1]);
		}
	}

	if ((reflec_flag))												// reflect on surface normal. Assume particle is on the surface
																	// prevent case of inner reflection in case particle still got in
	{
		double surface_normal[DIM];
		double vec_norm = 0.0;
		double scalar_prod = 0.0;
		// find vector towards target
		for(j = 0; j< DIM; j++)
		{
			surface_normal[j] = target_vec[j + DIM * target_nr_hit] - x[j];
		}
		for(j = 0; j< DIM; j++)
		{
			vec_norm += surface_normal[j] * surface_normal[j];
		}
		// point surface vector outwards and normalize
		vec_norm = sqrt(vec_norm);
		for(j = 0; j< DIM; j++)
		{
			surface_normal[j] = surface_normal[j]/vec_norm;
		}
		// find part of v in surface direction
		for(j = 0; j< DIM; j++)
		{
			scalar_prod += surface_normal[j] * v[j];
		}
		// reflect that part
		for(j = 0; j< DIM; j++)
		{
			v[j] = v[j] - 2.0 * scalar_prod * surface_normal[j];
		}

	}
}

void VVerlet_Step_hardsphere_reflect_all(const int N,  double  x[] , double v[], double a[], double *t, //mach einzelnen update Schritt nach Velocity-Verlet ohne Hinderniss
					void (*derivmethod) (double *y, double *ans, double t,int N)) {

	int i,j;
	SETNUMTHREADS
	double sum[DIM];						// variable stepsize maximum MAXSTEPSIZE
	double target_vec[20];
	int reflec_flag = 0; 					// set to 1of target if root found
	int target_nr_hit = 0;					// number of actual target hit
	double contact_stepsize;
	// upper search bound
	int number_targets;
	double lower_t_global[2*DIM]; 			// array containing contact times for all targets in cell
	int reflec_flag_global[2*DIM];			// global array tro indicate reflection on target
	for (i=0; i < 2*DIM; i++)
	{
		lower_t_global[i] = MAXSTEPSIZE;	// all times a equal or smaller stepsize, so search for minimum later
		reflec_flag_global[i] = 0;

	}
	// determine time to collsion --------------------------------------------------------------------------------------------
	contact_stepsize = MAXSTEPSIZE;
	Update_Lattice_Position(lattice_position, x, LATTICE_SPACING);
	if 	(	(number_targets =  TARGET_METHOD(lattice_position, target_vec, 	// Find Targets
		LATTICE_SPACING, ESCAPE_CELLS)\
		)>0
		)
	{
		#pragma omp parallel for
		for (int target_i = 0; target_i < number_targets; target_i++)
		{
			double lower_t; 						// lower search bound
			double upper_t;
			int signold; int signnew;
			lower_t = 0.0; 							// lower search bound
			upper_t = MAXSTEPSIZE;					// upper search bound
			double steps = 500.0;
			int root_flag = 1;						// set to one if root found and also one at fitst loop
			int root_loop = 0;						// number of current loop
			int root_loop_max = 3;
			while ((root_flag == 1)	&&	(root_loop < root_loop_max))
			{
				root_flag = 0;
				root_loop += 1;
			// get sign lower bound before loop
				double eval_value = 0.0;
				for(i = 0; i< DIM; i++)
				{
					double v_temp = v[i] + 0.5 * a[i] * lower_t;
					double x_temp = x[i] + v_temp * lower_t;
					eval_value += pow( (x_temp - target_vec[i + DIM * target_i]) , 2.0);		// f = sum (x-R)^2
				}
				eval_value = eval_value - pow(TARGET_LENGTH , 2.0);
				if (eval_value > 0)
				{
					signold = 1;
				}else
				{
					signold = -1;
				}
				double dt = (upper_t - lower_t)/ steps;
			// check all signs in intervall till first change
				for (double t_check = lower_t; t_check < upper_t; t_check += dt)
				{
					eval_value = 0.0;
					for(i = 0; i< DIM; i++)
					{
						double v_temp = v[i] + 0.5 * a[i] * t_check;
						double x_temp = x[i] + v_temp * t_check;
						eval_value += pow( (x_temp - target_vec[i + DIM * target_i]) , 2.0);		// f = sum (x-R)^2
					}
					eval_value = eval_value - pow(TARGET_LENGTH , 2.0);
					//printf("%2.2e\n", eval_value);
					if (eval_value > 0)
					{
						signnew = 1;
					}else
					{
						signnew = -1;
					}
					if (signnew != signold)				// at first sign change adjust intervall
					{
						lower_t = t_check - dt;
						steps = 10000.0;				// refine only 4 orders of magnitude per loop
						upper_t = t_check;
						root_flag = 1;					// redo loop to refine root
						lower_t_global[target_i] = lower_t;
						break;
					}
					signold = signnew;
				}
			}
		}
	}// collision detection end--------------------------------------------------------------------------------------------------
	double minimum_t = MAXSTEPSIZE;
	for (i=0; i < number_targets; i++)				//find smallest reflection time
	{
		if(lower_t_global[i] < minimum_t)
		{
			reflec_flag = 1;
			target_nr_hit = i;
			minimum_t = lower_t_global[i];
		}
	}
	if (reflec_flag ==1 )
	{
		if (TARGET_FLAG)				// check for first contact
		{
			time_first_contact += *t + minimum_t;
			TARGET_FLAG = 0;// no need to test again after 1st contact
		}
		contact_stepsize =  minimum_t;				// pick lower boud so particle is JUST outside the boundary - Else zero search converges badly
	}else
	{
		contact_stepsize = MAXSTEPSIZE;
	}

	#pragma omp parallel 							// Values for Bath particles parralized//
	{
		int ID = omp_get_thread_num();
		int MAX_THREADS = omp_get_num_threads();

		if (MAX_THREADS < DIM) 						// Catch case : not enough threads
		{	for (j= 0; j<DIM; j++)
			{
				sum[j] = 0.0f;						// prepare sum for heavy particle
			}
		}else if (ID < DIM)  						// heavy particle dimensions go to threads
		{
			sum[ID] = 0.0f;							// prepare sum for heavy particle
		}
		#pragma omp for
		for (i= 0 ; i < N/2; i++)					// Erstes Halbupdate ////
		{
			double v_temp = v[i] = v[i] + 0.5 * a[i] * contact_stepsize;		//vnext = v + a/2 *dt
			x[i] = x[i] + v_temp * contact_stepsize;			//xnext = x + vnext *dt
		}
		#pragma omp for 		// Values for Bath particles parralized
		for (i = 0 ; i < OSSZI; i++)
		{
			double coup = coupling[i];
			double om = ommega[i];
			for(j = 0; j < DIM; j++)
			{		// Ableitung /
				double a_temp = a[DIM + j + i*DIM] = - om * om * x[DIM + j + i*DIM]\
											  + coup * x[j]; 		// dp_q/dt = -w^2 * q+ gamma * x
				v[DIM +j + i*DIM] = v[DIM +j + i*DIM]\
										+ 0.5 * a_temp * contact_stepsize; // Zweites Halbupdate ////

			}
		}
		double private_sum[DIM];
		double private_c[DIM] = {0.0,0.0};			// Kahan correction
		double private_t[DIM] = {0.0,0.0};
		double private_y[DIM] = {0.0,0.0};
		for(j = 0; j < DIM; j++)					// Values for heavy particle
		{
			private_sum[j] = 0.0; 					// ech thread calcs partial sum
		}
		#pragma omp for 							// calc sum gamma_i *q_i - gamma_i^2/ommega_i^2 * x
		for (i = 0 ; i < OSSZI; i++)
		{
			double coup = coupling[i];
			double om = ommega[i];
			for (j=0; j<DIM; j++)
			{										// Kahan summation scheme
				private_y[j] = coup * x[DIM + j + i*DIM]\
						 - pow(coup,2.0)/pow(om,2.0) * x[j]  - private_c[j];
				private_t[j] = private_sum[j] + private_y[j];
				private_c[j] = (private_t[j] - private_sum[j]) - private_y[j];
				private_sum[j] = private_t[j];
			}
		}
        for(j=0; j<DIM; j++)
        {
        	#pragma omp atomic		//addiere threadsummen auf
            sum[j] =sum[j] + private_sum[j];
        }
	    #pragma omp barrier
	    if ( MAX_THREADS < DIM)  // Catch case : not enough threads
		{	for (j= 0; j<DIM; j++)
			{
				a[j] = sum[j] /mass ;   	//p-Ableitung
			}
		}else if (ID < DIM)  // heavy particle dimensions go to threads
		{
			a[ID] = sum[ID] /mass ;   	//p-Ableitung
		}
	}

	for(j = 0; j< DIM; j++)
	{
		v[j] = v[j] + 0.5 * a[j] * contact_stepsize;
	}
	*t = *t + contact_stepsize;
	double norm = 0.0;
	for(j = 0; j< DIM; j++)
	{
			norm  += pow(target_vec[j + DIM * target_nr_hit] - x[j],2.0);
	}
	norm = sqrt(norm);
	if (number_targets)
	{
		if(norm < TARGET_LENGTH * 0.98)									// catch case particle entered sphere, only nonempty cells
		{
			FLAG_DUMP_VALUE  = 1;
			printf("\nsphere entered at t=%2.2e, lattice point %d %d", *t, lattice_position[0], lattice_position[1]);
		}
	}

	if ((reflec_flag))												// reflect on surface normal. Assume particle is on the surface
																	// prevent case of inner reflection in case particle still got in
	{
		double surface_normal[DIM];
		double vec_norm = 0.0;
		double scalar_prod = 0.0;
		// find vector towards target
		for(j = 0; j< DIM; j++)
		{
			surface_normal[j] = target_vec[j + DIM * target_nr_hit] - x[j];
		}
		for(j = 0; j< DIM; j++)
		{
			vec_norm += surface_normal[j] * surface_normal[j];
		}
		// point surface vector outwards and normalize
		vec_norm = sqrt(vec_norm);
		for(j = 0; j< DIM; j++)
		{
			surface_normal[j] = surface_normal[j]/vec_norm;
		}
		// find part of v in surface direction
		for(j = 0; j< DIM; j++)
		{
			scalar_prod += surface_normal[j] * v[j];
		}
		// reflect that part for ALL
		#pragma omp parallel for
		for (i = 0; i < OSSZI; i++ )
		{
			for(j = 0; j< DIM; j++)
			{
				v[i*DIM + j] = v[i*DIM + j] - 2.0 * scalar_prod * surface_normal[j];
			}
		}
	}
}

void VVerlet_Step_deriv_Kupf(const int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t,
					void (*derivmethod) (double *y, double *ans, double t,int N)) {
	//mach einzelnen update Schritt nach Velocity-Verlet ohne Hinderniss
	// use fully mobile Kupferman bath that is coupling U = sum gamma_i (q_i-X)^2

	int i,j;
	SETNUMTHREADS
	double sum[DIM];
	#pragma omp parallel 			// Values for Bath particles parralized///////////////////
	{
		int ID = omp_get_thread_num();
		int MAX_THREADS = omp_get_num_threads();

		if (MAX_THREADS < DIM)  // Catch case : not enough threads
		{	for (j= 0; j<DIM; j++)
			{
				sum[j] = 0.0f;				// prepare sum for heavy particle
			}
		}else if (ID < DIM)  // heavy particle dimensions go to threads
		{
			sum[ID] = 0.0f;				// prepare sum for heavy particle
		}
		#pragma omp for
		for (i= 0 ; i < N/2; i++)	// Erstes Halbupdate ////
		{
			double v_temp = v[i] = v[i] + 0.5 * a[i] * MAXSTEPSIZE;		//vnext = v + a/2 *dt
			x[i] = x[i] + v_temp * MAXSTEPSIZE;			//xnext = x + vnext *dt
		}
		#pragma omp for 		// Values for Bath particles parralized
		for (i = 0 ; i < OSSZI; i++)
		{
			double coup = coupling[i];
			double om = ommega[i];
			for(j = 0; j < DIM; j++)
			{		// Ableitung /
				double a_temp = a[DIM + j + i*DIM] = -coup * (x[DIM + j + i*DIM] - x[j])/massq[i];
				v[DIM +j + i*DIM] = v[DIM +j + i*DIM]\
										+ 0.5 * a_temp * MAXSTEPSIZE; // Zweites Halbupdate ////

			}
		}
		double private_sum[DIM];
		double private_c[DIM] = {0.0,0.0};			// Kahan correction
		double private_t[DIM] = {0.0,0.0};
		double private_y[DIM] = {0.0,0.0};
		for(j = 0; j < DIM; j++)		// Values for heavy particle
		{
			private_sum[j] = 0.0; 			// JEder thread berechnet Partialsumme
		}
		#pragma omp for 			// berechne die Summe gamma_i *q_i - gamma_i^2/ommega_i^2 * x
		for (i = 0 ; i < OSSZI; i++)
		{
			double coup = coupling[i];
			double om = ommega[i];
			for (j=0; j<DIM; j++)
			{	// Kahan summation scheme
				private_y[j] = coup * (x[DIM + j + i*DIM] - x[j]) - private_c[j];
				private_t[j] = private_sum[j] + private_y[j];
				private_c[j] = (private_t[j] - private_sum[j]) - private_y[j];
				private_sum[j] = private_t[j];
			}
		}
        for(j=0; j<DIM; j++)
        {
        	#pragma omp atomic		//addiere threadsummen auf
            sum[j] =sum[j] + private_sum[j];
        }
	    #pragma omp barrier
	    if ( MAX_THREADS < DIM)  // Catch case : not enough threads
		{	for (j= 0; j<DIM; j++)
			{
				a[j] = sum[j] /mass ;   	//p-Ableitung
			}
		}else if (ID < DIM)  // heavy particle dimensions go to threads
		{
			a[ID] = sum[ID] /mass ;   	//p-Ableitung
		}
	}

	for(j = 0; j< DIM; j++)
	{
		v[j] = v[j] + 0.5 * a[j] * MAXSTEPSIZE;
		VIRTUAL_V[j] = v[j];
		VIRTUAL_X[j] = x[j];
	}
	*t = *t + MAXSTEPSIZE;
}

void VVerlet_Step_hardsphere_reflect_Kupf(const int N,  double  x[] , double v[], double a[], double *t, //mach einzelnen update Schritt nach Velocity-Verlet ohne Hinderniss
					void (*derivmethod) (double *y, double *ans, double t,int N)) {

	int i,j;
	SETNUMTHREADS
	double sum[DIM];						// variable stepsize maximum MAXSTEPSIZE
	double target_vec[20];
	int reflec_flag = 0; 					// set to 1of target if root found
	int target_nr_hit = 0;					// number of actual target hit
	double contact_stepsize = LASTSTEPSIZE;
	// upper search bound
	int number_targets;
	double lower_t_global[2*DIM]; 			// array containing contact times for all targets in cell
	int reflec_flag_global[2*DIM];			// global array tro indicate reflection on target
	for (i=0; i < 2*DIM; i++)
	{
		lower_t_global[i] = LASTSTEPSIZE;	// all times a equal or smaller stepsize, so search for minimum later
		reflec_flag_global[i] = 0;

	}
	// determine time to collsion --------------------------------------------------------------------------------------------
	contact_stepsize = LASTSTEPSIZE;
	if(VIRTUAL_FLAG)
	{
		Update_Lattice_Position(lattice_position, VIRTUAL_X, LATTICE_SPACING);
	}else
	{
		Update_Lattice_Position(lattice_position, x, LATTICE_SPACING);
	}

	if 	(	(number_targets =  TARGET_METHOD(lattice_position, target_vec, 	// Find Targets
		LATTICE_SPACING, ESCAPE_CELLS)\
		)>0
		)
	{
		double N_CHECK = 1000.0;
		#pragma omp parallel for
		for (int target_i = 0; target_i < number_targets; target_i++)	// check for all targets in increments of LASTSTEPSIZE/N_CHECK for contact
		{
			double dh = LASTSTEPSIZE /  N_CHECK;
			for (double h_check = 0.0; h_check < LASTSTEPSIZE; h_check += dh)
			{
				double r_check = 0.0;
				for (int j = 0; j < DIM; j++)
				{
					double v_temp = v[j] + 0.5 * a[j] * h_check;
					double x_temp = x[j] + v_temp * h_check;
					if(VIRTUAL_FLAG)
					{
						v_temp = VIRTUAL_V[j] + 0.5 * a[j] * h_check;
						x_temp = VIRTUAL_X[j] + v_temp * h_check;
					}
					r_check += pow( (x_temp - target_vec[j + target_i * DIM] ) ,2.0 );
				}
				r_check = sqrt(r_check);
				if (r_check < TARGET_LENGTH)
				{
					reflec_flag = 1;
					target_nr_hit = target_i;
					contact_stepsize = h_check - dh;
					if (TARGET_FLAG)				// check for first contact
					{
						time_first_contact += *t + contact_stepsize;
						TARGET_FLAG = 0;// no need to test again after 1st contact
					}
					break;
					break;
				}
			}
		}
	}
	//-------------------------------------------------------------------------------------------------------------------
	#pragma omp parallel 							// Values for Bath particles parralized//
	{
		int ID = omp_get_thread_num();
		int MAX_THREADS = omp_get_num_threads();

		if (MAX_THREADS < DIM) 						// Catch case : not enough threads
		{	for (j= 0; j<DIM; j++)
			{
				sum[j] = 0.0f;						// prepare sum for heavy particle
			}
		}else if (ID < DIM)  						// heavy particle dimensions go to threads
		{
			sum[ID] = 0.0f;							// prepare sum for heavy particle
		}
		#pragma omp for
		for (i= 0 ; i < N/2; i++)					// Erstes Halbupdate ////
		{
			double v_temp = v[i] = v[i] + 0.5 * a[i] * contact_stepsize;		//vnext = v + a/2 *dt
			x[i] = x[i] + v_temp * contact_stepsize;			//xnext = x + vnext *dt
		}
		if(VIRTUAL_FLAG)
		{
			#pragma omp for
			for (i= 0 ; i < DIM; i++)					// Erstes Halbupdate ////
			{
				double v_temp = VIRTUAL_V[i] = VIRTUAL_V[i] + 0.5 * a[i] * contact_stepsize;		//vnext = v + a/2 *dt
				VIRTUAL_X[i] = VIRTUAL_X[i] + v_temp * contact_stepsize;			//xnext = x + vnext *dt
			}
		}
		#pragma omp for 		// Values for Bath particles parralized
		for (i = 0 ; i < OSSZI; i++)
		{
			double coup = coupling[i];
			double om = ommega[i];
			for(j = 0; j < DIM; j++)
			{		// Ableitung /
				double a_temp = a[DIM + j + i*DIM] = -coup * (x[DIM + j + i*DIM] - x[j])/massq[i]; 		// dp_q/dt = -w^2 * q+ gamma * x
				v[DIM +j + i*DIM] = v[DIM +j + i*DIM]\
										+ 0.5 * a_temp * contact_stepsize; // Zweites Halbupdate ////

			}
		}
		double private_sum[DIM];
		double private_c[DIM] = {0.0,0.0};			// Kahan correction
		double private_t[DIM] = {0.0,0.0};
		double private_y[DIM] = {0.0,0.0};
		for(j = 0; j < DIM; j++)					// Values for heavy particle
		{
			private_sum[j] = 0.0; 					// ech thread calcs partial sum
		}
		#pragma omp for 							// calc sum gamma_i *q_i - gamma_i^2/ommega_i^2 * x
		for (i = 0 ; i < OSSZI; i++)
		{
			double coup = coupling[i];
			double om = ommega[i];
			for (j=0; j<DIM; j++)
			{										// Kahan summation scheme
				private_y[j] = coup * (x[DIM + j + i*DIM] - x[j]) - private_c[j];
				private_t[j] = private_sum[j] + private_y[j];
				private_c[j] = (private_t[j] - private_sum[j]) - private_y[j];
				private_sum[j] = private_t[j];
			}
		}
        for(j=0; j<DIM; j++)
        {
        	#pragma omp atomic		//addiere threadsummen auf
            sum[j] =sum[j] + private_sum[j];
        }
	    #pragma omp barrier
	    if ( MAX_THREADS < DIM)  // Catch case : not enough threads
		{	for (j= 0; j<DIM; j++)
			{
				a[j] = sum[j] /mass ;   	//p-Ableitung
			}
		}else if (ID < DIM)  // heavy particle dimensions go to threads
		{
			a[ID] = sum[ID] /mass ;   	//p-Ableitung
		}
	}

	for(j = 0; j< DIM; j++)
	{
		v[j] = v[j] + 0.5 * a[j] * contact_stepsize;
	}
	if(VIRTUAL_FLAG)
	{
		for(j = 0; j< DIM; j++)
		{
			VIRTUAL_V[j] = VIRTUAL_V[j] + 0.5 * a[j] * contact_stepsize;
		}
	}
	*t = *t + contact_stepsize;
	double norm = 0.0;
	for(j = 0; j< DIM; j++)
	{
		if(VIRTUAL_FLAG)
		{
			norm  += pow(target_vec[j + DIM * target_nr_hit] - VIRTUAL_X[j],2.0);
		}else
		{
			norm  += pow(target_vec[j + DIM * target_nr_hit] - x[j],2.0);
		}
	}
	norm = sqrt(norm);
	if (number_targets)
	{
		if(norm < TARGET_LENGTH * 0.98)									// catch case particle entered sphere, only nonempty cells
		{
			FLAG_DUMP_VALUE  = 1;
			printf("\nsphere entered at t=%2.2e, lattice point %d %d", *t, lattice_position[0], lattice_position[1]);
		}
	}

	if ((reflec_flag))												// reflect on surface normal. Assume particle is on the surface
																	// prevent case of inner reflection in case particle still got in
	{
		double surface_normal[DIM];
		double vec_norm = 0.0;
		double scalar_prod = 0.0;
		// find vector towards target
		for(j = 0; j< DIM; j++)
		{

			if(VIRTUAL_FLAG)
			{
				surface_normal[j] = target_vec[j + DIM * target_nr_hit] - VIRTUAL_X[j];
			}else
			{
				surface_normal[j] = target_vec[j + DIM * target_nr_hit] - x[j];
			}
		}
		for(j = 0; j< DIM; j++)
		{
			vec_norm += surface_normal[j] * surface_normal[j];
		}
		// point surface vector outwards and normalize
		vec_norm = sqrt(vec_norm);
		for(j = 0; j< DIM; j++)
		{
			surface_normal[j] = surface_normal[j]/vec_norm;
		}
		// find part of v in surface direction
		for(j = 0; j< DIM; j++)
		{
			if(VIRTUAL_FLAG)
			{
				scalar_prod += surface_normal[j] * VIRTUAL_V[j];
			}else
			{
				scalar_prod += surface_normal[j] * v[j];
			}
		}
		// reflect that part
		for(j = 0; j< DIM; j++)
		{
			if(VIRTUAL_FLAG)
			{
				VIRTUAL_V[j] = VIRTUAL_V[j] - 2.0 * scalar_prod * surface_normal[j];
			}else
			{
				v[j] = v[j] - 2.0 * scalar_prod * surface_normal[j];
			}
		}
	}
	if ((reflec_flag))
	{
		LASTSTEPSIZE = LASTSTEPSIZE - contact_stepsize; 		// do remaining step time
	}
	else
	{
		LASTSTEPSIZE = MAXSTEPSIZE;
	}
}

void VVerlet_Step_deriv_Box(const int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t,
					void (*derivmethod) (double *y, double *ans, double t,int N)) {
	//macht einzelnen update schritt in Box der Länge Lattice_Spacting
	// use fully mobile Kupferman bath that is coupling U = sum gamma_i (q_i-X)^2

	int i,j;
	SETNUMTHREADS
	double sum = 0.0;
	double contact_stepsize = LASTSTEPSIZE;
	int reflec_flag = 0;
	// find contact time ------------------------------------
	double aa,b,c;
	aa=b=c=0;

	int REFLEC_WALL = 1;
	if (REFLEC_WALL)
	{											// debugging purposes, free evolution if false
		for (int i_bound = -1; i_bound < 2; i_bound += 2)
		{
			double root1 = 0.0; double  root2 = 0.0;
			if(VIRTUAL_FLAG)
			{
				c = VIRTUAL_X[0] + i_bound * LATTICE_SPACING/2.0;
				b = VIRTUAL_V[0];
				aa = 0.5 * a[0];
			}else
			{
				c = x[0] + i_bound * LATTICE_SPACING/2.0;
				b = v[0];
				aa = 0.5 * a[0];
			}

			if(b >= 0)													//  minimize rounoff err by case split
			{
				root1 = (-b - sqrt(b*b - 4.0 * aa * c)) / (2.0 * aa);
				root2 = 2.0 * c /(-b - sqrt(b*b - 4.0 * aa * c) );
			}else
			{
				root1 =  2.0 * c /(-b + sqrt(b*b - 4.0 * aa * c) );
				root2 = (-b + sqrt(b*b - 4.0 * aa * c)) / (2.0 * aa);
			}
			if( (root1 > 0) && (root1 < contact_stepsize) )
			{
				contact_stepsize = 0.995 * root1;
				reflec_flag = 1;
			}
			if( (root2 > 0) && (root2 < contact_stepsize) )
			{
				contact_stepsize = 0.995 * root2;
				reflec_flag = 1;
			}
		}
	}
	// contact time end ||||||||||||||||||||||||||||||||||||<

	#pragma omp parallel for
	for (i= 0 ; i < N/2; i++)	// Erstes Halbupdate ////
	{
		double v_temp = v[i] = v[i] + 0.5 * a[i] * contact_stepsize;		//vnext = v + a/2 *dt
		x[i] = x[i] + v_temp * contact_stepsize;							//xnext = x + vnext *dt

	}
	if(VIRTUAL_FLAG)
	{
		double v_temp = VIRTUAL_V[0] = VIRTUAL_V[0] + 0.5 * a[0] * contact_stepsize;	// virtual particle
		VIRTUAL_X[0] = VIRTUAL_X[0] + v_temp * contact_stepsize;
	}
	for (i = 0 ; i < OSSZI; i++)
	{
		double coup = coupling[i];
		double om = ommega[i];
		double a_temp = a[1 + i] = -coup * (x[1 + i] - x[0])/massq[i];
		v[1 + i] = v[1 + i]\
					+ 0.5 * a_temp * contact_stepsize; // Zweites Halbupdate ////
	}
	#pragma omp parallel for reduction(+:sum) 			// berechne die Summe gamma_i *q_i - gamma_i^2/ommega_i^2 * x
	for (i = 0 ; i < OSSZI; i++)
	{
		double coup = coupling[i];
		sum += coup * (x[1 + i] - x[0]);
	}
	a[0] = sum /mass;
	v[0] = v[0] + 0.5 * a[0] * contact_stepsize;
	if(VIRTUAL_FLAG)
	{
		VIRTUAL_V[0] = VIRTUAL_V[0] + 0.5 * a[0] * contact_stepsize;
	}
	*t = *t + contact_stepsize;
	if(reflec_flag)
	{
		LASTSTEPSIZE = LASTSTEPSIZE - contact_stepsize;
		if(REFLECTALL_FLAG)								// reflect whole bath if flag is set
		{
			for (i = 0 ; i < OSSZI; i++)
			{
				double v_rel = v[1 + i] - v[0];
				v[1 + i] = - v[0] + v_rel;						//keep relative velocity after reflection of heavy particle the same
			}
		}
		if(VIRTUAL_FLAG)
		{
			VIRTUAL_V[0] = -VIRTUAL_V[0];										// reflect heavy
		}else
		{
			v[0] = -v[0];										// reflect heavy
		}
	}
	else
	{
		LASTSTEPSIZE = MAXSTEPSIZE;
	}
}

void VVerlet_Step_deriv_Box_Ramp(const int N,  double  x[N/2] , double v[N/2], double a[N/2], double *t,
					void (*derivmethod) (double *y, double *ans, double t,int N)) {
	//macht einzelnen update schritt in Box der Länge Lattice_Spacting
	// use fully mobile Kupferman bath that is coupling U = sum gamma_i (q_i-X)^2

	int i,j;
	SETNUMTHREADS
	double sum = 0.0;
	double contact_stepsize = LASTSTEPSIZE;
	int reflec_flag = 0;

	#pragma omp parallel for
	for (i= 0 ; i < N/2; i++)	// Erstes Halbupdate ////
	{
		double v_temp = v[i] = v[i] + 0.5 * a[i] * contact_stepsize;		//vnext = v + a/2 *dt
		x[i] = x[i] + v_temp * contact_stepsize;							//xnext = x + vnext *dt
	}
	for (i = 0 ; i < OSSZI; i++)
	{
		double coup = coupling[i];
		double om = ommega[i];
		double a_temp = a[1 + i] = -coup * (x[1 + i] - x[0])/massq[i];
		v[1 + i] = v[1 + i]\
					+ 0.5 * a_temp * contact_stepsize; // Zweites Halbupdate ////
	}
	#pragma omp parallel for reduction(+:sum) 			// berechne die Summe gamma_i *q_i - gamma_i^2/ommega_i^2 * x
	for (i = 0 ; i < OSSZI; i++)
	{
		double coup = coupling[i];
		sum += coup * (x[1 + i] - x[0]);
	}
	a[0] = sum /mass;

	if(x[0] < -(LATTICE_SPACING/2.0 - L_BALL) )			//left wall;
	{
		double height_ramp = 2*KBOLTZ*TEMP;
		a[0] += height_ramp/L_BALL/mass;
	}
	if(x[0] > (LATTICE_SPACING/2.0 - L_BALL) )			//right wall;
	{
		double height_ramp = 2*KBOLTZ*TEMP;
		a[0] += -height_ramp/L_BALL/mass;
	}
	v[0] = v[0] + 0.5 * a[0] * contact_stepsize;

	*t = *t + contact_stepsize;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <nlopt.h>
#include "mpi.h"

int numTasks, mpiRank;

#include "myFunc.h"

void writeParam(char *filename, int nParam, double *x)
{
  FILE *fn=fopen(filename,"w");
  if(!fn) {
    fprintf(stderr,"Cannot open %s\n",filename);
    exit(1);
  }

  uint32_t n=nParam; // ensure size is uint32_t
  fwrite(&n,sizeof(uint32_t), 1, fn);
  for(int i=0; i < nParam; i++) {
    float tmp=x[i];
    fwrite(&tmp,sizeof(float), 1, fn);
  }
  fclose(fn);
}

int masterOP=1;
void startClient(void * restrict my_func_data)
{
  int op;
  double xFromMPI[N_PARAM];
  double partialError,sum;

  for(;;) { // loop until the master says I am done - then exit
    MPI_Bcast(&op, 1, MPI_INT, 0, MPI_COMM_WORLD); // receive the op code
    if(op==0) { // we are done, normal exit
      break;
    }
    MPI_Bcast(xFromMPI, N_PARAM, MPI_DOUBLE, 0, MPI_COMM_WORLD); // receive the parameters
    partialError = objFunc(N_PARAM,  xFromMPI, NULL, my_func_data);
    MPI_Reduce(&partialError, &sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  }
}

double mpiObjFunc(unsigned n,  const double * restrict x, double * restrict grad,
		  void * restrict my_func_data)
{
  int op;
  double partialError, totalError=0.;

  MPI_Bcast(&masterOP, 1, MPI_INT, 0, MPI_COMM_WORLD); // Send the master op code
  MPI_Bcast((void*) x, N_PARAM, MPI_DOUBLE, 0, MPI_COMM_WORLD); // Send the parameters
  partialError = objFunc(N_PARAM,  x, NULL, my_func_data);
  MPI_Reduce(&partialError, &totalError, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); // get the totalError

  return(totalError);
}

int main(int argc, char* argv[])
{
  nlopt_opt opt;
  userData_t uData = {0};
  FILE *fout_mpi;

  if(argc < 4) {
    fprintf(stderr,"Use: datafile paramFile mpiTimingFile\n");
    return -1;
  }

  int ret = MPI_Init(&argc,&argv);
  if (ret != MPI_SUCCESS) {
    printf ("Error in MPI_Init()!\n");
    MPI_Abort(MPI_COMM_WORLD, ret);
  }

  MPI_Comm_size(MPI_COMM_WORLD,&numTasks);
  MPI_Comm_rank(MPI_COMM_WORLD,&mpiRank);

#if defined(NO_IO_EXAMPLES)
#pragma message "TEST case: Clients use random memory values"
    init_noIO(NO_IO_EXAMPLES,&uData);
#else
#pragma message "TEST case: all clients perform file reads"
  { // for simplicity, append the mpiRank to the data filename
    char filename[256];
    sprintf(filename,"%s.%d",argv[1],mpiRank);
    //fprintf(stderr,"Loading %s into coprocessor %d\n", filename, MIC_DEV);
    init(filename,&uData);
  }
#endif

  if(mpiRank > 0) {
    startClient(&uData);
  } else { // Master code

    { // open MPI results file for this size run
      char buf[256];
      sprintf(buf, "%s.%04d.txt",argv[3],numTasks);
      fout_mpi = fopen(buf,"w");
      if(!fout_mpi) {
	fprintf(stderr,"Cannot open file %s\n",buf);
	return -1;
      }
    }

    printf ("Number of tasks= %d My rank= %d, number clients %d\n", numTasks,mpiRank,numTasks);
#ifdef MPI_NUM_COPROC_PER_NODE
    printf ("Number of coprocessors per node= %d\n", MPI_NUM_COPROC_PER_NODE);
#endif
    printf("myFunc %s\n", desc);
    printf("nExamples %d\n", uData.nExamples);
    printf("Number Parameters %d\n", N_PARAM);

    opt = nlopt_create(NLOPT_LN_PRAXIS, N_PARAM); // algorithm and dimensionality
    // NOTE: alternative optimization methods ...
    //opt = nlopt_create(NLOPT_LN_NEWUOA, N_PARAM);
    //opt = nlopt_create(NLOPT_LN_COBYLA, N_PARAM);
    //opt = nlopt_create(NLOPT_LN_BOBYQA, N_PARAM);
    //opt = nlopt_create(NLOPT_LN_AUGLAG, N_PARAM);

    nlopt_set_min_objective(opt, mpiObjFunc, (void*) &uData);
/*
#if defined(MAX_RUNTIME)
  fprintf(stderr,"Warning: performing a quick %d second test!\n", MAX_RUNTIME);
  nlopt_set_maxtime(opt, MAX_RUNTIME); // Use for running quick tests
#else
  fprintf(stderr,"MAX runtime %d seconds!\n", 120*60);
  nlopt_set_maxtime(opt, (120. * 60.)); // maximum runtime in seconds
#endif
*/
    double minf; /* the minimum objective value, upon return */

    __declspec(align(64)) double x[N_PARAM];
    for(int i=0; i < N_PARAM; i++) x[i] = 0.1*(random()/(double)RAND_MAX);

    double startTime=getTime();
    ret=nlopt_optimize(opt, x, &minf);
    printf("Optimization Time %g\n",getTime()-startTime);

    if (ret < 0) {
      printf("nlopt failed! ret %d\n", ret);
    } else {
      printf("found minimum %0.10g ret %d\n", minf,ret);
    }
    writeParam(argv[2],N_PARAM, x);

    nlopt_destroy(opt);

    masterOP = 0; // signal completion
    MPI_Bcast(&masterOP, 1, MPI_INT, 0, MPI_COMM_WORLD); // Send the master op code
    printf("----------- performance times for the Master ----------\n");
    fini(&uData);
  }

  int client_nExamples[numTasks];
  double client_timeObjFunc[numTasks];
  int client_countObjFunc[numTasks];
  double client_timeDataLoad[numTasks];
  double client_minTime[numTasks];
  double client_maxTime[numTasks];

  MPI_Gather(&uData.nExamples, 1, MPI_INT, client_nExamples, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Gather(&uData.countObjFunc, 1, MPI_INT, client_countObjFunc, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Gather(&uData.timeDataLoad, 1, MPI_DOUBLE, client_timeDataLoad, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Gather(&uData.timeObjFunc, 1, MPI_DOUBLE, client_timeObjFunc, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Gather(&uData.minTime, 1, MPI_DOUBLE, client_minTime, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Gather(&uData.maxTime, 1, MPI_DOUBLE, client_maxTime, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  if(mpiRank==0) {
    printf("----------- performance times for the MPI run ----------\n");
    printf("function: %s\n",desc);

    uint64_t totalExamples=0;
    for(int i=0; i < numTasks; i++) totalExamples += client_nExamples[i];
    printf("totalExamples %g\n",(double) totalExamples);

    printf("AveObjTime %g, countObjFunc %d, totalObjTime %g\n",
	   uData.timeObjFunc/uData.countObjFunc, uData.countObjFunc, uData.timeObjFunc);
#ifdef FLOP_ESTIMATE
    printf("Estimated flops in myFunc %d, estimated average TFlop/s %g, nClients %d\n", FLOP_ESTIMATE,
	   (((double)totalExamples * (double)FLOP_ESTIMATE)/(uData.timeObjFunc/uData.countObjFunc)/1.e12),
	   numTasks);
    printf("Estimated maximum TFlop/s %g, minimum TFLop/s %g\n",
	   (((double)totalExamples*(double)FLOP_ESTIMATE)/(uData.minTime)/1.e12),
	   (((double)totalExamples*(double)FLOP_ESTIMATE)/(uData.maxTime)/1.e12) );
#endif

    fprintf(fout_mpi, "nExamples countObjFunc timeDataLoad  timeObjFunc minObjFcnTime maxObjFcnTime\n");
    for(int i=0; i < numTasks; i++) {
      fprintf(fout_mpi, "%g %g %g %g %g %g\n",
	      (double) client_nExamples[i], (double) client_countObjFunc[i], client_timeDataLoad[i],
	      client_timeObjFunc[i], client_minTime[i], client_maxTime[i]);
    }
  }

  MPI_Finalize();

  return 0;
}

ext: c

max_stars_repo_name: puckbee/PhiBench

max_issues_repo_name: puckbee/PhiBench

max_forks_repo_name: puckbee/PhiBench


text:
/*
  Copyright [2020] [IBM Corporation]
  Licensed under the Apache License, Version 2.0 (the "License");
  you may not use this file except in compliance with the License.
  You may obtain a copy of the License at
  http://www.apache.org/licenses/LICENSE-2.0
  Unless required by applicable law or agreed to in writing, software
  distributed under the License is distributed on an "AS IS" BASIS,
  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  See the License for the specific language governing permissions and
  limitations under the License.
*/

#ifndef MCAS_CCPM_LOG_H__
#define MCAS_CCPM_LOG_H__

#include <ccpm/interfaces.h>
#include <gsl/pointers>
#include <cstring>
#include <cstddef>

namespace ccpm
{
struct block_header;

struct log
  : public ILog
{
  using heap_type = gsl::not_null<IHeap_expandable *>;
  using persister_type = gsl::not_null<ccpm::persister *>;
private:
  /*
   * The log needs to be stored persistently, and to use persistent storage.
   * Use persister_type for the former and heap_type for the latter.
   */
  persister_type _persister; // not owned
  heap_type _mr; // not owned
  /* The log needs a root */
  block_header *_root; // owned
  void clear_top();

  static constexpr size_t min_log_extend = std::size_t(65536U);

  void extend(std::size_t size);
public:
  explicit log(persister_type, heap_type mr_);

  log(const log &) = delete;
  log &operator=(const log &) = delete;
  /*
   * Restoration of a log from persistent memory is done by moving the
   * log to itself with a move constructor. The persisted log has a vft
   * pointer, which is invalid when the log is rediscovered in persistent
   * memory. The move constructor creates valid vft pointer and (we hope)
   * preserves the rest of the log member data.
   */
  log(log &&) = default;
  log(log &&other, persister_type, heap_type);

  ~log();
  /*
   * Make a record of an old value
   */
  void add(void *begin, std::size_t size) override;
  /*
   * Make a record of an allocate
   */
  void allocated(void *&p, std::size_t size) override;
  /*
   * Make a record of a free
   */
  void freed(void *&p, std::size_t size) override;
  /*
   * commit and discard all log entries
   */
  void commit() override;
  /*
   * Restore all data areas added after initialization (or the most recent clear)
   * to the values present at their last add command.
   */
  void rollback() override;

  bool includes(const void *addr) const { return _mr->includes(addr); }
};

struct bad_alloc_log
	: public std::bad_alloc
{
	bad_alloc_log()
	{}
	const char *what() const noexcept override
	{
		return "bad undo log allocate";
	}
};

}
#endif

ext: h

max_stars_repo_name: omriarad/mcas

max_issues_repo_name: omriarad/mcas

max_forks_repo_name: omriarad/mcas


text:
// Authors: David Blei (blei@cs.princeton.edu)
//          Sean Gerrish (sgerrish@cs.princeton.edu)
//
// Copyright 2011 Sean Gerrish and David Blei
// All Rights Reserved.
//
// See the README for this package for details about modifying or
// distributing this software.

#ifndef LDASEQ_H
#define LDASEQ_H

#include <sys/stat.h>
#include <sys/types.h>

#include "gsl-wrappers.h"
#include "lda.h"

#define LDA_SEQ_EM_THRESH 1e-4
#define SAVE_LAG 10

/*
 * an lda sequence is a collection of simplex sequences for K topics
 * and an alpha vector
 *
 */

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <stdlib.h>
#include <assert.h>

#include "param.h"
#include "ss-lm.h"
#include "data.h"
#include "lda.h"

#define LDA_SEQ_EM_THRESHOLD 1e-5;

// lda sequence variational posterior distribution

// === allocation and initialization ===

inf_var* inf_var_alloc(int number_topics,
		       corpus_seq_t* corpus_seq);
void inf_var_free(inf_var* ptr);

// initialize lda sequence from lda model topics

void init_lda_seq_from_ss(lda_seq* model,
                          double topic_chain_variance,
                          double topic_obs_variance,
                          double alpha,
                          gsl_matrix* init_suffstats);

// === fitting ===

// infer a corpus with an lda-seq

double update_inf_var(lda_seq* seq,
		      const corpus_seq_t* data,
		      gsl_matrix** phi,
		      size_t t,
		      const char* root);
double update_inf_var_multiple(lda_seq* seq,
			       const corpus_seq_t* data,
			       gsl_matrix** phi,
			       size_t t,
			       const char* root);
void update_inf_reg(lda_seq* seq,
		    const corpus_seq_t* data,
		    gsl_matrix** phi,
		    size_t t,
		    const char* root);

double lda_seq_infer(lda_seq* model,
                     const corpus_seq_t* data,
                     gsl_matrix** suffstats,
                     gsl_matrix* gammas,
                     gsl_matrix* lhoods,
		     int iter,
		     const char* file_root);

// fit lda sequence from sufficient statistics

double fit_lda_seq(lda_seq* m,
                   const corpus_seq_t* data,
                   const corpus_seq_t* heldout,
                   const char* file_root);

void update_lda_seq_ss(int time,
                       const doc_t* doc,
                       const lda_post* post,
                       gsl_matrix** ss);

double fit_lda_seq_topics(lda_seq* model,
                          gsl_matrix** ss);

// === reading and writing ===

// read and write a lda sequence

void write_lda_seq(const lda_seq* m, const char* root);

lda_seq* read_lda_seq(const char* root, corpus_seq_t* data);

// write lda sequence sufficient statistics

void write_lda_seq_suffstats(lda_seq* m,
                             gsl_matrix** topic_ss,
                             const char* root);

// new lda sequence

lda_seq* new_lda_seq(corpus_seq_t* data,
		     int W,
		     int T,
		     int K);

void make_lda_from_seq_slice(lda* lda_m,
                             lda_seq* lda_seq_m,
                             int time);

#endif

max_stars_repo_name: iwangjian/dtm-lab

max_issues_repo_name: iwangjian/topic-extractor

max_forks_repo_name: iwangjian/topic-extractor


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

#include "Exceptions.h"

#include <gsl/gsl>
#include <iosfwd>

// ===========================================================================
//
// Test case structs
//
// ===========================================================================

using ArgSpan = ::gsl::span<const char*const>;

struct CmdLineParseTestCase
{
	template<::std::size_t N>
	CmdLineParseTestCase(const char*const(&args)[N]) noexcept
		: m_args(::gsl::make_span(args + 1, N - 1)) {}

	ArgSpan m_args;
};

struct CmdLineParseFailTestCase : public CmdLineParseTestCase
{
	template<::std::size_t N>
	CmdLineParseFailTestCase(const char*const(&args)[N], char const* pExMsgPattern) noexcept :
		CmdLineParseTestCase(args),
		m_pExMsgPattern(pExMsgPattern)
		{}
	bool doesExMatch(CmdLineError const& ex) const;

	char const*	m_pExMsgPattern;
};

::std::ostream& operator<<(::std::ostream& ostrm, CmdLineParseTestCase const& tc);

#endif // TESTUTIL_H_INCLUDED

max_stars_repo_name: IanEmmons/CmdLineUtil

max_issues_repo_name: IanEmmons/CmdLineUtil

max_forks_repo_name: IanEmmons/CmdLineUtil


text:
#include <blas.h>
#include <inttypes.h>
#include "mex.h"
#include "mx_blas.h"

/* Level 1 BLAS */
void
mx_swap(mxArray *mxx, mxArray *mxy)
{
    const ptrdiff_t incx = 1;
    const ptrdiff_t incy = 1;
    const ptrdiff_t n = (const ptrdiff_t)mxGetNumberOfElements(mxx);

    if (mxIsSingle(mxx)) {
        float *x = (float *)mxGetData(mxx);
        float *y = (float *)mxGetData(mxy);
        sswap(&n, x, &incx, y, &incy);

    } else if (mxIsDouble(mxx)) {
        double *x = (double *)mxGetData(mxx);
        double *y = (double *)mxGetData(mxy);
        dswap(&n, x, &incx, y, &incy);

    } else {
        mexErrMsgTxt("not implemented");
    }

    return;
} /* mx_swap() */

void
mx_scal(const double a, mxArray *mxx)
{
    const ptrdiff_t incx = 1;
    const ptrdiff_t n = (const ptrdiff_t)mxGetNumberOfElements(mxx);

    if (mxIsSingle(mxx)) {
        const float af = (const float)a;
        float *x = (float *)mxGetData(mxx);
        sscal(&n, &af, x, &incx);

    } else if (mxIsDouble(mxx)) {
        double *x = (double *)mxGetData(mxx);
        dscal(&n, &a, x, &incx);

    } else {
        mexErrMsgTxt("not implemented");
    }

    return;
} /* mx_scal() */

void
mx_copy(const mxArray *mxx, mxArray *mxy)
{
    const ptrdiff_t incx = 1;
    const ptrdiff_t incy = 1;
    const ptrdiff_t n = (const ptrdiff_t)mxGetNumberOfElements(mxx);

    if (mxIsSingle(mxx)) {
        float *y = (float *)mxGetData(mxy);
        const float *x = (const float *)mxGetData(mxx);
        scopy(&n, x, &incx, y, &incy);

    } else if (mxIsDouble(mxx)) {
        double *y = (double *)mxGetData(mxy);
        const double *x = (const double *)mxGetData(mxx);
        dcopy(&n, x, &incx, y, &incy);

    } else {
        mexErrMsgTxt("not implemented");
    }

    return;
} /* mx_copy() */

void
mx_axpy(const double a, const mxArray *mxx, mxArray *mxy)
{
    const ptrdiff_t incx = 1;
    const ptrdiff_t incy = 1;
    const ptrdiff_t n = (const ptrdiff_t)mxGetNumberOfElements(mxx);

    if (mxIsSingle(mxx)) {
        float *y = (float *)mxGetData(mxy);
        const float *x = (const float *)mxGetData(mxx);
        const float af = (const float)a;
        saxpy(&n, &af, x, &incx, y, &incy);

    } else if (mxIsDouble(mxx)) {
        double *y = (double *)mxGetData(mxy);
        const double *x = (const double *)mxGetData(mxx);
        daxpy(&n, &a, x, &incx, y, &incy);

    } else {
        mexErrMsgTxt("not implemented");
    }

    return;
} /* mx_axpy() */

double
mx_dot(const mxArray *mxx, const mxArray *mxy)
{
    double p = -1.0;

    const ptrdiff_t incx = 1;
    const ptrdiff_t incy = 1;
    const ptrdiff_t n = (const ptrdiff_t)mxGetNumberOfElements(mxx);

    if (mxIsSingle(mxx)) {
        const float *x = (const float *)mxGetData(mxx);
        const float *y = (const float *)mxGetData(mxy);
        p = (double)sdot(&n, x, &incx, y, &incy);

    } else if (mxIsDouble(mxx)) {
        const double *x = (const double *)mxGetData(mxx);
        const double *y = (const double *)mxGetData(mxy);
        p = ddot(&n, x, &incx, y, &incy);

    } else {
        mexErrMsgTxt("not implemented");
    }

    return p;
} /* mx_dot() */

double
mx_nrm2(const mxArray *mxx)
{
    double p = -1.0;

    const ptrdiff_t incx = 1;
    const ptrdiff_t n = (const ptrdiff_t)mxGetNumberOfElements(mxx);

    if (mxIsSingle(mxx)) {
        const float *x = (const float *)mxGetData(mxx);
        p = (double)snrm2(&n, x, &incx);

    } else if (mxIsDouble(mxx)) {
        const double *x = (const double *)mxGetData(mxx);
        p = dnrm2(&n, x, &incx);

    } else {
        mexErrMsgTxt("not implemented");
    }

    return p;
} /* mx_nrm2() */

double
mx_asum(const mxArray *mxx)
{
    double p = -1.0;

    const ptrdiff_t incx = 1;
    const ptrdiff_t n = (const ptrdiff_t)mxGetNumberOfElements(mxx);

    if (mxIsSingle(mxx)) {
        const float *x = (const float *)mxGetData(mxx);
        p = (double)sasum(&n, x, &incx);

    } else if (mxIsDouble(mxx)) {
        const double *x = (const double *)mxGetData(mxx);
        p = dasum(&n, x, &incx);

    } else {
        mexErrMsgTxt("not implemented");
    }

    return p;
} /* mx_asum() */


text:
/*
 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch

 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, or (at your option)
 any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
 02110-1301, USA.

 $Id$
*/

#include <config.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <locale.h>
#include <dirent.h>
#include <sys/time.h>
#include <gsl/gsl_rng.h>

#include "string_f.h"

unsigned long int random_seed()
{
 unsigned long int seed;
 FILE *devrandom;

 if ((devrandom = fopen("/dev/urandom","r")) == NULL) {
#ifdef HAVE_GETTIMEOFDAY
   struct timeval tv;
   gettimeofday(&tv, 0);
   seed = tv.tv_sec + tv.tv_usec;
#else
   seed = 0;
#endif
 } else {
   fread(&seed, sizeof(seed), 1, devrandom);
   fclose(devrandom);
 }

 return seed;
}

void FC_FUNC_(oct_printrecipe, OCT_PRINTRECIPE)
  (STR_F_TYPE _dir, STR_F_TYPE filename STR_ARG2)
{

#if HAVE_SCANDIR && HAVE_ALPHASORT
  char *lang, *tmp, dir[512];
  struct dirent **namelist;
  int ii, nn;
  gsl_rng *rng;

  /* get language */
  lang = getenv("LANG");
  if(lang == NULL) lang = "en";

  /* convert directory from Fortran to C string */
  TO_C_STR1(_dir, tmp);
  strcpy(dir, tmp);
  free(tmp);

  strcat(dir, "/recipes");

  /* check out if lang dir exists */
  nn = scandir(dir, &namelist, 0, alphasort);
  if (nn < 0){
    printf("Directory does not exist: %s", dir);
    return;
  }

  for(ii=0; ii<nn; ii++)
    if(strncmp(lang, namelist[ii]->d_name, 2) == 0){
      strcat(dir, "/");
      strcat(dir, namelist[ii]->d_name);
      break;
    }

  if(ii == nn)
    strcat(dir, "/en"); /* default */

  /* clean up */
  for(ii=0; ii<nn; ii++)
    free(namelist[ii]);
  free(namelist);

  /* now we read the recipes */
  nn = scandir(dir, &namelist, 0, alphasort);

  /* initialize random numbers */
  gsl_rng_env_setup();
  rng = gsl_rng_alloc(gsl_rng_default);
  gsl_rng_set(rng, random_seed());
  ii = gsl_rng_uniform_int(rng, nn - 2);
  gsl_rng_free(rng);

  strcat(dir, "/");
  strcat(dir, namelist[ii+2]->d_name); /* skip ./ and ../ */

  /* clean up again */
  for(ii=0; ii<nn; ii++)
    free(namelist[ii]);
  free(namelist);

  TO_F_STR2(dir, filename);

#else
  printf("Sorry, recipes cannot be printed unless scandir and alphasort are available with your C compiler.\n");
#endif
}

max_stars_repo_name: neelravi/octopus

max_issues_repo_name: neelravi/octopus

max_forks_repo_name: neelravi/octopus


text:
/*
Copyright (C) 2019-2020 JingWeiZhangHuai <jingweizhanghuai@163.com>
Licensed under the Apache License, Version 2.0; you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.
*/

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

#include <cblas.h>
#include "morn_tensor.h"

void ConvTensorToMatData(MTensor *tns,int bc,float *mdata,int knl_height,int knl_width,int y_stride,int x_stride);
void ConvMatDataToTensor(float *mdata,MTensor *tns,int bc,int knl_height,int knl_width,int y_stride,int x_stride);

struct TensorSampleConvPara
{
    MLayer *prev;

    int knl_num;
    int knl_height;
    int knl_width;

    int x_stride;
    int y_stride;

    float sample_ratio;

    int res_valid;

    float rate;
    float decay;
    float momentum;
};
void *mTensorSampleConvPara(MSheet *ini,char *name)
{
    struct TensorSampleConvPara *para = (struct TensorSampleConvPara *)mMalloc(sizeof(struct TensorSampleConvPara));

    char *value = mINIRead(ini,name,"prev");
    para->prev = mNetworkLayer(ini,value);
    mException((para->prev == NULL),EXIT,"invalid prev");

    para->res_valid = (strcmp("Input",mLayerType(para->prev))!=0);

    value = mINIRead(ini,name,"knl_num");
    if(value != NULL) para->knl_num= atoi(value);else para->knl_num= 1;

    value = mINIRead(ini,name,"knl_height");
    if(value != NULL) para->knl_height= atoi(value);else para->knl_height= 1;

    value = mINIRead(ini,name,"knl_width");
    if(value != NULL) para->knl_width= atoi(value);else para->knl_width= 1;

    value = mINIRead(ini,name,"x_stride");
    if(value != NULL) para->x_stride= atoi(value);else para->x_stride= 1;

    value = mINIRead(ini,name,"y_stride");
    if(value != NULL) para->y_stride= atoi(value);else para->y_stride= 1;

    value = mINIRead(ini,name,"sample_ratio");
    if(value != NULL) para->sample_ratio = atof(value);else para->sample_ratio = 0.5;

    value = mINIRead(ini,name,"rate");
    if(value != NULL) para->rate = atof(value);
    else
    {
        value = mINIRead(ini,"para","rate");
        if(value != NULL) para->rate = atof(value);
        else              para->rate = 0.001;
    }

    value = mINIRead(ini,name,"decay");
    if(value != NULL) para->decay = atof(value);
    else
    {
        value = mINIRead(ini,"para","decay");
        if(value != NULL) para->decay = atof(value);
        else              para->decay = 0.01;
    }
    mException((para->decay<0.0f)||(para->decay>=1.0f),EXIT,"invalid para decay");

    value = mINIRead(ini,name,"momentum");
    if(value != NULL) para->momentum = atof(value);
    else
    {
        value = mINIRead(ini,"para","momentum");
        if(value != NULL) para->momentum = atof(value);
        else              para->momentum = 0.9;
    }

    return para;
}

struct HandleTensorSampleConv
{
    float *mat;
    float *data;
    float *kernel;
    float *update;
    int   *locate;
};
void endTensorSampleConv(void *info)
{
    struct HandleTensorSampleConv *handle = (struct HandleTensorSampleConv *)info;
    if(handle->mat   != NULL) mFree(handle->mat   );
    if(handle->data  != NULL) mFree(handle->data  );
    if(handle->kernel!= NULL) mFree(handle->kernel);
    if(handle->update!= NULL) mFree(handle->update);
    if(handle->locate!= NULL) mFree(handle->locate);
}
#define HASH_TensorSampleConv 0x8a31af2
void TensorSampleConvSet(MLayer *layer)
{
    if(layer->state != DFLT) return;

    struct TensorSampleConvPara *para = (struct TensorSampleConvPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out= layer->tns;

    MHandle *hdl=mHandle(out,TensorSampleConv);
    struct HandleTensorSampleConv *handle = (struct HandleTensorSampleConv *)(hdl->handle);

    int out_height= in->height/para->y_stride;
    int out_width = in->width /para->x_stride;
    int mheight = (out_height*out_width);
    int mwidth = (para->knl_height*para->knl_width*in->channel)*para->sample_ratio;

    int data_size = para->knl_num*(mwidth+1);

    mTensorRedefine(out,in->batch,para->knl_num,out_height,out_width,NULL);
    if(morn_network_flag == MORN_TRAIN)
    {
        if(INVALID_TENSOR(res)) mTensorRedefine(res,in->batch,in->channel,in->height,in->width,in->data);
        else                    mTensorRedefine(res,in->batch,in->channel,in->height,in->width,NULL);

        if(handle->update != NULL) mFree(handle->update);
        handle->update =(float *)mMalloc(data_size*sizeof(float));
        memset(handle->update,0,data_size*sizeof(float));
    }

    if(handle->kernel !=NULL) mFree(handle->kernel);
    handle->kernel = (float *)mMalloc(data_size*sizeof(float));

    if(morn_network_parafile==NULL)
    {
        float scale = sqrt(2.0f/mwidth);
        for(int i=0;i<data_size;i++)
            handle->kernel[i] = mNormalRand(0.0f,1.0f)*scale;
    }
    else
    {
        mNetworkParaRead(layer,"kernel",handle->kernel,data_size*sizeof(float));
    }

    int matwidth = para->knl_height*para->knl_width*in->channel+1;
    if(handle->mat!=NULL) mFree(handle->mat);
    handle->mat = (float *)mMalloc(mheight*matwidth*sizeof(float));

    if(handle->data!=NULL) mFree(handle->data);
    handle->data= (float *)mMalloc(mheight*(mwidth+1)*sizeof(float));

    if(handle->locate!=NULL) mFree(handle->locate);
    handle->locate = (int *)mMalloc(mwidth*para->knl_num*sizeof(int));
    if(morn_network_parafile==NULL)
    {
        char *flag = (char *)mMalloc(para->knl_height*para->knl_width*in->channel*sizeof(char));
        for(int j=0;j<para->knl_num;j++)
        {
            memset(flag,1,para->knl_height*para->knl_width*in->channel*sizeof(char));
            int *locate = handle->locate+j*mwidth;

            for(int i=0;i<mwidth;)
            {
                locate[i]=mRand(0,para->knl_height*para->knl_width*in->channel);
                if(flag[i]==0) continue;
                flag[i]=0; i++;
            }

            mAscSortS32(locate,NULL,locate,NULL,mwidth);
        }
        mFree(flag);
    }
    else
    {
        mNetworkParaRead(layer,"locate",handle->locate,mwidth*para->knl_num*sizeof(int));
    }

    hdl->valid = 1;
}

void mTensorSampleConvForward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("SampleConv",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorSampleConvPara *para = (struct TensorSampleConvPara *)(layer->para);

    MTensor *in = para->prev->tns;
    MTensor *out=layer->tns;

    TensorSampleConvSet(layer);

    MHandle *hdl=mHandle(out,TensorSampleConv);
    struct HandleTensorSampleConv *handle = (struct HandleTensorSampleConv *)(hdl->handle);

    int mheight = (out->height*out->width);
    int mwidth = (para->knl_height*para->knl_width*in->channel)*para->sample_ratio;
    int matwidth = para->knl_height*para->knl_width*in->channel+1;

    float *mat_data = handle->mat;
    float *in_data = handle->data;

    for(int b=0;b<in->batch;b++)
    {
        ConvTensorToMatData(in,b,mat_data,para->knl_height,para->knl_width,para->y_stride,para->x_stride);

        for(int k=0;k<para->knl_num;k++)
        {
            int   *locate = handle->locate+k* mwidth;
            float *kernel = handle->kernel+k*(mwidth+1);
            float *out_data = out->data[b]+k* mheight;
            for(int j=0;j<mheight;j++)
            {
                for(int i=0;i<mwidth;i++) in_data[j*mwidth+i]=mat_data[j*matwidth+locate[i]];
                in_data[j*mwidth+mwidth]=1.0f;
            }

            cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasTrans,
                        1,mheight,mwidth+1,
                        1.0f,
                          kernel,mwidth+1,
                         in_data,mwidth+1,
                   0.0f,out_data,mheight);
        }
    }

    layer->state = MORN_FORWARD;
}

void mTensorSampleConvBackward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("SampleConv",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorSampleConvPara *para = (struct TensorSampleConvPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out=layer->res;

    MHandle *hdl=mHandle(layer->tns,TensorSampleConv);
    struct HandleTensorSampleConv *handle = (struct HandleTensorSampleConv *)(hdl->handle);
    mException((hdl->valid == 0),EXIT,"no forward operate");

    int mheight = (out->height*out->width);
    int mwidth = (para->knl_height*para->knl_width*in->channel)*para->sample_ratio;
    int matwidth = para->knl_height*para->knl_width*in->channel+1;

    float *mat_data = handle->mat;
    float *in_data  = handle->data;
    float *res_data = handle->data;

    mNetworkParaWrite(layer,"locate",handle->locate,mwidth*para->knl_num*sizeof(int));
    mNetworkParaWrite(layer,"kernel",handle->kernel,para->knl_num*mwidth*sizeof(float));

    for(int b=0;b<out->batch;b++)
    {
        ConvTensorToMatData(in,b,mat_data,para->knl_height,para->knl_width,para->y_stride,para->x_stride);

        for(int k=0;k<para->knl_num;k++)
        {
            int   *locate = handle->locate+k* mwidth;
            float *update = handle->update+k*(mwidth+1);
            float *out_data = out->data[b]+k* mheight;
            for(int j=0;j<mheight;j++)
            {
                for(int i=0;i<mwidth;i++) in_data[j*mwidth+i]=mat_data[j*matwidth+locate[i]];
                in_data[j*mwidth+mwidth]=1.0f;
            }

            cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,
                        1,mwidth,mheight,
                        1.0f/mheight,
                        out_data,mheight,
                        in_data,mwidth,
                        (b==0)?para->momentum:1.0f,
                        update,mwidth);

        }
    }

    cblas_saxpby(para->knl_num*mwidth,
                 (0.0f-(para->rate/(float)(in->batch))),handle->update,1,
                 (1.0f-(para->decay*para->rate))       ,handle->kernel,1);

    if(para->res_valid==0) return;

    if(para->prev->state == MORN_FORWARD)
    {
        for(int b=0;b<res->batch;b++)
            memset(res->data[b],0,in->height*in->width*in->channel*sizeof(float));
        para->prev->state = MORN_BACKWARD;
    }

    for(int b=0;b<in->batch;b++)
    {
        memset(mat_data,0,matwidth*mheight*sizeof(float));

        for(int k=0;k<para->knl_num;k++)
        {
            int   *locate = handle->locate+k* mwidth;
            float *kernel = handle->kernel+k*(mwidth+1);
            float *out_data = out->data[b]+k* mheight;

            cblas_sgemm(CblasRowMajor,CblasTrans,CblasNoTrans,
                        mheight,mwidth,1,
                        1.0f,
                        out_data,mheight,
                          kernel,mwidth,
                    0.0,res_data,mwidth);

            for(int j=0;j<mheight;j++)
            {
                for(int i=0;i<mwidth;i++) mat_data[j*matwidth+locate[i]]+=res_data[j*mwidth+i];
            }
        }

        ConvMatDataToTensor(mat_data,res,b,para->knl_height,para->knl_width,para->y_stride,para->x_stride);
    }
}

void DirConvTensorToMatData(MTensor *tns,int bc,float *mdata,int knl_r,int knl_dir,int y_stride,int x_stride)
{
    int height = tns->height;
    int width  = tns->width;
    int channel= tns->channel;

    int out_width = width/x_stride;
    int out_height=height/y_stride;
    int mwidth = (knl_r+knl_r+1)*channel+1;
    int mheight= out_height*out_width;

    int sx=0,sy=0;
         if(knl_dir==0) {sx=0;sy= 1;}
    else if(knl_dir==1) {sx=1;sy= 1;}
    else if(knl_dir==2) {sx=1;sy= 0;}
    else if(knl_dir==3) {sx=1;sy=-1;}

    int x0=(width -(out_width -1)*x_stride)/2;
    int y0=(height-(out_height-1)*y_stride)/2;

    float *tdata = tns->data[bc];
    int tsize = tns->height*tns->width;

    int i,j,c;
    for(j=out_width;j<mheight;j++)
    {
        int n=y0+j/out_width*y_stride-knl_r*sy;
        int m=x0+j%out_width*x_stride-knl_r*sx;
        for(i=0;i<mwidth-1;i+=channel)
        {
            int h=n;if(h<0)h=0;else if(h>=height)h=height-1;
            int w=m;if(w<0)w=0;else if(w>= width)w= width-1;
            for(c=0;c<channel;c++) mdata[(j*mwidth)+i+c]=tdata[c*tsize+h*width+w];
            n=n+sy;m=m+sx;
        }
        mdata[(j*mwidth)+mwidth-1]=1.0f;
    }
}

void DirConvMatDataToTensor(float *mdata,MTensor *tns,int bc,int knl_r,int knl_dir,int y_stride,int x_stride)
{
    int height = tns->height;
    int width  = tns->width;
    int channel= tns->channel;

    int out_width = width/x_stride;
    int out_height=height/y_stride;
    int mwidth = (knl_r+knl_r+1)*channel+1;
    int mheight= out_height*out_width;

    int sx=0,sy=0;
         if(knl_dir==0) {sx=0;sy= 1;}
    else if(knl_dir==1) {sx=1;sy= 1;}
    else if(knl_dir==2) {sx=1;sy= 0;}
    else if(knl_dir==3) {sx=1;sy=-1;}

    int x0=(width -(out_width -1)*x_stride)/2;
    int y0=(height-(out_height-1)*y_stride)/2;

    float *tdata = tns->data[bc];
    int tsize = tns->height*tns->width;

    int i,j,c;
    for(j=0;j<mheight;j++)
    {
        int n=y0+j/out_width*y_stride-knl_r*sy;
        int m=x0+j%out_width*x_stride-knl_r*sx;
        for(i=0;i<mwidth-1;i+=channel)
        {
            int h=n;if(h<0)h=0;else if(h>=height)h=height-1;
            int w=m;if(w<0)w=0;else if(w>= width)w= width-1;
            for(c=0;c<channel;c++) tdata[c*tsize+h*width+w]+=mdata[(j*mwidth)+i+c];
            n=n+sy;m=m+sx;
        }
    }
}

struct TensorDirConvPara
{
    MLayer *prev;

    int knl_num;
    int knl_r;

    int dir;

    int x_stride;
    int y_stride;

    int res_valid;

    float rate;
    float decay;
    float momentum;
};
void *mTensorDirConvPara(MSheet *ini,char *name)
{
    struct TensorDirConvPara *para = (struct TensorDirConvPara *)mMalloc(sizeof(struct TensorDirConvPara));

    char *value = mINIRead(ini,name,"prev");
    para->prev = mNetworkLayer(ini,value);
    mException((para->prev == NULL),EXIT,"invalid prev");

    para->res_valid = (strcmp("Input",mLayerType(para->prev))!=0);

    value = mINIRead(ini,name,"knl_num");
    if(value != NULL) para->knl_num= atoi(value);else para->knl_num= 1;

    value = mINIRead(ini,name,"knl_r");
    if(value != NULL) para->knl_r= atoi(value);else para->knl_r= 1;

    value = mINIRead(ini,name,"dir");
    mException((value==NULL),EXIT,"invalid convolution direction");
         if((strcmp(value, "u-d" )==0)||(strcmp(value, "d-u" )==0)) para->dir=0;
    else if((strcmp(value,"lu-rd")==0)||(strcmp(value,"rd-lu")==0)) para->dir=1;
    else if((strcmp(value, "l-r" )==0)||(strcmp(value, "r-l" )==0)) para->dir=2;
    else if((strcmp(value,"ru-ld")==0)||(strcmp(value,"ld-ru")==0)) para->dir=3;
    else mException(1,EXIT,"invalid convolution direction");

    value = mINIRead(ini,name,"x_stride");
    if(value != NULL) para->x_stride= atoi(value);else para->x_stride= 1;

    value = mINIRead(ini,name,"y_stride");
    if(value != NULL) para->y_stride= atoi(value);else para->y_stride= 1;

    value = mINIRead(ini,name,"rate");
    if(value != NULL) para->rate = atof(value);
    else
    {
        value = mINIRead(ini,"para","rate");
        if(value != NULL) para->rate = atof(value);
        else              para->rate = 0.001;
    }

    value = mINIRead(ini,name,"decay");
    if(value != NULL) para->decay = atof(value);
    else
    {
        value = mINIRead(ini,"para","decay");
        if(value != NULL) para->decay = atof(value);
        else              para->decay = 0.01;
    }
    mException((para->decay<0.0f)||(para->decay>=1.0f),EXIT,"invalid para decay");

    value = mINIRead(ini,name,"momentum");
    if(value != NULL) para->momentum = atof(value);
    else
    {
        value = mINIRead(ini,"para","momentum");
        if(value != NULL) para->momentum = atof(value);
        else              para->momentum = 0.9;
    }

    return para;
}

struct HandleTensorDirConv
{
    float *mat;
    float *kernel;
    float *update;
};
void endTensorDirConv(void *info)
{
    struct HandleTensorDirConv *handle = (struct HandleTensorDirConv *)info;
    if(handle->mat   != NULL) mFree(handle->mat);
    if(handle->kernel!= NULL) mFree(handle->kernel);
    if(handle->update!= NULL) mFree(handle->update);
}
#define HASH_TensorDirConv 0x9e2674b3
void TensorDirConvSet(MLayer *layer)
{
    if(layer->state != DFLT) return;

    struct TensorDirConvPara *para = (struct TensorDirConvPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out= layer->tns;

    MHandle *hdl=mHandle(out,TensorDirConv);
    struct HandleTensorDirConv *handle = (struct HandleTensorDirConv *)(hdl->handle);

    int out_height= in->height/para->y_stride;
    int out_width = in->width /para->x_stride;
    int mheight = (out_height*out_width);
    int mwidth = (para->knl_r+para->knl_r+1)*in->channel+1;
    int data_size = para->knl_num*mwidth;

    mTensorRedefine(out,in->batch,para->knl_num,out_height,out_width,NULL);
    if(morn_network_flag == MORN_TRAIN)
    {
        if(INVALID_TENSOR(res)) mTensorRedefine(res,in->batch,in->channel,in->height,in->width,in->data);
        else                    mTensorRedefine(res,in->batch,in->channel,in->height,in->width,NULL);

        if(handle->update != NULL) mFree(handle->update);
        handle->update =(float *)mMalloc(data_size*sizeof(float));
        memset(handle->update,0,data_size*sizeof(float));
    }

    if(handle->kernel !=NULL) mFree(handle->kernel);
    handle->kernel = (float *)mMalloc(data_size*sizeof(float));

    if(morn_network_parafile==NULL)
    {
        float scale = sqrt(2.0f/mwidth);
        for(int i=0;i<data_size;i++)
            handle->kernel[i] = mNormalRand(0.0f,1.0f)*scale;
    }
    else
    {
        mNetworkParaRead(layer,"kernel",handle->kernel,data_size*sizeof(float));
    }

    if(handle->mat!=NULL) mFree(handle->mat);
    handle->mat = (float *)mMalloc(mheight*mwidth*sizeof(float));

    hdl->valid = 1;
}

void mTensorDirConvForward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("DirConv",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorDirConvPara *para = (struct TensorDirConvPara *)(layer->para);

    MTensor *in = para->prev->tns;
    MTensor *out=layer->tns;

    TensorDirConvSet(layer);

    MHandle *hdl=mHandle(out,TensorDirConv);
    struct HandleTensorDirConv *handle = (struct HandleTensorDirConv *)(hdl->handle);

    int mheight = (out->height*out->width);
    int mwidth = (para->knl_r+para->knl_r+1)*in->channel+1;

    float *kernel_data= handle->kernel;
    float *in_data = handle->mat;

    for(int b=0;b<in->batch;b++)
    {
        DirConvTensorToMatData(in,b,in_data,para->knl_r,para->dir,para->y_stride,para->x_stride);
        float *out_data = out->data[b];

        cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasTrans,
                    para->knl_num,mheight,mwidth,
                    1.0f,
                    kernel_data,mwidth,
                        in_data,mwidth,
               0.0f,   out_data,mheight);
    }

    layer->state = MORN_FORWARD;
}

void mTensorDirConvBackward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("DirConv",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorDirConvPara *para = (struct TensorDirConvPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out=layer->res;

    MHandle *hdl=mHandle(layer->tns,TensorDirConv);
    struct HandleTensorDirConv *handle = (struct HandleTensorDirConv *)(hdl->handle);
    mException((hdl->valid == 0),EXIT,"no forward operate");

    int mheight = (out->height*out->width);
    int mwidth = (para->knl_r+para->knl_r+1)*in->channel+1;

    float *kernel_data= handle->kernel;
    float *update_data= handle->update;
    float *    in_data= handle->mat;
    float *   res_data= handle->mat;

    mNetworkParaWrite(layer,"kernel",kernel_data,para->knl_num*mwidth*sizeof(float));

    for(int b=0;b<out->batch;b++)
    {
        DirConvTensorToMatData(in,b,in_data,para->knl_r,para->dir,para->y_stride,para->x_stride);
        float *out_data = out->data[b];

        cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,
                    para->knl_num,mwidth,mheight,
                    1.0f/mheight,
                       out_data,mheight,
                        in_data,mwidth,
                    (b==0)?para->momentum:1.0f,
                    update_data,mwidth);
    }

    cblas_saxpby(para->knl_num*mwidth,
                 (0.0f-(para->rate/(float)(in->batch))),update_data,1,
                 (1.0f-(para->decay*para->rate))       ,kernel_data,1);

    if(para->res_valid==0) return;

    if(para->prev->state == MORN_FORWARD)
    {
        for(int b=0;b<res->batch;b++)
            memset(res->data[b],0,in->height*in->width*in->channel*sizeof(float));
        para->prev->state = MORN_BACKWARD;
    }

    for(int b=0;b<in->batch;b++)
    {
        float *out_data = out->data[b];

        cblas_sgemm(CblasRowMajor,CblasTrans,CblasNoTrans,
                    mheight,mwidth,para->knl_num,
                    1.0f,
                       out_data,mheight,
                    kernel_data,mwidth,
                0.0,   res_data,mwidth);

        DirConvMatDataToTensor(res_data,res,b,para->knl_r,para->dir,para->y_stride,para->x_stride);
    }
}

max_stars_repo_name: dianjixz/Morn

max_issues_repo_name: dianjixz/Morn

max_forks_repo_name: dianjixz/Morn


text:
/*
 * This file is part of the CN24 semantic segmentation software,
 * copyright (C) 2015 Clemens-Alexander Brust (ikosa dot de at gmail dot com).
 *
 * For licensing information, see the LICENSE file included with this project.
 */
/**
 * @file MKLHelper.h
 * @brief Contains macros to help with various BLAS implementations.
 *
 * @author Clemens-Alexander Brust (ikosa dot de at gmail dot com)
 */

#ifndef CONV_MKLHELPER_H
#define CONV_MKLHELPER_H
#ifdef BUILD_BLAS

/**
 * Here we include our favorite BLAS library and define important functions
 * so that they are independent of the type
 */

#ifdef BLAS_ACCELERATE
#include <Accelerate/Accelerate.h>
#define INNERGEMM cblas_sgemm
#define INNERGEMV cblas_sgemv

#endif

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

#define INNERGEMM cblas_sgemm
#define INNERGEMV cblas_sgemv

#endif

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

#define INNERGEMM cblas_sgemm
#define INNERGEMV cblas_sgemv

#endif

#ifdef BLAS_MKL

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

#define INNERGEMM cblas_sgemm
#define INNERGEMV cblas_sgemv

#endif

#else

#include "Log.h"

#define INNERGEMM(...) { FATAL("GEMM function called without BLAS support!"); }
#define INNERGEMV(...) { FATAL("GEMV function called without BLAS support!"); }

#endif

#endif


text:
/*
 # This file is part of the Astrometry.net suite.
 # Licensed under a 3-clause BSD style license - see LICENSE
 */

/**
 Runs the verification procedure in stand-alone mode.
 */
#include "os-features.h"
#include "matchfile.h"
#include "matchobj.h"
#include "index.h"
#include "xylist.h"
#include "rdlist.h"
#include "log.h"
#include "errors.h"
#include "mathutil.h"

#include "verify.h"
#include "verify2.h"

#define SIGN(x) (((x) >= 0) ? (1) : (-1))

static const char* OPTIONS = "hvi:m:f:r:p";

static void print_help(const char* progname) {
    printf("Usage:   %s\n"
           "   -m <match-file>\n"
           "   -f <xylist-file>\n"
           "  (    -i <index-file>\n"
           "   OR  -r <index-rdls>\n"
           "  )\n"
           "   [-v]: verbose\n"
           "\n", progname);
}

int Npaths = 0;

static void explore_path(il** reflists, dl** problists, int i, int NT, int NR,
                         int* theta, double* logprobs,
                         anbool* refused, int mu,
                         double distractor, double logbg) {
    int j;
    double logprob;
    FILE* f = stderr;
    double logd = log(distractor + (1.0-distractor)*mu / (double)NR) + logbg;

    if (i == NT) {

        /*
         fprintf(f, "allpaths.append(array([");
         for (j=0; j<NT; j++)
         fprintf(f, "%g,", logprobs[j]);
         fprintf(f, "]))\n");
         */
        fprintf(f, "alllogprobs.append(%g)\n", logprobs[i-1]);
        Npaths++;
        return;
    }

    if (i == 0)
        logprob = 0.0;
    else
        logprob = logprobs[i-1];

    for (j=0; reflists[i] && j<il_size(reflists[i]); j++) {
        int refi;
        refi = il_get(reflists[i], j);
        if (refused[refi])
            continue;

        logprobs[i] = logprob + dl_get(problists[i], j) - logbg;
        theta[i] = refi;
        //fprintf(f, "plot([%i, %i], [%g, %g], 'r-')\n", i, i+1, logprob, logprobs[i]);
        //fprintf(f, "pathsx.append(%i)\npathsy.append(%g)\n", i+1, logprobs[i]);
        fprintf(f, "pathsx.append([%i, %i])\npathsy.append([%g, %g])\n", i, i+1, logprob, logprobs[i]);
        fprintf(f, "pathst.append(%i)\n", refi);
        refused[refi] = TRUE;
        explore_path(reflists, problists, i+1, NT, NR, theta, logprobs,
                     refused, mu+1, distractor, logbg);
        refused[refi] = FALSE;
    }

    logprobs[i] = logprob + logd - logbg;
    theta[i] = -1;
    //fprintf(f, "plot([%i, %i], [%g, %g], 'r-')\n", i, i+1, logprob, logprobs[i]);
    //fprintf(f, "pathsx.append(%i)\npathsy.append(%g)\n", i+1, logprobs[i]);
    fprintf(f, "pathsx.append([%i, %i])\npathsy.append([%g, %g])\n", i, i+1, logprob, logprobs[i]);
    fprintf(f, "pathst.append(%i)\n", -1);
    explore_path(reflists, problists, i+1, NT, NR, theta, logprobs,
                 refused, mu, distractor, logbg);
}

static void add_radial_and_tangential_correction(const double* in,
                                                 double r, double t,
                                                 const double* qc,
                                                 double* out,
                                                 int N) {
    int i;
    for (i=0; i<N; i++) {
        double rdir[2];
        // radial vector
        rdir[0] = in[2*i+0] - qc[0];
        rdir[1] = in[2*i+1] - qc[1];
        out[2*i+0] = in[2*i+0] - r * rdir[0] + t * rdir[1];
        out[2*i+1] = in[2*i+1] - r * rdir[1] - t * rdir[0];
    }
}

#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_vector_double.h>
#include "gslutils.h"

void find_cd_correction(const double* testxy, const double* sigma2s, int NT,
                        const int* theta, const double* refxy, int NR,
                        const double* crpix) {
    gsl_matrix *A;
    gsl_vector *B1, *B2, *X1, *X2;
    int M, N;
    int mu;
    int i;

    /*
     solve min(|B - A*X|^2) for X

     B: (refxy - crpix)_{x,y} / sigma
     X: CD matrix elements (a,b) and (c,d)
     A: (testxy - crpix)_{x,y} / sigma
     */

    mu = 0;
    for (i=0; i<NT; i++)
        if (theta[i] >= 0)
            mu++;
    // number of samples
    M = mu;
    // number of coefficients
    N = 2;

    A = gsl_matrix_alloc(M, N);
    B1 = gsl_vector_alloc(M);
    B2 = gsl_vector_alloc(M);

    mu = 0;
    for (i=0; i<NT; i++) {
        double w;
        if (theta[i] < 0)
            continue;
        w = 1.0 / sqrt(sigma2s[i]);
        gsl_matrix_set(A, mu, 0, w * (testxy[2*i]   - crpix[0]));
        gsl_matrix_set(A, mu, 1, w * (testxy[2*i+1] - crpix[1]));
        gsl_vector_set(B1, mu, w * (refxy[2*theta[i]]   - crpix[0]));
        gsl_vector_set(B2, mu, w * (refxy[2*theta[i]+1] - crpix[1]));
        mu++;
    }

    if (gslutils_solve_leastsquares_v(A, 2, B1, &X1, NULL, B2, &X2, NULL)) {
        ERROR("Failed to solve CD matrix correction\n");
        return;
    }

    logmsg("CD matrix:\n");
    logmsg("  %g\n", gsl_vector_get(X1, 0));
    logmsg("  %g\n", gsl_vector_get(X1, 1));
    logmsg("  %g\n", gsl_vector_get(X2, 0));
    logmsg("  %g\n", gsl_vector_get(X2, 1));

    gsl_matrix_free(A);
    gsl_vector_free(B1);
    gsl_vector_free(B2);
    gsl_vector_free(X1);
    gsl_vector_free(X2);
}

int main(int argc, char** args) {
    int argchar;
    int loglvl = LOG_MSG;
    char* indexfn = NULL;
    char* matchfn = NULL;
    char* xyfn = NULL;
    char* rdfn = NULL;

    index_t* index = NULL;
    matchfile* mf;
    MatchObj* mo;
    verify_field_t* vf;
    starxy_t* fieldxy;
    xylist_t* xyls;
    rdlist_t* rdls;

    double pix2 = 1.0;
    double distractors = 0.25;
    double fieldW=0, fieldH=0;
    double logbail = log(1e-100);
    double logkeep = log(1e12);
    double logaccept = LARGE_VAL;
    anbool growvariance = TRUE;
    anbool fake = FALSE;
    double logodds;
    anbool do_paths = FALSE;

    while ((argchar = getopt(argc, args, OPTIONS)) != -1)
        switch (argchar) {
        case 'p':
            do_paths = TRUE;
            break;
        case 'r':
            rdfn = optarg;
            break;
        case 'f':
            xyfn = optarg;
            break;
        case 'i':
            indexfn = optarg;
            break;
        case 'm':
            matchfn = optarg;
            break;
        case 'h':
            print_help(args[0]);
            exit(0);
        case 'v':
            loglvl++;
            break;
        }

    log_init(loglvl);

    if (!(indexfn || rdfn) || !matchfn || !xyfn) {
        logerr("You must specify (index (-i) or index rdls (-r)) and matchfile (-m) and field xylist (-f).\n");
        print_help(args[0]);
        exit(-1);
    }

    mf = matchfile_open(matchfn);
    if (!mf) {
        ERROR("Failed to read match file %s", matchfn);
        exit(-1);
    }

    xyls = xylist_open(xyfn);
    if (!xyls) {
        ERROR("Failed to open xylist %s", xyfn);
        exit(-1);
    }
    // don't need these...
    xylist_set_include_flux(xyls, FALSE);
    xylist_set_include_background(xyls, FALSE);
    fieldW = xylist_get_imagew(xyls);
    fieldH = xylist_get_imageh(xyls);

    logmsg("Field W,H = %g, %g\n", fieldW, fieldH);

    mo = matchfile_read_match(mf);
    if (!mo) {
        ERROR("Failed to read object from match file.");
        exit(-1);
    }
    mo->wcstan.imagew = fieldW;
    mo->wcstan.imageh = fieldH;

    fieldxy = xylist_read_field(xyls, NULL);
    if (!fieldxy) {
        ERROR("Failed to read a field from xylist %s", xyfn);
        exit(-1);
    }

    if (indexfn) {
        index = index_load(indexfn, 0);
        if (!index) {
            ERROR("Failed to open index %s", indexfn);
            exit(-1);
        }

        pix2 += square(index->meta.index_jitter / mo->scale);

    } else {
        double indexjitter;

        rdls = rdlist_open(rdfn);
        if (!rdls) {
            ERROR("Failed to open rdlist %s", rdfn);
            exit(-1);
        }

        // HACK
        indexjitter = 1.0; // arcsec.
        pix2 += square(indexjitter / mo->scale);
    }

    logmsg("Pixel jitter: %g pix\n", sqrt(pix2));

    vf = verify_field_preprocess(fieldxy);

    if (index) {
        mo->logodds = 0.0;
        mo->dimquads = index_get_quad_dim(index);
        verify_hit(index->starkd, index->meta.cutnside,
                   mo, NULL, vf,
                   pix2, distractors, fieldW, fieldH,
                   logbail, logkeep, logaccept, growvariance,
                   fake);

        logodds = mo->logodds;

        index_close(index);

    } else {
        int cutnside;
        int cutnsweeps;
        int indexid;
        int uni_nw, uni_nh;
        int* perm;
        double* testxy;
        double* refxy;
        int i, j, NT, NR;
        double* sigma2s = NULL;
        rd_t* rd;
        double effA;
        double qc[2], Q2;

        // -get reference stars
        rd = rdlist_read_field(rdls, NULL);
        if (!rd) {
            ERROR("Failed to read rdls field");
            exit(-1);
        }
        NR = rd_n(rd);
        refxy = malloc(2 * NR * sizeof(double));
        for (i=0; i<NR; i++) {
            double ra, dec;
            ra  = rd_getra (rd, i);
            dec = rd_getdec(rd, i);
            if (!tan_radec2pixelxy(&(mo->wcstan), ra, dec, refxy + 2*i, refxy + 2*i + 1)) {
                ERROR("rdls point projects to wrong side of sphere!");
                exit(-1);
            }
        }
        // -remove the ref star closest to each quad star.
        for (i=0; i<mo->dimquads; i++) {
            double qxy[2];
            int besti = -1;
            double bestd2 = LARGE_VAL;
            if (!tan_xyzarr2pixelxy(&(mo->wcstan), mo->quadxyz + 3*i, qxy, qxy+1)) {
                ERROR("quad star projects to wrong side of sphere!");
                exit(-1);
            }
            logmsg("Ref quad star %i is at (%.1f, %.1f)\n", i, qxy[0], qxy[1]);
            for (j=0; j<NR; j++) {
                double d2 = distsq(qxy, refxy + 2*j, 2);
                if (d2 < bestd2) {
                    bestd2 = d2;
                    besti = j;
                }
            }
            logmsg("Ref star %i is closest: (%.1f, %.1f)\n", besti, refxy[2*besti+0], refxy[2*besti+1]);
            // remove it!
            memmove(refxy + 2*besti, refxy + 2*(besti + 1),
                    2*(NR - besti - 1) * sizeof(double));
            NR--;
        }
        logmsg("Reference stars: %i\n", NR);

        indexid = mo->indexid;
        if (index_get_missing_cut_params(indexid, &cutnside, &cutnsweeps, NULL, NULL, NULL)) {
            ERROR("Failed to get index cut parameters for index id %i", indexid);
            exit(-1);
        }

        verify_get_quad_center(vf, mo, qc, &Q2);

        verify_apply_ror(refxy, NULL, &NR, cutnside, mo,
                         vf, pix2, distractors, fieldW, fieldH,
                         growvariance, fake,
                         &testxy, &sigma2s, &NT, &perm, &effA, &uni_nw, &uni_nh);

        /*{
         double d = distractors;
         // Predicted optimal number of reference stars:
         int mmax = (int)round(exp(log(effA*(1-d)/(2*M_PI*pix2)) + d*log(d)/(1-d) + (M_PI*Q2 + 1)/effA * log(M_PI*Q2 / (M_PI*Q2 + effA))));
         logmsg("mmax = %i\n", mmax);
         logmsg("first term: %g\n", effA*(1-d)/(2*M_PI*pix2));
         logmsg("second term: %g\n", exp(d*log(d) / (1-d)));
         logmsg("third term: %g\n", exp((M_PI*Q2 + 1)/effA * log(M_PI*Q2 / (M_PI*Q2 + effA))));

         // Predicted number of reference stars to allow the
         // accept threshold to be reached.
         double t1 = d*log(d) + (1-d)*(log(effA*(1-d)/(2*M_PI*pix2)) + (M_PI*Q2 + 1)/effA * log(M_PI*Q2 / (M_PI*Q2 + effA)));
         for (i=1; i<1000000; i++) {
         double logM = i*t1 - i*(1-d)*log(i);
         if (logM > logkeep) {
         logmsg("m = %i: M = %g\n", i, exp(logM));
         break;
         }
         }
         //NR = MIN(NR, 2 * i + 10);
         //logmsg("Setting NR to %i\n", NR);
         }*/

        FILE* f = stderr;

        fprintf(f, "distractor = %g\nNR=%i\nNT=%i\n", distractors, NR, NT);
        fprintf(f, "W=%i\nH=%i\n", (int)fieldW, (int)fieldH);
        fprintf(f, "effA=%g\n", effA);
        fprintf(f, "sig2=%g\n", pix2);

        fprintf(f, "quadxy = array([");
        for (i=0; i<mo->dimquads; i++)
            fprintf(f, "[%g,%g],", mo->quadpix[2*i+0], mo->quadpix[2*i+1]);
        fprintf(f, "])\n");

        fprintf(f, "testxy = array([");
        for (i=0; i<NT; i++)
            fprintf(f, "[%g,%g],", testxy[2*i+0], testxy[2*i+1]);
        fprintf(f, "])\n");

        fprintf(f, "sigmas = array([");
        for (i=0; i<NT; i++)
            fprintf(f, "%g,", sqrt(sigma2s[i]));
        fprintf(f, "])\n");

        double* rs2 = verify_compute_sigma2s_arr(refxy, NR, qc, Q2, pix2, !fake);
        fprintf(f, "refsigmas = array([");
        for (i=0; i<NR; i++)
            fprintf(f, "%g,", sqrt(rs2[i]));
        fprintf(f, "])\n");
        free(rs2);

        fprintf(f, "refxy = array([");
        for (i=0; i<NR; i++)
            fprintf(f, "[%g,%g],", refxy[2*i+0], refxy[2*i+1]);
        fprintf(f, "])\n");

        fprintf(f, "cutx = array([");
        for (i=0; i<=uni_nw; i++)
            fprintf(f, "%g,", i * fieldW / (float)uni_nw);
        fprintf(f, "])\n");

        fprintf(f, "cuty = array([");
        for (i=0; i<=uni_nh; i++)
            fprintf(f, "%g,", i * fieldH / (float)uni_nh);
        fprintf(f, "])\n");

        double* all_logodds;
        int* theta;
        int besti;
        double worst;

        logodds = verify_star_lists(refxy, NR, testxy, sigma2s, NT,
                                    effA, distractors, logbail, logaccept,
                                    &besti, &all_logodds, &theta, &worst);

        fprintf(f, "besti = %i\n", besti);

        fprintf(f, "worstlogodds = %g\n", worst);

        fprintf(f, "logodds = array([");
        for (i=0; i<NT; i++)
            fprintf(f, "%g,", all_logodds[i]);
        fprintf(f, "])\n");

        fprintf(f, "theta = array([");
        for (i=0; i<NT; i++)
            fprintf(f, "%i,", theta[i]);
        fprintf(f, "])\n");

        // compare observed sigmas to expected...
        fprintf(f, "obssigmas=array([");
        for (i=0; i<NT; i++) {
            double d2, r2;
            if (theta[i] < 0)
                continue;
            d2 = distsq(testxy + 2*i, refxy + 2*theta[i], 2);
            r2 = distsq(testxy + 2*i, qc, 2);
            fprintf(f, "[%g,%g,%g],", sigma2s[i], d2, r2/Q2);
        }
        fprintf(f, "])\n");

        find_cd_correction(testxy, sigma2s, NT, theta, refxy, NR, qc);

        {
            // introduce known radial and tangential terms and see if we can recover them...
            //add_radial_and_tangential_correction(testxy, -0.01, -0.01, qc, testxy, NT);

            // What is the ML correction to rotation and scaling?
            // (shear, translation?  distortion?)
            // -> may need all matches, not just nearest neighbour, to
            //    do this correctly.
            double racc = 0, tacc = 0;
            int mu = 0;

            for (i=0; i<NT; i++) {
                double dxy[2];
                double rdir[2];
                double R2, ddotr;
                double dr[2];
                double dt[2];
                double fr, ft;

                if (theta[i] == -1)
                    continue;
                mu++;

                // jitter vector
                dxy[0] = testxy[2*i+0] - refxy[2*theta[i]+0];
                dxy[1] = testxy[2*i+1] - refxy[2*theta[i]+1];
                // radial vector (this should perhaps be to the ref star, not test)
                rdir[0] = testxy[2*i+0] - qc[0];
                rdir[1] = testxy[2*i+1] - qc[1];
                //
                R2 = rdir[0]*rdir[0] + rdir[1]*rdir[1];
                ddotr = (dxy[0]*rdir[0] + dxy[1]*rdir[1]);
                // jitter vector projected onto radial vector.
                dr[0] = ddotr * rdir[0] / R2;
                dr[1] = ddotr * rdir[1] / R2;
                // tangential
                dt[0] = dxy[0] - dr[0];
                dt[1] = dxy[1] - dr[1];

                assert(fabs(dr[0] + dt[0] - dxy[0]) < 1e-10);
                assert(fabs(dr[1] + dt[1] - dxy[1]) < 1e-10);

                // fractional change in radial, tangential components.
                fr = SIGN(ddotr) * sqrt((dr[0]*dr[0] + dr[1]*dr[1]) / R2);
                ft = SIGN(rdir[0]*dt[1] - rdir[1]*dt[0]) * sqrt((dt[0]*dt[0] + dt[1]*dt[1]) / R2);

                racc += fr;
                tacc += ft;
            }
            racc /= (double)mu;
            tacc /= (double)mu;

            logmsg("Radial correction: %g\n", racc);
            logmsg("Tangential correction: %g\n", tacc);

            logmsg("Log-odds: %g\n", logodds);

            // Rotate and scale the test stars...
            double* t2xy = malloc(NT * 2 * sizeof(double));
            add_radial_and_tangential_correction(testxy, racc, tacc, qc, t2xy, NT);
            double logodds2 = verify_star_lists(refxy, NR, t2xy, sigma2s, NT,
                                                effA, distractors, logbail, logaccept,
                                                NULL, NULL, NULL, NULL);
            logmsg("Log-odds 2: %g\n", logodds2);

            fprintf(f, "t2xy = array([");
            for (i=0; i<NT; i++)
                fprintf(f, "[%g,%g],", t2xy[2*i+0], t2xy[2*i+1]);
            fprintf(f, "])\n");

            free(t2xy);

        }

        if (do_paths) {
            il** reflist;
            dl** problist;

            NT = besti+1;

            logmsg("Finding all matches...\n");

            verify_get_all_matches(refxy, NR, testxy, sigma2s, NT,
                                   effA, distractors, 5.0, 0.5,
                                   &reflist, &problist);

            /*
             --reflist contains one list per test star, containing the
             indices of reference stars within nsigma and within
             limit of the distractor rate.

             -a "regular" conflict occurs when one reference star
             appears in more than one list.

             -a "ref" conflict occurs when a list has more than one
             element in it.

             -each star has some "clique" of stars that it can
             interact with (transitive connectivity of the 'nearby'
             graph).  These will usually be small, but might not
             be...  We can compute analytically the sums over small,
             simple groups, but more complex ones will be very hairy.

             */

            double np;
            np = 1.0;
            for (i=0; i<NT; i++) {
                if (!reflist[i])
                    continue;
                np *= (1.0 + il_size(reflist[i]));
            }
            logmsg("Number of paths: about %g\n", np);

            fprintf(f, "allpaths=[]\n");
            fprintf(f, "clf()\n");
            fprintf(f, "alllogprobs = []\n");
            fprintf(f, "pathsx = []\n");
            fprintf(f, "pathsy = []\n");
            fprintf(f, "pathst = []\n");

            int theta[NT];
            double logprobs[NT];
            anbool refused[NR];
            for (i=0; i<NR; i++)
                refused[i] = FALSE;

            Npaths = 0;
            logmsg("Finding all paths...\n");
            explore_path(reflist, problist, 0, NT, NR, theta, logprobs, refused, 0, distractors, log(1.0/effA));
            logmsg("Number of paths: %i\n", Npaths);

            fprintf(f, "pathsx = array(pathsx)\npathsy = array(pathsy)\n");

            //fprintf(f, "axis([0, %i, -100, 100])\n", NT);

            for (i=0; i<NT; i++) {
                il_free(reflist[i]);
                dl_free(problist[i]);
            }
            free(reflist);
            free(problist);
        }

        free(theta);
        free(all_logodds);
        free(sigma2s);
        free(testxy);
        free(refxy);

        rd_free(rd);
        rdlist_close(rdls);
    }

    logmsg("Logodds: %g\n", logodds);
    logmsg("Odds: %g\n", logodds);

    verify_field_free(vf);
    starxy_free(fieldxy);

    xylist_close(xyls);
    matchfile_close(mf);

    return 0;
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet