text:
/*
  NAME:
     calc_splitnmerge
  PURPOSE:
     calculates the split and merge hierarchy after an proj_EM convergence
  CALLING SEQUENCE:
     calc_splitnmerge(struct datapoint * data,int N,
     struct gaussian * gaussians, int K, gsl_matrix * qij,
     int * snmhierarchy){
  INPUT:
     data      - the data
     N         - number of data points
     gaussians - model gaussians
     K         - number of gaussians
     qij       - matrix of log(posterior likelihoods)
     bs        - bs from previous EM estimate
  OUTPUT:
     snmhierarchy - the hierarchy, first row has the highest prioriry,
                    goes down from there
  REVISION HISTORY:
     2008-09-21 - Written Bovy
*/
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <proj_gauss_mixtures.h>

void calc_splitnmerge(struct datapoint * data,int N,
		      struct gaussian * gaussians, int K,
		      gsl_matrix * qij, int * snmhierarchy){
  gsl_matrix * tempqij = gsl_matrix_alloc(N,K);
  gsl_matrix_memcpy(tempqij,qij);
  gsl_matrix * Jmerge = gsl_matrix_alloc(K,K);
  gsl_matrix_set_all(Jmerge,-1.);
  int kk1, kk2, kk, ii,maxsnm= K*(K-1)*(K-2)/2;
  int d = (gaussians->VV)->size1;//dim of mm
  double temp1,temp2,temp;
  for (kk1 = 0; kk1 != K; ++kk1)
    for (kk2 = kk1+1; kk2 != K; ++kk2){
      temp = 0.;
      for (ii=0; ii != N; ++ii){
	//make them all exps
	temp1 = exp(gsl_matrix_get(qij,ii,kk1));
	temp2 = exp(gsl_matrix_get(qij,ii,kk2));
	temp += temp1*temp2;
      }
      gsl_matrix_set(Jmerge,kk1,kk2,temp);
    }

  //Then calculate Jsplit
  gsl_vector * Jsplit = gsl_vector_alloc(K);
  gsl_vector * Jsplit_temp = gsl_vector_alloc(K);
  gsl_vector_set_all(Jsplit,-1.);
  //if there is missing data, fill in the missing data
  struct missingdatapoint{
    gsl_vector *ww;
  };
  struct missingdatapoint * missingdata;
  missingdata = (struct missingdatapoint *) malloc(N * sizeof (struct missingdatapoint) );
  for (ii=0; ii != N; ++ii){
    missingdata->ww = gsl_vector_alloc(d);
    ++missingdata;
  }
  missingdata -= N;

  gsl_matrix * tempRR,* tempVV;
  gsl_vector * tempSS,* tempwork;
  gsl_vector * expectedww = gsl_vector_alloc(d);
  //gsl_vector_view tempUcol;
  double lambda;
  int di,signum;
  for (ii=0; ii != N; ++ii){
    //First check whether there is any missing data
    if ((data->ww)->size == d){
      gsl_vector_memcpy(missingdata->ww,data->ww);
      ++missingdata;
      ++data;
      continue;
    }

    /*AS IT STANDS THE MISSING DATA PART IS *NOT* IMPLEMENTED CORRECTLY:
      A CORRECT IMPLEMENTATION NEEDS THE NULL SPACE OF THE PROJECTION
      MATRIX WHICH CAN BE FOUND FROM THE FULL SINGULAR VALUE DECOMPOSITIIN,
      UNFORTUNATELY GSL DOES NOT COMPUTE THE FULL SVD, BUT ONLY THE THIN SVD.
      LAPACK MIGHT DO, BUT MIGHT NOT BE INSTALLED (?) AND THIS MIGHT BE HARD TO IMPLEMENT.

      INDICATED BELOW ARE THE SECTION THAT WOULD HAVE TO BE FIXED TO MAKE THIS WORK
     */

    //calculate expectation, for this we need to calculate the bbijs (EXACTLY THE SAME AS IN PROJ_EM, SHOULD WRITE GENERAL FUNCTION TO DO THIS)
    gsl_vector_set_zero(expectedww);
    for (kk = 0; kk != K; ++kk){
      //prepare...
      di = (data->SS)->size1;
      p = gsl_permutation_alloc (di);
      wminusRm = gsl_vector_alloc (di);
      gsl_vector_memcpy(wminusRm,data->ww);
      TinvwminusRm = gsl_vector_alloc (di);
      Tij = gsl_matrix_alloc(di,di);
      gsl_matrix_memcpy(Tij,data->SS);
      Tij_inv = gsl_matrix_alloc(di,di);
      VRT = gsl_matrix_alloc(d,di);
      VRTTinv = gsl_matrix_alloc(d,di);
      Rtrans = gsl_matrix_alloc(d,di);
      //Calculate Tij
      gsl_matrix_transpose_memcpy(Rtrans,data->RR);
      gsl_blas_dsymm(CblasLeft,CblasUpper,1.0,gaussians->VV,Rtrans,0.0,VRT);//Only the upper right part of VV is calculated
      gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1.0,data->RR,VRT,1.0,Tij);//This is Tij
      //Calculate LU decomp of Tij and Tij inverse
      gsl_linalg_LU_decomp(Tij,p,&signum);
      gsl_linalg_LU_invert(Tij,p,Tij_inv);
      //Calculate Tijinv*(w-Rm)
      gsl_blas_dgemv(CblasNoTrans,-1.0,data->RR,gaussians->mm,1.0,wminusRm);
      gsl_blas_dsymv(CblasUpper,1.0,Tij_inv,wminusRm,0.0,TinvwminusRm);
      //Now calculate bij and Bij
      gsl_vector_memcpy(bs->bbij,gaussians->mm);
      gsl_blas_dgemv(CblasNoTrans,1.0,VRT,TinvwminusRm,1.0,bs->bbij);
      //..and add the result to expectedww
      gsl_vector_scale(bs->bbij,exp(gsl_matrix_get(qij,ii,kk)));
      gsl_vector_add(expectedww,bs->bbij);
      //Clean up
      gsl_permutation_free (p);
      gsl_vector_free(wminusRm);
      gsl_vector_free(TinvwminusRm);
      gsl_matrix_free(Tij);
      gsl_matrix_free(Tij_inv);
      gsl_matrix_free(VRT);
      gsl_matrix_free(VRTTinv);
      gsl_matrix_free(Rtrans);
      ++gaussians;
    }
    gaussians -= K;
    //if missing, fill in the missing data
    tempRR = gsl_matrix_alloc((data->RR)->size2,(data->RR)->size1);//will hold the transpose of RR
    //tempVV = gsl_matrix_alloc((data->RR)->size1,(data->RR)->size1);
    //tempSS = gsl_vector_alloc((data->RR)->size1);
    //tempwork = gsl_vector_alloc((data->RR)->size1);
    gsl_matrix_transpose_memcpy(tempRR,data->RR);
    gsl_blas_dgemv(CblasNoTrans,1.,tempRR,data->ww,0.,missingdata->ww);
    //gsl_linalg_SV_decomp(tempRR,tempVV,tempSS,tempwork);
    //compute the missing data THIS PART IS NOT IMPLEMENTED CORRECTLY
    //for (kk = 0; kk != d-(data->ww)->size; ++kk){
      //tempUcol = gsl_matrix_column(tempRR,d-1-kk);
    //gsl_blas_ddot(&(tempUcol.vector),expectedww,&lambda);
    //gsl_vector_scale(&(tempUcol.vector),lambda);
    //gsl_vector_add(missingdata->ww,&(tempUcol.vector));
    //}
    ++missingdata;
    ++data;
    //free
    gsl_matrix_free(tempRR);
    //gsl_matrix_free(tempVV);
    //gsl_vector_free(tempSS);
    //gsl_vector_free(tempwork);
  }
  data -= N;
  missingdata -= N;

  //then for every gaussian, calculate the KL divergence between the local data density and the l-th gaussian
  double tempsplit;
  p = gsl_permutation_alloc (d);
  tempVV= gsl_matrix_alloc(d,d);
  tempRR= gsl_matrix_alloc(d,d);
  tempSS= gsl_vector_alloc(d);
  tempwork= gsl_vector_alloc(d);
  for (kk = 0; kk != K; ++kk){
    //calculate qil/ql factors
    normalize_row(tempqij,kk,false,true,0.);
    //calculate inverse of V and det(V)
    gsl_matrix_memcpy(tempVV,gaussians->VV);
    gsl_linalg_LU_decomp(tempVV,p,&signum);
    gsl_linalg_LU_invert(tempVV,p,tempRR);//tempRR now has the inverse of VV
    tempsplit = d * halflogtwopi + 0.5 * gsl_linalg_LU_lndet(tempVV);
    for (ii=0; ii != N; ++ii){
      if (exp(gsl_matrix_get(tempqij,ii,kk)) == 0.){
	++missingdata;
	continue;
      }
      tempsplit += gsl_matrix_get(tempqij,ii,kk) * exp(gsl_matrix_get(tempqij,ii,kk));
      gsl_vector_memcpy(tempSS,gaussians->mm);
      gsl_vector_scale(tempSS,-1.);
      gsl_vector_add(tempSS,missingdata->ww);
      gsl_blas_dgemv(CblasNoTrans,1.0,tempRR,tempSS,0.,tempwork);
      gsl_blas_ddot(tempSS,tempwork,&lambda);
      tempsplit += 0.5 * exp(gsl_matrix_get(tempqij,ii,kk)) * lambda;
      ++missingdata;
    }
    gsl_vector_set(Jsplit,kk,tempsplit);
    //printf("Jsplit for gaussian %i = %f\n",kk,tempsplit);
    missingdata -= N;
    ++gaussians;
  }
  gaussians -= K;

  //free
  gsl_permutation_free(p);
  gsl_matrix_free(tempRR);
  gsl_matrix_free(tempVV);
  gsl_vector_free(tempSS);
  gsl_vector_free(tempwork);

  //and put everything in the hierarchy
  size_t maxj, maxk, maxl;
  for (kk1 = 0; kk1 != maxsnm; kk1 += (K-2)){
    gsl_matrix_max_index(Jmerge,&maxj,&maxk);
    gsl_vector_memcpy(Jsplit_temp,Jsplit);
    gsl_vector_set(Jsplit_temp,maxj,-1.);
    gsl_vector_set(Jsplit_temp,maxk,-1.);
    for (kk2=0; kk2 != K-2; ++kk2){
      maxl = gsl_vector_max_index(Jsplit_temp);
      gsl_vector_set(Jsplit_temp,maxl,-1.);
      *(snmhierarchy++)= maxj;
      *(snmhierarchy++)= maxk;
      *(snmhierarchy++)= maxl;
      //printf("j = %i, k = %i, l = %i\n",(int)maxj,(int)maxk,(int)maxl);
    }
    //then set it to zero and find the next
    gsl_matrix_set(Jmerge,maxj,maxk,-1.);
  }
  snmhierarchy -= 3*maxsnm;

  //clean up
  gsl_matrix_free(Jmerge);
  gsl_vector_free(Jsplit);
  gsl_vector_free(Jsplit_temp);

  return ;
}


text:
/* -*- linux-c -*- */
/* fewbody_ks.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 <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>
#include "fewbody.h"

/* the dot product of two K-S vectors */
double fb_ks_dot(double x[4], double y[4])
{
	return(x[0] * y[0] + x[1] * y[1] + x[2] * y[2] + x[3] * y[3]);
}

/* the modulus of a K-S vector */
double fb_ks_mod(double x[4])
{
	return(sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3]));
}

/* calculate Q given q */
void fb_calc_Q(double q[4], double Q[4])
{
	double u[4];

	u[0] = sqrt(0.5 * (fb_mod(q) + fabs(q[0])));
	u[1] = 0.5 * q[1] / u[0];
	u[2] = 0.5 * q[2] / u[0];
	u[3] = 0.0;

	if (q[0] > 0.0) {
		Q[0] = u[0];
		Q[1] = u[1];
		Q[2] = u[2];
		Q[3] = u[3];
	} else {
		Q[0] = u[1];
		Q[1] = u[0];
		Q[2] = u[3];
		Q[3] = u[2];
	}
}

/* calculate the K-S matrix */
void fb_calc_ksmat(double Q[4], double Qmat[4][4])
{
	Qmat[0][0] = Q[0];
	Qmat[0][1] = -Q[1];
	Qmat[0][2] = -Q[2];
	Qmat[0][3] = Q[3];

	Qmat[1][0] = Q[1];
	Qmat[1][1] = Q[0];
	Qmat[1][2] = -Q[3];
	Qmat[1][3] = -Q[2];

	Qmat[2][0] = Q[2];
	Qmat[2][1] = Q[3];
	Qmat[2][2] = Q[0];
	Qmat[2][3] = Q[1];

	Qmat[3][0] = Q[3];
	Qmat[3][1] = -Q[2];
	Qmat[3][2] = Q[1];
	Qmat[3][3] = -Q[0];
}

/* calculate the a matrix */
void fb_calc_amat(double **a, int nstar, int kstar)
{
	int i, j, k;

	/* first zero out the matrix */
	for (i=0; i<nstar; i++) {
		for (k=0; k<kstar; k++) {
			a[i][k] = 0.0;
		}
	}

	/* then set the non-zero components */
	k = -1;
	for (i=0; i<nstar-1; i++) {
		for (j=i+1; j<nstar; j++) {
			k++;
			a[i][k] = 1.0;
			a[j][k] = -1.0;
		}
	}
}

