text:
/* -*- linux-c -*- */
/* fewbody.c

   Copyright (C) 2002-2004 John M. Fregeau

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

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

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/times.h>
#include <unistd.h>
#include <getopt.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>
#include "fewbody.h"

int fb_debug = 0;

fb_ret_t fewbody(fb_input_t input, fb_hier_t *hier, double *t)
{
	int i, j, k, status, done=0, forceclassify=0, restart, restep;
	long clk_tck;
	double s, slast, sstop=FB_SSTOP, tout, h=FB_H, *y, texpand, tnew, R[3];
	double Ei, E, Lint[3], Li[3], L[3], DeltaL[3];
	double s2, s2prev=GSL_POSINF, s2prevprev=GSL_POSINF, s2minprev=GSL_POSINF, s2max=0.0, s2min;
	struct tms firsttimebuf, currtimebuf;
	clock_t firstclock, currclock;
	fb_hier_t phier;
	fb_ret_t retval;
	fb_nonks_params_t nonks_params;
	fb_ks_params_t ks_params;
	char string1[FB_MAX_STRING_LENGTH], string2[FB_MAX_STRING_LENGTH], logentry[FB_MAX_LOGENTRY_LENGTH];
	const gsl_odeiv_step_type *ode_type=gsl_odeiv_step_rk8pd;
	gsl_odeiv_step *ode_step;
	gsl_odeiv_control *ode_control;
	gsl_odeiv_evolve *ode_evolve;
	gsl_odeiv_system ode_sys;

	/* initialize a few things */
	fb_init_hier(hier);
	retval.iclassify = 0;
	retval.Rmin = FB_RMIN;
	retval.Rmin_i = -1;
	retval.Rmin_j = -1;
	retval.Nosc = 0;
	strncpy(logentry, input.firstlogentry, FB_MAX_LOGENTRY_LENGTH);

	/* set up the perturbation tree, initially flat */
	phier.nstarinit = hier->nstar;
	phier.nstar = hier->nstar;
	fb_malloc_hier(&phier);
	fb_init_hier(&phier);
	for (i=0; i<phier.nstar; i++) {
		fb_objcpy(&(phier.hier[phier.hi[1]+i]), &(hier->hier[hier->hi[1]+i]));
	}

	/* initialize GSL integration routine */
	if (input.ks) {
		ode_step = gsl_odeiv_step_alloc(ode_type, 8 * (hier->nstar * (hier->nstar - 1) / 2) + 1);
		ode_control = gsl_odeiv_control_y_new(input.absacc, input.relacc);
		ode_evolve = gsl_odeiv_evolve_alloc(8 * (hier->nstar * (hier->nstar - 1) / 2) + 1);
		ode_sys.function = fb_ks_func;
		ode_sys.jacobian = NULL;
		ode_sys.dimension = 8 * (hier->nstar * (hier->nstar - 1) / 2) + 1;
		ode_sys.params = &ks_params;
	} else {
		ode_step = gsl_odeiv_step_alloc(ode_type, 6 * hier->nstar);
		ode_control = gsl_odeiv_control_y_new(input.absacc, input.relacc);
		ode_evolve = gsl_odeiv_evolve_alloc(6 * hier->nstar);
		ode_sys.function = fb_nonks_func;
		ode_sys.jacobian = fb_nonks_jac;
		ode_sys.dimension = 6 * hier->nstar;
		ode_sys.params = &nonks_params;
	}

	/* set parameters for integrator */
	if (input.ks) {
		ks_params.nstar = hier->nstar;
		ks_params.kstar = ks_params.nstar*(ks_params.nstar-1)/2;
		fb_malloc_ks_params(&ks_params);
		fb_init_ks_params(&ks_params, *hier);
	} else {
		nonks_params.nstar = hier->nstar;
		fb_malloc_nonks_params(&nonks_params);
		fb_init_nonks_params(&nonks_params, *hier);
	}

	/* set the initial conditions in y_i */
	if (input.ks) {
		y = fb_malloc_vector(8*ks_params.kstar+1);
		y[0] = *t;
		fb_euclidean_to_ks(phier.obj, y, ks_params.nstar, ks_params.kstar);
		s = 0.0;
	} else {
		y = fb_malloc_vector(6*nonks_params.nstar);
		fb_euclidean_to_nonks(phier.obj, y, nonks_params.nstar);
		s = *t;
	}

	/* store the initial energy and angular momentum */
	Ei = fb_petot(&(hier->hier[hier->hi[1]]), hier->nstar) + fb_ketot(&(hier->hier[hier->hi[1]]), hier->nstar) +
		fb_einttot(&(hier->hier[hier->hi[1]]), hier->nstar);
	fb_angmom(&(hier->hier[hier->hi[1]]), hier->nstar, Li);
	fb_angmomint(&(hier->hier[hier->hi[1]]), hier->nstar, Lint);
	for (i=0; i<3; i++) {
		Li[i] += Lint[i];
	}

	/* integrate along */
	fb_dprintf("fewbody: integrating...\n");

	done = 0;
	retval.count = 0;
	tout = *t;
	texpand = 0.0;
	clk_tck = sysconf(_SC_CLK_TCK);
	firstclock = times(&firsttimebuf);
	retval.tcpu = 0.0;
	while (*t < input.tstop && retval.tcpu < input.tcpustop && !done) {
		/* take one step */
		slast = s;
		status = gsl_odeiv_evolve_apply(ode_evolve, ode_control, ode_step, &ode_sys, &s, sstop, &h, y);
		if (status != GSL_SUCCESS) {
			break;
		}

		/* set objects' positions and velocities in phier */
		if (input.ks) {
			tnew = y[0];
			fb_ks_to_euclidean(y, phier.obj, ks_params.nstar, ks_params.kstar);
		} else {
			tnew = s;
			fb_nonks_to_euclidean(y, phier.obj, nonks_params.nstar);
		}

		/* restart: do we need to restart the integrator?
		   restep: do we need to repeat the last integration step?
		   forceclassify: do we need to force a call to the classify() routine? */
		restart = 0;
		restep = 0;
		forceclassify = 0;

		/* see if we need to expand or collapse the perturbation hierarchy */
		if (fb_expand(&phier, tnew, input.tidaltol)) {
			texpand = tnew;
			s = slast;
			restart = 1;
			restep = 1;
			for (i=0; i<hier->nstar; i++) {
				for (k=0; k<3; k++) {
					phier.hier[phier.hi[1]+i].x[k] = hier->hier[hier->hi[1]+i].x[k];
					phier.hier[phier.hi[1]+i].v[k] = hier->hier[hier->hi[1]+i].v[k];
				}
			}
			fb_elkcirt(&phier, *t);
		} else if (tnew >= texpand && fb_collapse(&phier, tnew, input.tidaltol)) {
			*t = tnew;
			/* if there is only one object, then it's stable---force classify() */
			if (phier.nobj == 1) {
				forceclassify = 1;
			} else {
				restart = 1;
			}
		} else {
			*t = tnew;
		}

		/* if we're not repeating the previous integration step, then do physics */
		if (!restep) {
			/* trickle down so updated information is in hier */
			fb_trickle(&phier, *t);
			for (i=0; i<hier->nstar; i++) {
				for (k=0; k<3; k++) {
					hier->hier[hier->hi[1]+i].x[k] = phier.hier[phier.hi[1]+i].x[k];
					hier->hier[hier->hi[1]+i].v[k] = phier.hier[phier.hi[1]+i].v[k];
				}
			}

			/* update Rmin (closest approach) */
			for (i=0; i<hier->nstar-1; i++) {
				for (j=i+1; j<hier->nstar; j++) {
					for (k=0; k<3; k++) {
						R[k] = hier->hier[hier->hi[1]+i].x[k] - hier->hier[hier->hi[1]+j].x[k];
					}
					if (fb_mod(R) < retval.Rmin) {
						retval.Rmin = fb_mod(R);
						retval.Rmin_i = i;
						retval.Rmin_j = j;
					}
				}
			}

			/* do physical collisions */
			if (fb_collide(hier, input.fexp)) {
				/* initialize phier to a flat tree */
				phier.nstar = hier->nstar;
				fb_init_hier(&phier);
				for (i=0; i<phier.nstar; i++) {
					fb_objcpy(&(phier.hier[phier.hi[1]+i]), &(hier->hier[hier->hi[1]+i]));
				}

				/* if there is only one object, then it's stable---force classify() */
				if (phier.nobj == 1) {
					restart = 0;
					forceclassify = 1;
				} else {
					restart = 1;
				}
			}

			/* check for oscillations in s^2 */
			/* first calculate s^2 */
			s2 = 0.0;
			for (i=0; i<hier->nstar-1; i++) {
				for (j=i+1; j<hier->nstar; j++) {
					for (k=0; k<3; k++) {
						s2 += fb_sqr(hier->hier[hier->hi[1]+i].x[k] - hier->hier[hier->hi[1]+j].x[k]);
					}
				}
			}
			/* local min */
			if (s2 > s2prev && s2prev < s2prevprev) {
				s2min = s2prev;
				/* increment Nosc if there is a valid oscillation */
				if (s2max >= 2.0*s2minprev && s2max >= 2.0*s2min) {
					retval.Nosc++;
					fb_dprintf("fewbody: oscillation in s^2 detected: s2max/s2minprev=%g s2max/s2min=%g: Nosc=%d\n", s2max/s2minprev, s2max/s2min, retval.Nosc);
					s2minprev = s2min;
				} else {
					fb_dprintf("fewbody: oscillation in s^2 not big enough: s2max/s2minprev=%g s2max/s2min=%g\n", s2max/s2minprev, s2max/s2min);
					/* that pesky first min... */
					if (s2max == 0.0) {
						s2minprev = s2min;
					}
				}
			}
			/* local max: just set value and wait for next min */
			if (s2 < s2prev && s2prev > s2prevprev) {
				s2max = s2prev;
			}
			/* set previous values */
			s2prevprev = s2prev;
			s2prev = s2;

			/* see if we're done */
			if (retval.count % input.ncount == 0 || forceclassify) {
				status = fb_classify(hier, *t, input.tidaltol);
				retval.iclassify++;
				fb_dprintf("fewbody: current status:  t=%.6g  %s  (%s)\n",
					   *t, fb_sprint_hier(*hier, string1),
					   fb_sprint_hier_hr(*hier, string2));
				snprintf(&(logentry[strlen(logentry)]), FB_MAX_LOGENTRY_LENGTH-strlen(logentry),
					 "  current status:  t=%.6g  %s  (%s)\n", *t, fb_sprint_hier(*hier, string1),
					 fb_sprint_hier_hr(*hier, string2));
				if (status) {
					done = 1;
				}
			}

			/* print stuff if necessary */
			if (input.Dflag == 1 && (*t >= tout || done)) {
				tout = *t + input.dt;
				fb_print_story(&(hier->hier[hier->hi[1]]), hier->nstar, *t, logentry);
			}
		}

		/* restart integrator if necessary */
		if (restart) {
			fb_dprintf("fewbody: restarting integrator: nobj=%d count=%ld\n", phier.nobj, retval.count);
			fb_free_vector(y);
			if (input.ks) {
				fb_free_ks_params(ks_params);
				ks_params.nstar = phier.nobj;
				ks_params.kstar = ks_params.nstar*(ks_params.nstar-1)/2;
				fb_malloc_ks_params(&ks_params);
				fb_init_ks_params(&ks_params, phier);

				y = fb_malloc_vector(8*ks_params.kstar+1);
				y[0] = *t;
				fb_euclidean_to_ks(phier.obj, y, ks_params.nstar, ks_params.kstar);
			} else {
				fb_free_nonks_params(nonks_params);
				nonks_params.nstar = phier.nobj;
				fb_malloc_nonks_params(&nonks_params);
				fb_init_nonks_params(&nonks_params, phier);

				y = fb_malloc_vector(6*nonks_params.nstar);
				fb_euclidean_to_nonks(phier.obj, y, nonks_params.nstar);
			}

			/* and re-allocate integrator */
			gsl_odeiv_evolve_free(ode_evolve);
			gsl_odeiv_control_free(ode_control);
			gsl_odeiv_step_free(ode_step);

			/* re-initialize integrator */
			if (input.ks) {
				ode_step = gsl_odeiv_step_alloc(ode_type, 8*ks_params.kstar+1);
				ode_control = gsl_odeiv_control_y_new(input.absacc, input.relacc);
				ode_evolve = gsl_odeiv_evolve_alloc(8*ks_params.kstar+1);
				ode_sys.dimension = 8*ks_params.kstar+1;
			} else {
				ode_step = gsl_odeiv_step_alloc(ode_type, 6*nonks_params.nstar);
				ode_control = gsl_odeiv_control_y_new(input.absacc, input.relacc);
				ode_evolve = gsl_odeiv_evolve_alloc(6*nonks_params.nstar);
				ode_sys.dimension = 6*nonks_params.nstar;
			}
		}

		/* update variables that change on every integration step */
		retval.count++;
		currclock = times(&currtimebuf);
		retval.tcpu = ((double) (currtimebuf.tms_utime + currtimebuf.tms_stime - firsttimebuf.tms_utime - firsttimebuf.tms_stime))/((double) clk_tck);
	}

	/* do final classification */
	retval.retval = fb_classify(hier, *t, input.tidaltol);
	retval.iclassify++;
	fb_dprintf("fewbody: current status:  t=%.6g  %s  (%s)\n",
		   *t, fb_sprint_hier(*hier, string1),
		   fb_sprint_hier_hr(*hier, string2));
	snprintf(&(logentry[strlen(logentry)]), FB_MAX_LOGENTRY_LENGTH-strlen(logentry),
		 "  current status:  t=%.6g  %s  (%s)\n", *t, fb_sprint_hier(*hier, string1),
		 fb_sprint_hier_hr(*hier, string2));

	/* print final story */
	if (input.Dflag == 1) {
		fb_print_story(&(hier->hier[hier->hi[1]]), hier->nstar, *t, logentry);
	}

	fb_dprintf("fewbody: final: phier.nobj = %d\n", phier.nobj);

	E = fb_petot(&(hier->hier[hier->hi[1]]), hier->nstar) + fb_ketot(&(hier->hier[hier->hi[1]]), hier->nstar) +
		fb_einttot(&(hier->hier[hier->hi[1]]), hier->nstar);
	fb_angmom(&(hier->hier[hier->hi[1]]), hier->nstar, L);
	fb_angmomint(&(hier->hier[hier->hi[1]]), hier->nstar, Lint);
	for (i=0; i<3; i++) {
		L[i] += Lint[i];
		DeltaL[i] = L[i] - Li[i];
	}

	/* free GSL stuff */
	gsl_odeiv_evolve_free(ode_evolve);
	gsl_odeiv_control_free(ode_control);
	gsl_odeiv_step_free(ode_step);

	/* free our own stuff */
	fb_free_vector(y);
	fb_free_hier(phier);

	if (input.ks) {
		fb_free_ks_params(ks_params);
	} else {
		fb_free_nonks_params(nonks_params);
	}

	/* done! */
	retval.DeltaE = E-Ei;
	retval.DeltaEfrac = E/Ei-1.0;
	retval.DeltaL = fb_mod(DeltaL);
	retval.DeltaLfrac = fb_mod(DeltaL)/fb_mod(Li);
	return(retval);
}

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos


text:
/*
 * MotorControl.h
 *
 *  Created on: 02/01/2014
 *      Author: james
 */

#ifndef MOTORCONTROL_H_
#define MOTORCONTROL_H_
#include "Arduino.h"
#include "MotorControl/Motor.h"
#include "MotorControl/Encoder.h"
#include "MathHelper/ComponentsFromVector.h"
#include "MathHelper/Vector.h"
#include "MotorControl/SpeedPWMStruct.h"
#include <gsl_math.h>

#define oneonroot3 0.5773502691896258
#define onethird 0.3333333333333333

class MotorControl {
public:
	MotorControl(Motor parMotors[], Encoder parEncoders[], SpeedPWMStruct speedpwmstruct);
	MotorControl();
	float getDesiredIndividualSpeed(int motor_num, float rotation_rate);
	void setOverallDesiredVector(float magnitude, float angl);
	float getNewAimFromFeedback(float desiredSpeed, float actualSpeed);
	int getPWMFromSpeed(float speed);
	void updateMotors();
	virtual ~MotorControl();
	void setMotorPWM(int motorNum, int, boolean FOR)
	static const int clicksperrev;
	static const int millispersec;
	Motor motors[3];
	Encoder encoders[3];
	float getIndividualSpeed(int encoder_num);
private:
	ComponentsFromVector motorSpeedComponents[3];
	Vector overallDesiredVector;
	SpeedPWMStruct speedpwmstruct;
	static const float mtrconst;
};

#endif /* MOTORCONTROL_H_ */

ext: h

max_stars_repo_name: JamesMBartlett/RobotSoccer

max_issues_repo_name: JamesMBartlett/RobotSoccer

max_forks_repo_name: JamesMBartlett/RobotSoccer


text:
#include <assert.h>
#include <complex.h>
#include <gsl/gsl_const_mksa.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <malloc.h>
#include <math.h>
#include <nlopt.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "tsvread.h"

typedef struct {
  gsl_spline *spline;
  gsl_interp_accel *acc;
  int flags;
  double d;
} f_param_t;

double opt_func(unsigned n, const double *x, double *grad, void *f_data) {
  f_param_t *fparam = (f_param_t *)f_data;
  gsl_spline *spline_nu = fparam->spline;
  gsl_interp_accel *acc_nu = fparam->acc;
  int flags = fparam->flags;
  double d = fparam->d;
  return gsl_spline_eval(spline_nu, x[0], acc_nu);
}

int main(int argc, char **argv) {
  /* open file for reading data */
  tsv_t *tsv = tsv_fread(argv[1]);

  if (tsv->columns != 2) {
    fprintf(stderr, "input file must be 2 colums tab seperated\n");
    return -1;
  }

  double *x_nu = (double *)malloc(tsv->rows * sizeof(double));
  double *y_nu = (double *)malloc(tsv->rows * sizeof(double));
  assert(x_nu != NULL);
  assert(y_nu != NULL);

  unsigned int i;
  for (i = 0; i < tsv->rows; i++) {
    x_nu[i] = tsv->data[i][0];
    y_nu[i] = tsv->data[i][1];
  }

  gsl_interp_accel *acc_nu = gsl_interp_accel_alloc();
  gsl_spline *spline_nu = gsl_spline_alloc(gsl_interp_cspline, tsv->rows);
  gsl_spline_init(spline_nu, x_nu, y_nu, tsv->rows);

  double x[1];  /* initial guess */
  double lb[1]; /* lower bound */
  double ub[1]; /* upper bound */
  double minf = 0;
  lb[0] = x_nu[0];
  ub[0] = x_nu[tsv->rows - 1];
  x[0] = x_nu[tsv->rows / 2];

  f_param_t *fparam = (f_param_t *)malloc(1 * sizeof(f_param_t));
  fparam->spline = spline_nu;
  fparam->acc = acc_nu;
  fparam->flags = 0;
  void *fvparam;
  fvparam = fparam;

  nlopt_opt opt;
  opt = nlopt_create(NLOPT_LN_NELDERMEAD, 1);
  nlopt_set_lower_bounds(opt, lb);
  nlopt_set_upper_bounds(opt, ub);
  nlopt_set_max_objective(opt, opt_func, fvparam);
  nlopt_set_xtol_rel(opt, 1e-16);
  // nlopt_set_ftol_abs(opt,1e-6);
  // nlopt_set_maxtime(opt,4);

  int err;
  double max_nu;
  err = nlopt_optimize(opt, x, &minf);
  if (err < 0) {
    fprintf(stderr, "nlopt failed! %d\n", err);
    fprintf(stderr, "lower: %g\n", lb[0]);
    fprintf(stderr, "upper: %g\n", ub[0]);
    fprintf(stderr, "guess: %g\n", x[0]);
    if (lb[0] > ub[0]) {
      fprintf(stderr, "lower bigger than upper\n");
    }
  } else {
    max_nu = gsl_spline_eval(spline_nu, x[0], acc_nu);
    printf("Maximum: %.20g %.20g \n", x[0], max_nu);
  }

  nlopt_destroy(opt);
  tsv_free(tsv);
  free(x_nu);
  free(y_nu);

  return 0;
}

ext: c

max_stars_repo_name: AaronWebster/directdebye

max_issues_repo_name: AaronWebster/directdebye

max_forks_repo_name: AaronWebster/directdebye


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

#include "../inc/knnring.h"

knnresult kNN(double *X, double *Y, int n, int m, int d, int k){

	knnresult k_nearest;

	//initalize the result's variable
	k_nearest.nidx  = (int *)malloc(m*k*sizeof(int));
	k_nearest.ndist = (double *)malloc(m*k*sizeof(double));
	k_nearest.m = m;
	k_nearest.k = k;

  //compute the distance Matrix D
	double *D = compute_D(X,Y,n,m,d,k);

  //select k-nearest neighbors for every point in QUERY
	k_select(D,k_nearest.ndist, k_nearest.nidx,n,m,k);
	return k_nearest;
}

double *compute_D(double *X, double *Y, int n, int m, int d, int k){

	//compute distances using
	//D = (X*X)ee_t - 2XY_t + ee_t(Y*Y)_t
	double *tempD   = (double *)malloc(n*m*sizeof(double));
	double *e_1     = (double *)malloc(d*m*sizeof(double));
	double *XX      = (double *)malloc(n*d*sizeof(double));

	for(int i=0;i<d*m;i++)
		e_1[i]=1;

	for(int i=0;i<n*d;i++)
		XX[i] = X[i]*X[i];

	cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n,m,d,1,XX,d,e_1,m,0,tempD,m);

	free(e_1);
	free(XX);

	double *e_2 = (double *)malloc(n*d*sizeof(double));
	double *YY = (double *)malloc(m*d*sizeof(double));

	for(int i=0;i<n*d;i++)
		e_2[i]=1;

	for(int i=0;i<m*d;i++)
		YY[i] = Y[i]*Y[i];

	cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, n,m,d,-2,X,d,Y,d,1,tempD,m);
	cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, n,m,d,1,e_2,d,YY,d,1,tempD,m);

	for(int i=0;i<n*m;i++)
		tempD[i] = sqrt(fabs(tempD[i]));

	free(e_2);
	free(YY);

	return tempD;
}

void k_select(double *D, double *ndist, int *nidx,int n, int m, int k){
	for(int j=0;j<m;j++){
		for(int i=0;i<k;i++){
			ndist[j*k+i] = D[j+i*m];
			nidx[j*k+i]  = i;
		}
		//sort the initial k neighbors
		quick_sort(ndist,nidx,j*k,k*(j+1)-1);
	}

	//search the rest points to find it there are closer neighbors
	for(int j=0;j<m;j++){
		for(int i=k;i<n;i++){
			if(D[j+i*m] < ndist[k*(j+1)-1]){
				ndist[k*(j+1)-1] = D[j+i*m];
				nidx[k*(j+1)-1] = i;
				quick_sort(ndist,nidx,j*k,(j+1)*k-1);
			}
		}
	}
}

void quick_sort(double *k_dist, int *k_neigh, int low, int high){
	if(low<high){

		//actuall sorting an element
		int pi = partition(k_dist,k_neigh,low,high);

		//recursivelly call quick_sort
		quick_sort(k_dist,k_neigh,low,pi-1);
		quick_sort(k_dist,k_neigh,pi+1,high);
	}
}

int partition(double *k_dist, int *k_neigh,int low, int high){

	double temp_dist;
	int temp_point;
	int checker = 0;

	double pivot = k_dist[high];
	//printf("the pivot is %lf \n",pivot);
	int i = (low-1);
	for(int j=low;j<=high-1;j++){
		if(k_dist[j]<=pivot){
			/*swap the two elements*/

			i=i+1;
			temp_dist  = k_dist[i];
			temp_point = k_neigh[i];

			k_dist[i]  = k_dist[j];
			k_neigh[i] = k_neigh[j];

			k_dist[j]  = temp_dist;
			k_neigh[j] = temp_point;
			checker = 1;
		}
	}

	/* swap the pivot */
	i=i+1;
	temp_dist  = k_dist[i];
	temp_point = k_neigh[i];

	k_dist[i]  = k_dist[high];
	k_neigh[i] = k_neigh[high];

	k_dist[high]  = temp_dist;
	k_neigh[high] = temp_point;

	return i;
}

void write_to_file(char *str1, int *nidx, double *ndist, int amount, int k){
  FILE *fp;
  //file for indexes
  fp = fopen(str1,"w+");

  int counter = 1;
  for(int i=0;i<amount;i++){
    fprintf(fp,"%d)%d -- %lf\n",counter,nidx[i],ndist[i]);
    counter ++;
		if((i+1)%k == 0){
      fprintf(fp,"k nearest neighbors for next point \n");
      counter = 1;
    }
  }
	fclose(fp);
}


text:
#pragma once

#include <gsl\gsl>
#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <map>
#include "DrawableGameComponent.h"
#include "FullScreenRenderTarget.h"
#include "GaussianBlur.h"

namespace Rendering
{
	class DiffuseLightingDemo;

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

		std::shared_ptr<DiffuseLightingDemo> DiffuseLighting() const;

		float BlurAmount() const;
		void SetBlurAmount(float blurAmount);

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

	private:
		std::shared_ptr<DiffuseLightingDemo> mDiffuseLightingDemo;
		Library::FullScreenRenderTarget mRenderTarget;
		Library::GaussianBlur mGaussianBlur;
	};
}


text:
#ifndef FUNCTION_ANY_FUNCTION_H
#define FUNCTION_ANY_FUNCTION_H

#include "wasm_base.h"
#include "wasm_value.h"
#include "function/CFunction.h"
#include <gsl/gsl>

