text:
/*System includes*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

/*GSL includes*/
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_multimin.h>
#include <pthread.h>

/*User includes*/
#include "ERarefaction.h"

static char *usage[] = {"ERarefaction - Determines an emprical rarefraction curve by resampling\n",
                        "Required parameters:\n",
			"   -in   filename      parameter file  \n",
                        "Optional:\n",
			"   -sample int         how frequently to sample\n",
			"   -v                  verbose\n"};

static int  nLines   = 6;

static int  verbose  = FALSE;

double expectedDiversity(int n, t_Data *ptData)
{
  int i = 0;
  double dDenom = gsl_sf_lnchoose(ptData->nL, n);
  double dSum = 0.0;

  for(i = 0; i < ptData->nNA; i++){
    int nA = ptData->aanAbund[i][0];
    if(ptData->nL - nA >= n){

      double dNumer = gsl_sf_lnchoose(ptData->nL - nA, n);

      dSum += ((double) ptData->aanAbund[i][1])*exp(dNumer - dDenom);
    }

  }

  return ((double) ptData->nS) - dSum;
}

int main(int argc, char* argv[])
{
  int  i = 0, j = 0, r = 1;
  t_Params tParams;
  t_Data   tData;

  /*get command line params*/
  getCommandLineParams(&tParams, argc, argv);

  /*read in abundance distribution*/
  readAbundanceData(tParams.szInputFile, &tData);

  for(r = 1; r < tData.nL; r++){
    if(r % tParams.nSample == 0){
      printf("%d %f\n",r, expectedDiversity(r, &tData));
    }
  }

  freeAbundanceData(&tData);

  exit(EXIT_SUCCESS);
}

void writeUsage(FILE* ofp)
{
  int i = 0;
  char *line;

  for(i = 0; i < nLines; i++){
    line = usage[i];
    fputs(line,ofp);
  }
}

char *extractParameter(int argc, char **argv, char *param,int when)
{
  int i = 0;

  while((i < argc) && (strcmp(param,argv[i]))){
    i++;
  }

  if(i < argc - 1){
    return(argv[i + 1]);
  }

  if((i == argc - 1) && (when == OPTION)){
    return "";
  }

  if(when == ALWAYS){
    fprintf(stdout,"Can't find asked option %s\n",param);
  }

  return (char *) NULL;
}

void getCommandLineParams(t_Params *ptParams,int argc,char *argv[])
{
  char *szTemp = NULL;
  char *cError = NULL;

  /*get parameter file name*/
  ptParams->szInputFile  = extractParameter(argc,argv, INPUT_FILE,ALWAYS);
  if(ptParams->szInputFile == NULL)
    goto error;

  /*get out file stub*/
  szTemp  = extractParameter(argc,argv,SAMPLE,OPTION);
  if(szTemp != NULL){
    ptParams->nSample = strtol(szTemp,&cError,10);
    if(*cError != '\0'){
      goto error;
    }
  }
  else{
    ptParams->nSample = 1;
  }

  /*verbosity*/
  szTemp = extractParameter(argc, argv, VERBOSE, OPTION);
  if(szTemp != NULL){
    verbose = TRUE;
  }

  return;

 error:
  writeUsage(stdout);
  exit(EXIT_FAILURE);
}

void readAbundanceData(const char *szFile, t_Data *ptData)
{
  int **aanAbund = NULL;
  int  i = 0, nNA = 0, nA = 0, nC = 0;
  int  nL = 0, nS = 0;
  char szLine[MAX_LINE_LENGTH];
  FILE* ifp = NULL;

  ifp = fopen(szFile, "r");

  if(ifp){
    char* szTok   = NULL;
    char* pcError = NULL;

    fgets(szLine, MAX_LINE_LENGTH, ifp);

    szTok = strtok(szLine, DELIM);

    nNA = strtol(szTok,&pcError,10);
    if(*pcError != '\0'){
      goto formatError;
    }

    aanAbund = (int **) malloc(nNA*sizeof(int*));

    for(i = 0; i < nNA; i++){
      aanAbund[i] = (int *) malloc(sizeof(int)*2);

      fgets(szLine, MAX_LINE_LENGTH, ifp);

      szTok = strtok(szLine, DELIM);

      nA = strtol(szTok,&pcError,10);
      if(*pcError != '\0'){
	goto formatError;
      }

      szTok = strtok(NULL, DELIM);

      nC = strtol(szTok,&pcError,10);
      if(*pcError != '\0'){
	goto formatError;
      }

      nS += nC;
      nL += nC*nA;

      aanAbund[i][0]  = nA;
      aanAbund[i][1]  = nC;
    }
  }
  else{
    fprintf(stderr, "Failed to open abundance data file %s aborting\n", szFile);
    fflush(stderr);
    exit(EXIT_FAILURE);
  }

  ptData->nS          = nS;
  ptData->nL          = nL;
  ptData->aanAbund    = aanAbund;
  ptData->nNA         = nNA;
  return;

 formatError:
  fprintf(stderr, "Incorrectly formatted abundance data file\n");
  fflush(stderr);
  exit(EXIT_FAILURE);
}