/* calculate the T matrix */
void fb_calc_Tmat(double **a, double *m, double **T, int nstar, int kstar)
{
	int u, v, e;

	for (u=0; u<kstar; u++) {
		for (v=0; v<kstar; v++) {
			T[u][v] = 0.0;
			for (e=0; e<nstar; e++) {
				T[u][v] += a[e][u] * a[e][v] / m[e];
			}
			T[u][v] *= 0.5;
		}
	}
}

/* the derivatives function for the GSL ODE integrator */
int fb_ks_func(double s, const double *y, double *f, void *params)
{
	int k, l, m, kstar;
	double *M, **amat, **Tmat, Einit;
	double **Q, **P, **p, **A, **TP, **TQ, **UQ, *d;
	double Qmat[4][4], Pmat[4][4], Astar[4], T, U, val, L, H, G, GT, GU;

	/* set parameters */
	kstar = (*(fb_ks_params_t *) params).kstar;
	M = (*(fb_ks_params_t *) params).M;
	amat = (*(fb_ks_params_t *) params).amat;
	Tmat = (*(fb_ks_params_t *) params).Tmat;
	Einit = (*(fb_ks_params_t *) params).Einit;

	/* allocate memory */
	Q = fb_malloc_matrix(kstar, 4);
	P = fb_malloc_matrix(kstar, 4);
	p = fb_malloc_matrix(kstar, 4);
	A = fb_malloc_matrix(kstar, 4);
	TP = fb_malloc_matrix(kstar, 4);
	TQ = fb_malloc_matrix(kstar, 4);
	UQ = fb_malloc_matrix(kstar, 4);
	d = fb_malloc_vector(kstar);

	/* set Q_k, P_k, and p_k */
	for (k=0; k<kstar; k++) {
		/* read Q_k and P_k directly from y[] */
		for (l=0; l<4; l++) {
			Q[k][l] = y[k*8+1+l];
			P[k][l] = y[k*8+4+1+l];
		}

		/* calculate the KS matrix from Q_k */
		fb_calc_ksmat(Q[k], Qmat);

		/* and then calculate p_k */
		for (l=0; l<4; l++) {
			p[k][l] = 0.0;
			for (m=0; m<4; m++) {
				p[k][l] += Qmat[l][m] * P[k][m];
			}
			p[k][l] /= 2.0 * fb_ks_dot(Q[k], Q[k]);
		}
	}

	/* set intermediate variables */
	T = 0.0;
	U = 0.0;
	for (k=0; k<kstar; k++) {
		/* first calculate the K-S matrix for this value of k */
		fb_calc_ksmat(Q[k], Qmat);

		/* calculate A_k */
		for (m=0; m<4; m++) {
			A[k][m] = 0.0;
			for (l=0; l<kstar; l++) {
				A[k][m] += Tmat[k][l] * p[l][m];
			}
		}

		/* then d_k */
		d[k] = fb_ks_dot(A[k], p[k]);

		/* increment the sums for T and U */
		T += d[k];
		U += M[k] / fb_ks_dot(Q[k], Q[k]);

		/* then calculate the first partial derivative of T with respect to P_k */
		for (m=0; m<4; m++) {
			TP[k][m] = 0.0;
			for (l=0; l<4; l++) {
				TP[k][m] += Qmat[l][m] * A[k][l];
			}
			TP[k][m] /= fb_ks_dot(Q[k], Q[k]);
		}

		/* A^*_k */
		Astar[0] = A[k][0];
		Astar[1] = A[k][1];
		Astar[2] = A[k][2];
		Astar[3] = -A[k][3];

		/* calculate the K-S matrix for P_k */
		fb_calc_ksmat(P[k], Pmat);

		/* then calculate the first partial derivatives of T and U */
		for (l=0; l<4; l++) {
			val = 0.0;
			for (m=0; m<4; m++) {
				val += Pmat[m][l] * Astar[m];
			}
			TQ[k][l] = (val - 4.0 * d[k] * Q[k][l]) / fb_ks_dot(Q[k], Q[k]);
			UQ[k][l] = -2.0 * M[k] * Q[k][l] / fb_sqr(fb_ks_dot(Q[k], Q[k]));
		}
	}

	/* set the Lagrangian, Hamiltonian, and the initial energy */
	L = T + U;
	H = T - U;
	G = (H - Einit) / L;
	GT = (1.0 - G) / L;
	GU = -(1.0 + G) / L;

	/* set derivatives */
	f[0] = 1.0 / L;
	for (k=0; k<kstar; k++) {
		for (l=0; l<4; l++) {
			f[k*8+1+l] = GT * TP[k][l];
			f[k*8+4+1+l] = -GT * TQ[k][l] - GU * UQ[k][l];
		}
	}

	/* free memory */
	fb_free_matrix(Q);
	fb_free_matrix(P);
	fb_free_matrix(p);
	fb_free_matrix(A);
	fb_free_matrix(TP);
	fb_free_matrix(TQ);
	fb_free_matrix(UQ);
	fb_free_vector(d);

	/* all done */
	return(GSL_SUCCESS);
}

/* function to calculate the Einit parameter for the integrator */
/* the code here is a verbatim copy from the first part of fb_ks_func */
double fb_ks_Einit(const double *y, fb_ks_params_t params)
{
	int k, l, m;
	double **Q, **P, **p, **A, *d;
	double Qmat[4][4], T, U;

	/* allocate memory */
	Q = fb_malloc_matrix(params.kstar, 4);
	P = fb_malloc_matrix(params.kstar, 4);
	p = fb_malloc_matrix(params.kstar, 4);
	A = fb_malloc_matrix(params.kstar, 4);
	d = fb_malloc_vector(params.kstar);

	/* set Q_k, P_k, and p_k */
	for (k=0; k<params.kstar; k++) {
		/* read Q_k and P_k directly from y[] */
		for (l=0; l<4; l++) {
			Q[k][l] = y[k*8+1+l];
			P[k][l] = y[k*8+4+1+l];
		}

		/* calculate the KS matrix from Q_k */
		fb_calc_ksmat(Q[k], Qmat);

		/* and then calculate p_k */
		for (l=0; l<4; l++) {
			p[k][l] = 0.0;
			for (m=0; m<4; m++) {
				p[k][l] += Qmat[l][m] * P[k][m];
			}
			p[k][l] /= 2.0 * fb_ks_dot(Q[k], Q[k]);
		}
	}

	/* set intermediate variables */
	T = 0.0;
	U = 0.0;
	for (k=0; k<params.kstar; k++) {
		/* first calculate the K-S matrix for this value of k */
		fb_calc_ksmat(Q[k], Qmat);

		/* calculate A_k */
		for (m=0; m<4; m++) {
			A[k][m] = 0.0;
			for (l=0; l<params.kstar; l++) {
				A[k][m] += params.Tmat[k][l] * p[l][m];
			}
		}

		/* then d_k */
		d[k] = fb_ks_dot(A[k], p[k]);

		/* increment the sums for T and U */
		T += d[k];
		U += params.M[k] / fb_ks_dot(Q[k], Q[k]);
	}

	/* free memory */
	fb_free_matrix(Q);
	fb_free_matrix(P);
	fb_free_matrix(p);
	fb_free_matrix(A);
	fb_free_vector(d);

	return(T-U);
}

/* function to convert from Euclidean coordinates to K-S coordinates */
void fb_euclidean_to_ks(fb_obj_t **star, double *y, int nstar, int kstar)
{
	int i, j, k, l, m;
	double **q, **p, Q[4], P[4], Qmat[4][4];

	q = fb_malloc_matrix(kstar, 4);
	p = fb_malloc_matrix(kstar, 4);

	/* then calculate q_k and p_k */
	k = -1;
	for (i=0; i<nstar-1; i++) {
		for (j=i+1; j<nstar; j++) {
			k++;
			for (l=0; l<3; l++) {
				q[k][l] = star[i]->x[l] - star[j]->x[l];
				p[k][l] = (star[i]->m * star[i]->v[l] - star[j]->m * star[j]->v[l])/((double) nstar);
			}
			q[k][3] = 0.0;
			p[k][3] = 0.0;
		}
	}

	/* and then Q_k and P_k */
	for (k=0; k<kstar; k++) {
		/* calculate Q_k */
		fb_calc_Q(q[k], Q);

		/* then P_k */
		fb_calc_ksmat(Q, Qmat);
		for (l=0; l<4; l++) {
			P[l] = 0.0;
			for (m=0; m<4; m++) {
				P[l] += Qmat[m][l] * p[k][m];
			}
			P[l] *= 2.0;

			/* then set y_i */
			y[8*k+l+1] = Q[l];
			y[8*k+l+4+1] = P[l];
		}
	}

	fb_free_matrix(q);
	fb_free_matrix(p);
}

/* function to convert from K-S coordinates to Euclidean coordinates */
void fb_ks_to_euclidean(double *y, fb_obj_t **star, int nstar, int kstar)
{
	int i, j, k, l;
	double *m, mtot, **Q, **P, **q, **p, Qmat[4][4];

	/* allocate memory */
	m = fb_malloc_vector(nstar);
	Q = fb_malloc_matrix(kstar, 4);
	P = fb_malloc_matrix(kstar, 4);
	q = fb_malloc_matrix(kstar, 4);
	p = fb_malloc_matrix(kstar, 4);

	/* set the masses */
	mtot = 0.0;
	for (i=0; i<nstar; i++) {
		m[i] = star[i]->m;
		mtot += m[i];
	}

	/* set Q_k, P_k, and p_k */
	for (k=0; k<kstar; k++) {
		/* read Q_k and P_k directly from y[] */
		for (l=0; l<4; l++) {
			Q[k][l] = y[k*8+1+l];
			P[k][l] = y[k*8+4+1+l];
		}

		/* calculate the KS matrix from Q_k */
		fb_calc_ksmat(Q[k], Qmat);

		/* and then calculate q_k and p_k */
		for (l=0; l<4; l++) {
			q[k][l] = 0.0;
			p[k][l] = 0.0;
			for (j=0; j<4; j++) {
				q[k][l] += Qmat[l][j] * Q[k][j];
				p[k][l] += Qmat[l][j] * P[k][j];
			}
			p[k][l] /= 2.0 * fb_ks_dot(Q[k], Q[k]);
		}
	}

	/* set r_i and v_i */
	for (i=0; i<nstar; i++) {
		for (k=0; k<3; k++) {
			star[i]->x[k] = 0.0;
			star[i]->v[k] = 0.0;
		}

		for (j=i+1; j<nstar; j++) {
			for (k=0; k<3; k++) {
				star[i]->x[k] += m[j] * q[FB_KS_K(i, j, nstar)][k];
				star[i]->v[k] += p[FB_KS_K(i, j, nstar)][k];
			}
		}

		for (j=0; j<i; j++) {
			for (k=0; k<3; k++) {
				star[i]->x[k] += - m[j] * q[FB_KS_K(j, i, nstar)][k];
				star[i]->v[k] += - p[FB_KS_K(j, i, nstar)][k];
			}
		}

		for (k=0; k<3; k++) {
			star[i]->x[k] /= mtot;
			star[i]->v[k] /= m[i];
		}
	}

	/* free memory */
	fb_free_vector(m);
	fb_free_matrix(Q);
	fb_free_matrix(P);
	fb_free_matrix(q);
	fb_free_matrix(p);
}

ext: c

max_stars_repo_name: gnodvi/cosmos

max_issues_repo_name: gnodvi/cosmos

max_forks_repo_name: gnodvi/cosmos


