text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <builtin_types.h>
#include "intT.h"
#include <fenv.h>
#include <signal.h>

  //electron mass squared
  double m2 = 1.0;
  //electron charge
  double e = 0.30282212;
  //double m2=1.0;
  //double e=1.0;

double EV(double T, void * p, int* WLlist);

double Exact(double T, double B, int fermion)
//Return the exact integrand for constant fields
{
  double TB = e*T*B;
  double Exact;
  if( B < 1.0e-40)
	Exact = 0.0;
  else
    if(fermion == 1)
	Exact = exp(-m2*T)/(T*T*T)*(TB/tanh(TB)-1.0-1.0/3.0*TB*TB);
    else
	if (TB < 50.0)
		Exact = exp(-m2*T)/(T*T*T)*(TB/sinh(TB)-1.0+1.0/6.0*TB*TB);
	else
		Exact = exp(-m2*T)/(T*T*T)*(1.0/6.0*TB*TB-1.0);
  //printf("%f %e\n",T, exp(-T));
  return Exact;
}

double bumpd(double x)
//The bump function
{
  if (abs(x) < 1.0)
	return exp(-1.0/(1.0-x*x));
  else
	return 0.0;
}

double fixfluxB(double rho2, double lambda2)
//Magnetic field strength for the fixflux profile
{
  float l = sqrt(lambda2);
  float lmlmin = (l - lmin)/(tubedist - lmin);
  if(sqrt(rho2) > tubedist/2.0) printf("Warning: inappropriate magnetic field being used. rho=%f\n", sqrt(rho2));
  return  2.0/(lambda2*q2)*(1.0-0.75*lmlmin)*bumpd(2.0*sqrt(rho2/lambda2))
	+ 3.0/(tubedist*tubedist)*lmlmin;
}

double periodicB(double rho2, double lambda2)
//magnetic field for the periodic flux tube profile
{
  double expr, rho, lambda, flraf, f, ceilraf;
  double expral, exprapl;
  rho = sqrt(rho2);
  lambda = sqrt(lambda2);
  expr = exp(-rho/lambda);
  flraf = floor(rho/tubedist);
  ceilraf = ceil(rho/tubedist);
  if(rho < tubedist)
  {
	expral = 0.0;
	if((tubedist-rho)/lambda < 10.0)
		expral = exp(-(rho-tubedist)/lambda);
	f = expr/(tubedist*(1.0 + expr)*(1.0 + expr));
	if(rho > 0.0)
		f += expral/(rho*(1.0+expral)*(1.0+expral));
  }
  else
  {
        if(abs(fmod(rho,tubedist)) < 1.0e-6)
 	{
		f = flraf/(4.0*rho);
	}
	else
	{
	  expral = exp(-(rho-flraf*tubedist)/lambda);
	  exprapl = exp(-(rho-ceilraf*tubedist)/lambda);
          f = flraf*expral/(rho*(1.0+expral)*(1.0+expral)) +
		ceilraf*exprapl/(rho*(1.0+exprapl)*(1.0+exprapl));
        }
  }
  //f=1.0f/lambda2*exp(-1.0f*rho2/lambda2);
  //f=lambda2/8.0f;
  f=tubedist/(2.0*log(2.0)*lambda2)*f;
  //f=1.0;
  return f;
}

double MagField(double rho2, void* p)
//Compute the magnetic field.
{
  double B;
  double sinrho;
  struct Wparams *params = (struct Wparams *) p;
  switch(profile){
	case step:
	  if(rho2<params->l2)
		B=1.0/params->l2;
	  else
		B=0.0;
	  break;
	case smooth:
	  B=params->l2/((params->l2+rho2)*(params->l2+rho2));
	  break;
	case quadratic:
	  if(rho2<params->l2)
		B=2.0/params->l2*(1.0-rho2/params->l2);
	  else
		B=0.0;
	  break;
	case gaussian:
	  if(rho2<100.0*params->l2)
	  	B=1.0/(params->l2*exp(rho2/params->l2));
	  else
		B=0.0;
	  break;
	case periodic:
	  B=periodicB(rho2,params->l2);
	  break;
	case spline:
	  sinrho=sin(pi*sqrt(rho2)/tubedist);
	  B=1.0/(params->l2*exp(sinrho*sinrho/params->l2));
	  break;
	case fixflux:
	  B = fixfluxB(rho2,params->l2);
	  break;
  }
  B=2.0*params->F/e*B;

  return B;
}

double smTscal(double B, double T)
//Effective action for T<<1.0 for scalar QED
{
	double B2=e*e*B*B;
	double B4=B2*B2;
	double B6=B2*B4;
	double T2=T*T;
	double T3=T2*T;
	double T4=T3*T;
	double result;
	result = 7.0/360.0*B4*T*(1.0-m2*T)+(147.0*B4*m2*m2-31.0*B6)/15120.0*T3
		+(31.0*B6*m2-49.0*B4*m2*m2*m2)/15120.0*T4;
	return result;
}

double smTferm(double B, double T)
//Effective action for T<<1.0 for fermionic QED
{
	double B2=e*e*B*B;
	double B4=B2*B2;
	double B6=B2*B4;
	double T2=T*T;
	double T3=T2*T;
	double T4=T3*T;
	double result;
	result = 1.0/45.0*B4*T*(m2*T-1.0) + (2.0*B6/945.0 - B4*m2*m2/90.0)*T3
		+(7.0*B4*m2*m2*m2-4.0*B6*m2)/1890.0*T4;
	//result = 1.0/45.0*B4*T*(T-1.0) + (2.0*B6/945.0 - B4/90.0)*T3
	//	+(7.0*B4-4.0*B6)/1890.0*T4;
	return result;
}

double meanigrand(double* igrand, int ngroups, double* igranderr)
//compute the mean and std. err. for igrand
{
  double meanig=0.0;
  int i;
  *igranderr=0.0;
  for(i=0;i<ngroups;i++)
  {
	meanig+=igrand[i];
  }
  meanig/=ngroups;
  for(i=0; i < ngroups; i++)
  {
	*igranderr+=(igrand[i]-meanig)*(igrand[i]-meanig);
  }
  *igranderr = sqrt(*igranderr);
  *igranderr /= ngroups;
  return meanig;
}

double* Tfunc(double T, void* p, int* order)
//The integrand of the proper time, T, integral
{
  struct Wparams *params = (struct Wparams *) p;
  double rho2 = params->xcm.x*params->xcm.x + params->xcm.y*params->xcm.y;
  double smallT;
  double B;
  double TB;
  double* igrand = malloc(params->ng*sizeof(*igrand));
  int groupsize = 128; //loops per group
  int* WLlist = malloc(groupsize*sizeof(*WLlist));
  int i, j;
  double rho = params->xcm.x;
  double l2 = params->l2;
  double sin2rho = sin(2.0*pi*rho/tubedist);
  double denominator;
  double igmean; double igerr;
  if(!igrand || !WLlist) printf("error Tfunc():  malloc failed\n");

  //printf("func: T=%f\n",T);
  B=MagField(rho2,p);
  TB=e*T*B;
  //printf("T=%f f'/f''=%f\n",T,(params->l2+params->xcm.x*params->xcm.x));
  //determine a suitable definition for small T
  switch(profile)
  {
	case step:
	  smallT=abs(rho*rho-l2);
          break;
        case smooth:
	  smallT= 0.5*(l2+rho*rho);
          break;
        case quadratic:
	  smallT= 0.5*(l2+rho*rho);
          break;
        case gaussian:
	  smallT=l2;
	  break;
        case periodic:
	  smallT=0.5*(l2+rho*rho);
          break;
	case spline:
	  if(rho/tubedist>1.0e-8)
	  	smallT=tubedist*tubedist*l2/(pi*pi*sin2rho*sin2rho);
	  else
		smallT=1.0e12;
	  break;
	case fixflux:
	  //small if T << l2, or if the loop is far from any bump functions.
	  if(rho*rho<l2/4.0)
		smallT = 0.25*l2;
	  else
		smallT = 0.25*l2+(rho-sqrt(l2)/2.0)*(rho-sqrt(l2)/2.0);
  }
  for(i = 0; i < params->ng; i++)
  {
	//printf("func igrand[%d]=%f\n",i,igrand[i]);
	if(T < 0.1 && T < 0.01*smallT)
	//if(T<0.7)
  	{
		if(verbosity >= 2)
		  printf("using small T expression: %f %f %f\n", rho2, T, smallT);
		if(params->fermion == 1)
			igrand[i] = smTferm(B, T);
		else
			igrand[i] = smTscal(B, T);
  	}
        else if(T < 0.01*smallT)
	//use the exact expression
	{
		if(verbosity >= 2)
		  printf("using const. field expression %f %f %f\n", rho2, T, smallT);
		if(TB > 1.0e-6)
			igrand[i] = Exact(T, B, params->fermion);
		else
			igrand[i] = 0.0;
	}
	else if(T > 49.0)
	//for large T, the contribution is exponentially supressed.
	{
		if(verbosity >= 2)
                  printf("using large T expression %f %f %f\n", rho2, T);
		if(params->fermion==1)
		{
			igrand[i]=-1.0*exp(-m2*T)/(3.0*T)*e*e*B*B;
			//if(B>1e-6)
			//	igrand[i]=exp(-m2*T)/(T*T*T)*(TB/tanh(TB)-1.0-1.0/3.0*TB*TB);
			//else
			//	igrand[i]=0.0;
			//printf("igrand(%f)=%f\n",T,igrand[i]);
		}
		else
			igrand[i]=exp(-m2*T)/(6.0*T)*e*e*B*B;
	}
  	else
  	{
		for(j=0;j<groupsize;j++)
        	{
			//printf("func j=%d\n",j);
			WLlist[j]=order[j + i*groupsize];
        	}
		if(verbosity >= 3)
		  printf("EV= %e Exact = %e\n",EV(T,p, WLlist),T/sinh(T));
		if(params->fermion==1)
			igrand[i]= exp(-m2*T)/(T*T*T)*(EV(T, p, WLlist)-1.0-1.0/3.0*TB*TB);
		else
			igrand[i]= exp(-m2*T)/(T*T*T)*(EV(T, p, WLlist)-1.0+1.0/6.0*TB*TB);
		//igrand[i]=Exact(T,1.0);
		//return Exact(T,1.0);

  	}

  }
  if(T > 0.0f && verbosity >= 2){
	igmean = meanigrand(igrand, params->ng, &igerr);
	printf("WLvconstExpression: ");
  		printf("%e %e %e %e %e %e\n",
  		  rho2, T, Exact(T, B, params->fermion), smTscal(B, T), igmean, igerr);
  }
  free(WLlist);
  return igrand;
}

void signal_handler(int sig)
//Floating point signal handler
{
  printf( "Error: Floating point signal detected.\n");
  exit(1);
}

void EVMean(double *EV, float4 *Wsscal_h, float4 *Wsferm_h, int n, int *WL, double T, int fermion)
//Compute the expectation value and error of the Wilson loops
//WL is a list representing a group of wordline indices
//to compute the expectation value of
{
  int i, WLi;
  double EVWlNb2, EVWlN, EVWlNb4;
  double SEWlNb2, SEWlN, SEWlNb4;
  double *EVpart1;
  double *EVpart2;
  double *EVpart4;
  double Nby4part;

  //Turn on floating point exceptions:
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW);

  // Set the signal handler:
  signal(SIGFPE, signal_handler);

  *EV = 0.0;
  EVWlNb2 = 0.0; EVWlN = 0.0; EVWlNb4 = 0.0;
  //allocate arrays to store the wilson loop of each worldline
  EVpart1 = (double *) malloc(n*sizeof(*EVpart1));
  EVpart2 = (double *) malloc(n*sizeof(*EVpart2));
  EVpart4 = (double *) malloc(n*sizeof(*EVpart4));
  if(!EVpart1 || !EVpart2 || !EVpart4) printf("EVonly(): malloc failed\n");
  //determine the sum of the contributions
  for(i = 0; i < n; i++){
	WLi=WL[i];
	//Compute the scalar part
	Nby4part = Wsscal_h[WLi].y + Wsscal_h[WLi].z;
	EVpart1[i] = cos(0.5*Wsscal_h[WLi].x+0.25*(Nby4part));
	EVpart2[i] = cos(Wsscal_h[WLi].x);
	EVpart4[i] = cos(0.5*(Nby4part));
	if(fermion == 1)
	{
		if(isinf(Wsferm_h[WLi].x) !=0 ||isinf(Wsferm_h[WLi].y) !=0 ||isinf(Wsferm_h[WLi].z) !=0
		  || isnan(Wsferm_h[WLi].x) !=0 || isnan(Wsferm_h[WLi].y) !=0 ||isnan(Wsferm_h[WLi].z) !=0 )
		{
			printf("Warning: WSferm is infinite. Worldline=%d\n",WLi);
		}
		//Compute the fermion part
		//printf("Wsferm_h[%d].x=%f\n",WLi,Wsferm_h[WLi].x);
		//printf("Wsferm_h[%d].y=%f\n",WLi,Wsferm_h[WLi].y);
		//printf("Wsferm_h[%d].z=%f\n",WLi,Wsferm_h[WLi].z);
		Nby4part = Wsferm_h[WLi].y+Wsferm_h[WLi].z;
		if(abs(Nby4part) > 600.0 | abs(Wsferm_h[WLi].x) > 600.0)
		{
			printf("Warning: Large cosh argument.  T=%f, WLi=%d\n",T,WLi);
			EVpart1[i] *= 1.0e10;
			EVpart2[i] *= 1.0e10;
			EVpart4[i] *= 1.0e10;
		}
                else
		{
			EVpart1[i] *= cosh(0.5*(double)Wsferm_h[WLi].x+0.25*((double)Nby4part));
			EVpart2[i] *= cosh((double)Wsferm_h[WLi].x);
			EVpart4[i] *= cosh(0.5*((double)Nby4part));
			//EVpart1[i]*=cosh(T);
			//EVpart2[i]*=cosh(T);
			//EVpart4[i]*=cosh(T);
		}
		if(isnan(EVpart1[i]) != 0 || isnan(EVpart2[i]) != 0 || isnan(EVpart4[i]) != 0
			|| isinf(EVpart1[i]) != 0 || isinf(EVpart2[i]) != 0 || isinf(EVpart4[i]) != 0)
			printf("Warning: Infinity detected: WL# %d\n",WLi);
	}
	EVWlNb4 += EVpart4[i];
	EVWlNb2 += EVpart2[i];
	EVWlN   += EVpart1[i];
        //printf("EVpart1[%d]=%f\n",i,EVpart1[i]);
  }
  EVWlNb4 /=n;
  EVWlNb2 /=n;
  EVWlN  /=n;
  //Extrapolate to N=infinity
  *EV=8.0/3.0*EVWlN - 2.0*EVWlNb2 + 1.0/3.0*EVWlNb4;

  free(EVpart1);
  free(EVpart2);
  free(EVpart4);

}