int compare_doubles(const void* a, const void* b)
{
  double* arg1 = (double *) a;
  double* arg2 = (double *) b;
  if( *arg1 < *arg2 ) return -1;
  else if( *arg1 == *arg2 ) return 0;
  else return 1;
}

void freeAbundanceData(t_Data *ptData)
{
  int i = 0;

  for(i = 0; i < ptData->nNA; i++){
    free(ptData->aanAbund[i]);
  }
  free(ptData->aanAbund);
}

max_stars_repo_name: chrisquince/DiversityEstimates

max_issues_repo_name: chrisquince/DiversityEstimates

max_forks_repo_name: chrisquince/DiversityEstimates


text:
#ifndef COLLECTDATA_H
#define COLLECTDATA_H

#include <cstdlib>
#include <ctime>
#include <vector>

#include <yarp/sig/all.h>
#include <yarp/os/all.h>
#include <yarp/dev/all.h>
#include <yarp/math/Math.h>

#include <iCub/iKin/iKinFwd.h>
#include <iCub/ctrl/math.h>

#include <gsl/gsl_math.h>

#include "collectData_IDL.h"

class collectData : public yarp::os::RFModule, public collectData_IDL
{
protected:
    yarp::os::RpcServer             rpcPort;
    double                          period;
    std::string                     name;
    yarp::os::Stamp                 ts;

    yarp::os::RpcClient     rpcToMotorBabbling;
    yarp::os::RpcClient     rpcToLearning;

    bool    sendCmdBabbling();
    bool    sendCmdLearningPredict();
    bool    sendCmdLearningMove(const int &id, const std::string &partName);

    std::vector<int> reprTaxelsForearm_sim =
    {3, 15, 27, 39, 51, 75, 87, 183,
    207, 255, 291, 303, 315, 339, 351}; //23 taxels

    std::vector<int> reprTaxelsHandL =
    {3, 15, 27, 39, 51, 99, 101, 109, 122, 134};

    std::vector<int> reprTaxelsHandR =
    {3, 15, 39, 51, 101, 118, 137};

public:
    bool    configure(yarp::os::ResourceFinder &rf);
    bool    interruptModule();
    bool    close();
    bool    attach(yarp::os::RpcServer &source);
    double  getPeriod();
    bool    updateModule();

    //Thrift

    bool babble_then_reach(int32_t _id, const std::string &_part, int32_t _nb_object=1)
    {
        if (_nb_object>=1)
        {
            double delayBetweenReach = 15.0;
            bool ok=false;
            for (int i=0; i<_nb_object; i++)
            {
                ok = sendCmdBabbling();
                yarp::os::Time::delay(7);
                ok = ok & sendCmdLearningPredict();
                yarp::os::Time::delay(3);

                if (_id==-1)    // reach all taxel of part
                {
                    if (_part=="b" or _part=="both")
                    {
                        for (int i=0; i<reprTaxelsHandR.size(); i++)
                        {
                            ok = ok & sendCmdLearningMove(reprTaxelsHandR[i], "h");
                            yarp::os::Time::delay(delayBetweenReach);
                        }
                        for (int i=0; i<reprTaxelsForearm_sim.size(); i++)
                        {
                            ok = ok & sendCmdLearningMove(reprTaxelsForearm_sim[i], "a");
                            yarp::os::Time::delay(delayBetweenReach);
                        }
                    }
                    else if (_part=="h" or _part=="hand")
                    {
                        for (int i=0; i<reprTaxelsHandR.size(); i++)
                        {
                            ok = ok & sendCmdLearningMove(reprTaxelsHandR[i], "h");
                            yarp::os::Time::delay(delayBetweenReach);
                        }
                    }
                    else if (_part=="a" or _part=="arm")
                    {
                        for (int i=0; i<reprTaxelsForearm_sim.size(); i++)
                        {
                            ok = ok & sendCmdLearningMove(reprTaxelsForearm_sim[i], "a");
                            yarp::os::Time::delay(delayBetweenReach);
                        }
                    }
                }
                else
                {
                    ok = ok & sendCmdLearningMove(_id,_part);
                    yarp::os::Time::delay(delayBetweenReach);
                }

                if (!ok)
                    return ok;

            }
            return ok;
        }
        else
            return false;

    }

};

#endif // COLLECTDATA_H

ext: h

max_stars_repo_name: robotology/visuomotor-learning

max_issues_repo_name: robotology/visuomotor-learning

max_forks_repo_name: robotology/visuomotor-learning


text:
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

#include <assert.h>
#include <math.h>
#include <string.h>

#include "egsl.h"
#include "egsl_imp.h"
#define MAX_VALS 1024
#define MAX_CONTEXTS 1024

/*#define INVALID (val_from_context_and_index(1000,1000))*/

struct egsl_variable {
	gsl_matrix * gsl_m;
};

struct egsl_context {
	char name[256];
	int nallocated;
	int nvars;
	struct egsl_variable vars[MAX_VALS];
};

/* Current context */
int cid=0;
/* Maximum number of contexts */
int max_cid = 0;
static struct egsl_context egsl_contexts[MAX_CONTEXTS];