text:
/*
 *  iaf_cond_alpha_mc.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST 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.
 *
 *  NEST 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 NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef IAF_COND_ALPHA_MC_H
#define IAF_COND_ALPHA_MC_H

// Generated includes:
#include "config.h"

#ifdef HAVE_GSL

// C++ includes:
#include <vector>

// C includes:
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

// Includes from nestkernel:
#include "archiving_node.h"
#include "connection.h"
#include "event.h"
#include "nest_types.h"
#include "recordables_map.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"

// Includes from sli:
#include "dictdatum.h"
#include "name.h"

/* BeginDocumentation
Name: iaf_cond_alpha_mc - PROTOTYPE Multi-compartment conductance-based leaky
                          integrate-and-fire neuron model.

Description:
THIS MODEL IS A PROTOTYPE FOR ILLUSTRATION PURPOSES. IT IS NOT YET
FULLY TESTED. USE AT YOUR OWN PERIL!

iaf_cond_alpha_mc is an implementation of a multi-compartment spiking
neuron using IAF dynamics with conductance-based synapses. It serves
mainly to illustrate the implementation of multicompartment models in
NEST.

The model has three compartments: soma, proximal and distal dendrite,
labeled as s, p, and d, respectively. Compartments are connected through
passive conductances as follows

C_m.s d/dt V_m.s = ... - g_sp ( V_m.s - V_m.p )

C_m.p d/dt V_m.p = ... - g_sp ( V_m.p - V_m.s ) - g_pd ( V_m.p - V_m.d )

C_m.d d/dt V_m.d = ...                          - g_pd ( V_m.d - V_m.p )

A spike is fired when the somatic membrane potential exceeds threshold,
V_m.s >= V_th. After a spike, somatic membrane potential is clamped to
a reset potential, V_m.s == V_reset, for the refractory period. Dendritic
membrane potentials are not manipulated after a spike.

There is one excitatory and one inhibitory conductance-based synapse
onto each compartment, with alpha-function time course. The alpha
function is normalised such that an event of weight 1.0 results in a
peak current of 1 nS at t = tau_syn. Each compartment can also receive
current input from a current generator, and an external (rheobase)
current can be set for each compartment.

Synapses, including those for injection external currents, are addressed through
the receptor types given in the receptor_types entry of the state dictionary.
Note that in contrast to the single-compartment iaf_cond_alpha model, all
synaptic weights must be positive numbers!

Parameters:
The following parameters can be set in the status dictionary. Parameters
for each compartment are collected in a sub-dictionary; these sub-dictionaries
are called "soma", "proximal", and "distal", respectively. In the list below,
these parameters are marked with an asterisk.

V_m*         double - Membrane potential in mV
E_L*         double - Leak reversal potential in mV.
C_m*         double - Capacity of the membrane in pF
E_ex*        double - Excitatory reversal potential in mV.
E_in*        double - Inhibitory reversal potential in mV.
g_L*         double - Leak conductance in nS;
tau_syn_ex*  double - Rise time of the excitatory synaptic alpha function in ms.
tau_syn_in*  double - Rise time of the inhibitory synaptic alpha function in ms.
I_e*         double - Constant input current in pA.

g_sp         double - Conductance connecting soma and proximal dendrite, in nS.
g_pd         double - Conductance connecting proximal and distal dendrite, in
                      nS.
t_ref        double - Duration of refractory period in ms.
V_th         double - Spike threshold in mV.
V_reset      double - Reset potential of the membrane in mV.

Example:
See examples/nest/mc_neuron.py.

Remark:
This is a prototype for illustration which has undergone only limited testing.
Details of the implementation and user-interface will likely change.
USE AT YOUR OWN PERIL!

Sends: SpikeEvent

Receives: SpikeEvent, CurrentEvent, DataLoggingRequest

References:

Meffin, H., Burkitt, A. N., & Grayden, D. B. (2004). An analytical
model for the large, fluctuating synaptic conductance state typical of
neocortical neurons in vivo. J.  Comput. Neurosci., 16, 159-175.

Bernander, O ., Douglas, R. J., Martin, K. A. C., & Koch, C. (1991).
Synaptic background activity influences spatiotemporal integration in
single pyramidal cells.  Proc. Natl. Acad. Sci. USA, 88(24),
11569-11573.

Author: Plesser

SeeAlso: iaf_cond_alpha
*/

namespace nest
{
/**
 * Function computing right-hand side of ODE for GSL solver.
 * @note Must be declared here so we can befriend it in class.
 * @note Must have C-linkage for passing to GSL.
 * @note No point in declaring it inline, since it is called
 *       through a function pointer.
 */
extern "C" int
iaf_cond_alpha_mc_dynamics( double, const double*, double*, void* );

/**
 * @note All parameters that occur for both compartments
 *       and dendrite are stored as C arrays, with index 0 being soma.
 */
class iaf_cond_alpha_mc : public Archiving_Node
{

  // Boilerplate function declarations --------------------------------

public:
  iaf_cond_alpha_mc();
  iaf_cond_alpha_mc( const iaf_cond_alpha_mc& );
  ~iaf_cond_alpha_mc();

  /**
   * Import sets of overloaded virtual functions.
   * @see Technical Issues / Virtual Functions: Overriding, Overloading, and
   * Hiding
   */
  using Node::handle;
  using Node::handles_test_event;

  port send_test_event( Node&, rport, synindex, bool );

  void handle( SpikeEvent& );
  void handle( CurrentEvent& );
  void handle( DataLoggingRequest& );

  port handles_test_event( SpikeEvent&, rport );
  port handles_test_event( CurrentEvent&, rport );
  port handles_test_event( DataLoggingRequest&, rport );

  void get_status( DictionaryDatum& ) const;
  void set_status( const DictionaryDatum& );

private:
  void init_state_( const Node& proto );
  void init_buffers_();
  void calibrate();
  void update( Time const&, const long, const long );

  // Enumerations and constants specifying structure and properties ----

  //! Compartments, NCOMP is number
  enum Compartments_
  {
    SOMA = 0,
    PROX,
    DIST,
    NCOMP
  };

  /**
   * Minimal spike receptor type.
   * @note Start with 1 so we can forbid port 0 to avoid accidental
   *       creation of connections with no receptor type set.
   */
  static const port MIN_SPIKE_RECEPTOR = 1;

  /**
   * Spike receptors.
   */
  enum SpikeSynapseTypes
  {
    SOMA_EXC = MIN_SPIKE_RECEPTOR,
    SOMA_INH,
    PROX_EXC,
    PROX_INH,
    DIST_EXC,
    DIST_INH,
    SUP_SPIKE_RECEPTOR
  };

  static const size_t NUM_SPIKE_RECEPTORS =
    SUP_SPIKE_RECEPTOR - MIN_SPIKE_RECEPTOR;

  /**
   * Minimal current receptor type.
   *  @note Start with SUP_SPIKE_RECEPTOR to avoid any overlap and
   *        accidental mix-ups.
   */
  static const port MIN_CURR_RECEPTOR = SUP_SPIKE_RECEPTOR;

  /**
   * Current receptors.
   */
  enum CurrentSynapseTypes
  {
    I_SOMA = MIN_CURR_RECEPTOR,
    I_PROX,
    I_DIST,
    SUP_CURR_RECEPTOR
  };

  static const size_t NUM_CURR_RECEPTORS =
    SUP_CURR_RECEPTOR - MIN_CURR_RECEPTOR;

  // Friends --------------------------------------------------------

  friend int
  iaf_cond_alpha_mc_dynamics( double, const double*, double*, void* );

  friend class RecordablesMap< iaf_cond_alpha_mc >;
  friend class UniversalDataLogger< iaf_cond_alpha_mc >;

  // Parameters ------------------------------------------------------

  /**
   * Independent parameters of the model.
   * These parameters must be passed to the iteration function that
   * is passed to the GSL ODE solvers. Since the iteration function
   * is a C++ function with C linkage, the parameters can be stored
   * in a C++ struct with member functions, as long as we just pass
   * it by void* from C++ to C++ function. The struct must be public,
   * though, since the iteration function is a function with C-linkage,
   * whence it cannot be a member function of iaf_cond_alpha_mc.
   * @note One could achieve proper encapsulation by an extra level
   *       of indirection: Define the iteration function as a member
   *       function, plus an additional wrapper function with C linkage.
   *       Then pass a struct containing a pointer to the node and a
   *       pointer-to-member-function to the iteration function as void*
   *       to the wrapper function. The wrapper function can then invoke
   *       the iteration function on the node (Stroustrup, p 418). But
   *       this appears to involved, and the extra indirections cost.
   */
  struct Parameters_
  {
    double V_th;                //!< Threshold Potential in mV
    double V_reset;             //!< Reset Potential in mV
    double t_ref;               //!< Refractory period in ms
    double g_conn[ NCOMP - 1 ]; //!< Conductances connecting compartments,
                                //!< in nS
    double g_L[ NCOMP ];        //!< Leak Conductance in nS
    double C_m[ NCOMP ];        //!< Membrane Capacitance in pF
    double E_ex[ NCOMP ];       //!< Excitatory reversal Potential in mV
    double E_in[ NCOMP ];       //!< Inhibitory reversal Potential in mV
    double E_L[ NCOMP ]; //!< Leak reversal Potential (aka resting potential)
                         //!< in mV
    double tau_synE[ NCOMP ]; //!< Synaptic Time Constant Excitatory Synapse
                              //!< in ms
    double tau_synI[ NCOMP ]; //!< Synaptic Time Constant for Inhibitory
                              //!< Synapse in ms
    double I_e[ NCOMP ];      //!< Constant Current in pA

    Parameters_();                     //!< Sets default parameter values
    Parameters_( const Parameters_& ); //!< needed to copy C-arrays
    Parameters_& operator=( const Parameters_& ); //!< needed to copy C-arrays

    void get( DictionaryDatum& ) const; //!< Store current values in dictionary
    void set( const DictionaryDatum& ); //!< Set values from dicitonary
  };

  // State variables  ------------------------------------------------------

  /**
   * State variables of the model.
   * @note Copy constructor and assignment operator required because
   *       of C-style array.
   */
public:
  struct State_
  {

    /**
     * Elements of state vector.
     * For the multicompartmental case here, these are offset values.
     * The state variables are stored in contiguous blocks for each
     * compartment, beginning with the soma.
     */
    enum StateVecElems_
    {
      V_M = 0,
      DG_EXC,
      G_EXC,
      DG_INH,
      G_INH,
      STATE_VEC_COMPS
    };

    //! total size of state vector
    static const size_t STATE_VEC_SIZE = STATE_VEC_COMPS * NCOMP;

    //! neuron state, must be C-array for GSL solver
    double y_[ STATE_VEC_SIZE ];
    int r_; //!< number of refractory steps remaining

    State_( const Parameters_& ); //!< Default initialization
    State_( const State_& );
    State_& operator=( const State_& );

    void get( DictionaryDatum& ) const;
    void set( const DictionaryDatum&, const Parameters_& );

    /**
     * Compute linear index into state array from compartment and element.
     * @param comp compartment index
     * @param elem elemet index
     * @note compartment argument is not of type Compartments_, since looping
     *       over enumerations does not work.
     */
    static size_t
    idx( size_t comp, StateVecElems_ elem )
    {
      return comp * STATE_VEC_COMPS + elem;
    }
  };

private:
  // Internal buffers --------------------------------------------------------

  /**
   * Buffers of the model.
   */
  struct Buffers_
  {
    Buffers_( iaf_cond_alpha_mc& ); //!<Sets buffer pointers to 0
    //! Sets buffer pointers to 0
    Buffers_( const Buffers_&, iaf_cond_alpha_mc& );

    //! Logger for all analog data
    UniversalDataLogger< iaf_cond_alpha_mc > logger_;

    /** buffers and sums up incoming spikes/currents
     *  @note Using STL vectors here to ensure initialization.
     */
    std::vector< RingBuffer > spikes_;
    std::vector< RingBuffer > currents_;

    /** GSL ODE stuff */
    gsl_odeiv_step* s_;    //!< stepping function
    gsl_odeiv_control* c_; //!< adaptive stepsize control function
    gsl_odeiv_evolve* e_;  //!< evolution function
    gsl_odeiv_system sys_; //!< struct describing system

    // IntergrationStep_ should be reset with the neuron on ResetNetwork,
    // but remain unchanged during calibration. Since it is initialized with
    // step_, and the resolution cannot change after nodes have been created,
    // it is safe to place both here.
    double step_;            //!< step size in ms
    double IntegrationStep_; //!< current integration time step, updated by GSL

    /**
     * Input currents injected by CurrentEvent.
     * This variable is used to transport the current applied into the
     * _dynamics function computing the derivative of the state vector.
     * It must be a part of Buffers_, since it is initialized once before
     * the first simulation, but not modified before later Simulate calls.
     */
    double I_stim_[ NCOMP ]; //!< External Stimulus in pA
  };

  // Internal variables ---------------------------------------------

  /**
   * Internal variables of the model.
   */
  struct Variables_
  {
    /** initial value to normalise excitatory synaptic conductance */
    double PSConInit_E_[ NCOMP ];

    /** initial value to normalise inhibitory synaptic conductance */
    double PSConInit_I_[ NCOMP ];

    int RefractoryCounts_;
  };

  // Access functions for UniversalDataLogger -------------------------------

  /**
   * Read out state vector elements, used by UniversalDataLogger
   * First template argument is component "name", second compartment "name".
   */
  template < State_::StateVecElems_ elem, Compartments_ comp >
  double
  get_y_elem_() const
  {
    return S_.y_[ S_.idx( comp, elem ) ];
  }

  //! Read out number of refractory steps, used by UniversalDataLogger
  double
  get_r_() const
  {
    return Time::get_resolution().get_ms() * S_.r_;
  }

  // Data members ----------------------------------------------------

  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;

  //! Table of compartment names
  static std::vector< Name > comp_names_;

  //! Dictionary of receptor types, leads to seg fault on exit, see #328
  // static DictionaryDatum receptor_dict_;