max_stars_repo_name: QEDan/Worldline-Flux-Tube-Lattice

max_issues_repo_name: QEDan/Worldline-Flux-Tube-Lattice

max_forks_repo_name: QEDan/Worldline-Flux-Tube-Lattice


text:
#ifndef CTETRA_NUMSTATES_H
#define CTETRA_NUMSTATES_H

#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include "input.h"
#include "submesh.h"
#include "ecache.h"
#include "tetra.h"

double NumStates(double E, EnergyCache *Ecache);

double NumStatesContrib(double E, double E1, double E2, double E3, double E4, double num_tetra);

#endif // CTETRA_NUMSTATES_H

ext: h

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra


text:
#ifndef OPENMC_TALLIES_FILTER_ENERGY_H
#define OPENMC_TALLIES_FILTER_ENERGY_H

#include <gsl/gsl>

#include "openmc/tallies/filter.h"
#include "openmc/vector.h"

namespace openmc {

//==============================================================================
//! Bins the incident neutron energy.
//==============================================================================

class EnergyFilter : public Filter {
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors

  ~EnergyFilter() = default;

  //----------------------------------------------------------------------------
  // Methods

  std::string type() const override { return "energy"; }

  void from_xml(pugi::xml_node node) override;

  void get_all_bins(const Particle& p, TallyEstimator estimator,
    FilterMatch& match) const override;

  void to_statepoint(hid_t filter_group) const override;

  std::string text_label(int bin) const override;

  //----------------------------------------------------------------------------
  // Accessors

  const vector<double>& bins() const { return bins_; }
  void set_bins(gsl::span<const double> bins);

  bool matches_transport_groups() const { return matches_transport_groups_; }

protected:
  //----------------------------------------------------------------------------
  // Data members

  vector<double> bins_;

  //! True if transport group number can be used directly to get bin number
  bool matches_transport_groups_ {false};
};

//==============================================================================
//! Bins the outgoing neutron energy.
//!
//! Only scattering events use the get_all_bins functionality.  Nu-fission
//! tallies manually iterate over the filter bins.
//==============================================================================

class EnergyoutFilter : public EnergyFilter {
public:
  //----------------------------------------------------------------------------
  // Methods

  std::string type() const override { return "energyout"; }

  void get_all_bins(const Particle& p, TallyEstimator estimator,
    FilterMatch& match) const override;

  std::string text_label(int bin) const override;
};

} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_ENERGY_H


text:
///
/// @file
///
/// @author Angelika Schwarz (angies@cs.umu.se), Umeå University
/// @author Mirko Myllykoski (mirkom@cs.umu.se), Umeå University
///
/// @internal LICENSE
///
/// Copyright (c) 2019-2020, Umeå Universitet
///
/// Redistribution and use in source and binary forms, with or without
/// modification, are permitted provided that the following conditions are met:
///
/// 1. Redistributions of source code must retain the above copyright notice,
///    this list of conditions and the following disclaimer.
///
/// 2. Redistributions in binary form must reproduce the above copyright notice,
///    this list of conditions and the following disclaimer in the documentation
///    and/or other materials provided with the distribution.
///
/// 3. Neither the name of the copyright holder nor the names of its
///    contributors may be used to endorse or promote products derived from this
///    software without specific prior written permission.
///
/// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
/// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
/// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
/// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
/// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
/// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
/// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
/// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
/// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
/// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
/// POSSIBILITY OF SUCH DAMAGE.
///

#include <starneig_config.h>
#include <starneig/configuration.h>
#include <starneig/error.h>
#include "core.h"
#include "typedefs.h"
#include "robust.h"
#include "cpu.h"
#include "partition.h"
#include "../../common/common.h"
#include "../../common/node_internal.h"
#include "../../common/matrix.h"
#include <starneig/sep_sm.h>
#include <cblas.h>
#include <stdlib.h>
#include <starpu.h>
#include <math.h>
#include <float.h>

static int search_multiplicities(int n, double *lambda, int *lambda_type)
{
    int multiplicity = 0;
    for (int i = 0; i < n; i++) {
        if (lambda_type[i] == 0) { // REAL
            double real = lambda[i];
            for (int j = i+1; j < n; j++) {
                if (lambda_type[j] == 0 && lambda[j] == real) {
                    multiplicity = 1;
                    break;
                }
            }
        }
        else { // CMPLX
            double real = lambda[i];
            double imag = lambda[i+1];
            for (int j = i+2; j < n; j++) {
                if (lambda_type[j] == 1 && lambda_type[j+1]) {
                    if (lambda[j] == real && lambda[j+1] == imag) {
                        multiplicity = 1;
                        break;
                    }

                    j++;
                }
            }
            i++;
        }
    }

    return multiplicity;
}