int egsl_first_time = 1;

int egsl_total_allocations = 0;
int egsl_cache_hits = 0;

void egsl_error(void) {
	/* TODO: better handling of errors */

	egsl_print_stats();

	assert(0);
}

val assemble_val(int cid, int index, gsl_matrix*m) {
	val v;
		v.cid=cid;
		v.index=index;
		v.gslm = m;
	return v;
}

int its_context(val v) {
	return v.cid;
}

int its_var_index(val v) {
	return v.index;
}

#if 0
inline void check_valid_val(val v) { int i = v.cid; v.cid=i;}
#else
void check_valid_val(val v) {
	int context = its_context(v);
	if(context>cid) {
		fprintf(stderr, "Val is from invalid context (%d>%d)\n",context,cid);
		egsl_error();
	}
	int var_index = its_var_index(v);
	if(var_index >= egsl_contexts[context].nvars) {
		fprintf(stderr, "Val is invalid  (%d>%d)\n",var_index,
			egsl_contexts[context].nvars);
		egsl_error();
	}
}
#endif

gsl_matrix * egsl_gslm(val v) {
	check_valid_val(v);
	return v.gslm;
}

void egsl_push() { egsl_push_named("unnamed context"); }
void egsl_pop() { egsl_pop_named("unnamed context"); }

void egsl_push_named(const char*name) {
	if(egsl_first_time) {
		int c;
		for(c=0;c<MAX_CONTEXTS;c++) {
			egsl_contexts[c].nallocated = 0;
			egsl_contexts[c].nvars = 0;
			sprintf(egsl_contexts[c].name, "not yet used");
		}
		egsl_first_time  = 0;
	}
	cid++;

	if(cid >= MAX_CONTEXTS) {
		fprintf(stderr, "egsl: maximum number of contexts reached \n");
		egsl_print_stats();
		assert(0);
	}

	if(max_cid < cid) max_cid = cid;

	if(name != 0)
		sprintf(egsl_contexts[cid].name, "%s", name);
	else
		sprintf(egsl_contexts[cid].name, "Unnamed context");
}

void egsl_pop_named(const char*name) {
	assert(cid>=0);/*, "No egsl_push before?"); */
	if(name != 0) {
		if(strcmp(name, egsl_contexts[cid].name)) {
			fprintf(stderr, "egsl: context mismatch. You want to pop '%s', you are still at '%s'\n",
				name, egsl_contexts[cid].name);
			egsl_print_stats();
			assert(0);
		}
	}

	egsl_contexts[cid].nvars = 0;
	sprintf(egsl_contexts[cid].name, "Popped");
	cid--;
}
/*
void egsl_pop_check(int assumed) {
	if(assumed != cid) {
		fprintf(stderr, "egsl: You think you are in context %d while you are %d. \n", assumed, cid);
		if(assumed < cid)
			fprintf(stderr, "     It seems you miss %d egsl_pop() somewhere.\n", - assumed + cid);
		else
			fprintf(stderr, "     It seems you did %d egsl_pop() more.\n",  + assumed - cid);
		egsl_print_stats();
	}
	assert(cid>=0);
	egsl_contexts[cid].nvars = 0;
	cid--;
}*/

void egsl_print_stats() {
	fprintf(stderr, "egsl: total allocations: %d   cache hits: %d\n",
		egsl_total_allocations, egsl_cache_hits);
/*	printf("egsl: sizeof(val) = %d\n",(int)sizeof(val)); */
	int c; for(c=0;c<=max_cid&&c<MAX_CONTEXTS;c++) {
	/*	printf("egsl: context #%d\n ",c); */
/*	 	if(0==egsl_contexts[c].nallocated) break; */
		fprintf(stderr, "egsl: context #%d allocations: %d active: %d name: '%s' \n",
			c,	egsl_contexts[c].nallocated, egsl_contexts[c].nvars, egsl_contexts[c].name);
	}
}

val egsl_alloc(size_t rows, size_t columns) {
	struct egsl_context*c = egsl_contexts+cid;

/*	if(cid<3)
	printf("Alloc cid=%d nvars=%d nalloc=%d\n",cid,c->nvars,c->nallocated); */

	if(c->nvars>=MAX_VALS) {
		fprintf(stderr,"Limit reached, in context %d, nvars is %d\n",cid,c->nvars);
		egsl_error();
	}
	int index = c->nvars;
	if(index<c->nallocated) {
		gsl_matrix*m = c->vars[index].gsl_m;
		if(m->size1 == rows && m->size2 == columns) {
			egsl_cache_hits++;
			c->nvars++;
			return assemble_val(cid,index,c->vars[index].gsl_m);
		} else {
			gsl_matrix_free(m);
			egsl_total_allocations++;
			c->vars[index].gsl_m = gsl_matrix_alloc((size_t)rows,(size_t)columns);
			c->nvars++;
			return assemble_val(cid,index,c->vars[index].gsl_m);
		}
	} else {
		egsl_total_allocations++;
		c->vars[index].gsl_m = gsl_matrix_alloc((size_t)rows,(size_t)columns);
		c->nvars++;
		c->nallocated++;
		return assemble_val(cid,index,c->vars[index].gsl_m);
	}
}