namespace wasm {

struct AnyFunction:
	private std::variant<
		std::monostate,
		gsl::not_null<WasmFunction* const*>,
		CFunction
	>
{
	using null_function_alt = std::monostate;
	using wasm_function_alt = gsl::not_null<WasmFunction* const*>;
	using c_function_alt = CFunction;
	using variant_type = std::variant<
		null_function_alt,
		wasm_function_alt,
		c_function_alt
	>;

	using variant_type::variant_type;

	AnyFunction() = default;

	AnyFunction(const AnyFunction&) = default;
	AnyFunction(AnyFunction&&) = default;

	AnyFunction& operator=(const AnyFunction&) = default;
	AnyFunction& operator=(AnyFunction&&) = default;

	void emplace_wasm_function(gsl::not_null<wasm_function_alt> func)
	{ as_variant.emplace<wasm_function_alt>(func); }

	void emplace_c_function(CFunction cfunc)
	{ as_variant().emplace<c_function_alt>(std::move(cfunc)); }

	void emplace_null_function()
	{ as_variant().emplace<null_function_type>(std::monostate{}); }

	bool is_wasm_function() const
	{ return std::holds_alternative<wasm_function_alt>(as_variant()); }

	bool is_c_function() const
	{ return std::holds_alternative<c_function_alt>(as_variant()); }

	bool is_null() const
	{ return std::holds_alternative<std::monostate>(as_variant()); }

	const WasmFunction* get_wasm_function() const
	{
		auto f = std::get<wasm_function_alt>(as_variant());
		return *f;
	}

	WasmFunction* get_wasm_function()
	{ return const_cast<WasmFunction*>(get_wasm_function()); }

	c_function_type get_c_function() const
	{ return std::get<c_function_alt>(as_variant()); }

private:
	variant_type& as_variant()
	{ return static_cast<variant_type&>(*this); }

	const variant_type& as_variant() const
	{ return static_cast<const variant_type&>(*this); }

};

} /* namespace wasm */

#endif /* FUNCTION_ANY_FUNCTION_H */

max_stars_repo_name: tvanslyke/wasm-cpp

max_issues_repo_name: tvanslyke/wasm-cpp

max_forks_repo_name: tvanslyke/wasm-cpp


text:
//This computes the DFT (discrete Fourier transformation) along dim of matrix X.
//This uses a CBLAS matrix multiplication by the DFT matrix.

//This can be useful for smaller transform sizes, especially odd-length ones.
//It is nice that the scaling can be included in the matrix.

//This can also be useful if using only a few of the output frequencies, as given by F.
//Especially useful (with modification) if can be combined with another matrix multiply.

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

#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int dft_cblas_s (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const size_t F, const int sc);
int dft_cblas_d (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const size_t F, const int sc);
int dft_cblas_c (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const size_t F, const int sc);
int dft_cblas_z (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const size_t F, const int sc);

int dft_cblas_s (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const size_t F, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in dft_cblas_s: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndft<Lx) { fprintf(stderr,"error in dft_cblas_s: ndft must be >= Lx (length of vecs in X)\n"); return 1; }
    if (ndft<F) { fprintf(stderr,"error in dft_cblas_s: ndft must be >= F (num output freqs)\n"); return 1; }

    //Scaling
    const float s = sc ? 1.0f/sqrtf((float)(2u*ndft)) : 1.0f;

    if (N==0u || F==0u) {}
    else if (ndft==1u)
    {
        for (size_t n=N; n>0u; --n, ++X) { *Y++ = *X; *Y++ = 0.0f; }
        Y -= 2u*N;
    }
    else
    {
        //Initialize DFT matrix multiply (real, imag parts separate)
        const size_t LF = Lx * F;
        const float P2_N = (float)(2.0*M_PI/(double)ndft);
        float *DFTr, *DFTi;
        DFTr = (float *)aligned_alloc(sizeof(float),LF*sizeof(float));
        DFTi = (float *)aligned_alloc(sizeof(float),LF*sizeof(float));
        if (!DFTr) { fprintf(stderr,"error in dft_cblas_s: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        if (!DFTi) { fprintf(stderr,"error in dft_cblas_s: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t l=0u; l<Lx; ++l, ++DFTr, ++DFTi) { *DFTr = s; *DFTi = 0.0f; }
        for (size_t f=1u; f<F; ++f)
        {
            for (size_t l=0u; l<Lx; ++l)
            {
                *DFTr++ = s * cosf(P2_N*(float)f*(float)l);
                *DFTi++ = -s * sinf(P2_N*(float)f*(float)l);
            }
        }
        DFTr -= LF; DFTi -= LF;

        if (Lx==N)
        {
            //Matrix multiply
            cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0f,DFTr,(int)Lx,X,1,0.0f,Y,2);
            cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0f,DFTi,(int)Lx,X,1,0.0f,&Y[1],2);

            //Enforce real Nyquist
            if (ndft%2u==0u && F>ndft/2u) { Y[ndft+1u] = 0.0f; }
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=0u; v<V; ++v, Y+=2u*F)
                {
                    //Matrix multiply
                    cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0f,DFTr,(int)Lx,X,1,0.0f,Y,2);
                    cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0f,DFTi,(int)Lx,X,1,0.0f,&Y[1],2);

                    //Enforce real Nyquist
                    if (ndft%2u==0u && F>ndft/2u) { Y[ndft+1u] = 0.0f; }
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B*(Lx-1u), Y+=2u*B*(F-1u))
                {
                    for (size_t b=B; b>0u; --b, ++X, Y+=2)
                    {
                        //Matrix multiply
                        cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0f,DFTr,(int)Lx,X,(int)K,0.0f,Y,2*(int)K);
                        cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0f,DFTi,(int)Lx,X,(int)K,0.0f,&Y[1],2*(int)K);

                        //Enforce real Nyquist
                        //if (ndft%2u==0u && F>ndft/2u) { Y[2u*K*(ndft/2u+1u)+1u] = 0.0f; }
                    }
                }
            }
        }
        free(DFTr); free(DFTi);
    }

    return 0;
}

int dft_cblas_d (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const size_t F, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in dft_cblas_d: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndft<Lx) { fprintf(stderr,"error in dft_cblas_d: ndft must be >= Lx (length of vecs in X)\n"); return 1; }
    if (ndft<F) { fprintf(stderr,"error in dft_cblas_d: ndft must be >= F (num output freqs)\n"); return 1; }

    //Scaling
    const double s = sc ? 1.0/sqrt((double)(2u*ndft)) : 1.0;

    if (N==0u || F==0u) {}
    else if (ndft==1u)
    {
        for (size_t n=N; n>0u; --n, ++X) { *Y++ = *X; *Y++ = 0.0; }
        Y -= 2u*N;
    }
    else
    {
        //Initialize DFT matrix multiply
        const size_t LF = Lx * F;
        const double P2_N = 2.0*M_PI/(double)ndft;
        double *DFTr, *DFTi;
        DFTr = (double *)aligned_alloc(sizeof(double),LF*sizeof(double));
        DFTi = (double *)aligned_alloc(sizeof(double),LF*sizeof(double));
        if (!DFTr) { fprintf(stderr,"error in dft_cblas_d: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        if (!DFTi) { fprintf(stderr,"error in dft_cblas_d: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t l=0u; l<Lx; ++l, ++DFTr, ++DFTi) { *DFTr = s; *DFTi = 0.0; }
        for (size_t f=1u; f<F; ++f)
        {
            for (size_t l=0u; l<Lx; ++l)
            {
                *DFTr++ = s * cos(P2_N*(double)f*(double)l);
                *DFTi++ = -s * sin(P2_N*(double)f*(double)l);
            }
        }
        DFTr -= LF; DFTi -= LF;

        if (Lx==N)
        {
            //Matrix multiply
            cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0,DFTr,(int)Lx,X,1,0.0,Y,2);
            cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0,DFTi,(int)Lx,X,1,0.0,&Y[1],2);

            //Enforce real Nyquist
            if (ndft%2u==0u && F>ndft/2u) { Y[ndft+1u] = 0.0f; }
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=0u; v<V; ++v, Y+=2u*F)
                {
                    //Matrix multiply
                    cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0,DFTr,(int)Lx,X,1,0.0,Y,2);
                    cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0,DFTi,(int)Lx,X,1,0.0,&Y[1],2);

                    //Enforce real Nyquist
                    if (ndft%2u==0u && F>ndft/2u) { Y[ndft+1u] = 0.0; }
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B*(Lx-1u), Y+=2u*B*(F-1u))
                {
                    for (size_t b=B; b>0u; --b, ++X, Y+=2)
                    {
                        //Matrix multiply
                        cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0,DFTr,(int)Lx,X,(int)K,0.0,Y,2*(int)K);
                        cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,1.0,DFTi,(int)Lx,X,(int)K,0.0,&Y[1],2*(int)K);

                        //Enforce real Nyquist
                        //if (ndft%2u==0u && F>ndft/2u) { Y[2u*K*(ndft/2u+1u)+1u] = 0.0; }
                    }
                }
            }
        }
        free(DFTr); free(DFTi);
    }

    return 0;
}

int dft_cblas_c (float *Y, const float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const size_t F, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in dft_cblas_c: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndft<Lx) { fprintf(stderr,"error in dft_cblas_c: ndft must be >= Lx (length of vecs in X)\n"); return 1; }
    if (ndft<F) { fprintf(stderr,"error in dft_cblas_c: ndft must be >= F (num output freqs)\n"); return 1; }

    //Scaling
    const float s = sc ? 1.0f/sqrtf((float)(2u*ndft)) : 1.0f;

    if (N==0u || F==0u) {}
    else if (ndft==1u)
    {
        for (size_t n=2u*N; n>0u; --n, ++X, ++Y) { *Y = *X; }
        Y -= 2u*N;
    }
    else
    {
        //Init DFT matrix multiply
        const float z[2] = {0.0f,0.0f}, o[2] = {1.0f,0.0f};
        const size_t LF = Lx * F;
        const float P2_N = (float)(2.0*M_PI/(double)ndft);
        float *DFT;
        DFT = (float *)aligned_alloc(sizeof(float),2u*LF*sizeof(float));
        if (!DFT) { fprintf(stderr,"error in dft_cblas_c: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t f=0u; f<F; ++f)
        {
            for (size_t l=0u; l<Lx; ++l)
            {
                *DFT++ = s * cosf(P2_N*(float)f*(float)l);
                *DFT++ = -s * sinf(P2_N*(float)f*(float)l);
            }
        }
        DFT -= 2u*LF;

        if (Lx==N)
        {
            //Matrix multiply
            cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,o,DFT,(int)Lx,X,1,z,Y,1);
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=0u; v<V; ++v, Y+=2u*F)
                {
                    //Matrix multiply
                    cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,o,DFT,(int)Lx,X,1,z,Y,1);
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=2u*B*(F-1u))
                {
                    for (size_t b=B; b>0u; --b, X+=2, Y+=2)
                    {
                        //Matrix multiply
                        cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,o,DFT,(int)Lx,X,(int)K,z,Y,(int)K);
                    }
                }
            }
        }
        free(DFT);
    }

    return 0;
}

int dft_cblas_z (double *Y, const double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const size_t ndft, const size_t F, const int sc)
{
    if (dim>3u) { fprintf(stderr,"error in dft_cblas_z: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t Lx = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    if (ndft<Lx) { fprintf(stderr,"error in dft_cblas_z: ndft must be >= Lx (length of vecs in X)\n"); return 1; }
    if (ndft<F) { fprintf(stderr,"error in dft_cblas_z: ndft must be >= F (num output freqs)\n"); return 1; }

    //Scaling
    const double s = sc ? 1.0/sqrt((double)(2u*ndft)) : 1.0;

    if (N==0u || F==0u) {}
    else if (ndft==1u)
    {
        for (size_t n=2u*N; n>0u; --n, ++X, ++Y) { *Y = *X; }
        Y -= 2u*N;
    }
    else
    {
        //Init DFT matrix multiply
        const double z[2] = {0.0,0.0}, o[2] = {1.0,0.0};
        const size_t LF = Lx * F;
        const double P2_N = 2.0*M_PI/(double)ndft;
        double *DFT;
        DFT = (double *)aligned_alloc(sizeof(double),2u*LF*sizeof(double));
        if (!DFT) { fprintf(stderr,"error in dft_cblas_z: problem with aligned_alloc. "); perror("aligned_alloc"); return 1; }
        for (size_t f=0u; f<F; ++f)
        {
            for (size_t l=0u; l<Lx; ++l)
            {
                *DFT++ = s * cos(P2_N*(double)f*(double)l);
                *DFT++ = -s * sin(P2_N*(double)f*(double)l);
            }
        }
        DFT -= 2u*LF;

        if (Lx==N)
        {
            //Matrix multiply
            cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,o,DFT,(int)Lx,X,1,z,Y,1);
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/Lx, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=0u; v<V; ++v, Y+=2u*F)
                {
                    //Matrix multiply
                    cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,o,DFT,(int)Lx,X,1,z,Y,1);
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=2u*B*(Lx-1u), Y+=2u*B*(F-1u))
                {
                    for (size_t b=B; b>0u; --b, X+=2, Y+=2)
                    {
                        //Matrix multiply
                        cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)F,(int)Lx,o,DFT,(int)Lx,X,(int)K,z,Y,(int)K);
                    }
                }
            }
        }
        free(DFT);
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


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

/* Author:  G. Jungman */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_coupling.h>

#include "error.h"

inline
static
int locMax3(const int a, const int b, const int c)
{
  int d = GSL_MAX(a, b);
  return GSL_MAX(d, c);
}

inline
static
int locMin3(const int a, const int b, const int c)
{
  int d = GSL_MIN(a, b);
  return GSL_MIN(d, c);
}

inline
static
int locMin5(const int a, const int b, const int c, const int d, const int e)
{
  int f = GSL_MIN(a, b);
  int g = GSL_MIN(c, d);
  int h = GSL_MIN(f, g);
  return GSL_MIN(e, h);
}

/* See: [Thompson, Atlas for Computing Mathematical Functions] */