  //! Mapping of recordables names to access functions
  static RecordablesMap< iaf_cond_alpha_mc > recordablesMap_;
};

inline port
iaf_cond_alpha_mc::send_test_event( Node& target,
  rport receptor_type,
  synindex,
  bool )
{
  SpikeEvent e;
  e.set_sender( *this );
  return target.handles_test_event( e, receptor_type );
}

inline port
iaf_cond_alpha_mc::handles_test_event( SpikeEvent&, rport receptor_type )
{
  if ( receptor_type < MIN_SPIKE_RECEPTOR
    || receptor_type >= SUP_SPIKE_RECEPTOR )
  {
    if ( receptor_type < 0 || receptor_type >= SUP_CURR_RECEPTOR )
    {
      throw UnknownReceptorType( receptor_type, get_name() );
    }
    else
    {
      throw IncompatibleReceptorType( receptor_type, get_name(), "SpikeEvent" );
    }
  }
  return receptor_type - MIN_SPIKE_RECEPTOR;
}

inline port
iaf_cond_alpha_mc::handles_test_event( CurrentEvent&, rport receptor_type )
{
  if ( receptor_type < MIN_CURR_RECEPTOR || receptor_type >= SUP_CURR_RECEPTOR )
  {
    if ( receptor_type >= 0 && receptor_type < MIN_CURR_RECEPTOR )
    {
      throw IncompatibleReceptorType(
        receptor_type, get_name(), "CurrentEvent" );
    }
    else
    {
      throw UnknownReceptorType( receptor_type, get_name() );
    }
  }
  return receptor_type - MIN_CURR_RECEPTOR;
}

inline port
iaf_cond_alpha_mc::handles_test_event( DataLoggingRequest& dlr,
  rport receptor_type )
{
  if ( receptor_type != 0 )
  {
    if ( receptor_type < 0 || receptor_type >= SUP_CURR_RECEPTOR )
    {
      throw UnknownReceptorType( receptor_type, get_name() );
    }
    else
    {
      throw IncompatibleReceptorType(
        receptor_type, get_name(), "DataLoggingRequest" );
    }
  }
  return B_.logger_.connect_logging_device( dlr, recordablesMap_ );
}

inline void
iaf_cond_alpha_mc::get_status( DictionaryDatum& d ) const
{
  P_.get( d );
  S_.get( d );
  Archiving_Node::get_status( d );

  ( *d )[ names::recordables ] = recordablesMap_.get_list();

  /**
   * @todo dictionary construction should be done only once for
   * static member in default c'tor, but this leads to
   * a seg fault on exit, see #328
   */
  DictionaryDatum receptor_dict_ = new Dictionary();
  ( *receptor_dict_ )[ names::soma_exc ] = SOMA_EXC;
  ( *receptor_dict_ )[ names::soma_inh ] = SOMA_INH;
  ( *receptor_dict_ )[ names::soma_curr ] = I_SOMA;

  ( *receptor_dict_ )[ names::proximal_exc ] = PROX_EXC;
  ( *receptor_dict_ )[ names::proximal_inh ] = PROX_INH;
  ( *receptor_dict_ )[ names::proximal_curr ] = I_PROX;

  ( *receptor_dict_ )[ names::distal_exc ] = DIST_EXC;
  ( *receptor_dict_ )[ names::distal_inh ] = DIST_INH;
  ( *receptor_dict_ )[ names::distal_curr ] = I_DIST;

  ( *d )[ names::receptor_types ] = receptor_dict_;
}

inline void
iaf_cond_alpha_mc::set_status( const DictionaryDatum& d )
{
  Parameters_ ptmp = P_; // temporary copy in case of errors
  ptmp.set( d );         // throws if BadProperty
  State_ stmp = S_;      // temporary copy in case of errors
  stmp.set( d, ptmp );   // throws if BadProperty

  // We now know that (ptmp, stmp) are consistent. We do not
  // write them back to (P_, S_) before we are also sure that
  // the properties to be set in the parent class are internally
  // consistent.
  Archiving_Node::set_status( d );

  // if we get here, temporaries contain consistent set of properties
  P_ = ptmp;
  S_ = stmp;
}

} // namespace

#endif // HAVE_GSL
#endif // IAF_COND_ALPHA_MC_H

ext: h

max_stars_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster

max_issues_repo_name: zlchai/SNN-simulator-on-PYNQcluster

max_forks_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster


text:
/*!
 \file fasttransforms.h
 \brief fasttransforms.h is the main header file for FastTransforms.
*/
#ifndef FASTTRANSFORMS_H
#define FASTTRANSFORMS_H

#include <cblas.h>
#include <fftw3.h>

typedef double ft_complex[2];

#ifdef _OPENMP
    #include <omp.h>
    #define FT_GET_THREAD_NUM() omp_get_thread_num()
    #define FT_GET_NUM_THREADS() omp_get_num_threads()
    #define FT_GET_MAX_THREADS() omp_get_max_threads()
    #define FT_SET_NUM_THREADS(x) omp_set_num_threads(x)
#else
    #define FT_GET_THREAD_NUM() 0
    #define FT_GET_NUM_THREADS() 1
    #define FT_GET_MAX_THREADS() 1
    #define FT_SET_NUM_THREADS(x)
#endif

#define FT_CONCAT(prefix, name, suffix) prefix ## name ## suffix

void ft_horner(const int n, const double * c, const int incc, const int m, double * x, double * f);
void ft_hornerf(const int n, const float * c, const int incc, const int m, float * x, float * f);

void ft_clenshaw(const int n, const double * c, const int incc, const int m, double * x, double * f);
void ft_clenshawf(const int n, const float * c, const int incc, const int m, float * x, float * f);

void ft_orthogonal_polynomial_clenshaw(const int n, const double * c, const int incc, const double * A, const double * B, const double * C, const int m, double * x, double * phi0, double * f);
void ft_orthogonal_polynomial_clenshawf(const int n, const float * c, const int incc, const float * A, const float * B, const float * C, const int m, float * x, float * phi0, float * f);

void ft_eigen_eval(const int n, const double * c, const int incc, const double * A, const double * B, const double * C, const int m, double * x, const int sign, double * f);
void ft_eigen_evalf(const int n, const float * c, const int incc, const float * A, const float * B, const float * C, const int m, float * x, const int sign, float * f);
void ft_eigen_evall(const int n, const long double * c, const int incc, const long double * A, const long double * B, const long double * C, const int m, long double * x, const int sign, long double * f);
#if defined(FT_QUADMATH)
    #include <quadmath.h>
    typedef __float128 quadruple;
    void ft_eigen_evalq(const int n, const quadruple * c, const int incc, const quadruple * A, const quadruple * B, const quadruple * C, const int m, quadruple * x, const int sign, quadruple * f);
#endif

#define FT_SN (1U << 0)
#define FT_CN (1U << 1)
#define FT_DN (1U << 2)

#include "tdc.h"

/*!
  \brief Pre-compute a factorization of the connection coefficients between Legendre and Chebyshev polynomials in double precision so that ft_bfmv converts between expansions:
  \f[
  \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Legendre}} P_\ell(x) = \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Chebyshev}} T_\ell(x).
  \f]
  `normleg` and `normcheb` govern the normalizations, either standard ( == 0) or orthonormalized ( == 1).\n
  See also \ref ft_plan_legendre_to_chebyshevf, \ref ft_plan_legendre_to_chebyshevl, and \ref ft_mpfr_plan_legendre_to_chebyshev.
*/
ft_tb_eigen_FMM * ft_plan_legendre_to_chebyshev(const int normleg, const int normcheb, const int n);
/*!
  \brief Pre-compute a factorization of the connection coefficients between Chebyshev and Legendre polynomials in double precision so that ft_bfmv converts between expansions:
  \f[
  \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Chebyshev}} T_\ell(x) = \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Legendre}} P_\ell(x).
  \f]
  `normcheb` and `normleg` govern the normalizations, either standard ( == 0) or orthonormalized ( == 1).\n
  See also \ref ft_plan_chebyshev_to_legendref, \ref ft_plan_chebyshev_to_legendrel, and \ref ft_mpfr_plan_chebyshev_to_legendre.
*/
ft_tb_eigen_FMM * ft_plan_chebyshev_to_legendre(const int normcheb, const int normleg, const int n);
/*!
  \brief Pre-compute a factorization of the connection coefficients between ultraspherical polynomials in double precision so that ft_bfmv converts between expansions:
  \f[
  \sum_{\ell=0}^{n-1} c_\ell^{(1)} C_\ell^{(\lambda)}(x) = \sum_{\ell=0}^{n-1} c_\ell^{(2)} C_\ell^{(\mu)}(x).
  \f]
  `norm1` and `norm2` govern the normalizations, either standard ( == 0) or orthonormalized ( == 1).\n
  See also \ref ft_plan_ultraspherical_to_ultrasphericalf, \ref ft_plan_ultraspherical_to_ultrasphericall, and \ref ft_mpfr_plan_ultraspherical_to_ultraspherical.
*/
ft_tb_eigen_FMM * ft_plan_ultraspherical_to_ultraspherical(const int norm1, const int norm2, const int n, const double lambda, const double mu);
/*!
  \brief Pre-compute a factorization of the connection coefficients between Jacobi polynomials in double precision so that ft_bfmv converts between expansions:
  \f[
  \sum_{\ell=0}^{n-1} c_\ell^{(1)} P_\ell^{(\alpha,\beta)}(x) = \sum_{\ell=0}^{n-1} c_\ell^{(2)} P_\ell^{(\gamma,\delta)}(x).
  \f]
  `norm1` and `norm2` govern the normalizations, either standard ( == 0) or orthonormalized ( == 1).\n
  See also \ref ft_plan_jacobi_to_jacobif, \ref ft_plan_jacobi_to_jacobil, and \ref ft_mpfr_plan_jacobi_to_jacobi.
*/
ft_tb_eigen_FMM * ft_plan_jacobi_to_jacobi(const int norm1, const int norm2, const int n, const double alpha, const double beta, const double gamma, const double delta);
/*!
  \brief Pre-compute a factorization of the connection coefficients between Laguerre polynomials in double precision so that ft_bfmv converts between expansions:
  \f[
  \sum_{\ell=0}^{n-1} c_\ell^{(1)} L_\ell^{(\alpha)}(x) = \sum_{\ell=0}^{n-1} c_\ell^{(2)} L_\ell^{(\beta)}(x).
  \f]
  `norm1` and `norm2` govern the normalizations, either standard ( == 0) or orthonormalized ( == 1).\n
  See also \ref ft_plan_laguerre_to_laguerref, \ref ft_plan_laguerre_to_laguerrel, and \ref ft_mpfr_plan_laguerre_to_laguerre.
*/
ft_tb_eigen_FMM * ft_plan_laguerre_to_laguerre(const int norm1, const int norm2, const int n, const double alpha, const double beta);
/*!
  \brief Pre-compute a factorization of the connection coefficients between Jacobi and ultraspherical polynomials in double precision so that ft_bfmv converts between expansions:
  \f[
  \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Jacobi}} P_\ell^{(\alpha,\beta)}(x) = \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{ultraspherical}} C_\ell^{(\lambda)}(x).
  \f]
  `normjac` and `normultra` govern the normalizations, either standard ( == 0) or orthonormalized ( == 1).\n
  See also \ref ft_plan_jacobi_to_ultrasphericalf, \ref ft_plan_jacobi_to_ultrasphericall, and \ref ft_mpfr_plan_jacobi_to_ultraspherical.
*/
ft_tb_eigen_FMM * ft_plan_jacobi_to_ultraspherical(const int normjac, const int normultra, const int n, const double alpha, const double beta, const double lambda);
/*!
  \brief Pre-compute a factorization of the connection coefficients between ultraspherical and Jacobi polynomials in double precision so that ft_bfmv converts between expansions:
  \f[
  \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{ultraspherical}} C_\ell^{(\lambda)}(x) = \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Jacobi}} P_\ell^{(\alpha,\beta)}(x).
  \f]
  `normultra` and `normjac` govern the normalizations, either standard ( == 0) or orthonormalized ( == 1).\n
  See also \ref ft_plan_ultraspherical_to_jacobif, \ref ft_plan_ultraspherical_to_jacobil, and \ref ft_mpfr_plan_ultraspherical_to_jacobi.
*/
ft_tb_eigen_FMM * ft_plan_ultraspherical_to_jacobi(const int normultra, const int normjac, const int n, const double lambda, const double alpha, const double beta);
/*!
  \brief Pre-compute a factorization of the connection coefficients between Jacobi and Chebyshev polynomials in double precision so that ft_bfmv converts between expansions:
  \f[
  \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Jacobi}} P_\ell^{(\alpha,\beta)}(x) = \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Chebyshev}} T_\ell(x).
  \f]
  `normjac` and `normcheb` govern the normalizations, either standard ( == 0) or orthonormalized ( == 1).\n
  See also \ref ft_plan_jacobi_to_chebyshevf, \ref ft_plan_jacobi_to_chebyshevl, and \ref ft_mpfr_plan_jacobi_to_chebyshev.
*/
ft_tb_eigen_FMM * ft_plan_jacobi_to_chebyshev(const int normjac, const int normcheb, const int n, const double alpha, const double beta);
/*!
  \brief Pre-compute a factorization of the connection coefficients between Chebyshev and Jacobi polynomials in double precision so that ft_bfmv converts between expansions:
  \f[
  \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Chebyshev}} T_\ell(x) = \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Jacobi}} P_\ell^{(\alpha,\beta)}(x).
  \f]
  `normcheb` and `normjac` govern the normalizations, either standard ( == 0) or orthonormalized ( == 1).\n
  See also \ref ft_plan_chebyshev_to_jacobif, \ref ft_plan_chebyshev_to_jacobil, and \ref ft_mpfr_plan_chebyshev_to_jacobi.
*/
ft_tb_eigen_FMM * ft_plan_chebyshev_to_jacobi(const int normcheb, const int normjac, const int n, const double alpha, const double beta);
/*!
  \brief Pre-compute a factorization of the connection coefficients between ultraspherical and Chebyshev polynomials in double precision so that ft_bfmv converts between expansions:
  \f[
  \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{ultraspherical}} C_\ell^{(\lambda)}(x) = \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Chebyshev}} T_\ell(x).
  \f]
  `normultra` and `normcheb` govern the normalizations, either standard ( == 0) or orthonormalized ( == 1).\n
  See also \ref ft_plan_ultraspherical_to_chebyshevf, \ref ft_plan_ultraspherical_to_chebyshevl, and \ref ft_mpfr_plan_ultraspherical_to_chebyshev.
*/
ft_tb_eigen_FMM * ft_plan_ultraspherical_to_chebyshev(const int normultra, const int normcheb, const int n, const double lambda);
/*!
  \brief Pre-compute a factorization of the connection coefficients between Chebyshev and ultraspherical polynomials in double precision so that ft_bfmv converts between expansions:
  \f[
  \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{Chebyshev}} T_\ell(x) = \sum_{\ell=0}^{n-1} c_\ell^{\mathrm{ultraspherical}} C_\ell^{(\lambda)}(x).
  \f]
  `normcheb` and `normultra` govern the normalizations, either standard ( == 0) or orthonormalized ( == 1).\n
  See also \ref ft_plan_chebyshev_to_ultrasphericalf, \ref ft_plan_chebyshev_to_ultrasphericall, and \ref ft_mpfr_plan_chebyshev_to_ultraspherical.
*/
ft_tb_eigen_FMM * ft_plan_chebyshev_to_ultraspherical(const int normcheb, const int normultra, const int n, const double lambda);