val egsl_alloc_in_context(int context, size_t rows, size_t columns) {
	struct egsl_context*c = egsl_contexts+context;

	if(c->nvars>=MAX_VALS) {
		fprintf(stderr,"Limit reached, in context %d, nvars is %d\n",context,c->nvars);
		egsl_error();
	}
	int index = c->nvars;
	if(index<c->nallocated) {
		gsl_matrix*m = c->vars[index].gsl_m;
		if(m->size1 == rows && m->size2 == columns) {
			egsl_cache_hits++;
			c->nvars++;
			return assemble_val(context,index,c->vars[index].gsl_m);
		} else {
			gsl_matrix_free(m);
			egsl_total_allocations++;
			c->vars[index].gsl_m = gsl_matrix_alloc((size_t)rows,(size_t)columns);
			c->nvars++;
			return assemble_val(context,index,c->vars[index].gsl_m);
		}
	} else {
		egsl_total_allocations++;
		c->vars[index].gsl_m = gsl_matrix_alloc((size_t)rows,(size_t)columns);
		c->nvars++;
		c->nallocated++;
		return assemble_val(context,index,c->vars[index].gsl_m);
	}
}

/** Creates a copy of v in the previous context. */
val egsl_promote(val v) {
	if(cid==0) {
		egsl_error();
	}

	gsl_matrix * m = egsl_gslm(v);
	val v2 = egsl_alloc_in_context(cid-1, m->size1, m->size2);
	gsl_matrix * m2 = egsl_gslm(v2);
	gsl_matrix_memcpy(m2, m);
	return v2;
}

void egsl_expect_size(val v, size_t rows, size_t cols) {
	gsl_matrix * m = egsl_gslm(v);

	int bad = (rows && (m->size1!=rows)) || (cols && (m->size2!=cols));
	if(bad) {
		fprintf(stderr, "Matrix size is %d,%d while I expect %d,%d",
			(int)m->size1,(int)m->size2,(int)rows,(int)cols);

		egsl_error();
	}
}

void egsl_print(const char*str, val v) {
	gsl_matrix * m = egsl_gslm(v);
	size_t i,j;
	int context = its_context(v);
	int var_index = its_var_index(v);
	fprintf(stderr, "%s =  (%d x %d)  context=%d index=%d\n",
		str,(int)m->size1,(int)m->size2,  context, var_index);

	for(i=0;i<m->size1;i++) {
		if(i==0)
			fprintf(stderr, " [ ");
		else
			fprintf(stderr, "   ");

		for(j=0;j<m->size2;j++)
			fprintf(stderr, "%f ", gsl_matrix_get(m,i,j));

		if(i==m->size1-1)
			fprintf(stderr, "] \n");
		else
			fprintf(stderr, "; \n");
	}
}

double* egsl_atmp(val v, size_t i, size_t j) {
	gsl_matrix * m = egsl_gslm(v);
	return gsl_matrix_ptr(m,(size_t)i,(size_t)j);
}

double egsl_norm(val v1){
	egsl_expect_size(v1, 0, 1);
	double n=0;
	size_t i;
	gsl_matrix * m = egsl_gslm(v1);
	for(i=0;i<m->size1;i++) {
		double v = gsl_matrix_get(m,i,0);
		n += v * v;
	}
	return sqrt(n);
}

double egsl_atv(val v1,  size_t i){
	return *egsl_atmp(v1, i, 0);
}

double egsl_atm(val v1, size_t i, size_t j){
	return *egsl_atmp(v1, i, j);
}

void egsl_free(void){
	int c;
	for(c=0;c<=max_cid;c++) {
		for(int i=egsl_contexts[c].nvars; i<egsl_contexts[c].nallocated; i++){
			gsl_matrix_free(egsl_contexts[c].vars[i].gsl_m);
		}
	egsl_contexts[c].nallocated = egsl_contexts[c].nvars;
	}
}

ext: c


text:
/**
 * @file DG_Element_2d.h
 * @Synopsis  This is the header file for the class which defines the DG_Element_2d class. This class contains all the
 * declarations which are necessary for defining a general 2d element. This class defines all the necessities for a
 * structured 2d DG element.
 * @author Kaushik Kulkarni
 * @version 1.0
 * @date 2017-02-18
 */

/* Copyright (c) 2016, Shivasubramanian Gopalakrishnan and Kaushik Kulkarni.
  All rights reserved.

   Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

     -Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
     -Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

     THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS  **AS IS** AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
     */

#ifndef DG_ELEMENT_2D_H_
#define DG_ELEMENT_2D_H_

#include <map>
#include <string>
#include <vector>
#include <functional>
#include <iostream>
#include <cblas.h>

using namespace std;

class DG_Element_2d {
private:

    int N; /// This represents the order of polynomial interpolation which is done while solving the DG equations.
    double x_start, x_end, y_start, y_end; /// Since, this is a structured rectangular element. These variables define the position and size of the element.
    double* massMatrix;
    double* derivativeMatrix_x; /// This is the matrix which is laid out in 1-D, this would help us to find t    he $\frac{d}{dx}$ of any term.
    double* derivativeMatrix_y; /// This is the matrix which is laid out in 1-D, this would help us to find t    he $\frac{d}{dy}$ of any term.
    double* fluxMatrix_top; /// This is the flux matrix for the top edge.
    double* fluxMatrix_right; // The Flux Matrix for the right edge.
    double* fluxMatrix_bottom; /// This would be the flux term for the the bottom edge.
    double* fluxMatrix_left; /// The Flux matrix for the left edge.
    double* inverseMassMatrix; /// The inverse of the mass matrix is also created only once in the field function. And just passed to each element.

    DG_Element_2d* topNeighbor;
    DG_Element_2d* rightNeighbor;
    DG_Element_2d* leftNeighbor;
    DG_Element_2d* bottomNeighbor;

public:

    double *X, *Y;
    map<string, double*> variable; /// This is the map which would contain all the variable that are needed by the problem.

    map<string, double**> boundaryTop;
    map<string, double**> boundaryRight;
    map<string, double**> boundaryLeft;
    map<string, double**> boundaryBottom;

    map<string, double**> neighboringTop;
    map<string, double**> neighboringRight;
    map<string, double**> neighboringBottom;
    map<string, double**> neighboringLeft;

    vector<string> boundaryVariables; /// This is the variable which stores the name of all the variables whose boundary and neighboring points are stored.

    DG_Element_2d(int _N, double x1, double y1, double x2, double y2);
    void addVariable_withBoundary(string v);
    void addVariable_withoutBoundary(string v);
    void initializeVariable(string v, function<double(double, double)> f);
    void setNeighboringElement(char type, DG_Element_2d* neighbor );
    void setVariableNeighbors(string v);

    // Functions for various operations on the variables.
    void delByDelX(string v, string vDash, string fluxType, string fluxVariable);

    void delByDelY(string v, string vDash, string fluxType, string fluxVariable);

    // Functions to set the various operator matrices.
    void setMassMatrix(double* m);
    void setInverseMassMatrix(double* im);
    void setderivateMatrix_x(double* d);
    void setderivateMatrix_y(double* d);
    void setTopFluxMatrix(double* f);
    void setRightFluxMatrix(double* f);
    void setLeftFluxMatrix(double* f);
    void setBottomFluxMatrix(double* f);

    // Functions to apply linear operations on the variables.
    void axpy(double a, string x, string y);
    void scal(double a, string x);
    void setFunctionsForVariables(string x, string y, function<double(double, double)>, string z);
    void setFunctionsForVariables(string w, string x, string y, function<double(double, double, double)>, string z);

    // Functions to do various other operations on the elements.
    double l2Norm(string v1, string v2);

};

#endif

max_stars_repo_name: Krishna-Dutt/DG-MHD-Solver

max_issues_repo_name: Krishna-Dutt/DG-MHD-Solver

max_forks_repo_name: Krishna-Dutt/DG-MHD-Solver


text:
/* interpolation/gsl_spline.h
 *
 * Copyright (C) 2001, 2007 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 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_SPLINE_H__
#define __GSL_SPLINE_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif
#include <stdlib.h>
#include <gsl/gsl_interp.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/* general interpolation object */
typedef struct {
  gsl_interp * interp;
  double  * x;
  double  * y;
  size_t  size;
} gsl_spline;

GSL_FUN gsl_spline *
gsl_spline_alloc(const gsl_interp_type * T, size_t size);

GSL_FUN int
gsl_spline_init(gsl_spline * spline, const double xa[], const double ya[], size_t size);

GSL_FUN const char * gsl_spline_name(const gsl_spline * spline);
GSL_FUN unsigned int gsl_spline_min_size(const gsl_spline * spline);

GSL_FUN int
gsl_spline_eval_e(const gsl_spline * spline, double x,
                  gsl_interp_accel * a, double * y);

GSL_FUN double
gsl_spline_eval(const gsl_spline * spline, double x, gsl_interp_accel * a);

GSL_FUN int
gsl_spline_eval_deriv_e(const gsl_spline * spline,
                        double x,
                        gsl_interp_accel * a,
                        double * y);

GSL_FUN double
gsl_spline_eval_deriv(const gsl_spline * spline,
                      double x,
                      gsl_interp_accel * a);

GSL_FUN int
gsl_spline_eval_deriv2_e(const gsl_spline * spline,
                         double x,
                         gsl_interp_accel * a,
                         double * y);

GSL_FUN double
gsl_spline_eval_deriv2(const gsl_spline * spline,
                       double x,
                       gsl_interp_accel * a);

GSL_FUN int
gsl_spline_eval_integ_e(const gsl_spline * spline,
                        double a, double b,
                        gsl_interp_accel * acc,
                        double * y);

GSL_FUN double
gsl_spline_eval_integ(const gsl_spline * spline,
                      double a, double b,
                      gsl_interp_accel * acc);

GSL_FUN void
gsl_spline_free(gsl_spline * spline);

__END_DECLS