static starneig_error_t eigenvectors(
    struct starneig_eigenvectors_conf const *_conf,
    int n, int *selected,
    double *S, int ldS,
    double *Q, int ldQ,
    double *Y, int ldY)
{
#define S(i,j) S[(i) + (j) * (size_t)ldS]
#define Q(i,j) Q[(i) + (j) * (size_t)ldQ]
#define X(i,j) X[(i) + (j) * (size_t)ldX]
#define Y(i,j) Y[(i) + (j) * (size_t)ldY]

    // use default configuration if necessary
    struct starneig_eigenvectors_conf *conf;
    struct starneig_eigenvectors_conf local_conf;
    if (_conf == NULL)
        starneig_eigenvectors_init_conf(&local_conf);
    else
        local_conf = *_conf;
    conf = &local_conf;

    //
    // check mandatory arguments
    //

    if (selected == NULL) {
        starneig_error("Eigenvalue selection bitmap is NULL. Exiting...");
        return STARNEIG_INVALID_ARGUMENTS;
    }

    int num_selected = starneig_eigvec_std_count_selected(n, selected);
    if (num_selected == 0) {
        starneig_error("Eigenvalue selection bitmap does not have any "
                       "selected eigenvalues. Exiting...");
        return STARNEIG_INVALID_ARGUMENTS;
    }

    if (S == NULL) {
        starneig_error("Matrix S is NULL. Exiting...");
        return STARNEIG_INVALID_ARGUMENTS;
    }

    if (Q == NULL) {
        starneig_error("Matrix Q is NULL. Exiting...");
        return STARNEIG_INVALID_ARGUMENTS;
    }

    if (Y == NULL) {
        starneig_error("Eigenvector matrix is NULL. Exiting...");
        return STARNEIG_INVALID_ARGUMENTS;
    }

    //
    // check configuration
    //

    if (conf->tile_size == STARNEIG_EIGENVECTORS_DEFAULT_TILE_SIZE) {
        double select_ratio = (double) num_selected/n;
        conf->tile_size = MIN(MAX(240, sqrt(n)/sqrt(select_ratio)), 936);
        starneig_message("Setting tile size to %d.", conf->tile_size);
    }

    if (conf->tile_size <= 0) {
        starneig_error("Tile size is %d. Exiting...", conf->tile_size);
        return STARNEIG_INVALID_CONFIGURATION;
    }

    //
    // preprocess
    //

    // extract the eigenvalues and their type (1-by-1 or 2-by-2 block)
    double *lambda = (double *) malloc((size_t)n * sizeof(double));
    int *lambda_type = (int *) malloc((size_t)n * sizeof(int));
    lambda[n-1] = S(n-1,n-1);
    lambda_type[n-1] = 0;
    for (int i = 0; i < n - 1; i++) {
        if (S(i+1,i) != 0.0) {
            // A 2-by-2 block in canonical Schur form has the shape
            // [ S(i,i)      S(i,i+1)   ] = [ a b ]  or [ a -b ]
            // [ S(i+1,i)    S(i+1,i+1) ]   [-c a ]     [ c  a ].
            lambda_type[i] = 1;
            lambda_type[i+1] = 1;
            double real = S(i+1,i+1);
            double imag = sqrt(fabs(S(i+1,i)))*sqrt(fabs(S(i,i+1)));
            lambda[i] = real;
            lambda[i+1] = imag;
            i++;
        }
        else {
            // 1-by-1 block
            lambda_type[i] = 0;
            lambda[i] = S(i,i);
        }
    }

    //
    // sanity check
    //

    starneig_error_t ret = STARNEIG_SUCCESS;
    int illposed = search_multiplicities(n, lambda, lambda_type);
    if (illposed) {
        starneig_warning("Multiple eigenvalues detected.\n");
        ret = STARNEIG_CLOSE_EIGENVALUES;
    }

    //
    // overflow control
    //

    const double eps = DBL_EPSILON/2;
    const double smlnum = MAX(2*DBL_MIN, DBL_MIN*((double)n/eps));

    //
    // partition
    //

    int num_tiles = (n+conf->tile_size-1)/conf->tile_size;
    int *first_row = (int *) malloc((num_tiles+1)*sizeof(int));
    int *first_col = (int *) malloc((num_tiles+1)*sizeof(int));
    starneig_eigvec_std_partition(n, lambda_type, conf->tile_size, first_row);
    starneig_eigvec_std_partition_selected(n, first_row, selected, num_tiles, first_col);

    //
    // workspace
    //

    int ldX = ldY;
    double *X = (double *) malloc((size_t)ldX*num_selected*sizeof(double));

    size_t num_segments = (size_t) num_tiles*num_selected;

    double *Xnorms = (double *) malloc(num_segments*sizeof(double));
#define Xnorms(col, tilerow) Xnorms[(col) + (tilerow) * (size_t)num_selected]

    scaling_t *scales = (scaling_t*) malloc(num_segments*sizeof(scaling_t));
#define scales(col, tilerow) scales[(col) + (tilerow) * (size_t)num_selected]
    starneig_eigvec_std_init_scaling_factor(num_tiles*num_selected, scales);

    double *Snorms =
        (double *) malloc((size_t)num_tiles*num_tiles*sizeof(double));
#define Snorms(i,j) Snorms[(i) + (j) * (size_t)num_tiles]

    int *info = (int *) malloc((size_t)num_selected*sizeof(int));

    // Copy all selected eigenvalue types to a compact memory representation.
    int *selected_lambda_type = (int *) malloc((size_t)num_selected*sizeof(int));
    int idx = 0;
    for (int i = 0; i < n; i++) {
        if (selected[i]) {
            selected_lambda_type[idx] = lambda_type[i];
            idx++;
        }
    }

    //
    // register
    //

    starpu_data_handle_t **S_tiles;
    starpu_data_handle_t **Q_tiles;
    starpu_data_handle_t **X_tiles;
    starpu_data_handle_t **Y_tiles;
    starpu_data_handle_t *selected_tiles;
    starpu_data_handle_t *lambda_tiles;
    starpu_data_handle_t *lambda_type_tiles;
    starpu_data_handle_t *selected_lambda_type_tiles;
    starpu_data_handle_t **Xnorms_tiles;
    starpu_data_handle_t **scales_tiles;
    starpu_data_handle_t **S_tiles_norms;
    starpu_data_handle_t *info_tiles;
    S_tiles = malloc(num_tiles*sizeof(starpu_data_handle_t *));
    Q_tiles = malloc(num_tiles*sizeof(starpu_data_handle_t *));
    X_tiles = malloc(num_tiles*sizeof(starpu_data_handle_t *));
    Y_tiles = malloc(num_tiles*sizeof(starpu_data_handle_t *));
    selected_tiles = malloc(num_tiles*sizeof(starpu_data_handle_t));
    lambda_tiles = malloc(num_tiles*sizeof(starpu_data_handle_t));
    lambda_type_tiles = malloc(num_tiles*sizeof(starpu_data_handle_t));
    selected_lambda_type_tiles = malloc(num_tiles*sizeof(starpu_data_handle_t));
    Xnorms_tiles = malloc(num_tiles*sizeof(starpu_data_handle_t *));
    scales_tiles = malloc(num_tiles*sizeof(starpu_data_handle_t *));
    S_tiles_norms = malloc(num_tiles*sizeof(starpu_data_handle_t *));
    info_tiles = malloc(num_tiles*sizeof(starpu_data_handle_t));
    for (int i = 0; i < num_tiles; i++) {
        S_tiles[i] = malloc(num_tiles*sizeof(starpu_data_handle_t));
        Q_tiles[i] = malloc(num_tiles*sizeof(starpu_data_handle_t));
        X_tiles[i] = malloc(num_tiles*sizeof(starpu_data_handle_t));
        Y_tiles[i] = malloc(num_tiles*sizeof(starpu_data_handle_t));
        Xnorms_tiles[i] = malloc(num_tiles*sizeof(starpu_data_handle_t));
        scales_tiles[i] = malloc(num_tiles*sizeof(starpu_data_handle_t));
        S_tiles_norms[i] = malloc(num_tiles*sizeof(starpu_data_handle_t));

        starpu_vector_data_register(
            &selected_tiles[i],
            STARPU_MAIN_RAM,
            (uintptr_t)(&selected[first_row[i]]),
            first_row[i+1]-first_row[i],
            sizeof(int));

        starpu_vector_data_register(
            &lambda_tiles[i],
            STARPU_MAIN_RAM,
            (uintptr_t)(&lambda[first_row[i]]),
            first_row[i+1]-first_row[i],
            sizeof(double));

        starpu_vector_data_register(
            &lambda_type_tiles[i],
            STARPU_MAIN_RAM,
            (uintptr_t)(&lambda_type[first_row[i]]),
            first_row[i+1]-first_row[i],
            sizeof(int));

        starpu_vector_data_register(
            &selected_lambda_type_tiles[i],
            STARPU_MAIN_RAM,
            (uintptr_t)(&selected_lambda_type[first_col[i]]),
            first_col[i+1]-first_col[i],
            sizeof(int));

        starpu_vector_data_register(
            &info_tiles[i],
            STARPU_MAIN_RAM,
            (uintptr_t)(&info[first_col[i]]),
            first_col[i+1]-first_col[i],
            sizeof(int));

        for (int j = 0; j < num_tiles; j++) {
            if (i <= j) {
                starpu_matrix_data_register(
                    &S_tiles[i][j],
                    STARPU_MAIN_RAM,
                    (uintptr_t)(&S(first_row[i], first_row[j])),
                    ldS,
                    first_row[i+1]-first_row[i],
                    first_row[j+1]-first_row[j],
                    sizeof(double));

                starpu_variable_data_register(
                    &S_tiles_norms[i][j],
                    STARPU_MAIN_RAM,
                    (uintptr_t)(&Snorms(i,j)),
                    sizeof(double));

                starpu_matrix_data_register(
                    &X_tiles[i][j],
                    STARPU_MAIN_RAM,
                    (uintptr_t)(&X(first_row[i], first_col[j])),
                    ldX,
                    first_row[i+1]-first_row[i],
                    first_col[j+1]-first_col[j],
                    sizeof(double));

                starpu_vector_data_register(
                    &Xnorms_tiles[i][j],
                    STARPU_MAIN_RAM,
                    (uintptr_t)(&Xnorms(first_col[j],i)),
                    first_col[j+1]-first_col[j],
                    sizeof(double));

                starpu_vector_data_register(
                    &scales_tiles[i][j],
                    STARPU_MAIN_RAM,
                    (uintptr_t)(&scales(first_col[j],i)),
                    first_col[j+1]-first_col[j],
                    sizeof(scaling_t));
            }

            starpu_matrix_data_register(
                &Q_tiles[i][j],
                STARPU_MAIN_RAM,
                (uintptr_t)(&Q(first_row[i], first_row[j])),
                ldQ,
                first_row[i+1]-first_row[i],
                first_row[j+1]-first_row[j],
                sizeof(double));

            starpu_matrix_data_register(
                &Y_tiles[i][j],
                STARPU_MAIN_RAM,
                (uintptr_t)(&Y(first_row[i], first_col[j])),
                ldY,
                first_row[i+1]-first_row[i],
                first_col[j+1]-first_col[j],
                sizeof(double));
        }
    }

    //
    // insert tasks
    //

    starneig_eigvec_std_insert_backsolve_tasks(num_tiles,
        S_tiles, S_tiles_norms, lambda_tiles, lambda_type_tiles,
        X_tiles, scales_tiles, Xnorms_tiles, selected_tiles,
        selected_lambda_type_tiles, info_tiles, smlnum,
        STARPU_MAX_PRIO, STARPU_DEFAULT_PRIO);

    starpu_task_wait_for_all();

    starneig_eigvec_std_unify_scaling(num_tiles, first_row, first_col, scales, X, ldX,
        lambda_type, selected);

    starneig_eigvec_std_insert_backtransform_tasks(first_row, num_tiles,
        Q_tiles, X_tiles, Y_tiles);

    //
    // evaluate reliability
    //

    for (int i = 0; i < num_selected; i++) {
        if (info[i] != STARNEIG_SUCCESS) {
            starneig_warning("Eigenvector column X(:,%d) was perturbed and "
                             "cannot be trusted.", i);
            ret = STARNEIG_CLOSE_EIGENVALUES;
        }
    }

    starpu_task_wait_for_all();

    //
    // clean up
    //

    for (int i = 0; i < num_tiles; i++) {
        starpu_data_unregister(selected_tiles[i]);
        starpu_data_unregister(lambda_tiles[i]);
        starpu_data_unregister(lambda_type_tiles[i]);
        starpu_data_unregister(selected_lambda_type_tiles[i]);
        starpu_data_unregister(info_tiles[i]);
        for (int j = 0; j < num_tiles; j++) {
            if (i <= j) {
                starpu_data_unregister(S_tiles[i][j]);
                starpu_data_unregister(S_tiles_norms[i][j]);
                starpu_data_unregister(X_tiles[i][j]);
                starpu_data_unregister(Xnorms_tiles[i][j]);
                starpu_data_unregister(scales_tiles[i][j]);
            }
            starpu_data_unregister(Q_tiles[i][j]);
            starpu_data_unregister(Y_tiles[i][j]);
        }
        free(S_tiles[i]);
        free(Q_tiles[i]);
        free(X_tiles[i]);
        free(Y_tiles[i]);
        free(Xnorms_tiles[i]);
        free(scales_tiles[i]);
        free(S_tiles_norms[i]);
    }
    free(S_tiles);
    free(Q_tiles);
    free(X_tiles);
    free(Y_tiles);
    free(Xnorms_tiles);
    free(scales_tiles);
    free(selected_tiles);
    free(lambda_tiles);
    free(lambda_type_tiles);
    free(selected_lambda_type_tiles);
    free(info_tiles);
    free(X);
    free(Xnorms);
    free(Snorms);
    free(info);
    free(scales);
    free(lambda_type);
    free(selected_lambda_type);
    free(lambda);
    free(first_row);
    free(first_col);

#undef Xnorms
#undef Snorms
#undef scales
#undef S
#undef Q
#undef X
#undef Y

    return ret;
}

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

__attribute__ ((visibility ("default")))
void starneig_eigenvectors_init_conf(struct starneig_eigenvectors_conf *conf) {
    conf->tile_size = STARNEIG_EIGENVECTORS_DEFAULT_TILE_SIZE;
}

__attribute__ ((visibility ("default")))
int starneig_SEP_SM_Eigenvectors_expert(
    struct starneig_eigenvectors_conf *conf,
    int n,
    int selected[],
    double S[], int ldS,
    double Q[], int ldQ,
    double X[], int ldX)
{
    CHECK_INIT();

    starneig_node_set_blas_mode(STARNEIG_BLAS_MODE_SEQUENTIAL);
    starneig_node_set_mode(STARNEIG_MODE_SM);
    starneig_node_resume_starpu();

    starneig_error_t ret = eigenvectors(
        conf, n, selected, S, ldS, Q, ldQ, X, ldX);

    starpu_task_wait_for_all();
    starneig_node_pause_starpu();
    starneig_node_set_mode(STARNEIG_MODE_OFF);
    starneig_node_set_blas_mode(STARNEIG_BLAS_MODE_ORIGINAL);

    return ret;
}

__attribute__ ((visibility ("default")))
int starneig_SEP_SM_Eigenvectors(
    int n,
    int selected[],
    double S[], int ldS,
    double Q[], int ldQ,
    double X[], int ldX)
{
    CHECK_INIT();
    return starneig_SEP_SM_Eigenvectors_expert(
        NULL, n, selected, S, ldS, Q, ldQ, X, ldX);
}

ext: c

max_stars_repo_name: NLAFET/StarNEig

max_issues_repo_name: NLAFET/StarNEig

max_forks_repo_name: NLAFET/StarNEig


text:
/*****************************************************************\
           __
          / /
		 / /                     __  __
		/ /______    _______    / / / / ________   __       __
	   / ______  \  /_____  \  / / / / / _____  | / /      / /
	  / /      | / _______| / / / / / / /____/ / / /      / /
	 / /      / / / _____  / / / / / / _______/ / /      / /
	/ /      / / / /____/ / / / / / / |______  / |______/ /
   /_/      /_/ |________/ / / / /  \_______/  \_______  /
                          /_/ /_/                     / /
			                                         / /
		       High Level Game Framework            /_/

  ---------------------------------------------------------------

  Copyright (c) 2007-2011 - Rodrigo Braz Monteiro.
  This file is subject to the terms of halley_license.txt.

\*****************************************************************/

#pragma once

#include <string>
#include <sstream>
#include <halley/data_structures/vector.h>
#include <gsl/gsl_assert>
#include <iomanip>
#include <cstdint>

namespace Halley {

	typedef char Character;
	typedef wchar_t utf16type;
	typedef char32_t utf32type;
	typedef std::wstring StringUTF16;
	typedef std::u32string StringUTF32;

	// String class
	class String {
	public:
		const static size_t npos = size_t(-1);

		String();
		String(const char* utf8);
		String(const char* utf8,size_t bytes);
		String(const std::basic_string<Character>& str);
		String(const String& str) noexcept;
		String(String&& str) noexcept;

		explicit String(const wchar_t* utf16);
		explicit String(const StringUTF32 &utf32);
		explicit String(char character);
		explicit String(wchar_t character);
		explicit String(int character);
		explicit String(float number);
		explicit String(double number);

		String& operator=(const char* utf8);
		String& operator=(const std::basic_string<Character>& str);
		String& operator=(String&& str) noexcept;
		String& operator=(const String& str);

		operator std::string() const;

		bool isEmpty() const;
		size_t length() const;

		void setSize(size_t size);
		void truncate(size_t size);

		String& trim(bool fromRight);
		String& trimBoth();

		bool contains(const String& string) const;
		size_t find(String str) const;

		String replaceAll(const String& before, const String& after) const;
		String replaceOne(const String& before, const String& after) const;
		void shrink();

		String left(size_t n) const;
		String right(size_t n) const;
		String mid(size_t start,size_t count=npos) const;

		bool startsWith(const String& string,bool caseSensitive=true) const;
		bool endsWith(const String& string,bool caseSensitive=true) const;

		void writeText(const Character* src,size_t len,size_t &pos);
		void writeChar(const Character &src,size_t &pos);
		void writeNumber(Character *temp,int number,int pad,size_t &pos);

		bool isNumber() const;
		bool isInteger() const;

		String asciiLower() const;
		String asciiUpper() const;
		void asciiMakeUpper();
		void asciiMakeLower();
		bool asciiCompareNoCase(const Character *src) const;

		void appendCharacter(int unicode);

		// Convert a string to a number
		int toInteger() const;
		long long toInteger64() const;
		float toFloat() const;
		int subToInteger(size_t start,size_t end) const;

		// std::string methods
		const char* c_str() const;
		String substr(size_t pos, size_t len=npos) const;
		size_t find(Character character, size_t pos=0) const;
		size_t find(const char* str, size_t pos=0) const;
		size_t find_last_of(char character) const;
		size_t size() const;
		const char& operator[](size_t pos) const;
		char& operator[](size_t pos);

		static const Character* stringPtrTrim(Character *chr,size_t len,size_t startPos);
		static const Character* stringTrim(String &str,size_t startPos);

		// Number tidy up functions
		static String prettyFloat(String src);
		static String prettySize(long long bytes);

		// Unicode routines
		StringUTF16 getUTF16() const;
		StringUTF32 getUTF32() const;
		size_t getUTF32Len() const;

		// Static unicode routines
		static size_t getUTF8Len(const wchar_t *utf16);
		static size_t getUTF8Len(const StringUTF32 &utf32);
		static size_t getUTF16Len(const StringUTF32 &utf32);

		inline std::string& cppStr() { return str; }
		inline const std::string& cppStr() const { return str; }

		Vector<String> split(char delimiter) const;
		Vector<String> split(String delimiter) const;
		static String concatList(const Vector<String>& list, String separator);

		//////////

		String operator += (const String &p);
		String operator += (const char* p);
		String operator += (const wchar_t* p);
		String operator += (const double &p);
		String operator += (const int &p);
		String operator += (const Character &p);