ft_btb_eigen_FMM * ft_plan_associated_jacobi_to_jacobi(const int norm1, const int norm2, const int n, const int c, const double alpha, const double beta, const double gamma, const double delta);
ft_btb_eigen_FMM * ft_plan_associated_laguerre_to_laguerre(const int norm1, const int norm2, const int n, const int c, const double alpha, const double beta);
ft_btb_eigen_FMM * ft_plan_associated_hermite_to_hermite(const int norm1, const int norm2, const int n, const int c);

/// A single precision version of \ref ft_plan_legendre_to_chebyshev.
ft_tb_eigen_FMMf * ft_plan_legendre_to_chebyshevf(const int normleg, const int normcheb, const int n);
/// A single precision version of \ref ft_plan_chebyshev_to_legendre.
ft_tb_eigen_FMMf * ft_plan_chebyshev_to_legendref(const int normcheb, const int normleg, const int n);
/// A single precision version of \ref ft_plan_ultraspherical_to_ultraspherical.
ft_tb_eigen_FMMf * ft_plan_ultraspherical_to_ultrasphericalf(const int norm1, const int norm2, const int n, const float lambda, const float mu);
/// A single precision version of \ref ft_plan_jacobi_to_jacobi.
ft_tb_eigen_FMMf * ft_plan_jacobi_to_jacobif(const int norm1, const int norm2, const int n, const float alpha, const float beta, const float gamma, const float delta);
/// A single precision version of \ref ft_plan_laguerre_to_laguerre.
ft_tb_eigen_FMMf * ft_plan_laguerre_to_laguerref(const int norm1, const int norm2, const int n, const float alpha, const float beta);
/// A single precision version of \ref ft_plan_jacobi_to_ultraspherical.
ft_tb_eigen_FMMf * ft_plan_jacobi_to_ultrasphericalf(const int normjac, const int normultra, const int n, const float alpha, const float beta, const float lambda);
/// A single precision version of \ref ft_plan_ultraspherical_to_jacobi.
ft_tb_eigen_FMMf * ft_plan_ultraspherical_to_jacobif(const int normultra, const int normjac, const int n, const float lambda, const float alpha, const float beta);
/// A single precision version of \ref ft_plan_jacobi_to_chebyshev.
ft_tb_eigen_FMMf * ft_plan_jacobi_to_chebyshevf(const int normjac, const int normcheb, const int n, const float alpha, const float beta);
/// A single precision version of \ref ft_plan_chebyshev_to_jacobi.
ft_tb_eigen_FMMf * ft_plan_chebyshev_to_jacobif(const int normcheb, const int normjac, const int n, const float alpha, const float beta);
/// A single precision version of \ref ft_plan_ultraspherical_to_chebyshev.
ft_tb_eigen_FMMf * ft_plan_ultraspherical_to_chebyshevf(const int normultra, const int normcheb, const int n, const float lambda);
/// A single precision version of \ref ft_plan_chebyshev_to_ultraspherical.
ft_tb_eigen_FMMf * ft_plan_chebyshev_to_ultrasphericalf(const int normcheb, const int normultra, const int n, const float lambda);

ft_btb_eigen_FMMf * ft_plan_associated_jacobi_to_jacobif(const int norm1, const int norm2, const int n, const int c, const float alpha, const float beta, const float gamma, const float delta);
ft_btb_eigen_FMMf * ft_plan_associated_laguerre_to_laguerref(const int norm1, const int norm2, const int n, const int c, const float alpha, const float beta);
ft_btb_eigen_FMMf * ft_plan_associated_hermite_to_hermitef(const int norm1, const int norm2, const int n, const int c);

/// A long double precision version of \ref ft_plan_legendre_to_chebyshev.
ft_tb_eigen_FMMl * ft_plan_legendre_to_chebyshevl(const int normleg, const int normcheb, const int n);
/// A long double precision version of \ref ft_plan_chebyshev_to_legendre.
ft_tb_eigen_FMMl * ft_plan_chebyshev_to_legendrel(const int normcheb, const int normleg, const int n);
/// A long double precision version of \ref ft_plan_ultraspherical_to_ultraspherical.
ft_tb_eigen_FMMl * ft_plan_ultraspherical_to_ultrasphericall(const int norm1, const int norm2, const int n, const long double lambda, const long double mu);
/// A long double precision version of \ref ft_plan_jacobi_to_jacobi.
ft_tb_eigen_FMMl * ft_plan_jacobi_to_jacobil(const int norm1, const int norm2, const int n, const long double alpha, const long double beta, const long double gamma, const long double delta);
/// A long double precision version of \ref ft_plan_laguerre_to_laguerre.
ft_tb_eigen_FMMl * ft_plan_laguerre_to_laguerrel(const int norm1, const int norm2, const int n, const long double alpha, const long double beta);
/// A long double precision version of \ref ft_plan_jacobi_to_ultraspherical.
ft_tb_eigen_FMMl * ft_plan_jacobi_to_ultrasphericall(const int normjac, const int normultra, const int n, const long double alpha, const long double beta, const long double lambda);
/// A long double precision version of \ref ft_plan_ultraspherical_to_jacobi.
ft_tb_eigen_FMMl * ft_plan_ultraspherical_to_jacobil(const int normultra, const int normjac, const int n, const long double lambda, const long double alpha, const long double beta);
/// A long double precision version of \ref ft_plan_jacobi_to_chebyshev.
ft_tb_eigen_FMMl * ft_plan_jacobi_to_chebyshevl(const int normjac, const int normcheb, const int n, const long double alpha, const long double beta);
/// A long double precision version of \ref ft_plan_chebyshev_to_jacobi.
ft_tb_eigen_FMMl * ft_plan_chebyshev_to_jacobil(const int normcheb, const int normjac, const int n, const long double alpha, const long double beta);
/// A long double precision version of \ref ft_plan_ultraspherical_to_chebyshev.
ft_tb_eigen_FMMl * ft_plan_ultraspherical_to_chebyshevl(const int normultra, const int normcheb, const int n, const long double lambda);
/// A long double precision version of \ref ft_plan_chebyshev_to_ultraspherical.
ft_tb_eigen_FMMl * ft_plan_chebyshev_to_ultrasphericall(const int normcheb, const int normultra, const int n, const long double lambda);

ft_btb_eigen_FMMl * ft_plan_associated_jacobi_to_jacobil(const int norm1, const int norm2, const int n, const int c, const long double alpha, const long double beta, const long double gamma, const long double delta);
ft_btb_eigen_FMMl * ft_plan_associated_laguerre_to_laguerrel(const int norm1, const int norm2, const int n, const int c, const long double alpha, const long double beta);
ft_btb_eigen_FMMl * ft_plan_associated_hermite_to_hermitel(const int norm1, const int norm2, const int n, const int c);

#include <mpfr.h>

typedef struct {
    mpfr_t * data;
    int n;
    int b;
} ft_mpfr_triangular_banded;

void ft_mpfr_destroy_plan(mpfr_t * A, int n);
void ft_mpfr_trmv(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * x, mpfr_rnd_t rnd);
void ft_mpfr_trsv(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * x, mpfr_rnd_t rnd);
void ft_mpfr_trmm(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * B, int LDB, int N, mpfr_rnd_t rnd);
void ft_mpfr_trsm(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t * B, int LDB, int N, mpfr_rnd_t rnd);
// C -- Julia interoperability. Julia `BigFloat` does not have the same size as `mpfr_t`.
// So we give all Julia-owned data its own address, and dereference to retrieve the number.
void ft_mpfr_trmv_ptr(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t ** x, mpfr_rnd_t rnd);
void ft_mpfr_trsv_ptr(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t ** x, mpfr_rnd_t rnd);
void ft_mpfr_trmm_ptr(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t ** B, int LDB, int N, mpfr_rnd_t rnd);
void ft_mpfr_trsm_ptr(char TRANS, int n, mpfr_t * A, int LDA, mpfr_t ** B, int LDB, int N, mpfr_rnd_t rnd);