#endif /* __GSL_INTERP_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
/*
 * Copyright 2017 RIS.
 * Authored by Ron Dahlgren <ron@red83.net>
 */
#ifndef ENTITIES_H_
#define ENTITIES_H_

#include "globals.h"
#include "rendering/rendering.h"
#include "animations/transform_animation.h"
#include "utility/c8math.h"

#include <gsl/gsl_matrix.h>

/**
 * This structure is used to manage all of the data associated with
 * a basic renderable entity. Comments below indicate which pointers
 * are optional and which are required. After calling entities.entity_alloc,
 * the resulting entity definition will be a basically zeroed but sane
 * copy of an internal prototype. Add some geometry and go to town.
 */
typedef struct entity_info {
  /* Required - specifies the material to use for rendering */
  shader_info shader;
  /* Required - should assume the projection matrix has already been set */
  render_func renderer;
  /* The animation state. Set current_cycle to actually use it */
  animation_state animation_state;
  /* Provides a reusable workspace for math operations */
  math_context math_context;
  /* Layout of these matrix arrays are row0, row1, row2, row3 */
  double transform[16]; /* Used in rendering */
  double reference_frame[16]; /* Animation solution is within this frame of ref */

  /* (OPT) Should claim additional resources needed by the entity */
  void (* initialize)(struct entity_info *);
  /* (OPT) Should release any claimed resources */
  void (* terminate)(struct entity_info *);

  /* (OPT) callback will be called (if set) before update */
  void (* pre_update)(double deltaT, struct entity_info *);
  /* (OPT) callback will be called (if set) after update */
  void (* post_update)(double deltaT, struct entity_info *);
  /* (OPT) Any additional data that may be needed */
  void* data;
} entity_info;

/**
 * Entity management operations are grouper here.
 */
typedef struct {
  /* Create an entity */
  entity_info* (* const entity_alloc)();
  /* Release an entity */
  void (__attribute__((nonnull(1)))* const entity_free)(entity_info* info);
  /* Add an entity */
  void (__attribute__((nonnull(1)))* const add_entity)(entity_info* info);
  /* Remove an entity */
  void (__attribute__((nonnull(1)))* const remove_entity)(entity_info* info);
  /* Conduct an update pass */
  void (* const update)(double deltaT);

  sd_result (* const initialize)(void);
  sd_result (* const release_resources)(void);
  sd_result (* const terminate)(void);
} entities_struct;

extern entities_struct const entities;
#endif


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

#include <config.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_spblas.h>
#include <gsl/gsl_splinalg.h>

/*
create_random_sparse()
  Create a random sparse matrix with approximately
M*N*density non-zero entries

Inputs: M       - number of rows
        N       - number of columns
        density - sparse density \in [0,1]
                  0 = no non-zero entries
                  1 = all m*n entries are filled
        r       - random number generator

Return: pointer to sparse matrix in triplet format (must be freed by caller)

Notes:
1) non-zero matrix entries are uniformly distributed in [0,1]
*/

static gsl_spmatrix *
create_random_sparse(const size_t M, const size_t N, const double density,
                     const gsl_rng *r)
{
  size_t nnzwanted = (size_t) floor(M * N * GSL_MIN(density, 1.0));
  gsl_spmatrix *m = gsl_spmatrix_alloc_nzmax(M, N,
                                             nnzwanted,
                                             GSL_SPMATRIX_TRIPLET);
  size_t i;

  /* set diagonal entries to try to ensure non-singularity */
  for (i = 0; i < GSL_MIN(M, N); ++i)
    {
      double x = gsl_rng_uniform(r);
      gsl_spmatrix_set(m, i, i, x);
    }

  while (gsl_spmatrix_nnz(m) < nnzwanted)
    {
      /* generate a random row and column */
      size_t i = gsl_rng_uniform(r) * M;
      size_t j = gsl_rng_uniform(r) * N;

      /* generate random m_{ij} and add it */
      double x = gsl_rng_uniform(r);
      gsl_spmatrix_set(m, i, j, x);
    }

  return m;
} /* create_random_sparse() */

static void
create_random_vector(gsl_vector *v, const gsl_rng *r)
{
  size_t i;

  for (i = 0; i < v->size; ++i)
    {
      double x = gsl_rng_uniform(r);
      gsl_vector_set(v, i, x);
    }
} /* create_random_vector() */