		bool operator== (const String& rhp) const;
		bool operator!= (const String& rhp) const;
		bool operator< (const String& rhp) const;
		bool operator> (const String& rhp) const;
		bool operator<= (const String& rhp) const;
		bool operator>= (const String& rhp) const;

	private:
		Character* getCharPointer(size_t pos);
		static size_t UTF8toUTF16(const char *utf8,wchar_t *utf16);
		static size_t UTF16toUTF8(const wchar_t *utf16,char *utf8);
		static size_t UTF32toUTF8(const utf32type *utf32,char *utf8);

		std::string str;
	};

	String operator+ (const String& lhp, const String& rhp);
	std::ostream& operator<< (std::ostream& os, const String& rhp);
	std::istream& operator>> (std::istream& is, String& rhp);

	using StringArray = Vector<String>;

}

namespace std {
	template<>
	struct hash<Halley::String>
	{
		size_t operator()(const Halley::String& s) const
		{
			return std::hash<std::string>()(s.cppStr());
		}
	};
}

max_stars_repo_name: Healthire/halley

max_issues_repo_name: Healthire/halley

max_forks_repo_name: Healthire/halley


text:
/**
 * Last modified: Tue Aug 30 16:29:05 CST 2016
 * @author: Lie Yan
 * @email: robin.lie.yan@outlook.com
 *
 */

#pragma once

#include <gsl.h>
#include <iostream>

namespace grund {

template<int TypeTag, typename ValueType>
class Unit {
public:
  using intrinsic_type = ValueType;

  Unit() = default;
  Unit(const Unit &rhs) = default;
  Unit(Unit &&rhs) = default;
  Unit &operator=(const Unit &rhs) = default;
  Unit &operator=(Unit &&rhs)      = default;

  constexpr explicit Unit(const ValueType &v) : _value(v) {}
  constexpr explicit operator ValueType() const { return _value; }
  constexpr ValueType value() const { return _value; }

  // unary minus
  constexpr Unit operator-() const { return Unit(-_value); }

  constexpr Unit &operator+=(const Unit &rhs) {
    _value += rhs._value;
    return *this;
  }

  constexpr Unit operator+(const Unit &rhs) const {
    Unit ret(*this);
    return ret += rhs;
  }

  constexpr Unit &operator-=(const Unit &rhs) {
    _value -= rhs._value;
    return *this;
  }

  constexpr Unit operator-(const Unit &rhs) const {
    Unit ret(*this);
    return ret -= rhs;
  }

  constexpr bool operator==(const Unit &rhs) const {
    return _value == rhs._value;
  }

  // derived operator
  constexpr bool operator!=(const Unit &rhs) const {
    return not(*this == rhs);
  }

  constexpr bool operator<(const Unit &rhs) const {
    return _value < rhs._value;
  }

  // derived operator
  constexpr bool operator<=(const Unit &rhs) const {
    return not(rhs < *this);
  }
  // derived operator
  constexpr bool operator>(const Unit &rhs) const { return rhs < *this; }
  // derived operator
  constexpr bool operator>=(const Unit &rhs) const {
    return not(*this < rhs);
  }

  friend std::ostream &operator<<(std::ostream &out, const Unit &m) {
    out << ValueType(m);
    return out;
  }

  template<typename T>
  friend constexpr Unit operator*(const Unit &u, const T &v) {
    return Unit{u.value() * v};
  }

  // derived operator
  template<typename T>
  friend constexpr Unit operator*(const T &v, const Unit &u) {
    return operator*(u, v);
  }

  template<typename T>
  friend constexpr Unit operator/(const Unit &u, const T &v) {
    return Unit{u.value() / v};
  }

  template <typename U>
  U to() const { return U(_value); }

private:
  ValueType _value;
};
}


text:
/* population genetics tools- mainly for working with site frequency spectrum
/
/
/ Andrew Kern
*/

#include "stdio.h"
#include "math.h"
#include "popGenTools.h"
#include "numerical.h"
#include "assert.h"
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf.h>

/* function for sfs generation */
double my_f(double q, void * p){
  struct my_f_params * params = (struct my_f_params *) p;
  double beta = params->beta;
  int i = params->i;
  int n = params->n;

  return ((1.0 - exp(-2.0 * beta * (1.0 - q))) / (1.0 - exp(-2.0 * beta))) * (2.0 / (q * (1.0 - q))) * gsl_ran_binomial_pdf(i, q, n);
}

/* integral of f over allele freq */
double my_F(double beta, void * p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_f_params fParams;
  gsl_function f;
  gsl_integration_workspace * w = gsl_integration_workspace_alloc (100000);
  double result,error;

  if (beta == 0){
    gsl_integration_workspace_free(w);
    return 1.0 / params->i;
  }
  else{
    fParams.beta = beta;
    fParams.i = params->i;
    fParams.n = params->n;
    f.function = &(my_f);
    f.params = &fParams;
    gsl_integration_qags(&f,0.0,1.0,1e-5,1e-5, 100000,w, &result,&error);
    gsl_integration_workspace_free(w);
    return result;
  }
}

/* integral of f over allele freq using numerical recipes integration */
double my_F2(double beta, void * p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_f_params fParams;
  double result;

  if (beta == 0){
    return 1.0 / params->i;
  }
  else{
    fParams.beta = beta;
    fParams.i = params->i;
    fParams.n = params->n;
    result = d_qromb2(&my_f, 0.0, 1.0, &fParams);
    return result;
  }
}

/* this returns the prob of SNP freq i given n and beta */
double snpProb(double beta, void * p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_F_params tempParams;
  long double tot;
  int i;
  gsl_function F;

  /* add up total prob space */
  tot = 0.0;
  tempParams.n = params->n;

  for(i = 1; i < params->n ; i++){
    tempParams.i = i;
    F.function =  &my_F;
    F.params = &tempParams;
    tot += GSL_FN_EVAL(&F, beta);
  }

  /* reset parameters */
  F.function =  &my_F;
  F.params = params;

  return   GSL_FN_EVAL(&F, beta) / tot ;
}

/* trying to speedup the snpProb routine with a lookup of the denominators (in snpProb_params) */
double snpProbDenomLookup(double beta, void * p){
  struct my_snpProb_params * params = (struct my_snpProb_params *) p;
  struct my_F_params tempParams;
  gsl_function F;

 /* set parameters */
  tempParams.n = params->n;
  tempParams.i = params->i;
  F.function =  &my_F;
  F.params = &tempParams;

  return   GSL_FN_EVAL(&F, beta) / params->denom ;
}

/* sampleSize vector - returns a vector of 1s or 0s depending on sample sizes in the data */
gsl_vector *sampleSizeVector(struct snp data[], int snpNumber, int maxSampleSize){
  gsl_vector *bools;
  int i;

  //alloc and initialize bools vector
  bools = gsl_vector_alloc(maxSampleSize + 1);
  gsl_vector_set_zero(bools);

  //go through data, set each sampleSize value to 1
  for(i = 0; i < snpNumber; i++){
    gsl_vector_set(bools,data[i].n, 1);
  }
  return(bools);
}

/* returns a vector of snpProbDenominators from i = 2 to n, 1 indexed */
gsl_vector *makeSnpProbDenomVector(double beta, int maxSampleSize, gsl_vector *sampleSizeVector){
  gsl_vector *probs;
  int i, j, test;
  double tot;
  struct my_F_params tempParams;
  gsl_function F;

  probs = gsl_vector_alloc(maxSampleSize + 1);

  /* go from n = 2 to maxSampleSize, calc total prob space for each, put in vector probs */
  for(j = 2; j < maxSampleSize + 1; j++){
    //check if sampleSize is included
    test = gsl_vector_get(sampleSizeVector,j);
      if (test){
      tot = 0.0;
      tempParams.n = j;
      for(i = 1; i < j ; i++){
	tempParams.i = i;
	F.function =  &my_F;
	F.params = &tempParams;
	tot += GSL_FN_EVAL(&F, beta);
      }
      gsl_vector_set(probs, j, tot);
    }
  }
   return probs;
}

/* snpProbDenom - returns denom for snpProb */
double snpProbDenom(double beta, int sampleSize){
  int i;
  double tot;
  struct my_F_params tempParams;
  gsl_function F;

  tot = 0.0;
  tempParams.n = sampleSize;
  for(i = 1; i < sampleSize ; i++){
	tempParams.i = i;
	F.function =  &my_F;
	F.params = &tempParams;
	tot += GSL_FN_EVAL(&F, beta);
      }
  return(tot);
}

/* snpProbMatrix- returns a matrix of log snpProbs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
   variable sampleSizes (2 to max), columns represent freqs (1 to max -1) */
gsl_matrix *snpProbMatrix(double beta,  int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools){
  gsl_vector *denoms;
  gsl_matrix *probs;
  struct my_snpProb_params FParams;
  gsl_function f;
  int i, j;

  //alloc probs
  probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
  gsl_matrix_set_zero(probs);
  //get denoms
  denoms = makeSnpProbDenomVector(beta,maxSampleSize, sampleSizeVector);

  //go through sampleSizes
  for(i = 2; i <= maxSampleSize; i++){
    //have sample size i?
    if (gsl_vector_get(sampleSizeVector, i)){
      //go through freqs
      for(j = 1; j < maxSampleSize; j++){
	//have freq j?
	if(gsl_matrix_get(sfsBools, i, j)){
	  //calc prob
	  FParams.i = j;
	  FParams.n = i;
	  FParams.denom = gsl_vector_get(denoms,i);
	  f.function = &snpProbDenomLookup;
	  f.params = &FParams;
	  gsl_matrix_set(probs, i, j, log(GSL_FN_EVAL(&f, beta)));
	}
      }
    }
  }
  gsl_vector_free(denoms);
  return(probs);
}
/* snpProbMatrix- returns a matrix of  snpProbs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
   variable sampleSizes (2 to max), columns represent freqs (1 to max -1) */
gsl_matrix *snpProbMatrixNotLog(double beta,  int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools){
  gsl_vector *denoms;
  gsl_matrix *probs;
  struct my_snpProb_params FParams;
  gsl_function f;
  int i, j;

  //alloc probs
  probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
  gsl_matrix_set_zero(probs);
  //get denoms
  denoms = makeSnpProbDenomVector(beta,maxSampleSize, sampleSizeVector);

  //go through sampleSizes
  for(i = 2; i <= maxSampleSize; i++){
    //have sample size i?
    if (gsl_vector_get(sampleSizeVector, i)){
      //go through freqs
      for(j = 1; j < maxSampleSize; j++){
	//have freq j?
	if(gsl_matrix_get(sfsBools, i, j)){
	  //calc prob
	  FParams.i = j;
	  FParams.n = i;
	  FParams.denom = gsl_vector_get(denoms,i);
	  f.function = &snpProbDenomLookup;
	  f.params = &FParams;
	  gsl_matrix_set(probs, i, j, GSL_FN_EVAL(&f, beta));
	}
      }
    }
  }
  gsl_vector_free(denoms);
  return(probs);
}
/* snpProbMatrixNotLogFull- returns a matrix of  snpProbs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
   variable sampleSizes (2 to max), columns represent freqs (1 to max -1). All freqs/samplesizes included */
gsl_matrix *snpProbMatrixNotLogFull(double beta,  int maxSampleSize, gsl_vector *sampleSizeVector){
  gsl_vector *denoms;
  gsl_matrix *probs;
  struct my_snpProb_params FParams;
  gsl_function f;
  int i, j;

  //alloc probs
  probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
  gsl_matrix_set_zero(probs);

  //get denoms

  denoms = makeSnpProbDenomVector(beta,maxSampleSize, sampleSizeVector);

  //go through sampleSizes
  for(i = 2; i <= maxSampleSize; i++){
    //check for sampleSize
    if (gsl_vector_get(sampleSizeVector,i)){
      //go through freqs
      for(j = 1; j < maxSampleSize; j++){
	//calc prob
	FParams.i = j;
	FParams.n = i;
	FParams.denom = gsl_vector_get(denoms,i);
	f.function = &snpProbDenomLookup;
	f.params = &FParams;
	gsl_matrix_set(probs, i, j, GSL_FN_EVAL(&f, beta));
      }
    }
  }
  gsl_vector_free(denoms);
  return(probs);
}

/* snpProbVectorNotLog- returns a vector of log snpProbs at specified sampleSize columns represent freqs (1 to max -1) */
gsl_vector *snpProbVectorNotLog(double beta,  int sampleSize){
  double denom;
  gsl_vector *probs;
  struct my_snpProb_params FParams;
  gsl_function f;
  int i;

  //alloc probs
  probs = gsl_vector_alloc(sampleSize);
  gsl_vector_set_zero(probs);
  //get denoms
  denom = snpProbDenom(beta, sampleSize);

  //go through freqs
  for(i = 1; i < sampleSize; i++){
    FParams.i = i;
    FParams.n = sampleSize;
    FParams.denom = denom;
    f.function = &snpProbDenomLookup;
    f.params = &FParams;
    gsl_vector_set(probs, i, GSL_FN_EVAL(&f, beta));
  }
  return(probs);
}

/* likelihood function for SFS given the data; p here contains the data and the weights */
double my_lik(double beta, void * p){
	struct my_lik_params * params = (struct my_lik_params *) p;
	struct my_F_params FParams;
	gsl_function f;
	double lik;
	int i;

	lik = 0;
	for(i = 0; i < params->snpNumber; i++){
		FParams.i = params->data[i].i;
		FParams.n = params->data[i].n;
		f.function = &snpProb;
		f.params = &FParams;
		lik += log(GSL_FN_EVAL(&f, beta));
		}
	return -lik;
}