static
int
delta(int ta, int tb, int tc, gsl_sf_result * d)
{
  gsl_sf_result f1, f2, f3, f4;
  int status = 0;
  status += gsl_sf_fact_e((ta + tb - tc)/2, &f1);
  status += gsl_sf_fact_e((ta + tc - tb)/2, &f2);
  status += gsl_sf_fact_e((tb + tc - ta)/2, &f3);
  status += gsl_sf_fact_e((ta + tb + tc)/2 + 1, &f4);
  if(status != 0) {
    OVERFLOW_ERROR(d);
  }
  d->val = f1.val * f2.val * f3.val / f4.val;
  d->err = 4.0 * GSL_DBL_EPSILON * fabs(d->val);
  return GSL_SUCCESS;
}

static
int
triangle_selection_fails(int two_ja, int two_jb, int two_jc)
{
  return ((two_jb < abs(two_ja - two_jc)) || (two_jb > two_ja + two_jc));
}

static
int
m_selection_fails(int two_ja, int two_jb, int two_jc,
                  int two_ma, int two_mb, int two_mc)
{
  return (
         abs(two_ma) > two_ja
      || abs(two_mb) > two_jb
      || abs(two_mc) > two_jc
      || GSL_IS_ODD(two_ja + two_ma)
      || GSL_IS_ODD(two_jb + two_mb)
      || GSL_IS_ODD(two_jc + two_mc)
      || (two_ma + two_mb + two_mc) != 0
          );
}

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

int
gsl_sf_coupling_3j_e (int two_ja, int two_jb, int two_jc,
                      int two_ma, int two_mb, int two_mc,
                      gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(two_ja < 0 || two_jb < 0 || two_jc < 0) {
    DOMAIN_ERROR(result);
  }
  else if (   triangle_selection_fails(two_ja, two_jb, two_jc)
           || m_selection_fails(two_ja, two_jb, two_jc, two_ma, two_mb, two_mc)
     ) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    int jca  = (-two_ja + two_jb + two_jc) / 2,
        jcb  = ( two_ja - two_jb + two_jc) / 2,
        jcc  = ( two_ja + two_jb - two_jc) / 2,
        jmma = ( two_ja - two_ma) / 2,
        jmmb = ( two_jb - two_mb) / 2,
        jmmc = ( two_jc - two_mc) / 2,
        jpma = ( two_ja + two_ma) / 2,
        jpmb = ( two_jb + two_mb) / 2,
        jpmc = ( two_jc + two_mc) / 2,
        jsum = ( two_ja + two_jb + two_jc) / 2,
        kmin = locMax3 (0, jpmb - jmmc, jmma - jpmc),
        kmax = locMin3 (jcc, jmma, jpmb),
        k, sign = GSL_IS_ODD (kmin - jpma + jmmb) ? -1 : 1,
        status = 0;
    double sum_pos = 0.0, sum_neg = 0.0, norm, term;
    gsl_sf_result bc1, bc2, bc3, bcn1, bcn2, bcd1, bcd2, bcd3, bcd4;

    status += gsl_sf_choose_e (two_ja, jcc , &bcn1);
    status += gsl_sf_choose_e (two_jb, jcc , &bcn2);
    status += gsl_sf_choose_e (jsum+1, jcc , &bcd1);
    status += gsl_sf_choose_e (two_ja, jmma, &bcd2);
    status += gsl_sf_choose_e (two_jb, jmmb, &bcd3);
    status += gsl_sf_choose_e (two_jc, jpmc, &bcd4);

    if (status != 0) {
      OVERFLOW_ERROR (result);
    }

    norm = sqrt (bcn1.val * bcn2.val)
           / sqrt (bcd1.val * bcd2.val * bcd3.val * bcd4.val * ((double) two_jc + 1.0));

    for (k = kmin; k <= kmax; k++) {
      status += gsl_sf_choose_e (jcc, k, &bc1);
      status += gsl_sf_choose_e (jcb, jmma - k, &bc2);
      status += gsl_sf_choose_e (jca, jpmb - k, &bc3);

      if (status != 0) {
        OVERFLOW_ERROR (result);
      }

      term = bc1.val * bc2.val * bc3.val;

      if (sign < 0) {
        sum_neg += norm * term;
      } else {
        sum_pos += norm * term;
      }

      sign = -sign;
    }

    result->val  = sum_pos - sum_neg;
    result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
    result->err += 2.0 * GSL_DBL_EPSILON * (kmax - kmin) * fabs(result->val);

    return GSL_SUCCESS;
  }
}

int
gsl_sf_coupling_6j_INCORRECT_e(int two_ja, int two_jb, int two_jc,
                               int two_jd, int two_je, int two_jf,
                               gsl_sf_result * result)
{
  return gsl_sf_coupling_6j_e(two_ja, two_jb, two_je, two_jd, two_jc, two_jf, result);
}

int
gsl_sf_coupling_6j_e(int two_ja, int two_jb, int two_jc,
                     int two_jd, int two_je, int two_jf,
                     gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(   two_ja < 0 || two_jb < 0 || two_jc < 0
     || two_jd < 0 || two_je < 0 || two_jf < 0
     ) {
    DOMAIN_ERROR(result);
  }
  else if(   triangle_selection_fails(two_ja, two_jb, two_jc)
          || triangle_selection_fails(two_ja, two_je, two_jf)
          || triangle_selection_fails(two_jb, two_jd, two_jf)
          || triangle_selection_fails(two_je, two_jd, two_jc)
     ) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result n1;
    gsl_sf_result d1, d2, d3, d4, d5, d6;
    double norm;
    int tk, tkmin, tkmax;
    double phase;
    double sum_pos = 0.0;
    double sum_neg = 0.0;
    double sumsq_err = 0.0;
    int status = 0;
    status += delta(two_ja, two_jb, two_jc, &d1);
    status += delta(two_ja, two_je, two_jf, &d2);
    status += delta(two_jb, two_jd, two_jf, &d3);
    status += delta(two_je, two_jd, two_jc, &d4);
    if(status != GSL_SUCCESS) {
      OVERFLOW_ERROR(result);
    }
    norm = sqrt(d1.val) * sqrt(d2.val) * sqrt(d3.val) * sqrt(d4.val);

    tkmin = locMax3(0,
                   two_ja + two_jd - two_jc - two_jf,
                   two_jb + two_je - two_jc - two_jf);

    tkmax = locMin5(two_ja + two_jb + two_je + two_jd + 2,
                    two_ja + two_jb - two_jc,
                    two_je + two_jd - two_jc,
                    two_ja + two_je - two_jf,
                    two_jb + two_jd - two_jf);

    phase = GSL_IS_ODD((two_ja + two_jb + two_je + two_jd + tkmin)/2)
            ? -1.0
            :  1.0;

    for(tk=tkmin; tk<=tkmax; tk += 2) {
      double term;
      double term_err;
      gsl_sf_result den_1, den_2;
      gsl_sf_result d1_a, d1_b;
      status = 0;

      status += gsl_sf_fact_e((two_ja + two_jb + two_je + two_jd - tk)/2 + 1, &n1);
      status += gsl_sf_fact_e(tk/2, &d1_a);
      status += gsl_sf_fact_e((two_jc + two_jf - two_ja - two_jd + tk)/2, &d1_b);
      status += gsl_sf_fact_e((two_jc + two_jf - two_jb - two_je + tk)/2, &d2);
      status += gsl_sf_fact_e((two_ja + two_jb - two_jc - tk)/2, &d3);
      status += gsl_sf_fact_e((two_je + two_jd - two_jc - tk)/2, &d4);
      status += gsl_sf_fact_e((two_ja + two_je - two_jf - tk)/2, &d5);
      status += gsl_sf_fact_e((two_jb + two_jd - two_jf - tk)/2, &d6);

      if(status != GSL_SUCCESS) {
        OVERFLOW_ERROR(result);
      }

      d1.val = d1_a.val * d1_b.val;
      d1.err = d1_a.err * fabs(d1_b.val) + fabs(d1_a.val) * d1_b.err;

      den_1.val  = d1.val*d2.val*d3.val;
      den_1.err  = d1.err * fabs(d2.val*d3.val);
      den_1.err += d2.err * fabs(d1.val*d3.val);
      den_1.err += d3.err * fabs(d1.val*d2.val);

      den_2.val  = d4.val*d5.val*d6.val;
      den_2.err  = d4.err * fabs(d5.val*d6.val);
      den_2.err += d5.err * fabs(d4.val*d6.val);
      den_2.err += d6.err * fabs(d4.val*d5.val);

      term  = phase * n1.val / den_1.val / den_2.val;
      phase = -phase;
      term_err  = n1.err / fabs(den_1.val) / fabs(den_2.val);
      term_err += fabs(term / den_1.val) * den_1.err;
      term_err += fabs(term / den_2.val) * den_2.err;

      if(term >= 0.0) {
        sum_pos += norm*term;
      }
      else {
        sum_neg -= norm*term;
      }

      sumsq_err += norm*norm * term_err*term_err;
    }

    result->val  = sum_pos - sum_neg;
    result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
    result->err += sqrt(sumsq_err / (0.5*(tkmax-tkmin)+1.0));
    result->err += 2.0 * GSL_DBL_EPSILON * (tkmax - tkmin + 2.0) * fabs(result->val);

    return GSL_SUCCESS;
  }
}

int
gsl_sf_coupling_RacahW_e(int two_ja, int two_jb, int two_jc,
                         int two_jd, int two_je, int two_jf,
                         gsl_sf_result * result)
{
  int status = gsl_sf_coupling_6j_e(two_ja, two_jb, two_je, two_jd, two_jc, two_jf, result);
  int phase_sum = (two_ja + two_jb + two_jc + two_jd)/2;
  result->val *= ( GSL_IS_ODD(phase_sum) ? -1.0 : 1.0 );
  return status;
}