/*
test_poisson()
  Solve u''(x) = -pi^2 sin(pi*x), u(x) = sin(pi*x)
  epsrel is the relative error threshold with the exact solution
*/
static void
test_poisson(const size_t N, const double epsrel, const int compress)
{
  const gsl_splinalg_itersolve_type *T = gsl_splinalg_itersolve_gmres;
  const size_t n = N - 2;                     /* subtract 2 to exclude boundaries */
  const double h = 1.0 / (N - 1.0);           /* grid spacing */
  const double tol = 1.0e-9;
  const size_t max_iter = 10;
  size_t iter = 0;
  gsl_spmatrix *A = gsl_spmatrix_alloc(n ,n); /* triplet format */
  gsl_spmatrix *B;
  gsl_vector *b = gsl_vector_alloc(n);        /* right hand side vector */
  gsl_vector *u = gsl_vector_calloc(n);       /* solution vector, u0 = 0 */
  gsl_splinalg_itersolve *w = gsl_splinalg_itersolve_alloc(T, n, 0);
  const char *desc = gsl_splinalg_itersolve_name(w);
  size_t i;
  int status;

  /* construct the sparse matrix for the finite difference equation */

  /* first row of matrix */
  gsl_spmatrix_set(A, 0, 0, -2.0);
  gsl_spmatrix_set(A, 0, 1, 1.0);

  /* loop over interior grid points */
  for (i = 1; i < n - 1; ++i)
    {
      gsl_spmatrix_set(A, i, i + 1, 1.0);
      gsl_spmatrix_set(A, i, i, -2.0);
      gsl_spmatrix_set(A, i, i - 1, 1.0);
    }

  /* last row of matrix */
  gsl_spmatrix_set(A, n - 1, n - 1, -2.0);
  gsl_spmatrix_set(A, n - 1, n - 2, 1.0);

  /* scale by h^2 */
  gsl_spmatrix_scale(A, 1.0 / (h * h));

  /* construct right hand side vector */
  for (i = 0; i < n; ++i)
    {
      double xi = (i + 1) * h;
      double bi = -M_PI * M_PI * sin(M_PI * xi);
      gsl_vector_set(b, i, bi);
    }

  if (compress)
    B = gsl_spmatrix_compcol(A);
  else
    B = A;

  /* solve the system */
  do
    {
      status = gsl_splinalg_itersolve_iterate(B, b, tol, u, w);
    }
  while (status == GSL_CONTINUE && ++iter < max_iter);

  gsl_test(status, "%s poisson status s=%d N=%zu", desc, status, N);

  /* check solution against analytic */
  for (i = 0; i < n; ++i)
    {
      double xi = (i + 1) * h;
      double u_gsl = gsl_vector_get(u, i);
      double u_exact = sin(M_PI * xi);

      gsl_test_rel(u_gsl, u_exact, epsrel, "%s poisson N=%zu i=%zu",
                   desc, N, i);
    }

  /* check that the residual satisfies ||r|| <= tol*||b|| */
  {
    gsl_vector *r = gsl_vector_alloc(n);
    double normr, normb;

    gsl_vector_memcpy(r, b);
    gsl_spblas_dgemv(CblasNoTrans, -1.0, A, u, 1.0, r);

    normr = gsl_blas_dnrm2(r);
    normb = gsl_blas_dnrm2(b);

    status = (normr <= tol*normb) != 1;
    gsl_test(status, "%s poisson residual N=%zu normr=%.12e normb=%.12e",
             desc, N, normr, normb);

    gsl_vector_free(r);
  }

  gsl_splinalg_itersolve_free(w);
  gsl_spmatrix_free(A);
  gsl_vector_free(b);
  gsl_vector_free(u);

  if (compress)
    gsl_spmatrix_free(B);
} /* test_poisson() */

/*
test_toeplitz()
  Test Toeplitz matrix T = A + B + C where:
A = diag(a,-1)
B = diag(b)
C = diag(c,1)

rhs = ones(n,1)
*/

static void
test_toeplitz(const size_t N, const double a, const double b,
              const double c)
{
  int status;
  const double tol = 1.0e-10;
  const size_t max_iter = 10;
  const gsl_splinalg_itersolve_type *T = gsl_splinalg_itersolve_gmres;
  const char *desc;
  gsl_spmatrix *A;
  gsl_vector *rhs, *x;
  gsl_splinalg_itersolve *w;
  size_t i, iter = 0;

  if (N <= 1)
    return;

  A = gsl_spmatrix_alloc(N ,N);
  rhs = gsl_vector_alloc(N);
  x = gsl_vector_calloc(N);
  w = gsl_splinalg_itersolve_alloc(T, N, 0);
  desc = gsl_splinalg_itersolve_name(w);

  /* first row */
  gsl_spmatrix_set(A, 0, 0, b);
  gsl_spmatrix_set(A, 0, 1, c);

  /* interior rows */
  for (i = 1; i < N - 1; ++i)
    {
      gsl_spmatrix_set(A, i, i - 1, a);
      gsl_spmatrix_set(A, i, i, b);
      gsl_spmatrix_set(A, i, i + 1, c);
    }

  /* last row */
  gsl_spmatrix_set(A, N - 1, N - 2, a);
  gsl_spmatrix_set(A, N - 1, N - 1, b);

  /* set rhs vector */
  gsl_vector_set_all(rhs, 1.0);

  /* solve the system */
  do
    {
      status = gsl_splinalg_itersolve_iterate(A, rhs, tol, x, w);
    }
  while (status == GSL_CONTINUE && ++iter < max_iter);

  gsl_test(status, "%s toeplitz status s=%d N=%zu a=%f b=%f c=%f",
           desc, status, N, a, b, c);

  /* check that the residual satisfies ||r|| <= tol*||b|| */
  {
    gsl_vector *r = gsl_vector_alloc(N);
    double normr, normb;

    gsl_vector_memcpy(r, rhs);
    gsl_spblas_dgemv(CblasNoTrans, -1.0, A, x, 1.0, r);

    normr = gsl_blas_dnrm2(r);
    normb = gsl_blas_dnrm2(rhs);

    status = (normr <= tol*normb) != 1;
    gsl_test(status, "%s toeplitz residual N=%zu a=%f b=%f c=%f normr=%.12e normb=%.12e",
             desc, N, a, b, c, normr, normb);

    gsl_vector_free(r);
  }

  gsl_vector_free(x);
  gsl_vector_free(rhs);
  gsl_spmatrix_free(A);
  gsl_splinalg_itersolve_free(w);
} /* test_toeplitz() */