/* likelihood function for SFS where each SNP has independent selection coeff */
double sfsLikBetaVector(gsl_vector *betas, void * p){
	struct my_lik_params * params = (struct my_lik_params *) p;
	struct my_F_params FParams;
	gsl_function f;
	double lik;
	int i;

	lik = 0;
	for(i = 0; i < params->snpNumber; i++){
		FParams.i = params->data[i].i;
		FParams.n = params->data[i].n;
		f.function = &snpProb;
		f.params = &FParams;
		lik += log(GSL_FN_EVAL(&f, gsl_vector_get(betas, i)));
		}
	return -lik;
}

/* weighted likelihood function for SFS given the data; here contains the data and the weights */
double weightedLik(double beta, void * p){
  struct my_lik_params * params = (struct my_lik_params *) p;
  struct my_F_params FParams;
  gsl_function f;
  double lik;
  int i;

  lik = 0;
  for(i = 0; i < params->snpNumber; i++){
    FParams.i = params->data[i].i;
    FParams.n = params->data[i].n;
    f.function = &snpProb;
    f.params = &FParams;
    lik += gsl_vector_get(params->weights, i) * log(GSL_FN_EVAL(&f, beta));
  }
  return -lik;
}

/* weighted likelihood function for SFS given the data; here contains the data and the weights and uses the snpProbLookup routine */
double weightedLikLook(double beta, void * p){
  struct my_lik_params * params = (struct my_lik_params *) p;
  double lik;
  int i;
  gsl_matrix *probs;

  //make prob matrix (note this is in logs)
  probs = snpProbMatrix(beta, params->maxSampleSize, params->sampleSizeVector, params->sfsBools);
  //now go and tally likelihood
  lik = 0;
  for(i = 0; i < params->snpNumber; i++){
    lik += gsl_matrix_get(probs,params->data[i].n,params->data[i].i)  * gsl_vector_get(params->weights,i);
  }
  gsl_matrix_free(probs);
  return -lik;
}

double weightedLikLookCI(double beta, void * p){
  struct my_likCI_params *params = (struct my_likCI_params *) p;
  struct my_lik_params likParams;
  gsl_function l;
  double result;

  likParams.data = params->data;
  likParams.snpNumber = params->snpNumber;
  likParams.sampleSizeVector = params->sampleSizeVector;
  likParams.maxSampleSize = params->maxSampleSize;
  likParams.weights = params->weights;
  likParams.sfsBools = params->sfsBools;
  l.function = &weightedLikLook;
  l.params = &likParams;
  result = GSL_FN_EVAL(&l, beta) - (params->lMax) - (params->logUnits) ;
  //  printf("%f\t%f\t%f\n",params->lMax,GSL_FN_EVAL(&l, beta),result);
  return result;
}

/* wrapper for my_lik for use in mnbrak */
double likWrap(double beta){
	gsl_function f;

	f.function = &my_lik;
	f.params = NULL;
	return GSL_FN_EVAL(&f, beta);
}