/// A multi-precision version of \ref ft_plan_legendre_to_chebyshev that returns a dense array of connection coefficients.
mpfr_t * ft_mpfr_plan_legendre_to_chebyshev(const int normleg, const int normcheb, const int n, mpfr_prec_t prec, mpfr_rnd_t rnd);
/// A multi-precision version of \ref ft_plan_chebyshev_to_legendre that returns a dense array of connection coefficients.
mpfr_t * ft_mpfr_plan_chebyshev_to_legendre(const int normcheb, const int normleg, const int n, mpfr_prec_t prec, mpfr_rnd_t rnd);
/// A multi-precision version of \ref ft_plan_ultraspherical_to_ultraspherical that returns a dense array of connection coefficients.
mpfr_t * ft_mpfr_plan_ultraspherical_to_ultraspherical(const int norm1, const int norm2, const int n, mpfr_srcptr lambda, mpfr_srcptr mu, mpfr_prec_t prec, mpfr_rnd_t rnd);
/// A multi-precision version of \ref ft_plan_jacobi_to_jacobi that returns a dense array of connection coefficients.
mpfr_t * ft_mpfr_plan_jacobi_to_jacobi(const int norm1, const int norm2, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_srcptr gamma, mpfr_srcptr delta, mpfr_prec_t prec, mpfr_rnd_t rnd);
/// A multi-precision version of \ref ft_plan_laguerre_to_laguerre that returns a dense array of connection coefficients.
mpfr_t * ft_mpfr_plan_laguerre_to_laguerre(const int norm1, const int norm2, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
/// A multi-precision version of \ref ft_plan_jacobi_to_ultraspherical that returns a dense array of connection coefficients.
mpfr_t * ft_mpfr_plan_jacobi_to_ultraspherical(const int normjac, const int normultra, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_srcptr lambda, mpfr_prec_t prec, mpfr_rnd_t rnd);
/// A multi-precision version of \ref ft_plan_ultraspherical_to_jacobi that returns a dense array of connection coefficients.
mpfr_t * ft_mpfr_plan_ultraspherical_to_jacobi(const int normultra, const int normjac, const int n, mpfr_srcptr lambda, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
/// A multi-precision version of \ref ft_plan_jacobi_to_chebyshev that returns a dense array of connection coefficients.
mpfr_t * ft_mpfr_plan_jacobi_to_chebyshev(const int normjac, const int normcheb, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
/// A multi-precision version of \ref ft_plan_chebyshev_to_jacobi that returns a dense array of connection coefficients.
mpfr_t * ft_mpfr_plan_chebyshev_to_jacobi(const int normcheb, const int normjac, const int n, mpfr_srcptr alpha, mpfr_srcptr beta, mpfr_prec_t prec, mpfr_rnd_t rnd);
/// A multi-precision version of \ref ft_plan_ultraspherical_to_chebyshev that returns a dense array of connection coefficients.
mpfr_t * ft_mpfr_plan_ultraspherical_to_chebyshev(const int normultra, const int normcheb, const int n, mpfr_srcptr lambda, mpfr_prec_t prec, mpfr_rnd_t rnd);
/// A multi-precision version of \ref ft_plan_chebyshev_to_ultraspherical that returns a dense array of connection coefficients.
mpfr_t * ft_mpfr_plan_chebyshev_to_ultraspherical(const int normcheb, const int normultra, const int n, mpfr_srcptr lambda, mpfr_prec_t prec, mpfr_rnd_t rnd);

/// Set the number of OpenMP threads.
void ft_set_num_threads(const int n);

/// Data structure to store sines and cosines of Givens rotations.
typedef struct ft_rotation_plan_s ft_rotation_plan;

/// Destroy a \ref ft_rotation_plan.
void ft_destroy_rotation_plan(ft_rotation_plan * RP);

ft_rotation_plan * ft_plan_rotsphere(const int n);

/// Convert a single vector of spherical harmonic coefficients in A with stride S from order m2 down to order m1.
void ft_kernel_sph_hi2lo(const ft_rotation_plan * RP, const int m1, const int m2, double * A, const int S);
/// Convert a single vector of spherical harmonic coefficients in A with stride S from order m1 up to order m2.
void ft_kernel_sph_lo2hi(const ft_rotation_plan * RP, const int m1, const int m2, double * A, const int S);

ft_rotation_plan * ft_plan_rottriangle(const int n, const double alpha, const double beta, const double gamma);

/// Convert a single vector of triangular harmonic coefficients in A with stride S from order m2 down to order m1.
void ft_kernel_tri_hi2lo(const ft_rotation_plan * RP, const int m1, const int m2, double * A, const int S);
/// Convert a single vector of triangular harmonic coefficients in A with stride S from order m1 up to order m2.
void ft_kernel_tri_lo2hi(const ft_rotation_plan * RP, const int m1, const int m2, double * A, const int S);

ft_rotation_plan * ft_plan_rotdisk(const int n, const double alpha, const double beta);

/// Convert a single vector of disk harmonic coefficients in A with stride S from order m2 down to order m1.
void ft_kernel_disk_hi2lo(const ft_rotation_plan * RP, const int m1, const int m2, double * A, const int S);
/// Convert a single vector of disk harmonic coefficients in A with stride S from order m1 up to order m2.
void ft_kernel_disk_lo2hi(const ft_rotation_plan * RP, const int m1, const int m2, double * A, const int S);

ft_rotation_plan * ft_plan_rotrectdisk(const int n, const double beta);

/// Convert a single vector of rectangularized disk harmonic coefficients in A with stride S from order m2 down to order m1.
void ft_kernel_rectdisk_hi2lo(const ft_rotation_plan * RP, const int m1, const int m2, double * A, const int S);
/// Convert a single vector of rectangularized disk harmonic coefficients in A with stride S from order m1 up to order m2.
void ft_kernel_rectdisk_lo2hi(const ft_rotation_plan * RP, const int m1, const int m2, double * A, const int S);

void ft_kernel_tet_hi2lo(const ft_rotation_plan * RP, const int L, const int m, double * A);
void ft_kernel_tet_lo2hi(const ft_rotation_plan * RP, const int L, const int m, double * A);

/// Data structure to store sines and cosines of Givens rotations for spin-weighted harmonics.
typedef struct ft_spin_rotation_plan_s ft_spin_rotation_plan;

/// Destroy a \ref ft_spin_rotation_plan.
void ft_destroy_spin_rotation_plan(ft_spin_rotation_plan * SRP);

ft_spin_rotation_plan * ft_plan_rotspinsphere(const int n, const int s);

/// Convert a single vector of spin-weighted spherical harmonic coefficients in A with stride S from order m down to order m%2.
void ft_kernel_spinsph_hi2lo(const ft_spin_rotation_plan * SRP, const int m, ft_complex * A, const int S);
/// Convert a single vector of spin-weighted spherical harmonic coefficients in A with stride S from order m%2 up to order m.
void ft_kernel_spinsph_lo2hi(const ft_spin_rotation_plan * SRP, const int m, ft_complex * A, const int S);

void ft_execute_sph_hi2lo(const ft_rotation_plan * RP, double * A, double * B, const int M);
void ft_execute_sph_lo2hi(const ft_rotation_plan * RP, double * A, double * B, const int M);

void ft_execute_sphv_hi2lo(const ft_rotation_plan * RP, double * A, double * B, const int M);
void ft_execute_sphv_lo2hi(const ft_rotation_plan * RP, double * A, double * B, const int M);

void ft_execute_tri_hi2lo(const ft_rotation_plan * RP, double * A, double * B, const int M);
void ft_execute_tri_lo2hi(const ft_rotation_plan * RP, double * A, double * B, const int M);

void ft_execute_disk_hi2lo(const ft_rotation_plan * RP, double * A, double * B, const int M);
void ft_execute_disk_lo2hi(const ft_rotation_plan * RP, double * A, double * B, const int M);

void ft_execute_rectdisk_hi2lo(const ft_rotation_plan * RP, double * A, double * B, const int M);
void ft_execute_rectdisk_lo2hi(const ft_rotation_plan * RP, double * A, double * B, const int M);

void ft_execute_tet_hi2lo(const ft_rotation_plan * RP1, const ft_rotation_plan * RP2, double * A, const int L, const int M);
void ft_execute_tet_lo2hi(const ft_rotation_plan * RP1, const ft_rotation_plan * RP2, double * A, const int L, const int M);

void ft_execute_spinsph_hi2lo(const ft_spin_rotation_plan * SRP, ft_complex * A, ft_complex * B, const int M);
void ft_execute_spinsph_lo2hi(const ft_spin_rotation_plan * SRP, ft_complex * A, ft_complex * B, const int M);

/// Data structure to store \ref ft_rotation_plan "ft_rotation_plan"s, arrays to represent 1D orthogonal polynomial transforms and their inverses, and Greek parameters.
typedef struct {
    ft_rotation_plan ** RP;
    double * B;
    double ** P;
    double ** Pinv;
    double alpha;
    double beta;
    double gamma;
    double delta;
    int NRP;
    int NP;
} ft_harmonic_plan;

/// Destroy a \ref ft_harmonic_plan.
void ft_destroy_harmonic_plan(ft_harmonic_plan * P);

/// Plan a spherical harmonic transform.
ft_harmonic_plan * ft_plan_sph2fourier(const int n);

/// Transform a spherical harmonic expansion to a bivariate Fourier series.
void ft_execute_sph2fourier(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int M);
/// Transform a bivariate Fourier series to a spherical harmonic expansion.
void ft_execute_fourier2sph(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int M);

void ft_execute_sphv2fourier(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int M);
void ft_execute_fourier2sphv(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int M);

/// Plan a triangular harmonic transform.
ft_harmonic_plan * ft_plan_tri2cheb(const int n, const double alpha, const double beta, const double gamma);

/// Transform a triangular harmonic expansion to a bivariate Chebyshev series.
void ft_execute_tri2cheb(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int M);
/// Transform a bivariate Chebyshev series to a triangular harmonic expansion.
void ft_execute_cheb2tri(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int M);

/// Plan a disk harmonic transform.
ft_harmonic_plan * ft_plan_disk2cxf(const int n, const double alpha, const double beta);

/// Transform a disk harmonic expansion to a Chebyshev--Fourier series.
void ft_execute_disk2cxf(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int M);
/// Transform a Chebyshev--Fourier series to a disk harmonic expansion.
void ft_execute_cxf2disk(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int M);

/// Plan a rectangularized disk harmonic transform.
ft_harmonic_plan * ft_plan_rectdisk2cheb(const int n, const double beta);

/// Transform a rectangularized disk harmonic expansion to a Chebyshev series.
void ft_execute_rectdisk2cheb(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int M);
/// Transform a Chebyshev series to a rectangularized disk harmonic expansion.
void ft_execute_cheb2rectdisk(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int M);

/// Plan a tetrahedral harmonic transform.
ft_harmonic_plan * ft_plan_tet2cheb(const int n, const double alpha, const double beta, const double gamma, const double delta);

/// Transform a tetrahedral harmonic expansion to a trivariate Chebyshev series.
void ft_execute_tet2cheb(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int L, const int M);
/// Transform a trivariate Chebyshev series to a tetrahedral harmonic expansion.
void ft_execute_cheb2tet(const char TRANS, const ft_harmonic_plan * P, double * A, const int N, const int L, const int M);

/// Data structure to store a \ref ft_spin_rotation_plan, and various arrays to represent 1D orthogonal polynomial transforms.
typedef struct {
    ft_spin_rotation_plan * SRP;
    ft_complex * B;
    ft_complex * P1;
    ft_complex * P2;
    ft_complex * P1inv;
    ft_complex * P2inv;
    int s;
} ft_spin_harmonic_plan;

/// Destroy a \ref ft_spin_harmonic_plan.
void ft_destroy_spin_harmonic_plan(ft_spin_harmonic_plan * P);

/// Plan a spin-weighted spherical harmonic transform.
ft_spin_harmonic_plan * ft_plan_spinsph2fourier(const int n, const int s);

/// Transform a spin-weighted spherical harmonic expansion to a bivariate Fourier series.
void ft_execute_spinsph2fourier(const char TRANS, const ft_spin_harmonic_plan * P, ft_complex * A, const int N, const int M);
/// Transform a bivariate Fourier series to a spin-weighted spherical harmonic expansion.
void ft_execute_fourier2spinsph(const char TRANS, const ft_spin_harmonic_plan * P, ft_complex * A, const int N, const int M);

int ft_fftw_init_threads(void);
void ft_fftw_plan_with_nthreads(const int n);

typedef struct {
    fftw_plan plantheta1;
    fftw_plan plantheta2;
    fftw_plan plantheta3;
    fftw_plan plantheta4;
    fftw_plan planphi;
    double * Y;
} ft_sphere_fftw_plan;

/// Destroy a \ref ft_sphere_fftw_plan.
void ft_destroy_sphere_fftw_plan(ft_sphere_fftw_plan * P);

ft_sphere_fftw_plan * ft_plan_sph_with_kind(const int N, const int M, const fftw_r2r_kind kind[3][1]);
/// Plan FFTW synthesis on the sphere.
ft_sphere_fftw_plan * ft_plan_sph_synthesis(const int N, const int M);
/// Plan FFTW analysis on the sphere.
ft_sphere_fftw_plan * ft_plan_sph_analysis(const int N, const int M);
ft_sphere_fftw_plan * ft_plan_sphv_synthesis(const int N, const int M);
ft_sphere_fftw_plan * ft_plan_sphv_analysis(const int N, const int M);

/// Execute FFTW synthesis on the sphere.
void ft_execute_sph_synthesis(const char TRANS, const ft_sphere_fftw_plan * P, double * X, const int N, const int M);
/// Execute FFTW analysis on the sphere.
void ft_execute_sph_analysis(const char TRANS, const ft_sphere_fftw_plan * P, double * X, const int N, const int M);

void ft_execute_sphv_synthesis(const char TRANS, const ft_sphere_fftw_plan * P, double * X, const int N, const int M);
void ft_execute_sphv_analysis(const char TRANS, const ft_sphere_fftw_plan * P, double * X, const int N, const int M);

typedef struct {
    fftw_plan planxy;
} ft_triangle_fftw_plan;

/// Destroy a \ref ft_triangle_fftw_plan.
void ft_destroy_triangle_fftw_plan(ft_triangle_fftw_plan * P);

ft_triangle_fftw_plan * ft_plan_tri_with_kind(const int N, const int M, const fftw_r2r_kind kind0, const fftw_r2r_kind kind1);
/// Plan FFTW synthesis on the triangle.
ft_triangle_fftw_plan * ft_plan_tri_synthesis(const int N, const int M);
/// Plan FFTW analysis on the triangle.
ft_triangle_fftw_plan * ft_plan_tri_analysis(const int N, const int M);

/// Execute FFTW synthesis on the triangle.
void ft_execute_tri_synthesis(const char TRANS, const ft_triangle_fftw_plan * P, double * X, const int N, const int M);
/// Execute FFTW analysis on the triangle.
void ft_execute_tri_analysis(const char TRANS, const ft_triangle_fftw_plan * P, double * X, const int N, const int M);

typedef struct {
    fftw_plan planxyz;
} ft_tetrahedron_fftw_plan;

void ft_destroy_tetrahedron_fftw_plan(ft_tetrahedron_fftw_plan * P);

ft_tetrahedron_fftw_plan * ft_plan_tet_with_kind(const int N, const int L, const int M, const fftw_r2r_kind kind0, const fftw_r2r_kind kind1, const fftw_r2r_kind kind2);
ft_tetrahedron_fftw_plan * ft_plan_tet_synthesis(const int N, const int L, const int M);
ft_tetrahedron_fftw_plan * ft_plan_tet_analysis(const int N, const int L, const int M);

void ft_execute_tet_synthesis(const char TRANS, const ft_tetrahedron_fftw_plan * P, double * X, const int N, const int L, const int M);
void ft_execute_tet_analysis(const char TRANS, const ft_tetrahedron_fftw_plan * P, double * X, const int N, const int L, const int M);

typedef struct {
    fftw_plan planr1;
    fftw_plan planr2;
    fftw_plan planr3;
    fftw_plan planr4;
    fftw_plan plantheta;
    double * Y;
} ft_disk_fftw_plan;

/// Destroy a \ref ft_disk_fftw_plan.
void ft_destroy_disk_fftw_plan(ft_disk_fftw_plan * P);

ft_disk_fftw_plan * ft_plan_disk_with_kind(const int N, const int M, const fftw_r2r_kind kind[3][1]);
/// Plan FFTW synthesis on the disk.
ft_disk_fftw_plan * ft_plan_disk_synthesis(const int N, const int M);
/// Plan FFTW analysis on the disk.
ft_disk_fftw_plan * ft_plan_disk_analysis(const int N, const int M);

/// Execute FFTW synthesis on the disk.
void ft_execute_disk_synthesis(const char TRANS, const ft_disk_fftw_plan * P, double * X, const int N, const int M);
/// Execute FFTW analysis on the disk.
void ft_execute_disk_analysis(const char TRANS, const ft_disk_fftw_plan * P, double * X, const int N, const int M);

typedef struct {
    fftw_plan planx1;
    fftw_plan planx2;
    fftw_plan plany;
} ft_rectdisk_fftw_plan;

/// Destroy a \ref ft_rectdisk_fftw_plan.
void ft_destroy_rectdisk_fftw_plan(ft_rectdisk_fftw_plan * P);

ft_rectdisk_fftw_plan * ft_plan_rectdisk_with_kind(const int N, const int M, const fftw_r2r_kind kind[3][1]);
/// Plan FFTW synthesis on the rectangularized disk.
ft_rectdisk_fftw_plan * ft_plan_rectdisk_synthesis(const int N, const int M);
/// Plan FFTW analysis on the rectangularized disk.
ft_rectdisk_fftw_plan * ft_plan_rectdisk_analysis(const int N, const int M);

/// Execute FFTW synthesis on the rectangularized disk.
void ft_execute_rectdisk_synthesis(const char TRANS, const ft_rectdisk_fftw_plan * P, double * X, const int N, const int M);
/// Execute FFTW analysis on the rectangularized disk.
void ft_execute_rectdisk_analysis(const char TRANS, const ft_rectdisk_fftw_plan * P, double * X, const int N, const int M);

typedef struct {
    fftw_plan plantheta1;
    fftw_plan plantheta2;
    fftw_plan plantheta3;
    fftw_plan plantheta4;
    fftw_plan planphi;
    double * Y;
    int S;
} ft_spinsphere_fftw_plan;

/// Destroy a \ref ft_spinsphere_fftw_plan.
void ft_destroy_spinsphere_fftw_plan(ft_spinsphere_fftw_plan * P);

int ft_get_spin_spinsphere_fftw_plan(const ft_spinsphere_fftw_plan * P);

ft_spinsphere_fftw_plan * ft_plan_spinsph_with_kind(const int N, const int M, const int S, const fftw_r2r_kind kind[2][1], const int sign);
/// Plan FFTW synthesis on the sphere with spin.
ft_spinsphere_fftw_plan * ft_plan_spinsph_synthesis(const int N, const int M, const int S);
/// Plan FFTW analysis on the sphere with spin.
ft_spinsphere_fftw_plan * ft_plan_spinsph_analysis(const int N, const int M, const int S);

/// Execute FFTW synthesis on the sphere with spin.
void ft_execute_spinsph_synthesis(const char TRANS, const ft_spinsphere_fftw_plan * P, ft_complex * X, const int N, const int M);
/// Execute FFTW analysis on the sphere with spin.
void ft_execute_spinsph_analysis(const char TRANS, const ft_spinsphere_fftw_plan * P, ft_complex * X, const int N, const int M);

typedef struct {
    ft_banded ** B;
    ft_triangular_banded ** T;
    int n;
} ft_gradient_plan;

void ft_destroy_gradient_plan(ft_gradient_plan * P);

ft_gradient_plan * ft_plan_sph_gradient(const int n);

void ft_execute_sph_gradient(ft_gradient_plan * P, double * U, double * Ut, double * Up, const int N, const int M);
void ft_execute_sph_curl(ft_gradient_plan * P, double * U, double * Ut, double * Up, const int N, const int M);

typedef struct {
    ft_triangular_banded ** T;
    ft_banded_qr ** F;
    double * X;
    int n;
} ft_helmholtzhodge_plan;

void ft_destroy_helmholtzhodge_plan(ft_helmholtzhodge_plan * P);

ft_helmholtzhodge_plan * ft_plan_sph_helmholtzhodge(const int n);

void ft_execute_sph_helmholtzhodge(ft_helmholtzhodge_plan * P, double * U1, double * U2, double * V1, double * V2, const int N, const int M);

/*!
  \brief A static struct to store an orthogonal matrix \f$Q \in \mathbb{R}^{3\times3}\f$, such that \f$Q^\top Q = I\f$.
  \f$Q\f$ has column-major storage.
*/
typedef struct {
    double Q[9];
} ft_orthogonal_transformation;

/*!
  \brief Every orthogonal matrix \f$Q \in \mathbb{R}^{3\times3}\f$ can be decomposed as a product of \f$ZYZ\f$ Euler angles and, if necessary, a reflection \f$R\f$ about the \f$xy\f$-plane.
  \f[
  Q = ZYZR,
  \f]
  where the \f$z\f$-axis rotations are:
  \f[
  Z = \begin{pmatrix} c & -s & 0\\ s & c & 0\\ 0 & 0 & 1\end{pmatrix},
  \f]
  the \f$y\f$-axis rotation is:
  \f[
  Y = \begin{pmatrix} c & 0 & -s\\ 0 & 1 & 0\\ s & 0 & c\end{pmatrix},
  \f]
  and the potential reflection is:
  \f[
  R = \begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & \pm 1\end{pmatrix}.
  \f]
  The reflection is stored as an integer `sign` corresponding to the bottom right entry.
*/
typedef struct {
    double s[3];
    double c[3];
    int sign;
} ft_ZYZR;

/*!
  \brief A static struct to store a reflection about the plane \f$w\cdot x = 0\f$ in \f$\mathbb{R}^3\f$.
*/
typedef struct {
    double w[3];
} ft_reflection;

ft_ZYZR ft_create_ZYZR(ft_orthogonal_transformation Q);

void ft_apply_ZYZR(ft_ZYZR Q, ft_orthogonal_transformation * U);

void ft_apply_reflection(ft_reflection Q, ft_orthogonal_transformation * U);

void ft_execute_sph_polar_rotation(double * A, const int N, const int M, double s, double c);

void ft_execute_sph_polar_reflection(double * A, const int N, const int M);

typedef struct {
    ft_symmetric_tridiagonal_symmetric_eigen * F11;
    ft_symmetric_tridiagonal_symmetric_eigen * F21;
    ft_symmetric_tridiagonal_symmetric_eigen * F12;
    ft_symmetric_tridiagonal_symmetric_eigen * F22;
    int l;
} ft_partial_sph_isometry_plan;

typedef struct {
    ft_partial_sph_isometry_plan ** F;
    int n;
} ft_sph_isometry_plan;

void ft_destroy_partial_sph_isometry_plan(ft_partial_sph_isometry_plan * F);

void ft_destroy_sph_isometry_plan(ft_sph_isometry_plan * F);

ft_partial_sph_isometry_plan * ft_plan_partial_sph_isometry(const int l);

ft_sph_isometry_plan * ft_plan_sph_isometry(const int n);

void ft_execute_sph_yz_axis_exchange(ft_sph_isometry_plan * J, double * A, const int N, const int M);

void ft_execute_sph_isometry(ft_sph_isometry_plan * J, ft_ZYZR Q, double * A, const int N, const int M);

void ft_execute_sph_rotation(ft_sph_isometry_plan * J, const double alpha, const double beta, const double gamma, double * A, const int N, const int M);

void ft_execute_sph_reflection(ft_sph_isometry_plan * J, ft_reflection W, double * A, const int N, const int M);

void ft_execute_sph_orthogonal_transformation(ft_sph_isometry_plan * J, ft_orthogonal_transformation Q, double * A, const int N, const int M);

#endif // FASTTRANSFORMS_H

max_stars_repo_name: MikaelSlevinsky/FastTransforms

max_issues_repo_name: MikaelSlevinsky/FastTransforms

max_forks_repo_name: MikaelSlevinsky/FastTransforms


text:
/*
** False disovery rate (FDR)
**
** G.Lohmann, Jan 2017
*/

#include <viaio/VImage.h>
#include <viaio/Vlib.h>
#include <viaio/mu.h>
#include <viaio/option.h>

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

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_math.h>

/* get array of FDR values */
double *GetFdr(gsl_histogram *nullhist,gsl_histogram *realhist)
{
  size_t i=0;
  size_t nbins = gsl_histogram_bins (nullhist);

  /* cumulative distribution functions (CDF) */
  gsl_histogram_pdf *cdfz = gsl_histogram_pdf_alloc(nbins);
  gsl_histogram_pdf *cdf0 = gsl_histogram_pdf_alloc(nbins);
  gsl_histogram_pdf_init (cdfz,realhist);
  gsl_histogram_pdf_init (cdf0,nullhist);

  /* fdr array */
  double *fdr = (double *)VCalloc(nbins,sizeof(double));
  for (i=0; i<nbins; i++) {
    double Fz = cdfz->sum[i];
    double F0 = cdf0->sum[i];
    double xFdr = 1.0;
    if (Fz  < 1.0) xFdr = (1.0-F0)/(1.0-Fz);
    if (xFdr > 1.0) xFdr = 1.0;
    if (xFdr < 0.0) xFdr = 0.0;
    fdr[i] = xFdr;
  }

  /* enforce monotonicity */
  fdr[0] = 1.0;
  for (i=1; i<nbins; i++) {
    if (fdr[i] > fdr[i-1]) fdr[i] = fdr[i-1];
  }
  return fdr;
}

/* thresholding */
void ApplyFdrThreshold(VImage dest,double alpha)
{
  double u=0,beta=1.0-alpha;
  size_t i,n=0;
  VFloat *pp = VImageData(dest);
  for (i=0; i<VImageNPixels(dest); i++) {
    u = (*pp);
    if (u < beta) *pp = 0;
    else n++;
    pp++;
  }
  fprintf(stderr," Number of voxels at Fdr < %.4f:  %lu\n",alpha,n);
}

/* estimage false discovery rates */
void FDR(VImage src,VImage dest,gsl_histogram *nullhist,gsl_histogram *realhist,double alpha)
{
  size_t j;
  double u=0,tiny=0.002;
  double hmin = gsl_histogram_min (realhist);
  double hmax = gsl_histogram_max (realhist);

  /* compute Fdr */
  double *Fdr = GetFdr(nullhist,realhist);

  /* apply FDR */
  VFillImage(dest,VAllBands,0);
  int b,r,c;
  for (b=0; b<VImageNBands(dest); b++) {
    for (r=0; r<VImageNRows(dest); r++) {
      for (c=0; c<VImageNColumns(dest); c++) {

	u = VGetPixel(src,b,r,c);
	if (u > hmax) u = hmax-tiny;
	if (u < hmin) u = hmin+tiny;
	if (fabs(u) > 0) {
	  if (gsl_histogram_find (realhist,(double)u,&j) == GSL_SUCCESS) {
	    VPixel(dest,b,r,c,VFloat) = 1.0 - Fdr[j];
	  }
	}
      }
    }
  }

  /* apply threshold if needed */
  if (alpha < 0.99) ApplyFdrThreshold(dest,alpha);
  return;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
// Copyright (c) 2013 Vasili Baranau
// Distributed under the MIT software license
// See the accompanying file License.txt or http://opensource.org/licenses/MIT

#ifndef Generation_PackingGenerators_Headers_BaseConjugateGradientStep_h
#define Generation_PackingGenerators_Headers_BaseConjugateGradientStep_h

#ifdef GSL_AVAILABLE

#include <gsl/gsl_multimin.h>
#include "BasePackingStep.h"
namespace PackingServices { class IPairPotential; }
namespace PackingServices { class IEnergyService; }

namespace PackingGenerators
{
    // Implements a base conjugate gradient optimization step
    class BaseConjugateGradientStep : public BasePackingStep
    {
    protected:
        PackingServices::IEnergyService* energyService;
        PackingServices::IPairPotential* pairPotential;
        std::vector<Core::SpatialVector> particleForces;

        // Minimization parameters
        gsl_multimin_function_fdf functionToMinimize;
        gsl_vector* coordinates;
        gsl_multimin_fdfminimizer* minimizer;
        const gsl_multimin_fdfminimizer_type* algorithmType;

        static const Core::FLOAT_TYPE potentialNormalizer;
        static const Core::FLOAT_TYPE initialStep;
        static const Core::FLOAT_TYPE gradientOrthogonalityTolerance;

    public:
        void SetParticles(Model::Packing* particles);

        void SetContext(const Model::ModellingContext& context);

        ~BaseConjugateGradientStep();

        // Methods for integration with GSL conjugate gradient
        void UpdateParticleCoordinates(const gsl_vector& x);

        virtual Core::FLOAT_TYPE GetPotentialEnergy();

        virtual void FillParticleForces(std::vector<Core::SpatialVector>* particleForces);

        void FillEnergyGradient(gsl_vector* g);

        void ResetGeneration();

    protected:
        BaseConjugateGradientStep(PackingServices::GeometryService* geometryService,
                PackingServices::INeighborProvider* neighborProvider,
                PackingServices::MathService* mathService,
                PackingServices::IPairPotential* pairPotential,
                PackingServices::IEnergyService* energyService);

        void InitializeMinimization();

        virtual void ResetOuterDiameterRatio() = 0;

    private:
        void InitializeMinimizationAlgorithm();

        void FreeGenerationData();

        void FillCoordinates(const Model::Packing& particles, gsl_vector* coordinates);
    };

    double GetPotentialEnergyGlobal(const gsl_vector* x, void * params);

    void FillEnergyGradientGlobal(const gsl_vector* x, void* params, gsl_vector* g);

    void FillEnergyAndGradientGlobal(const gsl_vector* x, void* params, double* f, gsl_vector* g);
}

#endif // GSL

#endif /* Generation_PackingGenerators_Headers_BaseConjugateGradientStep_h */

max_stars_repo_name: MINATILO/packing-generation

max_issues_repo_name: MINATILO/packing-generation

max_forks_repo_name: MINATILO/packing-generation


text:
/**
 * author: Jochen K"upper
 * created: Jan 2002
 * file: pygsl/src/statistics/longmodule.c
 * $Id: shortmodule.c,v 1.7 2004/03/24 08:40:45 schnizer Exp $
 *
 * "
 */

#include <Python.h>
#include <pygsl/error_helpers.h>
#include <pygsl/block_helpers.h>
#include <gsl/gsl_statistics.h>

/* include real functions for different data-types */

#define STATMOD_APPEND_PY_TYPE(X) X ## Int
#define STATMOD_APPEND_PYC_TYPE(X) X ## SHORT
#define STATMOD_FUNC_EXT(X, Y) X ## _short ## Y
#define STATMOD_PY_AS_C PyInt_AsLong
#define STATMOD_C_TYPE short int
#include "functions.c"

/* initialization */

PyGSL_STATISTICS_INIT(short, "short")

/*
 * Local Variables:
 * mode: c
 * c-file-style: "Stroustrup"
 * End:
 */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/* hka chi-square stat estimation
based on original Hudson, Kreitman, Aguade 1987 implementation.

The parameters of the model are calculated by solving the system of equations numerically.

A. D. Kern  6/2005
*/

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "string.h"
#include "hka.h"
#include "ctype.h"

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

/* currently the maximum number of loci is at 10^7, and the maximum locus name is 80 characters long */
struct hkaData data[MAXLOCI];
int locusNumber = 0;
int lineageNumber = 1;

int main(int argc, char *argv[]){
  int i = 0;
  gsl_vector  *chiSquares, *expectedPolyA, *expectedD, *devs;
  gsl_vector  *params; //0-locusNumber-1 = 4Nu's, locusNumber = tHat, locusNumber + 1 = fHat
  double x;

  getData(argc, argv);
  if(lineageNumber > 1){
    params = gsl_vector_alloc(locusNumber + 2);
  }
  else{
    params = gsl_vector_alloc(locusNumber + 1);
  }

  //solve the system of equations
  solveSystem(params);

  //calculate chiSquare Stats
  chiSquares =gsl_vector_alloc(locusNumber);
  expectedPolyA = gsl_vector_alloc(locusNumber);
  expectedD = gsl_vector_alloc(locusNumber);
  devs =  gsl_vector_alloc(locusNumber);

  x =  calculateHKA(params, chiSquares, expectedPolyA, expectedD);

  //print some stuff
  if(lineageNumber > 1){
    printf("chiSquared = %f\nt_hat = %f\nf_hat = %f\n",x, gsl_vector_get(params,locusNumber), gsl_vector_get(params,locusNumber + 1));
  }
  else{
    printf("chiSquared = %f\nt_hat = %f\nf_hat = 1.0\n",x, gsl_vector_get(params,locusNumber));
  }
    for(i = 0; i < locusNumber; i++){
      printf("locus %s \ttheta:\t %f\texpSA:\t %f\tobsSA:\t %f\texpD:\t %f\tobsD:\t %f\t chiSquared:\t%f\n",data[i].name,   gsl_vector_get(params,i),gsl_vector_get(expectedPolyA,i), data[i].sa, gsl_vector_get(expectedD,i), data[i].d, gsl_vector_get(chiSquares,i));
  }

  gsl_vector_free(params);
  return(0);
}

void usage(){
  printf("usage:\nhka infile\n");
  exit(1);
}

/* parses the data file and reads in options (currently no options boss) */
void getData(int argc, char *argv[]){
  FILE *infile;
  int  l, na, nb, sa, sb, args, flag, i,lineCount;
  double d;
  char string[81];  //here's the character length limit
  char *sArray[30], line[500];

  if (argc < 2){
    usage();
  }
  else{
    infile = fopen(argv[1],"r");
    if (infile == NULL){
      fprintf(stderr,"Error opening infile!!!\n");
      exit(1);
    }
    flag = 0;
    args = 2;
    while(args < argc){
      switch(argv[args][1]){
      case 'm' : //missing data mode
	flag = 1;
      }
      args++;
    }
    if (flag){
      //missing data mode expect- na, d,lArray, sArray
      while(fgets(line, 100, infile)){
	lineCount = chopByWhite(line,sArray,30);
	//	data[locusNumber].name = (char) sArray[0];
	data[locusNumber].na = atoi(sArray[1]);
	data[locusNumber].d = (double) atof(sArray[2]);
	for(i = 0; i < data[locusNumber].na - 2; i++){
	  data[locusNumber].lArray[i] = atoi(sArray[i+3]);
	  data[locusNumber].sArray[i] = atoi(sArray[i+data[locusNumber].na - 1 + 3]);
	}

	printf("%d\n",data[locusNumber].sArray[2]);
      }
    }
    else{

      /* collect hka data from infile- expect length, na, nb, sa, sb, d */
      while(fscanf(infile, "%80s %d %d %d %d %d %lf", string,  &l, &na, &nb, &sa, &sb, &d) != EOF){
	strcpy(data[locusNumber].name, string);
	data[locusNumber].length = l;
	data[locusNumber].na = na;
	data[locusNumber].nb = nb;
	if (nb > 1){
	  lineageNumber = 2;
	}
	data[locusNumber].sa = (double) sa;
	data[locusNumber].sb = (double) sb;
	data[locusNumber].d = (double) d;
	data[locusNumber].lArray[0] = 666;
	locusNumber++;
      }
    }
    fclose(infile);
  }
}

/* harmSum returns the denominator portion of Watterson's estimator */
double harmSum(int n){
  int i;
  double sum = 0.0;

  for(i = 1; i < n; i++){
    sum += 1.0 /  i;
  }
  return(sum);
}

/* harmSumSquare is for calculating the variance of Watterson's estimator */
double harmSumSquare(int n){
  int i;
  double sum = 0.0;

  for(i = 1; i < n; i++){
    sum += 1.0 / ( i * i);
  }
  return(sum);
}

/*setFunction- this sets up the multiroot function for solving */
int setHKAFunction(const gsl_vector *x, void *p, gsl_vector *f){
  double sumS, trueS, sumD, trueD, yi, trueYi;
  double sumSb, trueSb;
  int i;

  /*set up all the equations. the first are the sum of segSites and sum of divergence respectively */
  if (lineageNumber == 1){
    sumS = yi = sumD = trueS = trueD =  0.0;
    for(i = 0; i < locusNumber; i++){
      sumS += data[i].length * harmSum(data[i].na) * gsl_vector_get(x,i);
      sumD += data[i].length * (gsl_vector_get(x,locusNumber) + 1) * gsl_vector_get(x,i);
      trueS += data[i].sa;
      trueD += data[i].d;
    }
    gsl_vector_set(f,0,fabs(trueS - sumS));
    gsl_vector_set(f,1,fabs(trueD - sumD));
    //Sa+D eqns
    for(i = 0; i < locusNumber - 1; i++){
      trueYi = data[i].sa + data[i].d;
      yi = gsl_vector_get(x,i) *						\
	((data[i].length * gsl_vector_get(x,locusNumber) +		\
	  data[i].length + (data[i].length * harmSum(data[i].na))));
      gsl_vector_set(f,i+2,fabs(trueYi - yi));
    }
  }
  else{
    sumS = sumSb = yi = sumD = trueS = trueSb = trueD =  0.0;
    for(i = 0; i < locusNumber; i++){
      sumS += data[i].length * harmSum(data[i].na) * gsl_vector_get(x,i);
      sumSb += data[i].length * harmSum(data[i].nb) * gsl_vector_get(x,i) * gsl_vector_get(x, locusNumber + 1);
      sumD += data[i].length * (gsl_vector_get(x,locusNumber) + ((1.0 + gsl_vector_get(x,locusNumber+1)) / 2.0)) * gsl_vector_get(x,i);
      trueS += data[i].sa;
      trueSb += data[i].sb;
      trueD += data[i].d;
    }
    gsl_vector_set(f,0,fabs(trueS - sumS));
    gsl_vector_set(f,1,fabs(trueSb - sumSb));
    gsl_vector_set(f,2,fabs(trueD - sumD));
    //Sa+Sb+D eqns
    for(i = 0; i < locusNumber - 1; i++){
      trueYi = data[i].sa + data[i].sb + data[i].d;
      yi = (data[i].length * (gsl_vector_get(x,locusNumber) + ((1.0 + gsl_vector_get(x,locusNumber+1)) / 2.0)) * gsl_vector_get(x,i)) \
	+ ( data[i].length * harmSum(data[i].na) * gsl_vector_get(x,i)) + \
	(data[i].length * harmSum(data[i].nb) * gsl_vector_get(x,i) * gsl_vector_get(x, locusNumber + 1));
      gsl_vector_set(f,i+3,fabs(trueYi - yi));
    }
  }
  return GSL_SUCCESS;
}

int solveSystem(gsl_vector *params){
  const gsl_multiroot_fsolver_type *T;
  gsl_multiroot_fsolver *s;
  gsl_vector *x;
  int status;
  size_t  iter = 0;
  void *point=NULL;
  size_t n;

  if (lineageNumber == 1){
    n = locusNumber + 1;
    x = gsl_vector_alloc(locusNumber + 1);
  }
  else{
    n = locusNumber + 2;
    x = gsl_vector_alloc(locusNumber + 2);
  }

  gsl_multiroot_function f = {&setHKAFunction, n, point};
  gsl_vector_set_all(x, 1);
  T = gsl_multiroot_fsolver_hybrids;
  s = gsl_multiroot_fsolver_alloc(T, n);
  gsl_multiroot_fsolver_set(s, &f, x);

  do{
    iter++;
    status = gsl_multiroot_fsolver_iterate(s);
    if(status){
      break;
    }
    status = gsl_multiroot_test_residual(s->f, 1e-7);
  }
  while (status == GSL_CONTINUE && iter < 100);

  gsl_vector_memcpy(params, s->x);
  gsl_vector_free(x);
  gsl_multiroot_fsolver_free (s);

  return 0;
}

/* calculateHKA- this is the main beast that calculates estimates of t_hat, f(the ratio of popn sizes),
 and theta's for each locus. it returns the chi-square statistic, and takes pointers to the other
 relavent quantities (t_hat, etc.). */
double calculateHKA(gsl_vector *params,gsl_vector *chiSquares, gsl_vector *expectedPolyA, gsl_vector *expectedD){
  double  chiSquared, chiSquared_i;
  double expSA, expSB, varSA, varSB, expD, varD;
  double hsA, hsB;
  double devA, devB, devD,fHat;
  int i;

  if(lineageNumber > 1){
    fHat = gsl_vector_get(params,locusNumber+1);
  }
  else{
    fHat = 1.0;
  }
  chiSquared = 0.0;
  for(i=0; i < locusNumber; i++){
    hsA = harmSum(data[i].na);
    hsB = harmSum(data[i].nb);
    //calculate expectations and variances at the ith locus
    expSA = gsl_vector_get(params,i) * data[i].length * hsA;
    varSA = expSA + ((gsl_vector_get(params,i) * data[i].length) * (gsl_vector_get(params,i) * data[i].length) * harmSumSquare(data[i].na));
    expSB =  gsl_vector_get(params,i) * fHat * data[i].length * hsB;
    varSB = expSB + ((gsl_vector_get(params,i) * fHat * data[i].length) * (gsl_vector_get(params,i) * fHat * data[i].length) * harmSumSquare(data[i].nb)); // equals zero if no poly data from speciesB
    expD = (gsl_vector_get(params,i) * data[i].length) * (gsl_vector_get(params,locusNumber) + (0.5 * (1.0 + fHat)));
    varD = expD + (((gsl_vector_get(params,i) * data[i].length) * 0.5 * (1.0 + fHat)) * ((gsl_vector_get(params,i) *  data[i].length) * 0.5 * (1.0 + fHat)));

    //calculate deviations for species_a_theta, species_b_theta, and divergence.
    // then add to chi-square stat
    devA = ((data[i].sa - expSA) *  (data[i].sa - expSA)) / varSA;
    devB = ((data[i].sb - expSB) *  (data[i].sb - expSB)) / varSB;  //careful here will not be number if no poly data for speciesB
    devD = ((data[i].d - expD) * (data[i].d - expD)) / varD;

    //is deviation for species B finite?
    if (isnan(devB)){
      chiSquared_i = devA + devD;
    }
    else{
      chiSquared_i = devA + devB + devD;
    }
    chiSquared += chiSquared_i;
    gsl_vector_set(chiSquares,i,chiSquared_i);
    gsl_vector_set(expectedPolyA, i, expSA);
    gsl_vector_set(expectedD, i, expD);
  }

  return(chiSquared);
}

/* took this out of Jim Kent's tree */
int chopByWhite(char *in, char *outArray[], int outSize)
/* Like chopString, but specialized for white space separators. */
{
int recordCount = 0;
char c;
for (;;)
    {
    if (outArray != NULL && recordCount >= outSize)
	break;

    /* Skip initial separators. */
    while (isspace(*in)) ++in;
    if (*in == 0)
	break;

    /* Store start of word and look for end of word. */
    if (outArray != NULL)
	outArray[recordCount] = in;
    recordCount += 1;
    for (;;)
        {
        if ((c = *in) == 0)
            break;
        if (isspace(c))
            break;
        ++in;
        }
    if (*in == 0)
	break;

    /* Tag end of word with zero. */
    if (outArray != NULL)
	*in = 0;
    /* And skip over the zero. */
    in += 1;
    }
return recordCount;
}

max_stars_repo_name: kern-lab/hka

max_issues_repo_name: kern-lab/hka

max_forks_repo_name: kern-lab/hka