int
gsl_sf_coupling_9j_e(int two_ja, int two_jb, int two_jc,
                     int two_jd, int two_je, int two_jf,
                     int two_jg, int two_jh, int two_ji,
                     gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(   two_ja < 0 || two_jb < 0 || two_jc < 0
     || two_jd < 0 || two_je < 0 || two_jf < 0
     || two_jg < 0 || two_jh < 0 || two_ji < 0
     ) {
    DOMAIN_ERROR(result);
  }
  else if(   triangle_selection_fails(two_ja, two_jb, two_jc)
          || triangle_selection_fails(two_jd, two_je, two_jf)
          || triangle_selection_fails(two_jg, two_jh, two_ji)
          || triangle_selection_fails(two_ja, two_jd, two_jg)
          || triangle_selection_fails(two_jb, two_je, two_jh)
          || triangle_selection_fails(two_jc, two_jf, two_ji)
     ) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    int tk;
    int tkmin = locMax3(abs(two_ja-two_ji), abs(two_jh-two_jd), abs(two_jb-two_jf));
    int tkmax = locMin3(two_ja + two_ji, two_jh + two_jd, two_jb + two_jf);
    double sum_pos = 0.0;
    double sum_neg = 0.0;
    double sumsq_err = 0.0;
    double phase;
    for(tk=tkmin; tk<=tkmax; tk += 2) {
      gsl_sf_result s1, s2, s3;
      double term;
      double term_err;
      int status = 0;

      status += gsl_sf_coupling_6j_e(two_ja, two_ji, tk,  two_jh, two_jd, two_jg,  &s1);
      status += gsl_sf_coupling_6j_e(two_jb, two_jf, tk,  two_jd, two_jh, two_je,  &s2);
      status += gsl_sf_coupling_6j_e(two_ja, two_ji, tk,  two_jf, two_jb, two_jc,  &s3);

      if(status != GSL_SUCCESS) {
        OVERFLOW_ERROR(result);
      }
      term = s1.val * s2.val * s3.val;
      term_err  = s1.err * fabs(s2.val*s3.val);
      term_err += s2.err * fabs(s1.val*s3.val);
      term_err += s3.err * fabs(s1.val*s2.val);

      if(term >= 0.0) {
        sum_pos += (tk + 1) * term;
      }
      else {
        sum_neg -= (tk + 1) * term;
      }

      sumsq_err += ((tk+1) * term_err) * ((tk+1) * term_err);
    }

    phase = GSL_IS_ODD(tkmin) ? -1.0 : 1.0;

    result->val  = phase * (sum_pos - sum_neg);
    result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
    result->err += sqrt(sumsq_err / (0.5*(tkmax-tkmin)+1.0));
    result->err += 2.0 * GSL_DBL_EPSILON * (tkmax-tkmin + 2.0) * fabs(result->val);

    return GSL_SUCCESS;
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_coupling_3j(int two_ja, int two_jb, int two_jc,
                          int two_ma, int two_mb, int two_mc)
{
  EVAL_RESULT(gsl_sf_coupling_3j_e(two_ja, two_jb, two_jc,
                                   two_ma, two_mb, two_mc,
                                   &result));
}

double gsl_sf_coupling_6j_INCORRECT(int two_ja, int two_jb, int two_jc,
                                    int two_jd, int two_je, int two_jf)
{
  EVAL_RESULT(gsl_sf_coupling_6j_INCORRECT_e(two_ja, two_jb, two_jc,
                                             two_jd, two_je, two_jf,
                                             &result));
}

double gsl_sf_coupling_6j(int two_ja, int two_jb, int two_jc,
                          int two_jd, int two_je, int two_jf)
{
  EVAL_RESULT(gsl_sf_coupling_6j_e(two_ja, two_jb, two_jc,
                                   two_jd, two_je, two_jf,
                                   &result));
}

double gsl_sf_coupling_RacahW(int two_ja, int two_jb, int two_jc,
                          int two_jd, int two_je, int two_jf)
{
  EVAL_RESULT(gsl_sf_coupling_RacahW_e(two_ja, two_jb, two_jc,
                                      two_jd, two_je, two_jf,
                                      &result));
}

double gsl_sf_coupling_9j(int two_ja, int two_jb, int two_jc,
                          int two_jd, int two_je, int two_jf,
                          int two_jg, int two_jh, int two_ji)
{
  EVAL_RESULT(gsl_sf_coupling_9j_e(two_ja, two_jb, two_jc,
                                   two_jd, two_je, two_jf,
                                   two_jg, two_jh, two_ji,
                                   &result));
}


text:
#pragma once

#include <chrono>
#include <unordered_map>
#include <vector>

#include <gsl/span>
#include <json.hpp>

#include "Quiver/Animation/AnimationId.h"
#include "Quiver/Animation/AnimationLibrary.h"
#include "Quiver/Animation/AnimatorId.h"
#include "Quiver/Animation/Rect.h"
#include "Quiver/Graphics/ViewBuffer.h"

namespace qvr {

struct AnimatorTarget {
	ViewBuffer views;

	AnimatorTarget() = default;

	AnimatorTarget(const AnimatorTarget&) = delete;
	AnimatorTarget(const AnimatorTarget&&) = delete;

	AnimatorTarget& operator=(AnimatorTarget&) = delete;
	AnimatorTarget& operator=(AnimatorTarget&&) = delete;
};

class AnimatorRepeatSetting
{
public:
	explicit AnimatorRepeatSetting(const int repeatCount) : m_RepeatCount(repeatCount) {}

	AnimatorRepeatSetting() = default;

	int GetRepeatCount() const { return m_RepeatCount; }

	static const AnimatorRepeatSetting Forever;
	static const AnimatorRepeatSetting Never;
	static const AnimatorRepeatSetting Once;
	static const AnimatorRepeatSetting Twice;

private:
	int m_RepeatCount = Forever.GetRepeatCount();
};

struct AnimatorStartSetting
{
	AnimatorStartSetting(AnimationId animationId)
		: m_AnimationId(animationId)
	{}

	AnimatorStartSetting(
		AnimationId animationId,
		AnimatorRepeatSetting repeatSetting)
		: m_AnimationId(animationId)
		, m_RepeatSetting(repeatSetting)
	{}

	AnimationId m_AnimationId;
	AnimatorRepeatSetting m_RepeatSetting;
};

inline bool operator==(const AnimatorRepeatSetting& a, const AnimatorRepeatSetting& b) {
	return a.GetRepeatCount() == b.GetRepeatCount();
}

inline bool operator!=(const AnimatorRepeatSetting& a, const AnimatorRepeatSetting& b) {
	return a.GetRepeatCount() != b.GetRepeatCount();
}

struct AnimationLibraryEditorData;

class AnimationData;

class AnimatorCollection {
public:
	auto GetAnimations() const -> const AnimationLibrary& {
		return animations;
	}

	AnimationId AddAnimation(const AnimationData& anim);
	AnimationId AddAnimation(
		const AnimationData& anim,
		const AnimationSourceInfo& animSourceInfo);

	bool RemoveAnimation(const AnimationId id);

	int GetReferenceCount(const AnimationId animation) const {
		if (animationReferenceCounts.count(animation) == 0) {
			return 0;
		}

		return animationReferenceCounts.at(animation);
	}

	AnimatorId Add(
		AnimatorTarget& target,
		const AnimatorStartSetting& startSetting);

	bool Remove(const AnimatorId id);

	bool Exists(const AnimatorId id) const {
		return animators.states.count(id) != 0;
	}

	int GetCount() const {
		return animators.states.size();
	}

	void AnimatorGui(const AnimatorId id);

	bool SetAnimation(
		const AnimatorId animatorId,
		const AnimatorStartSetting& animation,
		const bool clearQueue = true);

	bool SetTarget(
		const AnimatorId id,
		AnimatorTarget& newTarget);

	bool SetFrame(
		const AnimatorId id,
		const int index);

	bool QueueAnimation(
		const AnimatorId animatorId,
		const AnimatorStartSetting& pendingAnimation);

	bool ClearAnimationQueue(const AnimatorId id);

	unsigned GetFrame(const AnimatorId animatorId) const;

	AnimationId GetAnimation(const AnimatorId animatorId) const;

	void Animate(const Animation::TimeUnit ms);

private:
	struct AnimatorState {
		AnimatorState(
			const AnimationId animation,
			const int frameIndex,
			const int index,
			AnimatorTarget& target,
			const AnimatorRepeatSetting& repeat)
			: index(index)
			, currentFrame(frameIndex)
			, currentAnimation(animation)
			, target(&target)
			, repeatSetting(repeat)
			, repeatCount(0)
		{}

		AnimatorState() = default;

		int index;
		int currentFrame;
		int repeatCount;
		AnimationId currentAnimation;
		AnimatorRepeatSetting repeatSetting;
		AnimatorTarget* target;
		std::vector<AnimatorStartSetting> queuedAnimations;
	};

	struct AnimatorState_Hot {
		AnimatorState_Hot(
			const AnimatorId animatorId,
			const Animation::TimeUnit timeLeft)
			: animatorId(animatorId)
			, timeLeftInFrame(timeLeft)
		{}

		Animation::TimeUnit timeLeftInFrame;
		AnimatorId animatorId;
	};

	struct Animators {
		std::vector<AnimatorState_Hot> hotStates;
		std::unordered_map<AnimatorId, AnimatorState> states;
		AnimatorId GetNextAnimatorId();
	private:
		AnimatorId lastId = AnimatorId::Invalid;
	} animators;

	AnimationLibrary animations;

	std::unordered_map<AnimationId, unsigned> animationReferenceCounts;

	// friends:

	friend void GuiControls(
		AnimatorCollection& animators,
		AnimationLibraryEditorData& editorData);

	friend void to_json(nlohmann::json& j, const AnimatorCollection& animators);
	friend void from_json(const nlohmann::json& j, AnimatorCollection& animators);
};

void GuiControls(AnimatorCollection& animators, AnimationLibraryEditorData& editorData);

void to_json(nlohmann::json& j, const AnimatorCollection& animators);
void from_json(const nlohmann::json& j, AnimatorCollection& animators);

}

max_stars_repo_name: rachelnertia/Quarrel

max_issues_repo_name: rachelnertia/Quarrel

max_forks_repo_name: rachelnertia/Quiver


text:
#include <Exceptions.h>
#include <Types.h>

#ifndef GENERIC_API_NAMESPACE
	#error "Wrong usage of this header"
#endif

#ifndef GENERIC_API_ROUTINES_NAMESPACE
	#error "Wrong usage of this header"
#endif

#define ROUTINES_NAMESPACE namespace GENERIC_API_ROUTINES_NAMESPACE
#define BLAS_NAMESPACE namespace GENERIC_API_NAMESPACE

#ifndef GENERIC_API_DEFINE

namespace cl
{
	namespace routines
	{
		ROUTINES_NAMESPACE
		{
			template<MathDomain md>
			void Copy(MemoryBuffer&, const MemoryBuffer&)
			{
				throw NotImplementedException();
			}
		}
	}	 // namespace routines
}	 // namespace cl

#else

BLAS_NAMESPACE
{
	#ifdef CBLAS_H
		#undef CBLAS_H
	#endif

	#include <cblas.h>
}

namespace cl
{
	namespace routines
	{
		ROUTINES_NAMESPACE
		{
			template<MathDomain md>
			void Copy(MemoryBuffer & dest, const MemoryBuffer& source);

			template<>
			inline void Copy<MathDomain::Float>(MemoryBuffer & dest, const MemoryBuffer& source)
			{
				GENERIC_API_NAMESPACE::cblas_scopy(static_cast<int>(dest.size), reinterpret_cast<const float*>(source.pointer), 1, reinterpret_cast<float*>(dest.pointer), 1);
			}
			template<>
			inline void Copy<MathDomain::Double>(MemoryBuffer & dest, const MemoryBuffer& source)
			{
				GENERIC_API_NAMESPACE::cblas_dcopy(static_cast<int>(dest.size), reinterpret_cast<const double*>(source.pointer), 1, reinterpret_cast<double*>(dest.pointer), 1);
			}
		}
	}	 // namespace routines
}	 // namespace cl

#endif

#undef GENERIC_API_DEFINE
#undef GENERIC_API_NAMESPACE
#undef GENERIC_API_ROUTINES_NAMESPACE
#undef ROUTINES_NAMESPACE
#undef BLAS_NAMESPACE

max_stars_repo_name: pmontalb/CudaLight

max_issues_repo_name: pmontalb/CudaLight

max_forks_repo_name: pmontalb/CudaLight


text:
/* stable/stable_koutrouvelis.h
 *
 * Koutrouvelis method for parameter estimation of alpha-stable
 * distributions. Based on [1] and MATLAB code available in [2].
 *
 * [1] Koutrouvelis, I. A. An Iterative Procedure for the Estimation
 *     of the Parameters of Stable Laws Communications in Statistics
 *     - Simulation and Computation, 1981, 10, 17-28
 * [2] Mark S. Veillete. Alpha-Stable Distributions in MATLAB
 *     http://math.bu.edu/people/mveillet/research.html
 *
 * Copyright (C) 2013. Javier Royuela del Val
 *                     Federico Simmross Wattenberg
 *
 * 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; version 3 of the License.
 *
 * 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, see <http://www.gnu.org/licenses/>.
 *
 *
 *  Javier Royuela del Val.
 *  E.T.S.I. Telecomunicación
 *  Universidad de Valladolid
 *  Paseo de Belén 15, 47002 Valladolid, Spain.
 *  jroyval@lpi.tel.uva.es
 */
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit.h>
#include "stable.h"

inline double sign(double x) {
  return ((.0 < x) - (x < .0));
}

double var(const double * data, const int N) {
  double acum = .0, acum2 = .0;
  double v;

  int i;
  for(i=0;i<N;i++) {
    acum  += data[i];
    acum2 += data[i]*data[i];
  }

  v = (1.0/(N-1.0)) * (acum2 - (1.0/N)*acum*acum);
  return v;
}

int chooseK(double,int);

int chooseL(double,int);

void setcovYY(const double * t, int K, int N, double alpha, double beta, double gam, double **covYY);

void setcovZZ(const double * t, int K, int N, double alpha, double beta, double gam, double **covZZ);

double ecfRoot(const double * data, const int N);

void stable_samplecharfunc(const double* x, const unsigned int Nx,
             const double* t, const unsigned int Nt, gsl_complex *z);

gsl_complex stable_samplecharfunc_point(const double* x,
            const unsigned int N, double t);

int stable_fit_koutrouvelis(StableDist * dist, const double * data, const unsigned int N) {

  int maxiter = 10;
  double xTol = 0.01;

  double * s = NULL;
  gsl_complex * phi = NULL;
  gsl_complex * phi2 = NULL;

  int i;

  double alpha = dist->alpha;
  double beta  = dist->beta;
  double sigma = dist->sigma;
  double mu1   = dist->mu_1;

  double alphaold = alpha;
  double mu1old   = mu1;

  double diff     = .0;
  double diffbest = DBL_MAX;
  double diffmax  = .0;

  double alphabest = alpha;
  double betabest  = beta;
  double sigmabest = sigma;
  double mu1best   = mu1;

  double * t = NULL;
  double * w = NULL;
  double * y = NULL;
  double * p = NULL;

  double * t2 = NULL;
  double * w2 = NULL;
  double * y2 = NULL;
  double * p2 = NULL;

  gsl_matrix * X;
  gsl_matrix * covmat;
  gsl_vector * cvout;
  gsl_vector * ydat;
  gsl_vector * weights;

  int stat;
  double sumsq = .0;
  double c0 = .0, c1 = .0;
  double cov00 = .0, cov01 = .0, cov11 = .0;
  double sigmanew = 0;

  double **covYY;// **covZZ;
  double **covZZ;

  int iter = 0;
  int K,L;
  int row;
  double t_0=0;
  double step_u=0;
  double estshift;

  gsl_multifit_linear_workspace * linws;

  if (sigma == 0) {
    sigma = sqrt(var(data,N));
	stable_setparams(dist,alpha,beta,sigma,mu1,1);
  }

  s = (double *) malloc(N*sizeof(double));
  for (i=0;i<N;i++) {
    s[i] = (data[i]-mu1)/sigma;
  }

  covmat = gsl_matrix_alloc(2,2);
  cvout = gsl_vector_alloc(2);

  for (iter = 0; iter<maxiter; iter ++) {
    if (iter <= 1) {
      K = chooseK(alpha,N);

      t = (double *) realloc(t,K*sizeof(double));
      p = (double *) realloc(p,K*sizeof(double));
      w = (double *) realloc(w,K*sizeof(double));
      y = (double *) realloc(y,K*sizeof(double));

      phi = (gsl_complex*) realloc(phi,K*sizeof(gsl_complex));

      for(i=0;i<K;i++) {
        t[i]=((double)i+1.0)*M_PI/25.0;
        w[i]=log(fabs(t[i]));
      }
    }

    if (iter==1) {
      covYY = (double **)malloc(K*sizeof(double *));
      covYY[0] = (double*)malloc(K*K*sizeof(double));
      for (row=1;row<K;row++) {
        covYY[row] = covYY[0] + row*K;
      }
    }

    stable_samplecharfunc(s,N,t,K,phi);

    for(i=0;i<K;i++) {
      y[i] = log(-2.0*gsl_complex_logabs(phi[i]));
    }

    if (iter==0) { //Use ordinary least squares regression
 //     printf("ordls: ");fflush(stdout);
      stat = gsl_fit_linear (w, 1, y, 1, K, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
      alpha = c1;
      sigmanew = pow(exp(c0)/2.0,1.0/alpha);
      sigma = sigma*sigmanew;
//      printf("stat = %d, alpha=%f, sigma=%f\n",stat,alpha,sigma);
    }
    else { //Use weighted least squares regression
//      printf("wls: ");fflush(stdout);
      setcovYY(t, K, N, alpha, beta, 1.0, covYY);
      for (i=0;i<K;i++) {
        p[i] = 1.0/covYY[i][i];
      }
      stat = gsl_fit_wlinear(w, 1, p, 1, y, 1, K, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
      alpha = c1;
      sigmanew = pow(exp(c0)/2.0,1.0/alpha);
      sigma = sigma*sigmanew;
 //     printf("stat = %d, alpha=%f, sigma=%f\n",stat,alpha,sigma);
    }

    /*************** rescale data ******************/
//    printf("rescale\n");fflush(stdout);
    for (i=0;i<N;i++) {
      s[i] = s[i]/sigmanew;
    }
    if (alpha<0) alpha = 0;
    if (alpha>2) alpha = 2;
    if (beta<-1) beta = -1;
    if (beta> 1) beta = 1;
    if (sigma<0) sigma = 0;

    /********** refine beta and mu **************/
 //   printf("refbetamu\n");fflush(stdout);
    if (iter <= 1) {
      L = chooseL(alpha,N);
//      printf("alpha = %f, N = %d, L = %d\n",alpha,N,L);fflush(stdout);
      t_0 = ecfRoot(s,N);
      step_u = M_PI/50;
      if (t_0/L<step_u) step_u = t_0/L;

      t2 = (double *) realloc(t2,L*sizeof(double));
      p2 = (double *) realloc(p2,L*sizeof(double));
      w2 = (double *) realloc(w2,L*sizeof(double));
      y2 = (double *) realloc(y2,L*sizeof(double));

      phi2 = (gsl_complex*)realloc(phi2,L*sizeof(gsl_complex));

      for(i=0;i<L;i++) {
        t2[i]=((double)i+1.0)*step_u;
        w2[i]=sign(t2[i])*pow(fabs(t2[i]),alpha);
      }
    }
    if (iter==0){
//      printf("multimin alloc L = %d\n",L);fflush(stdout);
      linws = gsl_multifit_linear_alloc (L,2);
    }
    else if (iter ==1) {
      //Reserve covariance matrix mem
      covZZ = (double **)malloc(L*sizeof(double *));
      covZZ[0] = (double*)malloc(L*L*sizeof(double));
      for (row=0;row<L;row++) {
        covZZ[row] = covZZ[0] + row*L;
      }
      gsl_multifit_linear_free(linws);
      linws = gsl_multifit_linear_alloc (L,2);
    }

 //   printf("Second charfunc iter %d\n",iter);fflush(stdout);
    stable_samplecharfunc(s,N,t2,L,phi2);
 //   printf("done\n");

    for(i=0;i<L;i++) {
      y2[i] = gsl_complex_arg(phi2[i]);
    }

    X    = gsl_matrix_alloc(L,2);
    ydat = gsl_vector_alloc(L);

    for(row=0;row<L;row++) {
      gsl_matrix_set(X,row,0,t2[row]);
      gsl_matrix_set(X,row,1,w2[row]);
      gsl_vector_set(ydat,row,y2[row]);
    }

    if (iter==0) {
 //     printf("multi ordls: ");
      stat   = gsl_multifit_linear(X,ydat,cvout,covmat,&sumsq,linws);
 //     printf("stat = %d, beta=%f\n",stat,beta);
    }
    else {
//      printf("multi wls: ");
      setcovZZ(t2, L, N, alpha,  beta,  1.0, covZZ);
      weights = gsl_vector_alloc(L);
      for (i=0;i<L;i++) {
        p2[i] = 1.0/covZZ[i][i];
        gsl_vector_set(weights,i,p2[i]);
      }
      stat  = gsl_multifit_wlinear(X,weights,ydat,cvout,covmat,&sumsq,linws);
//      printf("stat = %d, beta=%f\n",stat,beta);
    }
    if (alpha>1.98 || alpha < 0.05 || fabs(alpha-1.0)<0.05 || isnan(alpha)) {
        beta = 0.0;
    } else {
        beta = gsl_vector_get(cvout,1)/tan(alpha*M_PI*0.5);
    }
    estshift = gsl_vector_get(cvout,0);
    mu1    = mu1 + sigma * estshift;

//    printf("rem shift\n");
    /*** remove estimated shift ***/
    for (i=0;i<N;i++) {
      s[i] = s[i] - estshift;
    }

    if (isnan(alpha) || isnan(beta) || isnan(sigma) || isnan(mu1))
    {
      iter++;
      break;
    }

//    printf("check conv");
    /** check for convergence or blow up**/
    diff = pow(alpha - alphaold,2) + pow(mu1 - mu1old,2);
    if (iter <= 2 && diff > diffmax) {
      diffmax = diff;
    }
//    else if (diff > 2.0*diffmax) {
#ifdef DEBUG
      printf("blow up on iter %d\n",iter);
#endif
//      break;
//    }

//    printf("check best\n");
    if (fabs(diff) < diffbest) {
      alphabest = alpha;
      betabest  = beta;
      sigmabest = sigma;
      mu1best   = mu1;

      diffbest = diff;

      if (diff < xTol) {
        gsl_matrix_free(X);
        gsl_vector_free(ydat);
        if (iter>0) gsl_vector_free(weights);
        iter++;
        break;
      }
    }

    alphaold = alpha;
    mu1old   = mu1;

    gsl_matrix_free(X);
    gsl_vector_free(ydat);
    if (iter>0) {
      gsl_vector_free(weights);
    }
  }

  if (maxiter > 0 && iter >= 0) {
      alpha = alphabest;
      beta = betabest;
      sigma = sigmabest;
      mu1   = mu1best;
  }

  if (alpha<0) alpha = 0;
  if (alpha>2) alpha = 2;
  if (beta<-1) beta = -1;
  if (beta> 1) beta = 1;
  if (sigma<0) sigma = 0;

  stable_setparams(dist,alpha,beta,sigma,mu1,1);

  free(t);
  free(p);
  free(w);
  free(y);
  free(phi);
  free(t2);
  free(p2);
  free(w2);
  free(y2);
  free(phi2);
  gsl_matrix_free(covmat);
  gsl_vector_free(cvout);
  free(s);

  if (iter>1) {
    free(covYY[0]);
    free(covYY);
    free(covZZ[0]);
    free(covZZ);
  }

  gsl_multifit_linear_free(linws);

  //printf(" iter %d diff %f a %f b %f s %f m %f\n",iter,diff,alpha,beta,sigma,dist->mu_0);

  return 0;
}  // end stable_fit_koutrouvelis

int chooseK(double alpha, int N) {
  double a[]={1.9,1.5,1.3,1.1,.9,.7,.5,.3};
  int n[] = {200, 800, 1600};

  double Kmat[8][3]={{ 9, 9, 9},
                     {11,11,11},
                     {22,16,14},
                     {24,18,15},
                     {28,22,18},
                     {30,24,20},
                     {86,68,56},
                     {134,124,118}};

  int i,j;
  double xi;
  double xj;
  double Kp;

  if (alpha < 0.3) alpha = 0.3;
  if (alpha > 1.9) alpha = 1.9;
  if (N<200) N=200;
  if (N>1600) N=1600;

  i=1;
  j=1;
  while (i<7 && a[i]>=alpha) {
    ++i;
  }
  while (j<2 && n[j]<=N) {
    ++j;
  }

  xi = 1.0 - (alpha-a[i])/(a[i-1]-a[i]);
  xj = 1.0 - (double)(n[j]-N)/(double)(n[j]-n[j-1]);

//  printf("i,j = %d, %d; xi,xj = %f,  %f\n",i,j,xi,xj);

  /* bilinear interpolation */
  Kp =       xj * (xi * Kmat[i][ j ] + (1.0-xi) * Kmat[i-1][ j ] ) +
       (1.0-xj) * (xi * Kmat[i][j-1] + (1.0-xi) * Kmat[i-1][j-1] );

  Kp = floor(Kp+.5);

  return (int)Kp;
}

int chooseL(double alpha, int N) {
  double a[]={1.9,1.5,1.1,.9,.7,.5,.3};
  int n[] = {200, 800, 1600};

  double Lmat[7][3]={{ 9,10,11},
                     {12,14,15},
                     {16,18,17},
                     {14,14,14},
                     {24,16,16},
                     {40,38,36},
                     {70,68,66}};

  int i,j;
  double xi, xj;
  double Lp;

  if (alpha < 0.3) alpha = 0.3;
  if (alpha > 1.9) alpha = 1.9;
  if (N<200) N=200;
  if (N>1600) N=1600;

  i=1;
  j=1;
  while (i<6 && a[i]>=alpha) {
    ++i;
  }
  while (j<2 && n[j]<=N) {
    ++j;
  }

  xi = 1.0 - (alpha-a[i])/(a[i-1]-a[i]);
  xj = 1.0 - (double)(n[j]-N)/(double)(n[j]-n[j-1]);

//  printf("i,j = %d, %d; xi,xj = %f,  %f\n",i,j,xi,xj);
  /* bilinear interpolation */
  Lp =       xj * (xi * Lmat[i][ j ] + (1.0-xi) * Lmat[i-1][ j ] ) +
       (1.0-xj) * (xi * Lmat[i][j-1] + (1.0-xi) * Lmat[i-1][j-1] );

  Lp = floor(Lp + .5);

  return (int)Lp;
}

void setcovYY(const double * t, int K, int N, double alpha, double beta, double gam, double **covYY) {
  double w = tan(alpha*M_PI_2);
  double calpha = pow(gam,alpha);
  int    j = 0, k = 0;

  double * talpha = (double *)malloc(K*sizeof(double));

  for (j=0;j<K;j++) {
    talpha[j]=pow(fabs(t[j]),alpha);
  }

  double tj,tja,tk,tka;
  double tjmtka,tjptka,stj,stk;
  double A,B,D,E;//F,G,H;
  for (j=0;j<K;j++) {
    tj = t[j];
    tja = talpha[j];
    stj = sign(tj);
    for (k=0;k<K;k++) {
      tk  = t[k];
      tka = talpha[k];
      stk = sign(tk);
      tjmtka = pow(fabs(tj-tk),alpha);
      tjptka = pow(fabs(tj+tk),alpha);

      A =  calpha * (tja + tka - tjmtka);
      B =  calpha * beta * (-tja * stj * w + tka * stk * w + tjmtka * sign(tj - tk) * w );
      D =  calpha * (tja + tka - tjptka);
      E =  calpha * beta * ( tja * stj * w + tja * stk * w - tjptka * sign(tj + tk) * w );
//      F =  calpha * (tja + tka);
//      G = -calpha * tjmtka;
//      H = -calpha * tjptka;

//      printf("[j][k] = [%d][%d]",j,k);fflush(stdout);
      covYY[j][k] = (exp(A)*cos(B)+exp(D)*cos(E)-2.0) / ( 2.0 * N * pow(gam,2.0*alpha)*pow(fabs(tj*tk),alpha) );
//      covZZ[j][k] =  exp(F)*( exp(G)*cos(B)-exp(H)*cos(E) ) /(2.0*N);
    }
  }

  free(talpha);

  return;
}
void setcovZZ(const double * t, int K, int N, double alpha, double beta, double gam, double **covZZ) {
  double w = tan(alpha*M_PI_2);
  double calpha = pow(gam,alpha);
  int    j = 0, k = 0;

  double * talpha = (double *)malloc(K*sizeof(double));

  for (j=0;j<K;j++) {
    talpha[j]=pow(fabs(t[j]),alpha);
  }

  double tj,tja,tk,tka;
  double tjmtka,tjptka,stj,stk;
  double B,E,F,G,H;
  for (j=0;j<K;j++) {
    tj = t[j];
    tja = talpha[j];
    stj = sign(tj);
    for (k=0;k<K;k++) {
      tk  = t[k];
      tka = talpha[k];
      stk = sign(tk);
      tjmtka = pow(fabs(tj-tk),alpha);
      tjptka = pow(fabs(tj+tk),alpha);

//      A =  calpha * (tja + tka - tjmtka);
      B =  calpha * beta * (-tja * stj * w + tka * stk * w + tjmtka * sign(tj - tk) * w );
//      D =  calpha * (tja + tka - tjptka);
      E =  calpha * beta * ( tja * stj * w + tja * stk * w - tjptka * sign(tj + tk) * w );
      F =  calpha * (tja + tka);
      G = -calpha * tjmtka;
      H = -calpha * tjptka;

//      covYY[j][k] = (exp(A)*cos(B)+exp(D)*cos(E)-2.0) / ( 2.0 * N * pow(gam,2.0*alpha)*pow(fabs(tj*tk),alpha) );
      covZZ[j][k] =  exp(F)*( exp(G)*cos(B)-exp(H)*cos(E) ) /(2.0*N);
    }
  }

  free(talpha);

  return;
}

double ecfRoot(const double * data, int N) {

  double m=0;
  int i;
  for(i=0;i<N;i++) {
    m += fabs(data[i]);
  }
  m = m/N;

  double t = 0;
  gsl_complex val = stable_samplecharfunc_point(data, N, t);
  int iter=0;
  while (iter<10000 && fabs(GSL_REAL(val))>1e-3) {
    t += GSL_REAL(val)/m;
    val = stable_samplecharfunc_point(data, N, t);
  }

  return t;
}

gsl_complex stable_samplecharfunc_point(const double* x,
            const unsigned int N, double t)
{
  unsigned int i;
  double zr=0.0,zi=0.0;
  gsl_complex z;

  for(i=0;i<N;i++)
   {
      zr+=cos(t*x[i]);
      zi+=sin(t*x[i]);
    }

  GSL_SET_COMPLEX(&z,zr/N,zi/N);

  return z;
}

void stable_samplecharfunc(const double* x, const unsigned int Nx,
             const double* t, const unsigned int Nt, gsl_complex *z)
{
  unsigned int ix,it;
  double zr=0.0,zi=0.0,w;

  for(it=0;it<Nt;it++)
   {
     w = t[it];
     zr = .0;
     zi = .0;
     for(ix=0;ix<Nx;ix++)
       {
         zr+=cos(w*x[ix]);
         zi+=sin(w*x[ix]);
       }
//     if (isnan(zr*zi) || isinf(zr*zi)) {
//        printf("Bad cfuncpoint: z[%d] = z(%f) = %f +i%f",it,w,zr,zi);
//       }
     GSL_SET_COMPLEX(&z[it],zr/Nx,zi/Nx);
   }
  return;
}

max_stars_repo_name: TantraLabs/gonum

max_issues_repo_name: TantraLabs/gonum

max_forks_repo_name: TantraLabs/gonum


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

typedef struct Vortex
{
    double x, y;
    double vorticity;

    int active;
} Vortex;

/* possibly should be replaced with one function? */
double Qfield_x(double x, double y)
{
    double r = x*x+y*y;
    return -y/(2*M_PI*fmax(r,1/1024/1024));
//    return -y/(2*M_PI*fmax(r,1/1024));
}

double Qfield_y(double x, double y)
{
    double r = x*x+y*y;
    return  x/(2*M_PI*fmax(r,1/1024/1024));
}

double Pfield_x(double x, double y)
{
    double r = hypot(x,y);
    if(r == 0)
        return 1.;
    return x/r*exp(-r);
}

double Pfield_y(double x, double y)
{
    double r = hypot(x,y);
    if(r == 0)
        return 1.;
    return y/r*exp(-r);
}

double Mfield(double x, double y)
{
    double r2= x*x+y*y;
    double r = sqrt(r2);
    if(x==0 && y==0)
        return 0.;
    return (r+1.)*exp(-r)/r2;
}

int VortexInBody(double x, double y)
{
    if(x*x+y*y<1)
        return 1;
    else
        return 0;
}

int main(int argc, char **argv)
{
    int NumberOfPanels = 12, MaxVortexes = 16500, Niterations = 200,      i,j,k,m,n, iterator;
    double nu = 1/13, tau = 0.005/*, eps = 0.002*/, rho = 1.;
    double vx_inf = 1., vy_inf = 0.;

    double  PanelNodes [NumberOfPanels][2],
            PanelMids  [NumberOfPanels][2],
            PanelNorms [NumberOfPanels][2],
            PanelTaus  [NumberOfPanels][2],
            PanelLength[NumberOfPanels];

    double DeployPoints[NumberOfPanels][2];

    Vortex InFlow[MaxVortexes], NextInFlow[MaxVortexes];

    int ActiveVortexesInFLow=0;

    if(1)/* cylinder profile forming & deploys */
    {

        for(i=0;i<NumberOfPanels;i++)
        {
            PanelNodes[i][0] = cos(2.*M_PI/NumberOfPanels*i);
            PanelNodes[i][1] = sin(2.*M_PI/NumberOfPanels*i);

            DeployPoints[i][0] = 1.05*PanelNodes[i][0];
            DeployPoints[i][1] = 1.05*PanelNodes[i][1];

            PanelNorms[i][0] = cos(2.*M_PI/NumberOfPanels*(i+0.5));
            PanelNorms[i][1] = sin(2.*M_PI/NumberOfPanels*(i+0.5));

        }

        for(i=0;i<NumberOfPanels-1;i++)
        {
            PanelMids[i][0] = (PanelNodes[i][0]+PanelNodes[i+1][0])/2;
            PanelMids[i][1] = (PanelNodes[i][1]+PanelNodes[i+1][1])/2;

            PanelTaus[i][0] = (PanelNodes[i+1][0]-PanelNodes[i][0]);
            PanelTaus[i][1] = (PanelNodes[i+1][1]-PanelNodes[i][1]);
            PanelLength[i]  = sqrt(PanelTaus[i][0]*PanelTaus[i][0]+
                                   PanelTaus[i][1]*PanelTaus[i][1]);
            PanelTaus[i][0] /= PanelLength[i];
            PanelTaus[i][1] /= PanelLength[i];
//            PanelNorms[i][0] = -PanelTaus[i][1];
//            PanelNorms[i][1] = PanelTaus[i][0];
        }
        PanelMids[NumberOfPanels-1][0] = (PanelNodes[NumberOfPanels-1][0]+PanelNodes[0][0])/2;
        PanelMids[NumberOfPanels-1][1] = (PanelNodes[NumberOfPanels-1][1]+PanelNodes[0][1])/2;
        PanelTaus[NumberOfPanels-1][0] = (PanelNodes[0][0]-PanelNodes[NumberOfPanels-1][0]);
        PanelTaus[NumberOfPanels-1][1] = (PanelNodes[0][1]-PanelNodes[NumberOfPanels-1][1]);
        PanelLength[NumberOfPanels-1]  = hypot(PanelTaus[NumberOfPanels-1][0],PanelTaus[NumberOfPanels-1][1]);
        PanelTaus[NumberOfPanels-1][0] /= PanelLength[NumberOfPanels-1];
        PanelTaus[NumberOfPanels-1][1] /= PanelLength[NumberOfPanels-1];
//        PanelNorms[Np-1][0] = -PanelTaus[Np-1][1];
//        PanelNorms[Np-1][1] = PanelTaus[Np-1][0];
    }

    // defining and completing Vortex Birth Matrix
    gsl_matrix *OriginalVortexGenerationMatrix = gsl_matrix_alloc(NumberOfPanels+1,NumberOfPanels+1);
    gsl_matrix *VortexGenerationMatrix = gsl_matrix_alloc(NumberOfPanels+1,NumberOfPanels+1);

    for(i=0; i<NumberOfPanels; i++)
    {
        for(j=0; j<NumberOfPanels; j++)
        {
            gsl_matrix_set(VortexGenerationMatrix,i,j,
            PanelNorms[i][0]*Qfield_x(PanelMids[i][0]-DeployPoints[j][0],
                                      PanelMids[i][1]-DeployPoints[j][1])
           +PanelNorms[i][1]*Qfield_y(PanelMids[i][0]-DeployPoints[j][0],
                                      PanelMids[i][1]-DeployPoints[j][1]));
        }
    }
    for(i=0;i<NumberOfPanels;i++)
    {
        gsl_matrix_set(VortexGenerationMatrix,NumberOfPanels,i,1);
        gsl_matrix_set(VortexGenerationMatrix,i,NumberOfPanels,1);
    }
    gsl_matrix_set(VortexGenerationMatrix,NumberOfPanels,NumberOfPanels,0);
    gsl_matrix_memcpy(OriginalVortexGenerationMatrix,VortexGenerationMatrix);

    // recieve initial aproximation for vortexes on surface
    gsl_vector *VelocityOnInfProjections = gsl_vector_alloc(NumberOfPanels+1); //Velocity On Infinity
    gsl_vector *VelocityProjections = gsl_vector_alloc(NumberOfPanels+1); //Velocity from active vortexes
//    gsl_vector *velocities_on_surface = gsl_vector_alloc(Np+1); //full velocity on surface = VOI + VFV
    gsl_vector *NewVorticities = gsl_vector_alloc(NumberOfPanels+1);

    for(i=0;i<NumberOfPanels;i++)
    {
        gsl_vector_set(VelocityOnInfProjections,i,-(PanelNorms[i][0]*vx_inf+PanelNorms[i][1]*vy_inf));
    }
    gsl_vector_set(VelocityOnInfProjections,NumberOfPanels,0.);

    gsl_permutation * p = gsl_permutation_alloc (NumberOfPanels+1);
    gsl_linalg_LU_decomp(VortexGenerationMatrix,p,&k);
    gsl_linalg_LU_solve (VortexGenerationMatrix,p,VelocityOnInfProjections,NewVorticities);

    FILE *forces;
    forces = fopen("forces.dat","w");

    for(iterator=0; iterator<Niterations; iterator++)
    {
        double f_x = 0., f_y = 0./*, M_z = 0.*/;

        //add generated vortexes to flow
        for(i=0;i<NumberOfPanels;i++)
        {
            InFlow[ActiveVortexesInFLow+i].active = 1;
            InFlow[ActiveVortexesInFLow+i].vorticity = gsl_vector_get(NewVorticities,i);
            InFlow[ActiveVortexesInFLow+i].x = DeployPoints[i][0];
            InFlow[ActiveVortexesInFLow+i].y = DeployPoints[i][1];
        }
        ActiveVortexesInFLow += NumberOfPanels;

        // time iteration
        printf("active vortexes %d\n",ActiveVortexesInFLow);

        for(i=0;i<ActiveVortexesInFLow;i++)
        {
            while(InFlow[i].active!=1)i++;
            double eps = sqrt(4/3)*PanelLength[0]/1000;
            double  I0 = 2*M_PI*eps*eps,
                    I1 = 0.,
                    I2_x = 0., I2_y = 0.,
                    I3_x = 0., I3_y = 0.;

            double u_x = vx_inf;
            double u_y = vy_inf;

            //Epsilon update
            if(iterator>5 && 0)
            {
                double sq1=1., sq2=1., sq3=1.;
                for(j=0;j<ActiveVortexesInFLow;j++)
                {
                    double temp = (InFlow[i].x-InFlow[j].x)*(InFlow[i].x-InFlow[j].x)
                                 +(InFlow[i].y-InFlow[j].y)*(InFlow[i].y-InFlow[j].y);
                    if(i!=j && temp < sq1
                            && InFlow[j].active==1
                            && InFlow[i].vorticity*InFlow[j].vorticity>0
                       )
                    {
                        sq1 = temp;
                        sq2 = sq1;
                        sq3 = sq2;
                    }
                }
                eps = sqrt( (sq1+sq2+sq3)/3 );
            }

            //Vortex-caused part of velocity
            for(j=0;j<ActiveVortexesInFLow;j++)
            {
                double rx = InFlow[i].x-InFlow[j].x;
                double ry = InFlow[i].y-InFlow[j].y;
                if( InFlow[j].active==1 &&  i!=j /*&& fabs(rx) > 1/1024.*0.25 && fabs(ry) > 1/1024.*0.25*/ )
                {
                    u_x += InFlow[j].vorticity*Qfield_x(rx,ry);
                    u_y += InFlow[j].vorticity*Qfield_y(rx,ry);
                }
            }

            //Viscous component of velocity
            for(j=0;j<ActiveVortexesInFLow;j++)
            {
                if( InFlow[j].active==1 )
                {
                    I2_x += InFlow[j].vorticity*Pfield_x((InFlow[i].x-InFlow[j].x)/eps,(InFlow[i].y-InFlow[j].y)/eps)/eps;
                    I2_y += InFlow[j].vorticity*Pfield_y((InFlow[i].x-InFlow[j].x)/eps,(InFlow[i].y-InFlow[j].y)/eps)/eps;
                    I1  += InFlow[j].vorticity*exp(-hypot((InFlow[i].x-InFlow[j].x),(InFlow[i].y-InFlow[j].y))/eps);
                }
            }
            for(j=0;j<NumberOfPanels;j++)
            {
                double Rr = hypot((InFlow[i].x-PanelMids[j][0]),(InFlow[i].y-PanelMids[j][1]));
                I3_x -= PanelNorms[j][0]*exp(-Rr/eps)*PanelLength[j];
                I3_y -= PanelNorms[j][1]*exp(-Rr/eps)*PanelLength[j];
                I0   -= ((InFlow[i].x-PanelMids[j][0])*PanelNorms[j][0]
                        +(InFlow[i].y-PanelMids[j][1])*PanelNorms[j][1])
                        *PanelLength[j]
                        *Mfield((InFlow[i].x-PanelMids[j][0])/eps,
                        (InFlow[i].y-PanelMids[j][1])/eps);

            }

            //        printf("%f %f\n",u_x,u_y);
            u_x += nu*(-I2_x/I1+I3_x/I0);
            u_y += nu*(-I2_y/I1+I3_y/I0);
            //        printf("%f %f\n",u_x,u_y);

            if(tau*u_x < 1000 && tau*u_y < 1000 )
            {
                NextInFlow[i].x = InFlow[i].x + tau*u_x;
                NextInFlow[i].y = InFlow[i].y + tau*u_y;

//                int minRmid = 0;
//контроль протекания
                double Rmid = sqrt((PanelMids[0][0]-NextInFlow[i].x)*(PanelMids[0][0]-NextInFlow[i].x)+
                                   (PanelMids[0][1]-NextInFlow[i].y)*(PanelMids[0][1]-NextInFlow[i].y));
                double projectionDistOnMidVector = (PanelMids[0][0]*NextInFlow[i].x+PanelMids[0][1]*NextInFlow[i].y)/
                                                   (PanelMids[0][0]*PanelMids[0][0]+PanelMids[0][1]*PanelMids[0][1]);

                for(m=1;m<NumberOfPanels;m++)
                {
                    double temp = sqrt((PanelMids[m][0]-NextInFlow[i].x)*(PanelMids[m][0]-NextInFlow[i].x)+
                                       (PanelMids[m][1]-NextInFlow[i].y)*(PanelMids[m][1]-NextInFlow[i].y));
                    if(fmin(Rmid,temp)<Rmid)
                    {
                        Rmid = temp;
//                        minRmid = m;
                        projectionDistOnMidVector = (PanelMids[m][0]*NextInFlow[i].x+PanelMids[m][1]*NextInFlow[i].y)/
                                                    (PanelMids[m][0]*PanelMids[m][0]+PanelMids[m][1]*PanelMids[m][1]);
                    }
                }

//                if(!VortexInBody(NextInFlow[i].x,NextInFlow[i].y))
                if(projectionDistOnMidVector > 1.)
                {
                    NextInFlow[i].active = InFlow[i].active;
                    NextInFlow[i].vorticity = InFlow[i].vorticity;
                }
                else
                {
                    NextInFlow[i].active = 3;
                    NextInFlow[i].vorticity = InFlow[i].vorticity;
                }
            }
            else
            {
                NextInFlow[i].x = InFlow[i].x;
                NextInFlow[i].y = InFlow[i].y;
                NextInFlow[i].active = 2;
                NextInFlow[i].vorticity = InFlow[i].vorticity;
            }
//            if(NextInFlow[i].vorticity<0.01)
        }

        for(m = 0; m < NumberOfPanels; m++)
        {
            f_x += -rho*(gsl_vector_get(NewVorticities,m)/tau* PanelMids[m][1]);
            f_y +=  rho*(gsl_vector_get(NewVorticities,m)/tau* PanelMids[m][0]);
        }

        fprintf(forces,"%f %f %f\n",iterator*tau, f_x, f_y);

        //updating coords
        for(i=0;i<ActiveVortexesInFLow;i++)
        {
            InFlow[i] = NextInFlow[i];
        }

        //definition of speed from vortexes after iteration
        double VelocityProjInControlPoint[NumberOfPanels];
        for(j=0; j<NumberOfPanels; j++)
            VelocityProjInControlPoint[j] = PanelNorms[j][0]*vx_inf + PanelNorms[j][1]*vy_inf;

        for(i=0; i<ActiveVortexesInFLow; i++)
        {
            if(InFlow[i].active==1)
            {
                for(j=0; j<NumberOfPanels; j++)
                {
                    VelocityProjInControlPoint[j] +=
                            (PanelNorms[j][0]*Qfield_x((PanelMids[j][0]-InFlow[i].x),
                                                       (PanelMids[j][1]-InFlow[i].y))
                            +PanelNorms[j][1]*Qfield_y((PanelMids[j][0]-InFlow[i].x),
                                                       (PanelMids[j][1]-InFlow[i].y)))
                            *InFlow[i].vorticity;
                }
            }
        }

        for(j=0; j<NumberOfPanels; j++)
            gsl_vector_set(VelocityProjections,j,-VelocityProjInControlPoint[j]);
        gsl_vector_set(VelocityProjections,NumberOfPanels,0);

        gsl_matrix_memcpy(VortexGenerationMatrix, OriginalVortexGenerationMatrix);

        gsl_linalg_LU_decomp(VortexGenerationMatrix,p,&k);
        gsl_linalg_LU_solve (VortexGenerationMatrix,p,VelocityProjections,NewVorticities);

        if(iterator%10==0)
        {
            /*
            FILE *VortexesOut, *BodyOut, *averagedOut;
            VortexesOut = fopen("MVE_flow.dat","w");
            BodyOut =     fopen("MVE_body.dat","w");
            averagedOut = fopen("MVE_aver.dat","w");

            for(i=0;i<ActiveVortexesInFLow;i++)
            {
                if(InFlow[i].active==1)
                    fprintf(VortexesOut,"%f %f %f\n",InFlow[i].x,InFlow[i].y,InFlow[i].vorticity);
            }

            for(i=0; i<Np; i++)
                fprintf(BodyOut,"%f %f\n", DeployPoints[i][0], DeployPoints[i][1]);

            double averagedVorticity[101][101];

            for(i=0;i<ActiveVortexesInFLow;i++)
            {
                if(fabs(InFlow[i].x)<=10 && fabs(InFlow[i].y)<=10 && InFlow[i].active == 1)
                {
                    averagedVorticity[(int)(50+(InFlow[i].x)/0.2)][(int)(50+(InFlow[i].y)/0.2)] += InFlow[i].vorticity*0.04;
                }
            }

            for(i=0;i<101;i++)
            {
                for(j=0;j<101;j++)
                    fprintf(averagedOut,"%f %f %f\n", -10+i*0.2, -10+j*0.02, averagedVorticity[i][j]);

            }

            fclose(averagedOut);
            fclose(VortexesOut);
            fclose(BodyOut);
            int temp;
            system("screen -S gnuplot -X stuff \"replot^M\" ");
            scanf("%d",&temp);
            if(ActiveVortexesInFLow<16000)
                Niterations +=temp;
            else
                return 0;
//            system("screen -S gnuplot -X stuff \"plot \'MVE_flow.dat\' using 1:2, \'MVE_body.dat\' w line^M\" ");
*/
            char FileName[32];
            sprintf(FileName,"VelocityFiled%6.6d.dat",iterator);
            FILE *VelocityField;
            VelocityField = fopen(FileName,"w");
            fprintf(VelocityField,"TITLE = \"Flow Model\"\n");
            fprintf(VelocityField,"VARIABLES = \"x\", \"y\", \"vx\", \"vy\"\n");
            fprintf(VelocityField,"ZONE T=\"Frame 0\", I=%d, J=%d\n", 200, 100);
            for(i=0;i<200;i++)
            {
                for(j=0;j<100;j++)
                {
                    double vx=vx_inf, vy=vy_inf;
                    for(n=0;n<ActiveVortexesInFLow;n++)
                    {
                        if(InFlow[n].active==1)
                        {
                            vx += InFlow[n].vorticity*Qfield_x((0.04*i-2)-InFlow[n].x,
                                                               (0.04*j-2)-InFlow[n].y);
                            vy += InFlow[n].vorticity*Qfield_y((0.04*i-2)-InFlow[n].x,
                                                               (0.04*j-2)-InFlow[n].y);
                        }
                    }
                    fprintf(VelocityField,"%f %f %f %f\n",0.04*i-2, 0.04*j-2, vx, vy);
                }
            }
        }
    }
/*
    //the most interesting part: output of results
    FILE *VortexesOut, *BodyOut, *averagedOut;
    VortexesOut = fopen("MVE_flow.dat","w");
    BodyOut =     fopen("MVE_body.dat","w");
    averagedOut = fopen("MVE_aver.dat","w");

    for(i=0;i<ActiveVortexesInFLow;i++)
    {
        if(InFlow[i].active==1)
            fprintf(VortexesOut,"%f %f %f\n",InFlow[i].x,InFlow[i].y,InFlow[i].vorticity);
    }

    for(i=0; i<Np; i++)
        fprintf(BodyOut,"%f %f\n", DeployPoints[i][0], DeployPoints[i][1]);

    double averagedVorticity[101][101];

    for(i=0;i<ActiveVortexesInFLow;i++)
    {
        if(fabs(InFlow[i].x)<=10 && fabs(InFlow[i].y)<=10 && InFlow[i].active == 1)
        {
            averagedVorticity[(int)(50+(InFlow[i].x)/0.2)][(int)(50+(InFlow[i].y)/0.2)] += InFlow[i].vorticity*0.04;
        }
    }

    for(i=0;i<101;i++)
    {
        for(j=0;j<101;j++)
            fprintf(averagedOut,"%f %f %f\n", -10+i*0.2, -10+j*0.02, averagedVorticity[i][j]);

    }

    fclose(averagedOut);
    fclose(VortexesOut);
    fclose(BodyOut);
    */
    return 0;
}


text:
#pragma once

#include <memory>
#include <string>

#include <gsl/gsl>
#include <nonstd/optional.hpp>

#include "chainerx/backend.h"
#include "chainerx/context.h"
#include "chainerx/device.h"
#include "chainerx/kernel_registry.h"

namespace chainerx {
namespace cuda {

class CudaDevice;
class CudaBackend;

namespace cuda_internal {

// Creates a device instance.
// This function is meant to be used from the backend class. Never use it for other purpose.
// This is defined in cuda_internal namespace in order to make it a friend of CudaDevice
// class.
gsl::owner<CudaDevice*> CreateDevice(CudaBackend& backend, int index);

}  // namespace cuda_internal

class CudaBackend : public Backend {
public:
    static constexpr const char* kDefaultName = "cuda";
    static constexpr const size_t kCudnnDefaultMaxWorkspaceSize = 8 * 1024 * 1024;
    static constexpr const char* kCudnnMaxWorkspaceSizeEnvVarName = "CHAINERX_CUDNN_MAX_WORKSPACE_SIZE";

    using Backend::Backend;

    std::string GetName() const override;

    int GetDeviceCount() const override;

    bool SupportsTransfer(Device& src_device, Device& dst_device) override;

    // TODO(hvy): Move to CudaDevice.
    // Sets maximum cuDNN workspace size.
    // This value is shared across threads.
    void SetCudnnMaxWorkspaceSize(size_t max_workspace_size);

    // TODO(hvy): Move to CudaDevice.
    // Gets maximum cuDNN workspace size.
    size_t GetCudnnMaxWorkspaceSize();

    static KernelRegistry& GetGlobalKernelRegistry() {
        static gsl::owner<KernelRegistry*> global_kernel_registry = new KernelRegistry{};
        return *global_kernel_registry;
    }

protected:
    KernelRegistry& GetParentKernelRegistry() override { return GetGlobalKernelRegistry(); }

private:
    std::unique_ptr<Device> CreateDevice(int index) override;

    // TODO(hvy): Move to CudaDevice.
    nonstd::optional<size_t> cudnn_max_workspace_size_{};

    std::mutex mutex_;
};

}  // namespace cuda
}  // namespace chainerx

max_stars_repo_name: tkerola/chainer

max_issues_repo_name: tkerola/chainer

max_forks_repo_name: tkerola/chainer


text:
/**
 *
 * @file core_cpltmg_condex.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:47 2014
 *
 **/
#include <math.h>
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex32_t
 *
 *  CORE_cpltmg_condexq generates the Q used in condex matrix generation
 *
 *  See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-999898
 *  gallery('condex',n,4,100)
 *
 *  Returns a "counter-example" matrix to a condition estimator. It has order n
 *  and scalar parameter theta (default 100).
 *
 *  LAPACK (RCOND): It is the inverse of this matrix that is a counter-example.
 *
 *******************************************************************************
 *
 * @param[in] M
 *         The number of rows of the matrix Q used in condex generation. M >= 0.
 *
 * @param[in] N
 *         The number of columns of the matrix A to be generated. N >= 0.
 *
 * @param[out] Q
 *         On entry, the M-by-3 matrix to be initialized.
 *         On exit, the housholder reflectors required for condex generation.
 *
 * @param[in] LDQ
 *         The leading dimension of the matrix Q. LDQ >= max(1,M).
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cpltmg_condexq = PCORE_cpltmg_condexq
#define CORE_cpltmg_condexq PCORE_cpltmg_condexq
#endif
void CORE_cpltmg_condexq( int M, int N, PLASMA_Complex32_t *Q, int LDQ )
{
    PLASMA_Complex32_t tau[3];
    PLASMA_Complex32_t *tQ = Q;
    int i;

    /* First column is [ 1 ... 1 ] */
    for( i=0; i < M; i++, tQ++ )
        *tQ = (PLASMA_Complex32_t)1.0;

    /* Second column is [1 0 0 ... 0] */
    tQ = Q + LDQ;
    *tQ = (PLASMA_Complex32_t)1.;
    tQ++;
    memset( tQ, 0, (M-1) * sizeof(PLASMA_Complex32_t) );

    /* Third column is  (-1)^i * (1. + i / (N-1)) */
    tQ = Q + 2 * LDQ;
    for( i=0; i<M; i++, tQ++ )
        *tQ = (PLASMA_Complex32_t)( cpowf( -1.0, (float)i ) * (1.0 + (float)i/(N-1) ) );

    /* Generate orthogonal projector */
    LAPACKE_cgeqrf( LAPACK_COL_MAJOR, M, 3,    Q, LDQ, tau );
    LAPACKE_cungqr( LAPACK_COL_MAJOR, M, 3, 3, Q, LDQ, tau );

    return;
}