/* wrapper for my_lik for use in mnbrak */
double  wlikWrap(double beta, void * p){
	gsl_function f;
	struct my_lik_params * params = (struct my_lik_params *) p;

	f.function = &weightedLik;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

/* wrapper for my_lik for use in mnbrak */
double  wlikWrapLook(double beta, void * p){
	gsl_function f;
	struct my_lik_params * params = (struct my_lik_params *) p;

	f.function = &weightedLikLook;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

double  wlikCIWrapLook(double beta, void * p){
	gsl_function f;
	struct my_likCI_params * params = (struct my_likCI_params *) p;

	f.function = &weightedLikLookCI;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

/*error checking function for likelihood function */
void printProfile(double min, double max, void * p){
  gsl_function f;
  struct my_lik_params *params =  (struct my_lik_params *) p;
  int i;

  f.function = &weightedLikLook;
  f.params = params;
  for(i = min; i < max; i++){
    printf("%d\t%f\n",i, GSL_FN_EVAL(&f, i));
  }

}

/* get ML estimate of beta using likelihood function; BRENT method for gsl */
double ml_est(double * lik_beta_hat){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc;
	gsl_function L;

	/* get minimum bracket */
	ax = -10.0;
	bx = 0.0;
	cx = 10.0;

	mnbrak(&ax, &bx, &cx, &fa, &fb, &fc, &my_lik);

	/* check if there is a minimum, if not return inf */
	if (fa == fb || fc == fb || ax > cx){
		*lik_beta_hat = fa;
		return 666.0;
		}
	else{

		/* do mimization */

		L.function = &my_lik;
		//should have param line here but no params....

		T = gsl_min_fminimizer_brent;
		s = gsl_min_fminimizer_alloc(T);
		gsl_min_fminimizer_set(s, &L, bx, ax, cx);

		do{
			iter++;
			status = gsl_min_fminimizer_iterate(s);
			 m = gsl_min_fminimizer_x_minimum (s);
			  a = gsl_min_fminimizer_x_lower (s);
			  b = gsl_min_fminimizer_x_upper (s);

			  status = gsl_min_test_interval (a, b, 0.001, 0.0);

		}
	  while (status == GSL_CONTINUE && iter < max_iter);
	  gsl_min_fminimizer_free(s);

	  return m;
	}
}

/* get weighted ML estimate of beta using likelihood function; BRENT method for gsl */
double weighted_ml_est(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, min;
	gsl_function L;
	struct my_lik_params * params = (struct my_lik_params *) p;

	// get minimum bracket
	ax = -10.0;
	bx = 1.0;
	cx = 10;
	L.function = &weightedLik;
	L.params = params;
	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &wlikWrap, bx, params);

	//check if there is a minimum, if not return 666
	if (fa == fb || fc == fb || ax > cx){
	  min = fa;
	  if(fb < min){
	    min = fb;
	  }
	  if(fc < min){
	    min = fc;
	  }
	  *lik_beta_hat = min;

	  return 666.0;
	}
	else{
	  // do mimization

	  //initialize lik function
	  L.function = &weightedLik;
	  L.params = params;

	  //min routine
	  T = gsl_min_fminimizer_brent;
	  s = gsl_min_fminimizer_alloc(T);
	  gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	  do{
	    iter++;
	    status = gsl_min_fminimizer_iterate(s);
	    m = gsl_min_fminimizer_x_minimum (s);
	    a = gsl_min_fminimizer_x_lower (s);
	    b = gsl_min_fminimizer_x_upper (s);
	    status = gsl_min_test_interval (a, b, 0.001, 0.0);
	  }
	  while (status == GSL_CONTINUE && iter < max_iter);
	  *lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	  gsl_min_fminimizer_free(s);
	  return m;
	}
}

double weighted_ml_est_lookup(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, dummy;
	gsl_function L;
	struct my_lik_params * params = (struct my_lik_params *) p;

	// get minimum bracket
	ax = -20.0;
	bx = -1.0;
	cx = 20;
	L.function = &weightedLikLook;
	L.params = params;
	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &wlikWrapLook, bx, params);
	//swap bounds if needed
	if (cx < ax){
	  dummy = cx;
	  cx = ax;
	  ax = dummy;
	}

	// do mimization

	//initialize lik function
	L.function = &weightedLikLook;
	L.params = params;

	//min routine
	T = gsl_min_fminimizer_brent;
	s = gsl_min_fminimizer_alloc(T);
	gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	do{
	  iter++;
	  status = gsl_min_fminimizer_iterate(s);
	  m = gsl_min_fminimizer_x_minimum (s);
	  a = gsl_min_fminimizer_x_lower (s);
	  b = gsl_min_fminimizer_x_upper (s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	*lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	gsl_min_fminimizer_free(s);
	return m;
}

/* same as above but outputs 666 on minimization error */
double weighted_ml_est_lookup_errReport(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, dummy;
	gsl_function L;
	struct my_lik_params * params = (struct my_lik_params *) p;

	// get minimum bracket
	ax = -20.0;
	bx = -1.0;
	cx = 20;
	L.function = &weightedLikLook;
	L.params = params;
	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &wlikWrapLook, bx, params);
	//swap bounds if needed
	if (cx < ax){
	  dummy = cx;
	  cx = ax;
	  ax = dummy;
	}
	//check for error
	if (fb >= fc || fb >= fa){
	  return(666.0);
	}

	// do mimization

	//initialize lik function
	L.function = &weightedLikLook;
	L.params = params;

	//min routine
	T = gsl_min_fminimizer_brent;
	s = gsl_min_fminimizer_alloc(T);
	gsl_set_error_handler_off ();
	status = gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	if (status == GSL_EINVAL){
	  return(666);
	}
	do{
	  iter++;
	  status = gsl_min_fminimizer_iterate(s);
	  m = gsl_min_fminimizer_x_minimum (s);
	  a = gsl_min_fminimizer_x_lower (s);
	  b = gsl_min_fminimizer_x_upper (s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	*lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	gsl_min_fminimizer_free(s);
	return m;
}

/* get weighted ML lower CI beta using weightedLikLookCI method; BRENT method for gsl */
double weighted_ml_CILower_lookup(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_root_fsolver_type *T;
	gsl_root_fsolver *s;
	double m, a, b, ax, cx;
	gsl_function L;
	struct my_likCI_params * params = (struct my_likCI_params *) p;

	// get minimum bracket
	ax = -50.0;
	cx = params->beta_hat;

	// do root finding
	//initialize lik function
	L.function = &weightedLikLookCI;
	L.params = params;

	//min routine
	T = gsl_root_fsolver_brent;
	s = gsl_root_fsolver_alloc(T);
	gsl_root_fsolver_set(s, &L, ax, cx);
	do{
	  iter++;
	  status = gsl_root_fsolver_iterate(s);
	  m = gsl_root_fsolver_root(s);
	  a = gsl_root_fsolver_x_lower(s);
	  b = gsl_root_fsolver_x_upper(s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	gsl_root_fsolver_free(s);
	return m;
}

/* get weighted ML upper CI beta using weightedLikLookCI method; BRENT method for gsl */
double weighted_ml_CIUpper_lookup(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_root_fsolver_type *T;
	gsl_root_fsolver *s;
	double m, a, b, ax, cx;
	gsl_function L;
	struct my_likCI_params * params = (struct my_likCI_params *) p;

	// get minimum bracket
	cx = 100.0;
	ax = params->beta_hat;

	// do root finding
	//initialize lik function
	L.function = &weightedLikLookCI;
	L.params = params;

	//min routine
	T = gsl_root_fsolver_brent;
	s = gsl_root_fsolver_alloc(T);
	gsl_root_fsolver_set(s, &L, ax, cx);
	do{
	  iter++;
	  status = gsl_root_fsolver_iterate(s);
	  m = gsl_root_fsolver_root(s);
	  a = gsl_root_fsolver_x_lower(s);
	  b = gsl_root_fsolver_x_upper(s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	gsl_root_fsolver_free(s);
	return m;
}

/* summarizeSFS- returns a matrix of counts which summarize the SFS conditional
   upon sampleSize. rows = sampleSize, column = freqs */

gsl_matrix *summarizeSFS(int maxSampleSize, struct snp data[], int snpNumber){
  gsl_matrix *counts;
  int i;

  //alloc matrix of size maxSampleSize by maxSampleSize
  counts = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize + 1);
  gsl_matrix_set_zero(counts);
  assert(counts != NULL);

  gsl_matrix_set_zero(counts);

  //go through data, fill up matrix
  for(i = 0; i < snpNumber; i++){
    gsl_matrix_set(counts, data[i].n, data[i].i, gsl_matrix_get(counts, data[i].n, data[i].i) + 1);
  }
  return(counts);

}

/* summarizeSFS- returns a matrix of booleans which summarize the SFS conditional
   upon sampleSize. rows = sampleSize, column = freqs */

gsl_matrix *summarizeSFSBool(int maxSampleSize, struct snp data[], int snpNumber){
  gsl_matrix *counts;
  int i;

  //alloc matrix of size maxSampleSize by maxSampleSize
  counts = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize + 1);
  gsl_matrix_set_zero(counts);
  assert(counts != NULL);

  gsl_matrix_set_zero(counts);

  //go through data, fill up matrix
  for(i = 0; i < snpNumber; i++){
    gsl_matrix_set(counts, data[i].n, data[i].i, 1);
  }
  return(counts);

}

/* simulateSiteFreq- returns an int representing the frequency of
   of a site (unfolded), i, in sampleSize n, conditional on sel.
   coeff. beta , also needs a pointer to a gsl random number generator */

int simulateSiteFreq(int sampleSize, double alpha, void *r){
  int i;
  double sum, probs[sampleSize - 1], rand;
  gsl_function p;
  struct my_F_params simParams;
  gsl_rng * rn = (gsl_rng *) r;

  simParams.n = sampleSize;
  /* put probs in vector */
  for(i = 0; i < sampleSize - 1; i++){
    simParams.i = i + 1;
    p.function = &snpProb;
    p.params = &simParams;
    probs[i] = GSL_FN_EVAL(&p, alpha);
  }

  /* choose frequency */
  rand = gsl_rng_uniform(rn);
  sum = 0;
  for(i = 0; i < sampleSize - 1; i++){
    sum += probs[i];
    if (rand <= sum){
      return(i+1);
    }
  }
  return(666);
}

/* ascertainment stuff -- following Nielsen et al. 2004  */

/* simpleAscertain--  Nielsen et al. 2004 eqn 2 *, returns the prob of ascertainment
given a sample freq */
double simpleAscertain(double sampleFreq, void * p){
  struct my_F_params * params = (struct my_F_params *) p;
  double sampleSize, ascSize,tmpSampleFreq;
  double num1 = 0.0;
  double num2 = 0.0;

  sampleSize = params->n;
  ascSize = params->ascSize;
  tmpSampleFreq = (int) sampleFreq;

  //test for weirdness

  if (sampleFreq < ascSize){
    num1 = 0;
  }
  else{
    num1 = gsl_sf_choose(sampleFreq, ascSize);
  }

  if ((sampleSize -sampleFreq) < ascSize){
    num2 = 0;
  }
  else{
    num2 =  gsl_sf_choose(sampleSize - sampleFreq, ascSize);
  }

  return(1.0 - ((num1 + num2) / gsl_sf_choose(sampleSize, ascSize)));
}

/* outgroupAscertain-- equivalent to Nielsen et al. 2004 eqn 2,
 returns the prob of ascertainment given a sample freq, but
is for use when ascertainment was based on single outgroup comparisons
(e.g. divergence between reference genomes). specifically this
calculates the prob of sampling only ancestral alleles.
IMPORTANT- sampleFreq here is derived allele freq */
double outgroupAscertain(double sampleFreq, void * p){
  struct my_F_params * params = (struct my_F_params *) p;
  double sampleSize, ascSize;
  double num = 0.0;

  sampleSize = params->n;
  ascSize = params->ascSize;

  //test for weirdness
  if ((sampleSize - sampleFreq) < ascSize){
    num = 0;
  }
  else{
    num =  gsl_sf_choose((sampleSize - sampleFreq), ascSize);
  }

  return(num / gsl_sf_choose(sampleSize, ascSize));
}

/* probAscertainmentGivenModel-- Nielsen et al. 2005 eqn 8
   returns the prob of ascertainment given the model */
double probAscertainmentGivenModel(double beta, void *p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_F_params tempParams;

  long double tot, tmpres, tmpres2;
  int i=1;
  gsl_function F, G;

  /* add up total prob space */
  tot = 0.0;
  tempParams.n = params->n;
  tempParams.ascSize = params->ascSize;

  for(i = 1; i < params->n ; i++){
    tempParams.i = i;

    F.function =  &(snpProb);
    F.params = &tempParams;
    tmpres = 0.0;
    tmpres = GSL_FN_EVAL(&F, beta);
    G.function = &simpleAscertain;
    G.params = &tempParams;
    tmpres2 = GSL_FN_EVAL(&G, i) ;
    tot += tmpres * tmpres2;

  }

  return(tot) ;
}

/* probAscertainmentGivenModel-- Nielsen et al. 2005 eqn 8
   returns the prob of ascertainment given the model */
double probOutgroupAscertainmentGivenModel(double beta, void *p){
  struct my_F_params * params = (struct my_F_params *) p;
  struct my_F_params tempParams;

  long double tot, tmpres, tmpres2;
  int i=1;
  gsl_function F, G;

  /* add up total prob space */
  tot = 0.0;
  tempParams.n = params->n;
  tempParams.ascSize = params->ascSize;

  for(i = 1; i < params->n ; i++){
    tempParams.i = i;

    F.function =  &(snpProb);
    F.params = &tempParams;
    tmpres = 0.0;
    tmpres = GSL_FN_EVAL(&F, beta);
    G.function = &outgroupAscertain;
    G.params = &tempParams;
    tmpres2 = GSL_FN_EVAL(&G, i) ;
    tot += tmpres * tmpres2;

  }

  return(tot) ;
}

/* probAscertainmentGivenModelLookup-- Nielsen et al. 2005 eqn 8
   returns the prob of ascertainment given the model, lookup version
   requires snpProb matrix and ascProb matrix */
double probAscertainmentGivenModelLookup(double beta, int sampSize, gsl_matrix *snpProbs, gsl_matrix *ascProbs){

  long double tot;
  int i=1;

  /* add up total prob space */
  tot = 0.0;

  for(i = 1; i < sampSize ; i++){
    tot += gsl_matrix_get(snpProbs, sampSize, i) * gsl_matrix_get(ascProbs, sampSize, i);
  }
  return(tot) ;
}

/* probAscertainmentGivenModelHemiLookup-- same
   as above but expects vector not matrix of snpProbs */
double probAscertainmentGivenModelHemiLookup(double beta, int sampSize, gsl_vector *snpProbs, gsl_matrix *ascProbs){

  long double tot;
  int i=1;

  /* add up total prob space */
  tot = 0.0;

  for(i = 1; i < sampSize ; i++){
    tot += gsl_vector_get(snpProbs, i) * gsl_matrix_get(ascProbs, sampSize, i);
  }
  return(tot) ;
}

/* ascSize vector - returns a vector of 1s or 0s depending on ascertainment sample sizes in the data */
gsl_vector *ascSizeVector(struct snp data[], int snpNumber, int maxSampleSize){
  gsl_vector *bools;
  int i;

  //alloc and initialize bools vector
  bools = gsl_vector_alloc(maxSampleSize + 1);
  gsl_vector_set_zero(bools);

  //go through data, set each sampleSize value to 1
  for(i = 0; i < snpNumber; i++){
    gsl_vector_set(bools,data[i].ascSize, 1);
  }
  return(bools);
}

/* snpAscMatrix- returns a matrix of simpleAscertainment probs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
   variable sampleSizes (2 to max), columns represent freqs (1 to max -1). Currently this only supports a single ascertainment size  */
gsl_matrix *snpAscMatrix(int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools, \
		 int ascSize){
  gsl_matrix *probs;
  struct my_F_params FParams;
  gsl_function f;
  int i, j;

  //alloc probs
  probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
  gsl_matrix_set_zero(probs);
  //go through sampleSizes
  for(i = 2; i <= maxSampleSize; i++){
    //have sample size i?
    if (gsl_vector_get(sampleSizeVector, i)){
      //go through freqs
      for(j = 1; j < maxSampleSize; j++){
	//have freq j?
	if(gsl_matrix_get(sfsBools, i, j)){
	  //calc prob
	  FParams.i = j;
	  FParams.n = i;
	  FParams.ascSize = ascSize;
	  f.function = &simpleAscertain;
	  f.params = &FParams;
	  gsl_matrix_set(probs, i, j, GSL_FN_EVAL(&f, j));
	}
      }
    }
  }
  return(probs);
}

/* snpOutgroupAscMatrix- returns a matrix of outgroupAscertainment probs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
   variable sampleSizes (2 to max), columns represent freqs (1 to max -1). Currently this only supports a single ascertainment size  */
gsl_matrix *snpOutgroupAscMatrix(int maxSampleSize, gsl_vector *sampleSizeVector, int ascSize){
  gsl_matrix *probs;
  struct my_F_params FParams;
  gsl_function f;
  int i, j;

  //alloc probs
  probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
  gsl_matrix_set_zero(probs);
  //go through sampleSizes
  for(i = 2; i <= maxSampleSize; i++){
    //have sample size i?
    if (gsl_vector_get(sampleSizeVector, i)){
      //go through freqs
      for(j = 1; j < i; j++){
	//calc prob
	FParams.i = j;
	FParams.n = i;
	FParams.ascSize = ascSize;
	f.function = &outgroupAscertain;
	f.params = &FParams;
	gsl_matrix_set(probs, i, j, GSL_FN_EVAL(&f, j));
      }
    }
  }
  return(probs);
}

/*estimateAscSFS-- ML estimates of SFS, given ascertainment. returns
 a matrix of probs. matrix is maxSampleSize+1 by maxSampleSize. rows represent
variable sampleSizes (2 to max), columns represent freqs (1 to max -1). Currently this only supports a single ascertainment size */

gsl_matrix *estimateAscSFS(int maxSampleSize, gsl_vector *sampleSizeVector, gsl_matrix *sfsBools, \
			   int ascSize, gsl_matrix *sfsSummary){
  gsl_matrix *probs, *ascProbs;
    int i, j;
  double ssTmp, p_hat;
  double xij;

  //alloc probs
  probs = gsl_matrix_alloc(maxSampleSize + 1, maxSampleSize);
  gsl_matrix_set_zero(probs);
  //calculate ascProbs
  ascProbs = snpAscMatrix(maxSampleSize, sampleSizeVector, sfsBools, \
			  ascSize);

  //go through sampleSizes
  for(i = 2; i <= maxSampleSize; i++){
    //have sample size i?
    if (gsl_vector_get(sampleSizeVector, i)){
      //go through freqs to tally up denom at sampleSize
      ssTmp = 0.0;
      for(j = 1; j < maxSampleSize; j++){
	xij = gsl_matrix_get(sfsSummary, i, j);
	if (xij > 0){
	  ssTmp += xij / gsl_matrix_get(ascProbs, i, j);
	}
      }
      //go back through freqs to assign p_hats at sampleSize
      for(j = 1; j < maxSampleSize; j++){
	xij = gsl_matrix_get(sfsSummary, i, j);
	p_hat = (double) xij / gsl_matrix_get(ascProbs, i, j);
	gsl_matrix_set(probs, i, j, p_hat / ssTmp);
      }
    }
  }
  return(probs);
}

/* weighted likelihood function for SFS given the data and ascertainment; here contains the data and the weights and uses the snpProbLookup routine */
double weightedLikLookAsc(double beta, void * p){
  struct my_lik_params * params = (struct my_lik_params *) p;
  double lik, num;
  int i;
  gsl_matrix *probs, *ascProbs;
  gsl_vector *ascDenoms;

  //make prob matrix (note this is in logs)
  probs = snpProbMatrix(beta, params->maxSampleSize, params->sampleSizeVector, params->sfsBools);
  //make ascProb matrix (NOT in logs!)
  ascProbs = snpAscMatrix(params->maxSampleSize, params->sampleSizeVector, params->sfsBools, params->ascSize);

  //make vector of denoms
  ascDenoms = gsl_vector_alloc(params->maxSampleSize + 1);
  gsl_vector_set_zero(ascDenoms);
  for(i = params->ascSize; i <= params->maxSampleSize; i++){
     gsl_vector_set(ascDenoms,i,probAscertainmentGivenModelLookup(beta, i, probs, ascProbs));
  }

  //now go and tally likelihood
  lik = 0;
  num = 0;
  for(i = 0; i < params->snpNumber; i++){
    lik += gsl_matrix_get(probs,params->data[i].n,params->data[i].i)	\
      + log(gsl_matrix_get(ascProbs,params->data[i].n,params->data[i].i)) \
      - log(gsl_vector_get(ascDenoms, params->data[i].n));
     }
  gsl_matrix_free(probs);
  gsl_matrix_free(ascProbs);
  gsl_vector_free(ascDenoms);

  return -lik;
}

/* weighted likelihood function for SFS given the data and outgroup ascertainment; here contains the data and the weights and uses the snpProbLookup routine */
double weightedLikLookOutgroupAsc(double beta, void * p){
  struct my_lik_params * params = (struct my_lik_params *) p;
  double lik, num;
  int i;
  gsl_matrix *probs, *ascProbs;
  gsl_vector *ascDenoms;

  //make prob matrix
  probs = snpProbMatrixNotLogFull(beta, params->maxSampleSize, params->sampleSizeVector);
  //make ascProb matrix
  ascProbs = snpOutgroupAscMatrix(params->maxSampleSize, params->sampleSizeVector, params->ascSize);

  //make vector of denoms
  ascDenoms = gsl_vector_alloc(params->maxSampleSize + 1);
  gsl_vector_set_zero(ascDenoms);
  for(i = params->ascSize; i <= params->maxSampleSize; i++){
     gsl_vector_set(ascDenoms,i,probAscertainmentGivenModelLookup(beta, i, probs, ascProbs));
  }

  //now go and tally likelihood
  lik = 0;
  num = 0;
  for(i = 0; i < params->snpNumber; i++){
    num = (gsl_matrix_get(probs,params->data[i].n,params->data[i].i) * gsl_matrix_get(ascProbs,params->data[i].n,params->data[i].i)) \
      / gsl_vector_get(ascDenoms, params->data[i].n);
    lik += log(num);
    //printf("first: %f\t second: %f\t third: %f\n",gsl_matrix_get(probs,params->data[i].n,params->data[i].i),gsl_matrix_get(ascProbs,params->data[i].n,params->data[i].i),gsl_vector_get(ascDenoms, params->data[i].n));
    }

  gsl_matrix_free(probs);
  gsl_matrix_free(ascProbs);
  gsl_vector_free(ascDenoms);

  return -lik;
}

/* wrapper for my_lik for use in mnbrak */
double  wlikWrapLookAsc(double beta, void * p){
	gsl_function f;
	struct my_lik_params * params = (struct my_lik_params *) p;

	f.function = &weightedLikLookAsc;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

/* wrapper for my_lik for use in mnbrak */
double  wlikWrapLookOutgroupAsc(double beta, void * p){
	gsl_function f;
	struct my_lik_params * params = (struct my_lik_params *) p;

	f.function = &weightedLikLookOutgroupAsc;
	f.params = params;
	return GSL_FN_EVAL(&f, beta);
}

/* get weighted ML estimate of beta using likelihood function (lookup) conditional on ascertainment
; BRENT method for gsl */
double weighted_ml_est_lookup_asc(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, dummy;
	gsl_function L;
	struct my_lik_params * params = (struct my_lik_params *) p;

	// get minimum bracket
	ax = -20.0;
	bx = -1.0;
	cx = 20;
	L.function = &weightedLikLookAsc;
	L.params = params;
	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &wlikWrapLookAsc, bx, params);
	//swap bounds if needed
	if (cx < ax){
	  dummy = cx;
	  cx = ax;
	  ax = dummy;
	}
	// do mimization

	//initialize lik function
	L.function = &weightedLikLookAsc;
	L.params = params;

	//min routine
	T = gsl_min_fminimizer_brent;
	s = gsl_min_fminimizer_alloc(T);
	gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	do{
	  iter++;
	  status = gsl_min_fminimizer_iterate(s);
	  m = gsl_min_fminimizer_x_minimum (s);
	  a = gsl_min_fminimizer_x_lower (s);
	  b = gsl_min_fminimizer_x_upper (s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	*lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	gsl_min_fminimizer_free(s);
	return m;
}

/* get weighted ML estimate of beta using likelihood function (lookup) conditional on outgroup ascertainment
; BRENT method for gsl */
double weighted_ml_est_lookup_outgroup_asc(double * lik_beta_hat, void * p){
	int status;
	int iter = 0;
	int max_iter = 100;
	const gsl_min_fminimizer_type *T;
	gsl_min_fminimizer *s;
	double m, a, b, ax, bx, cx, fa, fb, fc, dummy;
	gsl_function L;
	struct my_lik_params * params = (struct my_lik_params *) p;

	// get minimum bracket
	ax = -20.0;
	bx = -1.0;
	cx = 20;
	L.function = &weightedLikLookOutgroupAsc;
	L.params = params;
	mnbrak2(&ax, &bx, &cx, &fa, &fb, &fc, &wlikWrapLookOutgroupAsc, bx, params);
	//swap bounds if needed
	if (cx < ax){
	  dummy = cx;
	  cx = ax;
	  ax = dummy;
	}
	// do mimization

	//initialize lik function
	L.function = &weightedLikLookOutgroupAsc;
	L.params = params;

	//min routine
	T = gsl_min_fminimizer_brent;
	s = gsl_min_fminimizer_alloc(T);
	gsl_min_fminimizer_set(s, &L, bx, ax, cx);
	do{
	  iter++;
	  status = gsl_min_fminimizer_iterate(s);
	  m = gsl_min_fminimizer_x_minimum (s);
	  a = gsl_min_fminimizer_x_lower (s);
	  b = gsl_min_fminimizer_x_upper (s);
	  status = gsl_min_test_interval (a, b, 0.001, 0.0);
	}
	while (status == GSL_CONTINUE && iter < max_iter);
	*lik_beta_hat = gsl_min_fminimizer_f_minimum(s);
	gsl_min_fminimizer_free(s);
	return m;
}

/* likelihood function for SFS where each SNP has
 independent selection coeff and outgroup ascertainment */
double sfsLikBetaVectorOutgroupAsc(gsl_vector *betas, void * p){
  struct my_lik_params * params = (struct my_lik_params *) p;
  double pSNP, pAsc, pAscModel,lik, tmp;
  int i;
  gsl_matrix *ascProbs;
  gsl_vector *snpProbs;

  lik = 0;
  //make ascProb matrix
  ascProbs = snpOutgroupAscMatrix(params->maxSampleSize, params->sampleSizeVector, \
				  params->ascSize);
  //go through snps
  for(i = 0; i < params->snpNumber; i++){
    //create vector of snpProbs, find pSNP_i
    snpProbs = snpProbVectorNotLog(gsl_vector_get(betas, i), params->data[i].n);
    pSNP = gsl_vector_get(snpProbs, params->data[i].i);

    //lookup pAsc_i
    pAsc = gsl_matrix_get(ascProbs, params->data[i].n, params->data[i].i);

    //get pAscModel
    pAscModel = probAscertainmentGivenModelHemiLookup(gsl_vector_get(betas, i),params->data[i].n, \
						      snpProbs, ascProbs);
    tmp = pSNP * pAsc / pAscModel;
    lik += log(tmp);
    gsl_vector_free(snpProbs);
  }
  gsl_matrix_free(ascProbs);
  return -lik;
}

max_stars_repo_name: andrewkern/segSiteHMM

max_issues_repo_name: andrewkern/segSiteHMM

max_forks_repo_name: andrewkern/segSiteHMM


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

/* Balance a general matrix by scaling the columns
 *
 * B =  A D
 *
 * where D is a diagonal matrix
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

#include "gsl_linalg.h"

int
gsl_linalg_balance_columns (gsl_matrix * A, gsl_vector * D)
{
  const size_t N = A->size2;
  size_t j;

  if (D->size != A->size2)
    {
      GSL_ERROR("length of D must match second dimension of A", GSL_EINVAL);
    }

  gsl_vector_set_all (D, 1.0);

  for (j = 0; j < N; j++)
    {
      gsl_vector_view A_j = gsl_matrix_column (A, j);

      double s = gsl_blas_dasum(&A_j.vector);

      double f = 1.0;

      if (s == 0.0)
        {
          gsl_vector_set (D, j, f);
          continue;
        }

      while (s > 1.0)
        {
          s /= 2.0;
          f *= 2.0;
        }

      while (s < 0.5)
        {
          s *= 2.0;
          f /= 2.0;
        }

      gsl_vector_set (D, j, f);

      if (f != 1.0)
        {
          gsl_blas_dscal(1.0/f, &A_j.vector);
        }
    }

  return GSL_SUCCESS;
}


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

#ifndef SRC_ADAPTOR_BLAS_H_
#define SRC_ADAPTOR_BLAS_H_
#include <assert.h>

#include <cblas.h>

namespace juzhen {

template<typename T>
void gemm(
    const int M, const int N,
    const int K, const T *A, const int lda, const T *B,
    const int ldb, T *c, const int ldc) {
  assert(0);  // always fails
}

template<> inline
void gemm<float>(
    const int M, const int N,
    const int K, const float *A, const int lda, const float *B,
    const int ldb, float *c, const int ldc) {
  cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, K, 1.,
              A, lda, B, ldb, 0., c, ldc);
}

template<> inline
void gemm<double>(
    const int M, const int N,
    const int K, const double *A, const int lda, const double *B,
    const int ldb, double *c, const int ldc) {
  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, K, 1.,
              A, lda, B, ldb, 0., c, ldc);
}

template<> inline
void gemm<CS>(
    const int M, const int N,
    const int K, const CS  *A, const int lda, const CS *B,
    const int ldb, CS  *c, const int ldc) {
  CS alpha(1., 0.);
  CS beta(0., 0.);
  cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, K, &alpha,
              A, lda, B, ldb, &beta, c, ldc);
}

template<> inline
void gemm<CD>(
    const int M, const int N,
    const int K, const CD  *A, const int lda, const CD  *B,
    const int ldb, CD  *c, const int ldc) {
  CD alpha(1., 0.);
  CD beta(0., 0.);
  cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, K, &alpha,
              A, lda, B, ldb, &beta, c, ldc);
}

template<typename T>
int geev(
    char nl, char nr, const int n,
    T *a, const int lda, CD * w, T *vl, const int ldvl,
    T *vr, const int ldvr) {
  assert(0);  // always fails
}

/*
 * Linear solver
 */
template<typename T>
int gesv(
    const int n, const int nrhs,
    T *a, const int lda, T *b, const int ldb) {
  assert(0);  // always fails
}

/*
 * Matrix inversion
 */
template<typename T>
int matrix_inverse(
    const int m, const int n, T *a, const int lda) {
  assert(0);  // always fails
}

/*
 * Matrix determinant
 */
template<typename T>
T matrix_determinant(const int m, T *a) {
  assert(0);  // always fails
}
}
#endif  // SRC_ADAPTOR_BLAS_H_

max_stars_repo_name: huichen/juzhen

max_issues_repo_name: huichen/juzhen

max_forks_repo_name: huichen/juzhen


text:
#include <nlopt.h>

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

/****************************************************************************/
/* test function from M. J. Box, "A new method of constrained optimization
   and a comparison with other methods," Computer J. 8 (1), 42-52 (1965) */

int testfuncs_verbose = 1;
int testfuncs_counter = 0;

static double testfuncs_status(int n, const double *x, double f)
{
    ++testfuncs_counter;
    if (testfuncs_verbose) {
        int i;
        printf("f_%d (%g", testfuncs_counter, x[0]);
        for (i = 1; i < n; ++i)
            printf(", %g", x[i]);
        printf(") = %g\n", f);
    }
    return f;
}

#define RETURN(f) return testfuncs_status(n, x, f);

static const double k1 = -145421.402, k2 = 2931.1506, k3 = -40.427932, k4 = 5106.192,
    k5 = 15711.36, k6 = -161622.577, k7 = 4176.15328, k8 = 2.8260078,
    k9 = 9200.476, k10 = 13160.295, k11 = -21686.9194, k12 = 123.56928,
    k13 = -21.1188894, k14 = 706.834, k15 = 2898.573, k16 = 28298.388,
    k17 = 60.81096, k18 = 31.242116, k19 = 329.574, k20 = -2882.082,
    k21 = 74095.3845, k22 = -306.262544, k23 = 16.243649, k24 = -3094.252,
    k25 = -5566.2628, k26 = -26237, k27 = 99, k28 = -0.42, k29 = 1300, k30 = 2100, k31 = 925548.252, k32 = -61968.8432, k33 = 23.3088196, k34 = -27096.648, k35 = -50843.766;

static double box(int n, const double *x, double *grad, void *data)
{
    double x1 = x[0], x2 = x[1], x3 = x[2], x4 = x[3], x5 = x[4];
    double b, x6, y1, y2, y3, y4, u;
    const double a0 = 9, a1 = 15, a2 = 50, a3 = 9.583, a4 = 20, a5 = 15, a6 = 6, a7 = 0.75;
    b = x2 + 0.01 * x3;
    x6 = (k1 + k2 * x2 + k3 * x3 + k4 * x4 + k5 * x5) * x1;
    y1 = k6 + k7 * x2 + k8 * x3 + k9 * x4 + k10 * x5;
    y2 = k11 + k12 * x2 + k13 * x3 + k14 * x4 + k15 * x5;
    y3 = k16 + k17 * x2 + k18 * x3 + k19 * x4 + k20 * x5;
    y4 = k21 + k22 * x2 + k23 * x3 + k24 * x4 + k25 * x5;
    u = a2 * y1 + a3 * y2 + a4 * y3 + a5 * y4 + 7840 * a6 - 100000 * a0 - 50800 * b * a7 + k31 + k32 * x2 + k33 * x3 + k34 * x4 + k35 * x5;
    if (grad) {
        int i;
        grad[0] = u + a1 * (k1 + k2 * x2 + k3 * x3 + k4 * x4 + k5 * x5);
        grad[1] = x1 * (a2 * k7 + a3 * k12 + a4 * k17 + a5 * k22 - 50800 * a7 + k32) + a1 * (k2 * x1);
        grad[2] = x1 * (a2 * k8 + a3 * k13 + a4 * k18 + a5 * k23 - 50800 * a7 * 0.01) + a1 * x1 * k3;
        grad[3] = x1 * (a2 * k9 + a3 * k14 + a4 * k19 + a5 * k24) + a1 * x1 * k4;
        grad[4] = x1 * (a2 * k10 + a3 * k15 + a4 * k20 + a5 * k25) + a1 * x1 * k5;
        for (i = 0; i < 5; ++i)
            grad[i] = -grad[i];
    }
    RETURN(-(u * x1 - 24345 + a1 * x6));
}

static double box_constraint(int n, const double *x, double *grad, void *data)
{
    int which_constraint = *((int *) data);
    double x1 = x[0], x2 = x[1], x3 = x[2], x4 = x[3], x5 = x[4];
    double x6, y1, y2, y3, x7, x8;
    int i;

    x6 = (k1 + k2 * x2 + k3 * x3 + k4 * x4 + k5 * x5) * x1;
    y1 = k6 + k7 * x2 + k8 * x3 + k9 * x4 + k10 * x5;
    y2 = k11 + k12 * x2 + k13 * x3 + k14 * x4 + k15 * x5;
    y3 = k16 + k17 * x2 + k18 * x3 + k19 * x4 + k20 * x5;
    x7 = (y1 + y2 + y3) * x1;
    x8 = (k26 + k27 * x2 + k28 * x3 + k29 * x4 + k30 * x5) * x1 + x6 + x7;

    if (grad) {
        grad[0] = grad[1] = grad[2] = grad[3] = grad[4] = 0;

        if (which_constraint != 2 && which_constraint != -2) {
            /* grad x6 */
            grad[0] += (k1 + k2 * x2 + k3 * x3 + k4 * x4 + k5 * x5);
            grad[1] += x1 * k2;
            grad[2] += x1 * k3;
            grad[3] += x1 * k4;
            grad[4] += x1 * k5;
        }
        if (which_constraint != 1 && which_constraint != -1) {
            /* grad x7 */
            grad[0] += (y1 + y2 + y3);
            grad[1] += x1 * (k7 + k12 + k17);
            grad[2] += x1 * (k8 + k13 + k18);
            grad[3] += x1 * (k9 + k14 + k19);
            grad[4] += x1 * (k10 + k15 + k20);
        }
        if (which_constraint == 3 || which_constraint == -3) {
            /* grad (x8 - x6 - x7) */
            grad[0] += k26 + k27 * x2 + k28 * x3 + k29 * x4 + k30 * x5;
            grad[1] += x1 * k27;
            grad[2] += x1 * k28;
            grad[3] += x1 * k29;
            grad[4] += x1 * k30;
        }
    }

    if (which_constraint == -1) {
        if (grad)
            for (i = 0; i < 5; ++i)
                grad[i] = -grad[i];
        return -x6;
    } else if (which_constraint == 1) {
        return x6 - 294000;
    } else if (which_constraint == -2) {
        if (grad)
            for (i = 0; i < 5; ++i)
                grad[i] = -grad[i];
        return -x7;
    } else if (which_constraint == 2) {
        return x7 - 294000;
    } else if (which_constraint == -3) {
        if (grad)
            for (i = 0; i < 5; ++i)
                grad[i] = -grad[i];
        return -x8;
    } else if (which_constraint == 3) {
        return x8 - 277200;
    }
    return 0;
}

int main(void)
{
    const double box_lb[5] = { 0, 1.2, 20, 9, 6.5 };
    const double box_ub[5] = { HUGE_VAL, 2.4, 60, 9.3, 7.0 };
    const double box_xmin[5] = { 4.53743, 2.4, 60, 9.3, 7.0 };

    int cdata[6] = { -1, 1, -2, 2, -3, 3 };
    double x[5] = { 2.52, 2, 37.5, 9.25, 6.8 };
    double minf;

    nlopt_minimize_constrained(NLOPT_LN_COBYLA, 5, box, NULL, 6, box_constraint, cdata, sizeof(int), box_lb, box_ub, x, &minf, -HUGE_VAL, 1e-10, 0, 0, NULL, 0, 0);
    printf("found f(%g,%g,%g,%g,%g) = %0.8f after %d iters\n", x[0], x[1], x[2], x[3], x[4], minf, testfuncs_counter);
    return 0;
}

max_stars_repo_name: rockandsalt/conan-center-index

max_issues_repo_name: rockandsalt/conan-center-index

max_forks_repo_name: rockandsalt/conan-center-index


text:
/*********************************************************************************/
/* Matrix product program with MPI on a virtual ring of processors               */
/* S. Vialle - October 2014                                                      */
/*********************************************************************************/

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

#include "main.h"
#include "Calcul.h"

/*-------------------------------------------------------------------------------*/
/* Sequential product of local matrixes (optimized seq product).                 */
/*-------------------------------------------------------------------------------*/
void OneLocalProduct(int step)
{
  int OffsetStepLigC;
  int i, j, k;

  // Set the number of OpenMP threads inside parallel regions
  omp_set_num_threads(NbThreads);

  // Compute the current step offset, in the MPI program, to access right C lines
  OffsetStepLigC = (Me * LOCAL_SIZE + step * LOCAL_SIZE) % SIZE;

  switch (KernelId) {

  // Kernel 0 : Optimized code implemented by an application developer
  case 0 :
    for (i = 0; i < LOCAL_SIZE; i++) {
      for (j = 0; j < LOCAL_SIZE; j++) {
        double accu[8];
        accu[0] = accu[1] = accu[2] = accu[3] =  accu[4] =  accu[5] =  accu[6] =  accu[7] = 0.0;
        for (k = 0; k < (SIZE/8)*8; k += 8) {
           accu[0] += A_Slice[i][k+0] * TB_Slice[j][k+0];
           accu[1] += A_Slice[i][k+1] * TB_Slice[j][k+1];
           accu[2] += A_Slice[i][k+2] * TB_Slice[j][k+2];
           accu[3] += A_Slice[i][k+3] * TB_Slice[j][k+3];
           accu[4] += A_Slice[i][k+4] * TB_Slice[j][k+4];
           accu[5] += A_Slice[i][k+5] * TB_Slice[j][k+5];
           accu[6] += A_Slice[i][k+6] * TB_Slice[j][k+6];
           accu[7] += A_Slice[i][k+7] * TB_Slice[j][k+7];
        }
        for (k = (SIZE/8)*8; k < SIZE; k++) {
           accu[0] += A_Slice[i][k] * TB_Slice[j][k];
        }
        C_Slice[i+OffsetStepLigC][j] = accu[0] + accu[1] + accu[2] + accu[3] +
                                       accu[4] + accu[5] + accu[6] + accu[7];
      }
    }
    break;

  // Kernel 1 : Very optimized computing kernel implemented in a HPC library
  case 1 :
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                LOCAL_SIZE, LOCAL_SIZE, SIZE,
                1.0, &A_Slice[0][0], SIZE,
                &B_Slice[0][0], LOCAL_SIZE,
                0.0, &C_Slice[OffsetStepLigC][0], LOCAL_SIZE);
    break;

  default :
    fprintf(stderr,"Error: kernel %d not implemented!\n",KernelId);
    exit(EXIT_FAILURE);
    break;
  }
}

max_stars_repo_name: pallamidessi/Courses

max_issues_repo_name: pallamidessi/Courses

max_forks_repo_name: pallamidessi/Courses


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_blas.h>

static int
hilbert_matrix(gsl_matrix * m)
{
  const size_t N = m->size1;
  const size_t M = m->size2;
  size_t i, j;

  for (i = 0; i < N; i++)
    {
      for (j = 0; j < M; j++)
        {
          gsl_matrix_set(m, i, j, 1.0/(i+j+1.0));
        }
    }

  return GSL_SUCCESS;
}

int
main()
{
  const size_t n = 10; /* number of observations */
  const size_t p = 8;  /* number of model parameters */
  size_t i;
  gsl_matrix *X = gsl_matrix_alloc(n, p);
  gsl_vector *y = gsl_vector_alloc(n);

  /* construct Hilbert matrix and rhs vector */
  hilbert_matrix(X);

  {
    double val = 1.0;
    for (i = 0; i < n; ++i)
      {
        gsl_vector_set(y, i, val);
        val *= -1.0;
      }
  }

  {
    const size_t npoints = 200;                   /* number of points on L-curve and GCV curve */
    gsl_multifit_linear_workspace *w =
      gsl_multifit_linear_alloc(n, p);
    gsl_vector *c = gsl_vector_alloc(p);          /* OLS solution */
    gsl_vector *c_lcurve = gsl_vector_alloc(p);   /* regularized solution (L-curve) */
    gsl_vector *c_gcv = gsl_vector_alloc(p);      /* regularized solution (GCV) */
    gsl_vector *reg_param = gsl_vector_alloc(npoints);
    gsl_vector *rho = gsl_vector_alloc(npoints);  /* residual norms */
    gsl_vector *eta = gsl_vector_alloc(npoints);  /* solution norms */
    gsl_vector *G = gsl_vector_alloc(npoints);    /* GCV function values */
    double lambda_l;                              /* optimal regularization parameter (L-curve) */
    double lambda_gcv;                            /* optimal regularization parameter (GCV) */
    double G_gcv;                                 /* G(lambda_gcv) */
    size_t reg_idx;                               /* index of optimal lambda */
    double rcond;                                 /* reciprocal condition number of X */
    double chisq, rnorm, snorm;

    /* compute SVD of X */
    gsl_multifit_linear_svd(X, w);

    rcond = gsl_multifit_linear_rcond(w);
    fprintf(stderr, "matrix condition number = %e\n", 1.0 / rcond);

    /* unregularized (standard) least squares fit, lambda = 0 */
    gsl_multifit_linear_solve(0.0, X, y, c, &rnorm, &snorm, w);
    chisq = pow(rnorm, 2.0);

    fprintf(stderr, "\n=== Unregularized fit ===\n");
    fprintf(stderr, "residual norm = %g\n", rnorm);
    fprintf(stderr, "solution norm = %g\n", snorm);
    fprintf(stderr, "chisq/dof = %g\n", chisq / (n - p));

    /* calculate L-curve and find its corner */
    gsl_multifit_linear_lcurve(y, reg_param, rho, eta, w);
    gsl_multifit_linear_lcorner(rho, eta, &reg_idx);

    /* store optimal regularization parameter */
    lambda_l = gsl_vector_get(reg_param, reg_idx);

    /* regularize with lambda_l */
    gsl_multifit_linear_solve(lambda_l, X, y, c_lcurve, &rnorm, &snorm, w);
    chisq = pow(rnorm, 2.0) + pow(lambda_l * snorm, 2.0);

    fprintf(stderr, "\n=== Regularized fit (L-curve) ===\n");
    fprintf(stderr, "optimal lambda: %g\n", lambda_l);
    fprintf(stderr, "residual norm = %g\n", rnorm);
    fprintf(stderr, "solution norm = %g\n", snorm);
    fprintf(stderr, "chisq/dof = %g\n", chisq / (n - p));

    /* calculate GCV curve and find its minimum */
    gsl_multifit_linear_gcv(y, reg_param, G, &lambda_gcv, &G_gcv, w);

    /* regularize with lambda_gcv */
    gsl_multifit_linear_solve(lambda_gcv, X, y, c_gcv, &rnorm, &snorm, w);
    chisq = pow(rnorm, 2.0) + pow(lambda_gcv * snorm, 2.0);

    fprintf(stderr, "\n=== Regularized fit (GCV) ===\n");
    fprintf(stderr, "optimal lambda: %g\n", lambda_gcv);
    fprintf(stderr, "residual norm = %g\n", rnorm);
    fprintf(stderr, "solution norm = %g\n", snorm);
    fprintf(stderr, "chisq/dof = %g\n", chisq / (n - p));

    /* output L-curve and GCV curve */
    for (i = 0; i < npoints; ++i)
      {
        printf("%e %e %e %e\n",
               gsl_vector_get(reg_param, i),
               gsl_vector_get(rho, i),
               gsl_vector_get(eta, i),
               gsl_vector_get(G, i));
      }

    /* output L-curve corner point */
    printf("\n\n%f %f\n",
           gsl_vector_get(rho, reg_idx),
           gsl_vector_get(eta, reg_idx));

    /* output GCV curve corner minimum */
    printf("\n\n%e %e\n",
           lambda_gcv,
           G_gcv);

    gsl_multifit_linear_free(w);
    gsl_vector_free(c);
    gsl_vector_free(c_lcurve);
    gsl_vector_free(reg_param);
    gsl_vector_free(rho);
    gsl_vector_free(eta);
    gsl_vector_free(G);
  }

  gsl_matrix_free(X);
  gsl_vector_free(y);

  return 0;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <config.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>

#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_complex_float.h>

#define BASE_DOUBLE
#include "templates_on.h"
#include "bitreverse.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "bitreverse.c"
#include "templates_off.h"
#undef  BASE_FLOAT

#include "factorize.c"

#define BASE_DOUBLE
#include "templates_on.h"
#include "c_init.c"
#include "c_main.c"
#include "c_pass_2.c"
#include "c_pass_3.c"
#include "c_pass_4.c"
#include "c_pass_5.c"
#include "c_pass_6.c"
#include "c_pass_7.c"
#include "c_pass_n.c"
#include "c_radix2.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "c_init.c"
#include "c_main.c"
#include "c_pass_2.c"
#include "c_pass_3.c"
#include "c_pass_4.c"
#include "c_pass_5.c"
#include "c_pass_6.c"
#include "c_pass_7.c"
#include "c_pass_n.c"
#include "c_radix2.c"
#include "templates_off.h"
#undef  BASE_FLOAT

#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_fft_halfcomplex_float.h>

#define BASE_DOUBLE
#include "templates_on.h"
#include "hc_init.c"
#include "hc_main.c"
#include "hc_pass_2.c"
#include "hc_pass_3.c"
#include "hc_pass_4.c"
#include "hc_pass_5.c"
#include "hc_pass_n.c"
#include "hc_radix2.c"
#include "hc_unpack.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "hc_init.c"
#include "hc_main.c"
#include "hc_pass_2.c"
#include "hc_pass_3.c"
#include "hc_pass_4.c"
#include "hc_pass_5.c"
#include "hc_pass_n.c"
#include "hc_radix2.c"
#include "hc_unpack.c"
#include "templates_off.h"
#undef  BASE_FLOAT

#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_real_float.h>

#define BASE_DOUBLE
#include "templates_on.h"
#include "real_init.c"
#include "real_main.c"
#include "real_pass_2.c"
#include "real_pass_3.c"
#include "real_pass_4.c"
#include "real_pass_5.c"
#include "real_pass_n.c"
#include "real_radix2.c"
#include "real_unpack.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "real_init.c"
#include "real_main.c"
#include "real_pass_2.c"
#include "real_pass_3.c"
#include "real_pass_4.c"
#include "real_pass_5.c"
#include "real_pass_n.c"
#include "real_radix2.c"
#include "real_unpack.c"
#include "templates_off.h"
#undef  BASE_FLOAT

max_stars_repo_name: manggoguy/parsec-modified