static void
test_random(const size_t N, const gsl_rng *r, const int compress)
{
  const gsl_splinalg_itersolve_type *T = gsl_splinalg_itersolve_gmres;
  const double tol = 1.0e-8;
  int status;
  gsl_spmatrix *A = create_random_sparse(N, N, 0.3, r);
  gsl_spmatrix *B;
  gsl_vector *b = gsl_vector_alloc(N);
  gsl_vector *x = gsl_vector_calloc(N);

  /* these random matrices require all N iterations to converge */
  gsl_splinalg_itersolve *w = gsl_splinalg_itersolve_alloc(T, N, N);

  const char *desc = gsl_splinalg_itersolve_name(w);

  create_random_vector(b, r);

  if (compress)
    B = gsl_spmatrix_compcol(A);
  else
    B = A;

  status = gsl_splinalg_itersolve_iterate(B, b, tol, x, w);
  gsl_test(status, "%s random status s=%d N=%zu", desc, status, N);

  /* check that the residual satisfies ||r|| <= tol*||b|| */
  {
    gsl_vector *res = gsl_vector_alloc(N);
    double normr, normb;

    gsl_vector_memcpy(res, b);
    gsl_spblas_dgemv(CblasNoTrans, -1.0, A, x, 1.0, res);

    normr = gsl_blas_dnrm2(res);
    normb = gsl_blas_dnrm2(b);

    status = (normr <= tol*normb) != 1;
    gsl_test(status, "%s random residual N=%zu normr=%.12e normb=%.12e",
             desc, N, normr, normb);

    gsl_vector_free(res);
  }

  gsl_spmatrix_free(A);
  gsl_vector_free(b);
  gsl_vector_free(x);
  gsl_splinalg_itersolve_free(w);

  if (compress)
    gsl_spmatrix_free(B);
} /* test_random() */

int
main()
{
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  size_t n;

  test_poisson(7, 1.0e-1, 0);
  test_poisson(7, 1.0e-1, 1);

  test_poisson(543, 1.0e-5, 0);
  test_poisson(543, 1.0e-5, 1);

  test_poisson(1000, 1.0e-6, 0);
  test_poisson(1000, 1.0e-6, 1);

  test_poisson(5000, 1.0e-7, 0);
  test_poisson(5000, 1.0e-7, 1);

  test_toeplitz(15, 0.01, 1.0, 0.01);
  test_toeplitz(15, 1.0, 1.0, 0.01);
  test_toeplitz(50, 1.0, 2.0, 0.01);
  test_toeplitz(1000, 0.5, 1.0, 0.01);

  for (n = 1; n <= 100; ++n)
    {
      test_random(n, r, 0);
      test_random(n, r, 1);
    }

  gsl_rng_free(r);

  exit (gsl_test_summary());
} /* main() */

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
#ifndef LIC_ASSEMBLY
#define LIC_ASSEMBLY

#include <vector>
#include <gsl/gsl>

#include "TypeDefinition.h"
#include "MethodDefinition.h"

namespace lic
{

class AssemblyLoader;

// We need to use it in the class definition below, hence the duplicity
constexpr size_t MetadataTablesCount = 45;

class Assembly
{
public:
    Assembly(size_t baseRva, std::unique_ptr<gsl::byte[]> &&data);
    ~Assembly();

    // Disable copying
    Assembly(const Assembly&) = delete;
    Assembly& operator=(const Assembly&) = delete;

    const gsl::span<TypeDefinition> Types();
    const gsl::span<MethodDefinition> Methods();

private:
    friend class AssemblyLoader;
    friend class TypeDefinition;
    friend class MethodDefinition;

    size_t baseRva;
    std::unique_ptr<gsl::byte[]> data;

    gsl::span<gsl::byte> metadataStream;
    gsl::span<gsl::byte> stringStream;
    gsl::span<gsl::byte> userStringStream;
    gsl::span<gsl::byte> guidStream;
    gsl::span<gsl::byte> blobStream;

    std::array<size_t, MetadataTablesCount> metadataRowSizes;
    std::array<size_t, MetadataTablesCount> metadataRowCounts;
    std::array<gsl::span<gsl::byte>, MetadataTablesCount> metadataTables;

    std::vector<TypeDefinition> types;
    std::vector<MethodDefinition> methods;
};

}

#endif // !LIC_ASSEMBLY

max_stars_repo_name: roberthusak/lic

max_issues_repo_name: roberthusak/lic

max_forks_repo_name: roberthusak/lic