text:
#pragma once

#include <gsl/string_span>

#include <fmt/format.h>

#include <string>
#include <algorithm>
#include <functional>
#include <vector>
#include <sstream>
#include <cctype>
#include <iosfwd>

// Thanks Stackoverflow...

std::string ucs2_to_local(const std::wstring &);
std::wstring local_to_ucs2(const std::string &);
std::string ucs2_to_utf8(const std::wstring &);
std::wstring utf8_to_ucs2(const std::string &);

// trim from start
inline std::string &ltrim(std::string &s) {
	s.erase(s.begin(), find_if(s.begin(), s.end(), [](auto c) { return !std::isspace(c); }));
	return s;
}

// trim from end
inline std::string &rtrim(std::string &s) {
	s.erase(find_if(s.rbegin(), s.rend(), [](auto c) { return !std::isspace(c); }).base(), s.end());
	return s;
}

// trim from both ends
inline std::string &trim(std::string &s) {
	return ltrim(rtrim(s));
}

inline std::vector<gsl::cstring_span<>>& split(gsl::cstring_span<> s,
	char delim,
	std::vector<gsl::cstring_span<>>& elems,
	bool trimItems = false,
	bool keepEmpty = false) {

	int pos = 0;
	while (pos < s.size()) {
		auto ch = s[pos++];

		// Empty item
		if (ch == delim) {
			if (keepEmpty) {
				elems.push_back({});
			}
			continue;
		}

		if (trimItems && isspace(ch)) {
			// Skip all chars that are considered whitespace
			continue;
		}

		auto start = pos - 1; // Start of item
		size_t count = 1; // Size of item

		// Seek past all of the item's content
		while (pos < s.size()) {
			ch = s[pos++];
			if (ch == delim) {
				break; // reached the end of the item
			}
			count++;
		}

		// Trim whitespace at the end of the string
		if (trimItems) {
			while (count > 0 && isspace(s[start + count - 1])) {
				count--;
			}
		}

		if (count != 0) {
			elems.push_back(s.subspan(start, count));
		} else if (keepEmpty) {
			elems.push_back({});
		}

	}

	return elems;
}

inline std::vector<gsl::cstring_span<>> split(gsl::cstring_span<> s, char delim, bool trim = false, bool keepEmpty = false) {
	std::vector<gsl::cstring_span<>> elems;
	split(s, delim, elems, trim, keepEmpty);
	return elems;
}

inline std::vector<std::string>& split(const std::string& s, char delim, std::vector<std::string>& elems, bool trimItems = false, bool keepEmpty = false) {
	std::stringstream ss(s);
	std::string item;
	while (getline(ss, item, delim)) {
		if (trimItems) {
			item = trim(item);
		}
		if (keepEmpty || !item.empty()) {
			elems.push_back(item);
		}
	}
	return elems;
}

inline std::vector<std::string> split(const std::string& s, char delim, bool trim = false, bool keepEmpty = false) {
	std::vector<std::string> elems;
	split(s, delim, elems, trim, keepEmpty);
	return elems;
}

inline std::string tolower(const std::string &s) {
	auto needsConversion = std::any_of(s.begin(), s.end(), [](char a) {
		return std::tolower(a) != a;
	});
	if (needsConversion) {
		std::string result = s;
		std::transform(result.begin(), result.end(), result.begin(), [](char ch) { return std::tolower((int)ch); });
		return result;
	} else {
		return s;
	}
}

inline std::string tounderscore(const std::string &s) {
	auto needsConversion = std::any_of(s.begin(), s.end(), [](char a) {
		return a == ' ';
	});
	if (needsConversion) {
		std::string result = s;
		std::transform(result.begin(), result.end(), result.begin(), [](char ch){
			if (ch == ' '){
				return '_';
			}
			return ch;
		});
		return result;
	}
	else {
		return s;
	}
}

inline std::string toupper(const std::string &s) {
	auto needsConversion = std::any_of(s.begin(), s.end(), [](char a) {
		return std::toupper(a) != a;
	});
	if (needsConversion) {
		std::string result = s;
		std::transform(result.begin(), result.end(), result.begin(), [](char ch) { return std::toupper((int)ch); });
		return result;
	}
	else {
		return s;
	}
}

// Nice to have operator for serializing vectors in the logger
namespace std {
	template<typename T>
	void format_arg(fmt::BasicFormatter<char> &f, const char *&format_str, const std::vector<T> &v) {
		using namespace std;

		f.writer().write("[");
		for (size_t i = 0; i < v.size(); ++i) {
			f.writer().write("{}", v[i]);
			if (i != v.size() - 1) {
				f.writer().write(", ");
			}
		}
		f.writer().write("]");
	}
}

inline bool endsWith(const std::string &str, const std::string &suffix) {
	if (suffix.length() > str.length()) {
		return false; // Short-circuit
	}

	return str.compare(str.length() - suffix.length(), suffix.length(), suffix) == 0;
}

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
/* integration/exponential.c
 *
 * Copyright (C) 2017 Konrad Griessinger, 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.
 */

/*
 * The code in this module is based on IQPACK, specifically the LGPL
 * implementation found in HERMITE_RULE:
 * https://people.sc.fsu.edu/~jburkardt/c_src/hermite_rule/hermite_rule.html
 */

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_gamma.h>

static int
exponential_check(const size_t n, const gsl_integration_fixed_params * params)
{
  (void) n;

  if (fabs(params->b - params->a) <= GSL_DBL_EPSILON)
    {
      GSL_ERROR("|b - a| too small", GSL_EDOM);
    }
  else if (params->a >= params->b)
    {
      GSL_ERROR("lower integration limit must be smaller than upper limit", GSL_EDOM);
    }
  else if (params->alpha <= -1.0)
    {
      GSL_ERROR("alpha must be > -1", GSL_EDOM);
    }
  else
    {
      return GSL_SUCCESS;
    }
}

static int
exponential_init(const size_t n, double * diag, double * subdiag, gsl_integration_fixed_params * params)
{
  size_t i;
  double a2i = params->alpha;

  /* construct the diagonal and subdiagonal elements of Jacobi matrix */
  for (i = 1; i <= n; i++)
    {
      diag[i-1] = 0.0;
      a2i += 2.0;
      subdiag[i-1] = (i + params->alpha*(i%2))/sqrt( a2i*a2i - 1.0 );
    }

  params->zemu = 2.0/(params->alpha + 1.0);
  params->shft = 0.5*(params->b + params->a);
  params->slp = 0.5*(params->b - params->a);
  params->al = params->alpha;
  params->be = 0.0;

  return GSL_SUCCESS;
}

static const gsl_integration_fixed_type exponential_type =
{
  exponential_check,
  exponential_init
};

const gsl_integration_fixed_type *gsl_integration_fixed_exponential = &exponential_type;

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef parser_496ce533_26fe_412b_8db6_e61b852be838_h
#define parser_496ce533_26fe_412b_8db6_e61b852be838_h

#include <rathen\config.h>
#include <rathen\basis.h>
#include <gslib\tree.h>

__rathen_begin__

class __gs_novtable ps_obj:
    public object
{
public:
    enum
    {
        pt_root,
        pt_block,
        pt_value,
        pt_expression,
        pt_operator,
        pt_function,
        pt_calling,
        pt_if_statement,
        pt_loop_statement,
        pt_go_statement,
        pt_stop_statement,
        pt_exit_statement,
        pt_variable,
        pt_constant,
    };

    enum
    {
        vt_string,
        vt_integer,
        vt_float,
    };

public:
    ps_obj() { _parent = 0; }
    ps_obj(ps_obj* parent) { _parent = parent; }
    ps_obj* get_parent() const { return _parent; }

public:
    virtual ~ps_obj() {}
    virtual bool is_anonymous() const = 0;

protected:
    ps_obj*         _parent;
};

struct ps_trait_pass {};
struct ps_trait_hold {};

template<uint tag, bool anony, class ph_trait, class bsc = ps_obj>
class ps_proto:
    public bsc
{
public:
    ps_proto();
    ps_proto(ps_obj* parent);
    bool is_holder() const override;
    uint get_tag() const override;
    bool is_anonymous() const override;
};

template<uint tag, bool anony, class bsc>
class ps_proto<tag, anony, ps_trait_pass, bsc>:
    public bsc
{
public:
    ps_proto() {}
    ps_proto(ps_obj* parent): bsc(parent) {}
    bool is_holder() const override { return false; }
    uint get_tag() const override { return tag; }
    bool is_anonymous() const override { return anony; }
};

template<uint tag, bool anony, class bsc>
class ps_proto<tag, anony, ps_trait_hold, bsc>:
    public bsc
{
public:
    ps_proto() {}
    ps_proto(ps_obj* parent): bsc(parent) {}
    bool is_holder() const override { return true; }
    uint get_tag() const override { return tag; }
    bool is_anonymous() const override { return anony; }

public:
    virtual void add_indexing(object* obj) { _indexing.add_value(obj); }
    virtual indexing& get_indexing() { return _indexing; }

protected:
    indexing        _indexing;

public:
    object* find_object(const gchar* name) { return _indexing.find_value(name); }
    const object* find_object(const gchar* name) const { return _indexing.find_value(name); }
};

template<class mc, class pc>
class ps_join:
    public mc,
    public pc
{
public:
    ps_join() {}
    ps_join(ps_obj* parent): mc(parent) {}
};

class ps_comp_value
{
public:
    ps_comp_value() { _type = 0, _global = false, _ref = false; }
    void set_type(type* ty) { _type = ty; }
    void set_global(bool g) { _global = g; }
    void set_ref(bool r) { _ref = r; }
    type* get_type() const { return _type; }
    bool is_complex() const { return _type->is_complex(); }
    bool is_global() const { return _global; }
    bool is_ref() const { return _ref; }

protected:
    type*           _type;
    bool            _global;
    bool            _ref;
};

class ps_comp_operator
{
public:
    ps_comp_operator() { _opr = 0; }
    void set_opr(const oprnode* opr) { _opr = opr; }
    void set_opr(const oprinfo& opr) { _opr = opr_manager::get_singleton_ptr()->get_operator(opr); }
    const oprnode* get_opr_info() const { return _opr; }
    bool is_unary() const { return _opr->unary; }
    int get_priority() const { return _opr->prior; }
    const gchar* get_opr_name() const { return _opr->opr; }

protected:
    const oprnode*  _opr;
};

class ps_comp_variable
{
public:
    ps_comp_variable() { _value = 0; }
    void set_value(ps_obj* value) { _value = value; }
    ps_obj* get_value() const { return _value; }

protected:
    ps_obj*         _value;
};

class ps_comp_constant
{
public:
    union vdif
    {
        int         i;
        real32      r;
        int64       i64;
        real        r64;
    };

public:
    ps_comp_constant() {}
    void set_string(const gchar* str) { _cvtag = ps_obj::vt_string, _vstr = str; }
    void set_integer(int i) { _cvtag = ps_obj::vt_integer, _value.i = i; }
    void set_float(real32 r) { _cvtag = ps_obj::vt_float, _value.r = r; }
    uint get_value_tag() const { return _cvtag; }
    const string& get_string() const { return _vstr; }
    int get_integer() const { return _value.i; }
    real32 get_float() const { return _value.r; }

protected:
    uint            _cvtag;
    string          _vstr;
    vdif            _value;
};

class ps_comp_function
{
public:
    ps_comp_function() { _rety = 0, _parent = 0, _retr = false; }
    void set_retv_type(type* ty) { _rety = ty; }
    void set_retv_ref(bool tr) { _retr = tr; }
    type* get_retv_type() const { return _rety; }
    bool get_retv_ref() const { return _retr; }
    void set_parent(ps_obj* obj) { _parent = obj; }
    ps_obj* get_parent() const { return _parent; }
    bool is_top_level() const { return _parent->get_tag() == ps_obj::pt_root; }

protected:
    type*           _rety;
    bool            _retr;
    ps_obj*         _parent;
};

class ps_comp_statement
{
public:
    ps_comp_statement() { _condition = 0, _block = 0; }
    void set_condition(ps_obj* condition) { _condition = condition; }
    void set_block(ps_obj* block) { _block = block; }
    ps_obj* get_condition() const { return _condition; }
    ps_obj* get_block() const { return _block; }

protected:
    ps_obj*         _condition;
    ps_obj*         _block;
};

class ps_comp_if_statement:
    public ps_comp_statement
{
public:
    ps_comp_if_statement() { _prevbr = _nextbr = _lastbr = 0; }
    void set_prev_branch(ps_obj* branch) { _prevbr = branch; }
    void set_next_branch(ps_obj* branch) { _nextbr = branch; }
    void set_last_branch(ps_obj* branch) { _lastbr = branch; }
    ps_obj* get_prev_branch() const { return _prevbr; }
    ps_obj* get_next_branch() const { return _nextbr; }
    ps_obj* get_last_branch() const { return _lastbr; }

protected:
    ps_obj*         _prevbr;
    ps_obj*         _nextbr;
    ps_obj*         _lastbr;
};

class ps_comp_loop_statement:
    public ps_comp_statement
{
};

class ps_comp_exit_statement
{
public:
    ps_comp_exit_statement() { _conclusion = 0; }
    void set_conclusion(ps_obj* conclusion) { _conclusion = conclusion; }
    ps_obj* get_conclusion() const { return _conclusion; }

protected:
    ps_obj*         _conclusion;
};

typedef ps_proto<ps_obj::pt_root, true, ps_trait_hold> ps_root;
typedef ps_proto<ps_obj::pt_block, true, ps_trait_hold> ps_block;
typedef ps_proto<ps_obj::pt_value, false, ps_trait_pass, ps_join<ps_obj, ps_comp_value> > ps_value;
typedef ps_proto<ps_obj::pt_expression, true, ps_trait_pass> ps_expression;
typedef ps_proto<ps_obj::pt_operator, true, ps_trait_pass, ps_join<ps_obj, ps_comp_operator> > ps_operator;
typedef ps_proto<ps_obj::pt_function, false, ps_trait_hold, ps_join<ps_obj, ps_comp_function> > ps_function;
typedef ps_proto<ps_obj::pt_calling, true, ps_trait_pass> ps_calling;
typedef ps_proto<ps_obj::pt_go_statement, true, ps_trait_pass> ps_go_statement;
typedef ps_proto<ps_obj::pt_stop_statement, true, ps_trait_pass> ps_stop_statement;
typedef ps_proto<ps_obj::pt_if_statement, true, ps_trait_pass, ps_join<ps_obj, ps_comp_if_statement> > ps_if_statement;
typedef ps_proto<ps_obj::pt_loop_statement, true, ps_trait_pass, ps_join<ps_obj, ps_comp_loop_statement> > ps_loop_statement;
typedef ps_proto<ps_obj::pt_exit_statement, true, ps_trait_pass, ps_join<ps_obj, ps_comp_exit_statement> > ps_exit_statement;
typedef ps_proto<ps_obj::pt_variable, false, ps_trait_pass, ps_join<ps_obj, ps_comp_variable> > ps_variable;
typedef ps_proto<ps_obj::pt_constant, false, ps_trait_pass, ps_join<ps_obj, ps_comp_constant> > ps_constant;

class parser:
    private tree<ps_obj, _treenode_wrapper<ps_obj> >
{
public:
    typedef _treenode_wrapper<ps_obj> wrapper;
    typedef ps_obj value;
    typedef tree<ps_obj, wrapper> superref;
    typedef superref::iterator iterator;
    typedef superref::const_iterator const_iterator;

public:
    parser() { insert<ps_root>(iterator()); }
    const gchar* get_name() const;
    void set_root(const gchar* name);
    template<class _cst>
    iterator insert(iterator i) { return superref::insert_tail<_cst>(i); }
    template<class _cst>
    iterator birth(iterator i) { return superref::birth_tail<_cst>(i); }

public:
    using superref::clear;
    using superref::erase;
    using superref::get_root;
    using superref::for_each;
    using superref::detach;
    using superref::attach;
    using superref::debug_check;

public:
    iterator find(iterator i, const gchar* name);
};

__rathen_end__

#endif

ext: h

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#ifndef _SPC_DRIZ_H
#define _SPC_DRIZ_H

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_vector.h>

double sgarea(double x1, double y1, double x2, double y2, int is, int js);
double boxer(int is,int js,double *x,double *y);

#endif

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#pragma once

#define NOMINMAX

// Windows
#include <windows.h>
#include <winrt\Windows.Foundation.h>

// Standard
#include <string>
#include <sstream>

#if defined(DEBUG) || defined(_DEBUG)
#define _CRTDBG_MAP_ALLOC
#include <stdlib.h>
#include <crtdbg.h>
#endif

// Guidelines Support Library
#include <gsl\gsl>

// DirectX
#include <d3d11_4.h>
#include <dxgi1_6.h>
#include <DirectXMath.h>


text:
#include <stdio.h>
#include <stdlib.h>

#include <gsl/gsl_spmatrix.h>

int
main()
{
  gsl_spmatrix *A = gsl_spmatrix_alloc(5, 4); /* triplet format */
  gsl_spmatrix *B, *C;
  size_t i, j;

  /* build the sparse matrix */
  gsl_spmatrix_set(A, 0, 2, 3.1);
  gsl_spmatrix_set(A, 0, 3, 4.6);
  gsl_spmatrix_set(A, 1, 0, 1.0);
  gsl_spmatrix_set(A, 1, 2, 7.2);
  gsl_spmatrix_set(A, 3, 0, 2.1);
  gsl_spmatrix_set(A, 3, 1, 2.9);
  gsl_spmatrix_set(A, 3, 3, 8.5);
  gsl_spmatrix_set(A, 4, 0, 4.1);

  printf("printing all matrix elements:\n");
  for (i = 0; i < 5; ++i)
    for (j = 0; j < 4; ++j)
      printf("A(%zu,%zu) = %g\n", i, j,
             gsl_spmatrix_get(A, i, j));

  /* print out elements in triplet format */
  printf("matrix in triplet format (i,j,Aij):\n");
  gsl_spmatrix_fprintf(stdout, A, "%.1f");

  /* convert to compressed column format */
  B = gsl_spmatrix_ccs(A);

  printf("matrix in compressed column format:\n");
  printf("i = [ ");
  for (i = 0; i < B->nz; ++i)
    printf("%d, ", B->i[i]);
  printf("]\n");

  printf("p = [ ");
  for (i = 0; i < B->size2 + 1; ++i)
    printf("%d, ", B->p[i]);
  printf("]\n");

  printf("d = [ ");
  for (i = 0; i < B->nz; ++i)
    printf("%g, ", B->data[i]);
  printf("]\n");

  /* convert to compressed row format */
  C = gsl_spmatrix_crs(A);

  printf("matrix in compressed row format:\n");
  printf("i = [ ");
  for (i = 0; i < C->nz; ++i)
    printf("%d, ", C->i[i]);
  printf("]\n");

  printf("p = [ ");
  for (i = 0; i < C->size1 + 1; ++i)
    printf("%d, ", C->p[i]);
  printf("]\n");

  printf("d = [ ");
  for (i = 0; i < C->nz; ++i)
    printf("%g, ", C->data[i]);
  printf("]\n");

  gsl_spmatrix_free(A);
  gsl_spmatrix_free(B);
  gsl_spmatrix_free(C);

  return 0;
}

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <stdio.h>
#include <math.h>  /*per sin i cos*/
#include <stdlib.h> /*pel malloc i calloc*/
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "memory.h"

// what data to save
int writeData; /* command line input */

// simulation run time
double tStop; /* command line input */

//capacitance
double c=1.;

// leak
double glk=1;
double vlk=-70;

// Na channel
double gna=37;
double vna=55;

// Na activation [m]
double thetam=30;
double sm=15;

// Na inactivation [h]
double thetah=-39;
double sh=3.1;
double tauh0=1;
double tauh1=500;
double thh=57;
double sigmah=3;

// K channel
double gk=45;
double vk=-80;

// K activation
double thetan=-32.;
double sn=-8.;
double taun0=1;
double taun1=100.;
double thn=80.;
double sigman=26.;

// A channel
double ga; /* command line input */

// A activation [a]
double taua   = 2;
double thetaa = -50;
double sigmaa = 20;

// A inactivation [b]
double taub=150;
double thetab=-70;
double sigmab=-6;

// kinetics rescaling
double phi=0.75;

// synapse
double vsyne = 0.;
double vsyni = -85.;
double gsyne, gsyni; /* command line input */
double rsyne, rsyni; /* command line input */
double betae = .2;
double betai = .18;
double te, ti, teNext,tiNext; /* synapse event times */

double v,m,b,n,a,se,si;
double minf,hinf,ninf,binf,ainf;
double tauh,taun,taua;
double ina,ik,ilk,ia,isyne,isyni;
double alphae;

//double twopi=8*atan(1.);
double twopi=3.1415926535;

int func (double t, const double x[], double dx[], void *params);

int jac (double t, const double y[], double *dfdy, double dfdt[], void *params);

int main(int argc, char *argv[]){

    double *x,*dx;
    double vlast,tspike;
    double si0;
    double rand;
    int ndim=6;
    int i;
    FILE *output;
    double tol=1.0e-5;
    double mu;

    //const gsl_odeiv_step_type * T  = gsl_odeiv_step_rk8pd;  // Runge-Kutta Prince-Dormand 8-9.
    //const gsl_odeiv_step_type * T  = gsl_odeiv_step_rk2imp; // Runge-Kutta 2 implicit.
    const gsl_odeiv_step_type * T  = gsl_odeiv_step_rk4imp; // Runge-Kutta 4 implicit.
    //const gsl_odeiv_step_type * T  = gsl_odeiv_step_bsimp; // Burlirsch-Stoer implicit.. Need jacobian!
    //const gsl_odeiv_step_type * T  = gsl_odeiv_step_gear1; // Gear 1 implicit.
    //const gsl_odeiv_step_type * T  = gsl_odeiv_step_gear2; // Gear 2 implicit.

    gsl_odeiv_step * s  = gsl_odeiv_step_alloc (T, ndim);
    gsl_odeiv_control * c  = gsl_odeiv_control_y_new (tol, 0.0);
    gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (ndim);

    gsl_odeiv_system sys = {func, jac, ndim, &mu};

    /* Random generator stuff */
    const gsl_rng_type *U;
    gsl_rng *r;
    gsl_rng_env_setup();
    U = gsl_rng_default;
    r = gsl_rng_alloc (U);

    /* Parameters From Command Line */
    writeData = atof(argv[1]);
    tStop     = atof(argv[2]);
    ga        = atof(argv[3]);
    gsyne     = atof(argv[4]);
    gsyni     = atof(argv[5]);
    rsyne     = atof(argv[6]);
    rsyni     = atof(argv[7]);

    double t, tf;
    double h = 1.e-3;
    int status;

    char name1[100];
    int number1,number2,number3;
    int ncount=0;
    int nite;
    double sum;

    if (writeData==1){  /* save all data */
      sprintf(name1, "data1cpt_ga%g_ge%g_gi%g_re%.1f_ri%g.txt",ga,gsyne, gsyni, rsyne,rsyni);}
    else if (writeData==2){ /* save spikes */
      sprintf(name1, "spike1cpt_ga%g_ge%g_gi%g_re%.1f_ri%g.txt",ga,gsyne, gsyni, rsyne,rsyni);}
    output=fopen(name1, "w");

   /*open the data file*/
    output=fopen(name1, "w");
    if (output==NULL){
    printf("Error in opening data file\n");
    exit(1);
    }

/*memory allocation*/
x  = new_vector(0, ndim-1);
dx = new_vector(0, ndim-1);

/*Initial point*/
double v0=-70;
x[0]=v0; /* v */
x[1]=1./(1.+exp((v0-thetan)/sn));     /* n */
x[2]=1/(1+exp(-(v0-thetaa)/sigmaa));  /* a */
x[3]=1/(1+exp(-(v0-thetab)/sigmab));  /* b */
x[4]=0; /* se */
x[5]=0; /* si */

/*integrate*/

t=0.;
nite=0.;
vlast=v0;
te = 0;
ti = -999.; tiNext = 1000./rsyni; si0 = 0;
te = -999.; rand =gsl_ran_flat(r, 0,1); teNext = -1000.*log(1-rand)/(rsyne);
if (teNext<tiNext){tf=teNext;}
else {tf=tiNext;}

while (t<tStop){

  while (t<tf){

   status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, tf, &h, x);

   if (status != GSL_SUCCESS)
    break;

  if (writeData==1){ /* save everything */
  fprintf(output," %.5f ",t);
  for (i=0;i<ndim;i++){
   fprintf(output," %.5f ",x[i]);
   }
  fprintf(output, "\n");
  }

  /* spike detection */
  tspike =0.;
  if ((x[0]>=-10)&&(vlast<-10)){
     ncount++;
     tspike=t;
     if (writeData==2){ /* save spikes */
        fprintf(output," %.5f %.5f",tspike, si0);  fprintf(output, "\n");
     }
  }

  vlast=x[0];

  nite++;
  }

  // generate next excitatory event
  if   (tf>=teNext){
    te = teNext;
    rand=gsl_ran_flat(r, 0,1); teNext+=-1000.*log(1-rand)/(rsyne);

    // for instantaneous rise se
    x[4]=1.; // ceiling at 1

    si0 =  x[5]; // si value at the excitatory event onset
    if (writeData==2){fprintf(output," %.5f %.5f ",-tf, x[5]);  fprintf(output, "\n");} /* save excitatory events */
  }

  // generate next inhibitory event
  if   (tf>=tiNext){tiNext+=1000./rsyni; x[5]=1.;} // ceiling at 1
  if (teNext<tiNext){tf=teNext;}
  else {tf=tiNext;}

  }

  if (writeData==2){fprintf(output," %.5f  %.5f ",tStop , tStop);  fprintf(output, "\n");}
  free_vector(x,0);
  free_vector(dx,0);

  fclose(output);

  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);

}

/* Vector field we want to integrate*/

int func (double t, const double x[], double dx[], void *params)
{
    double mu = *(double *)params;

   double m3,n4,a3;

//  general variables

v = x[0];
m = 1./(1.+exp(-(v+thetam)/sm));   // m—>minf instantaneous
n = x[1];
a = x[2] ; // 1/(1+exp(-(v-thetaa)/sigmaa)); // a->ainf instantaneous
b = x[3];
se= x[4];
si= x[5];

//  currents
m3=m*m*m;
n4=n*n*n*n;
a3=a*a*a;
ina= gna*m3*(1-n)*(v-vna);
ik = gk*n4*(v-vk);
ia = ga*a3*b*(v-vk);
ilk= glk*(v-vlk);
isyne = gsyne*se*(v-vsyne);
isyni = gsyni*si*(v-vsyni);

/* update voltage */
dx[0] = - (ilk + ina + ik + ia + isyni + isyne)/c;

/* update n gating */
ninf=1./(1.+exp((v-thetan)/sn));
taun=taun0+taun1/(1+exp((v+thn)/sigman));
dx[1]=phi*(ninf-n)/taun;

/* update a gating */
ainf = 1/(1+exp(-(v-thetaa)/sigmaa));
dx[2] = (ainf-a)/taua;

/* update b gating */
binf = 1./(1+exp(-(v-thetab)/sigmab));
dx[3]=(binf-b)/taub;

/* update excitation */
dx[4] = -betae*se; // instantaneous rise, exponential decay

/* update inhibition */
dx[5] = -betai*si;

  return GSL_SUCCESS;
}

/* The jacobian
 * In case we wanna use implict Burslish-Stoer methods
 */
int jac (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
    double mu = *(double *)params;
    gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);
    gsl_matrix * m = &dfdy_mat.matrix;

    // Si usamos Burlish-Stoer tenemos de dar el Jacobiano
    /*gsl_matrix_set (m, 0, 0, 0.0);
    gsl_matrix_set (m, 0, 1, 1.0);
    gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
    gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));*/

    // Si usamos otro método podemos dejar el jacobiano
    // todo a cero.
    gsl_matrix_set (m, 0, 0, 0.0);
    gsl_matrix_set (m, 0, 0, 0.0);
    gsl_matrix_set (m, 0, 0, 0.0);
    gsl_matrix_set (m, 0, 0, 0.0);
    dfdt[0] = 0.0;
    dfdt[1] = 0.0;

    return GSL_SUCCESS;
}

max_stars_repo_name: jgoldwyn/Gain-Control-With-IA

max_issues_repo_name: jgoldwyn/Gain-Control-With-IA

max_forks_repo_name: jgoldwyn/Gain-Control-With-IA


text:
/* Copyright (C) 2010-2019, The Regents of The University of Michigan.
 All rights reserved.

 This software was developed as part of the The Vulcan project in the Intelligent Robotics Lab
 under the direction of Benjamin Kuipers, kuipers@umich.edu. Use of this code is governed by an
 MIT-style License that can be found at "https://github.com/h2ssh/Vulcan".
*/

/**
* \file     mpepc_pose_follower.h
* \author   Jong Jin Park
*
* Declaration of MPEPCPoseFollower.
*/

#ifndef MPEPC_POSE_FOLLOWER_H
#define MPEPC_POSE_FOLLOWER_H

#include <mpepc/trajectory/params.h>
#include <mpepc/trajectory/robot_trajectory_info.h>
#include <core/pose.h>
#include <mpepc/control/control_law_coordinates.h>
#include <nlopt.h>
#include <vector>
#include <memory>

namespace vulcan
{

namespace mpepc
{

class  RobotSimulator;
class  TrajectoryEvaluator;
class  TaskManifold;
struct trajectory_planner_debug_info_t;

/**
* MPEPCPoseFollower optimizes velocity gain given a pose to converge to.
*
* The optimizer handles the simulate/evaluate loop internally and finds the optimal
* target along with the intermediate trajectories generated by the optimization
* algorithm. NLopt package is used for the numerical optimzation, and a number of
* different algorithms can be implemented.
*/
class MPEPCPoseFollower
{

public:

    /**
    * Constructor for MPEPCPoseFollower.
    *
    * \param    params          Parameters for the optimization
    */
    MPEPCPoseFollower(const mpepc_optimizer_params_t& params);

    /**
    * Destructor for MPEPCPoseFollower.
    */
    ~MPEPCPoseFollower(void);

    /**
    * Setup sets up the optmizer with the robot simulator and trajectory evaluator to be used.
    *
    * \param    simulator       RobotSimulator to use for the trajectory generation.
    * \param    evaluator       TrajectoryEvaluator for calculating the cost of a trajectory.
    */
    void setup(RobotSimulator& simulator, TrajectoryEvaluator& evaluator); // NOTE: This feels very unsafe. A better way?

    /**
    * setInitialGuesses adds a goal pose (if navigation task) and the previous optimum
    * to the initial guesses for local optimizer.
    *
    * \param    task                 TaskManifold that contains the goal.
    * \param    previousMotionTarget Previous output of the metric planner, i.e. the previous optimum.
    */
    void setInitialGuesses(const TaskManifold&    task,
                           const motion_target_t& previousMotionTarget);

    /**
    * run runs the optimizer and finds the best motion target by evaluating the expected
    * costs of an estimated trajectory generated from a candidate motion target.
    *
    * It returns the optimal target found, while the intermediate information including
    * all trajectories considered are stored in the provided debug info instance.
    *
    * \param    debugInfo        Place to store relevant information
    * \return   Optimal target found
    */
    motion_target_t runOptimizer(trajectory_planner_debug_info_t& debugInfo);

    /**
    * evaluateCost calculates the cost of a fixed-length trajectory represented by the given optimizer_target_t.
    * The length of the trajectory is specified in the given instance of the planner_data_t.
    *
    * \param    coords          Optimization variable, currently encoding a target pose in the robot frame and a velocity gain for the kinematic control law.
    * \return   Cost of the trajectory.
    */
    double evaluateCost(const double x[]);

    /**
    * evaluateGradient calculates the gradient around a particular point, by a simple forward differencing.
    *
    * \param    coords          Optimization variable, currently encoding a target pose in the robot frame and a velocity gain for the kinematic control law.
    * \param    cost            Evaluated cost at the given coordinate.
    * \param    gradient        Calculated gradient for each dimension (r, theta, delta, velocity) [output]
    */
    void evaluateGradient(const double x[], double cost, double gradient[]);

    bool haveFoundSolution(void) const { return haveFoundSolution_; };

//     /**
//     * evaluateSingleMotionTarget is a method which resets the data and info struct and evaluate a single optimizer target.
//     * The optimization_info_t will hold a single trajectory considered. TODO: Do i need this?
//     *
//     * \param candidateMotionTarget
//     * \param data
//     * \param info
//     */
//     double evaluateSingleMotionTarget(const motion_target_t& candidateMotionTarget, const planner_data_t& data, optimization_info_t& info);

private:

    // setting up the optimizer
    void setupNLOpt(void);

    // optimizers
    nlopt_opt globalOptimizer;
    nlopt_opt localOptimizer;

    // simulator and evaluator
    RobotSimulator*      simulator_;
    TrajectoryEvaluator* evaluator_;

    robot_trajectory_info_t trajectory;

    // debug info
    trajectory_planner_debug_info_t* optimizerInfo_;

    // inital robot state
    pose_t robotPose_;

    // target pose to converge to
    pose_t   targetPose_;
    motion_target_t candidateMotionTarget_;

    // indicator for finding a right solution
    bool haveFoundSolution_;

    // parameters
    mpepc_optimizer_params_t params_;
};

} // mpepc
} // vulcan

#endif // MPEPC_POSE_FOLLOWER_H


text:
#ifndef _EM_DEFS_H_
#define _EM_DEFS_H_ 1

#include <petsc.h>

#include <Eigen/Core>
#include <Eigen/Dense>

#include <complex>
#include <tuple>

typedef std::complex<double> Complex;

typedef std::tuple<int, int> TupleII;

const double PI = 3.14159265358979323846;
const double MU = 4 * PI * 1E-7;
const Complex II = Complex(0.0, 1.0);

const double RTOD = 180.0 / PI;
const double DTOR = PI / 180.0;

const double EPS = 1.0E-6;

const int TX_SIZE = 7;

typedef Eigen::Vector3d Point;
typedef Eigen::Vector3d VectorD;
typedef Eigen::Vector3cd VectorZ;

typedef Eigen::Matrix3d Tensor;

struct PETScBlockVector {
  Vec re, im;
};

enum AnisotropicForm { Isotropic = 0, Vertical = 1, Triaxial = 2, Arbitrary = 3 };

enum PolarType { XY_POLAR = -1, YX_POLAR = -2};

enum RefineStrategy { FixedNumber = 0, FixedFraction = 1 };

enum InnerPCType { Mixed = 0, AMS = 1, Direct = 2 };

enum DirectSolverType { MUMPS = 0, SUPERLUDIST = 1 };

enum ObsType {
  F_EX_RI = 111,
  F_EY_RI = 121,
  F_EZ_RI = 131,
  F_HX_RI = 141,
  F_HY_RI = 151,
  F_HZ_RI = 161,
  F_EX_AP = 112,
  F_EY_AP = 122,
  F_EZ_AP = 132,
  F_HX_AP = 142,
  F_HY_AP = 152,
  F_HZ_AP = 162,
  R_XY_AP = 212,
  R_YX_AP = 222,
  Z_XX_RI = 311,
  Z_XY_RI = 321,
  Z_YX_RI = 331,
  Z_YY_RI = 341,
  T_ZX_RI = 351,
  T_ZY_RI = 361,
  Z_XX_AP = 312,
  Z_XY_AP = 322,
  Z_YX_AP = 332,
  Z_YY_AP = 342
};

#define EM_ERR_USER -1

#endif

max_stars_repo_name: emfem/emfem

max_issues_repo_name: emfem/emfem

max_forks_repo_name: emfem/emfem


text:
/*
 *  vbHmmTs.c
 *  Model-specific core functions for VB-HMM-TS.
 *
 *  Created by OKAMOTO Kenji, SAKO Yasushi and RIKEN
 *  Copyright 2011-2015
 *  Cellular Informatics Laboratory, Advance Science Institute, RIKEN, Japan.
 *  All rights reserved.
 *
 *  Ver. 1.1.0
 *  Last modified on 2015.09.17
 */

#include "vbHmmTs.h"
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_psi.h>
#include <string.h>
#include "rand.h"

#ifdef _OPENMP
#include "omp.h"
#endif

//// Uncomment one/both of the following defenitions to activate constraint on I or/and K.
//#define  INTENSITY_CAP
//#define  TRANSITION_RATE_CAP
#ifdef INTENSITY_CAP
#define  maxIntensityRatio  10.0
#endif
#ifdef TRANSITION_RATE_CAP
#define  minPhotonNumPerState  10.0
#endif

#define  MAX(a,b)  ((a)>(b)?(a):(b))
#define  MIN(a,b)  ((a)<(b)?(a):(b))

//static int isGlobalAnalysis = 0;

void setFunctions_ts(){
    commonFunctions funcs;
    funcs.newModelParameters    = newModelParameters_ts;
    funcs.freeModelParameters   = freeModelParameters_ts;
    funcs.newModelStats         = newModelStats_ts;
    funcs.freeModelStats        = freeModelStats_ts;
    funcs.initializeVbHmm       = initializeVbHmm_ts;
    funcs.pTilde_z1             = pTilde_z1_ts;
    funcs.pTilde_zn_zn1         = pTilde_zn_zn1_ts;
    funcs.pTilde_xn_zn          = pTilde_xn_zn_ts;
    funcs.calcStatsVars         = calcStatsVars_ts;
    funcs.maximization          = maximization_ts;
    funcs.varLowerBound         = varLowerBound_ts;
    funcs.reorderParameters     = reorderParameters_ts;
    funcs.outputResults         = outputResults_ts;
    setFunctions( funcs );
}

//void setGFunctions_ts(){
//    gCommonFunctions funcs;
//    funcs.newModelParameters    = newModelParameters_ts;
//    funcs.freeModelParameters   = freeModelParameters_ts;
//    funcs.newModelStats         = newModelStats_ts;
//    funcs.freeModelStats        = freeModelStats_ts;
//    funcs.newModelStatsG        = newModelStatsG_ts;
//    funcs.freeModelStatsG       = freeModelStatsG_ts;
//    funcs.initializeVbHmmG      = initializeVbHmmG_ts;
//    funcs.pTilde_z1             = pTilde_z1_ts;
//    funcs.pTilde_zn_zn1         = pTilde_zn_zn1_ts;
//    funcs.pTilde_xn_zn          = pTilde_xn_zn_ts;
//    funcs.calcStatsVarsG        = calcStatsVarsG_ts;
//    funcs.maximizationG         = maximizationG_ts;
//    funcs.varLowerBoundG        = varLowerBoundG_ts;
//    funcs.reorderParametersG    = reorderParametersG_ts;
//    funcs.outputResultsG        = outputResultsG_ts;
//    setGFunctions( funcs );
//    isGlobalAnalysis = 1;
//}

void outputResults_ts( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    outputTsResults( xn, gv, iv, logFP );
}

//void outputResultsG_ts( xns, gv, ivs, logFP )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//FILE *logFP;
//{
//    outputTsResultsG( xns, gv, ivs, logFP );
//}

void *newModelParameters_ts( xn, sNo )
xnDataSet *xn;
int sNo;
{
    int i;
    tsParameters *p = (void*)malloc( sizeof(tsParameters) );

    p->uPiArr = (double*)malloc( sNo * sizeof(double) );
    p->sumUPi = 0.0;
    // hyperparameter for p( k(i,j) ) (i != j)
    p->uKMat = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->uKMat[i] = (double*)malloc( sNo * sizeof(double) );
    }
    p->sumUKArr = (double*)malloc( sNo * sizeof(double) );
    // hyperparameter for p( I(k) )
    p->aIArr = (double*)malloc( sNo * sizeof(double) );
    p->bIArr = (double*)malloc( sNo * sizeof(double) );

    p->avgPi = (double *)malloc( sNo * sizeof(double) );
    p->avgLnPi = (double *)malloc( sNo * sizeof(double) );
    p->avgI = (double *)malloc( sNo * sizeof(double) );
    p->avgLnI = (double *)malloc( sNo * sizeof(double) );
    p->avgK = (double **)malloc( sNo * sizeof(double*) );
    p->avgLnK = (double **)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ ){
        p->avgK[i] = (double *)malloc( sNo * sizeof(double) );
        p->avgLnK[i] = (double *)malloc( sNo * sizeof(double) );
    }
    p->avgLnKI = (double *)malloc( sNo * sizeof(double) );

    return p;
}

void freeModelParameters_ts( p, xn, sNo )
void **p;
xnDataSet *xn;
int sNo;
{
    tsParameters *gp = *p;
    int i;

    free( gp->uPiArr );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->uKMat[i] );
    }
    free( gp->uKMat );
    free( gp->sumUKArr );
    free( gp->aIArr );
    free( gp->bIArr );

    free( gp->avgPi );
    free( gp->avgLnPi );
    free( gp->avgI );
    free( gp->avgLnI );
    for( i = 0 ; i < sNo ; i++ ){
        free( gp->avgK[i] );
        free( gp->avgLnK[i] );
    }
    free( gp->avgK );
    free( gp->avgLnK );
    free( gp->avgLnKI );

    free( *p );
    *p = NULL;
}

void *newModelStats_ts( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
//    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        tsStats *s = (tsStats*)malloc( sizeof(tsStats) );

        int i;
        s->Ni = (double *)malloc( sNo * sizeof(double) );
        s->Ti = (double *)malloc( sNo * sizeof(double) );
        s->Nii = (double *)malloc( sNo * sizeof(double) );
        s->Nij = (double *)malloc( sNo * sizeof(double) );
        s->Mij = (double **)malloc( sNo * sizeof(double*) );
        for( i = 0 ; i < sNo ; i++ )
        {   s->Mij[i] = (double *)malloc( sNo * sizeof(double) );   }

        return s;

//    } else {
//
//        return NULL;
//
//    }
}

void freeModelStats_ts( s, xn, gv, iv )
void **s;
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
//    if( isGlobalAnalysis == 0 ){
        int sNo = gv->sNo;
        tsStats *gs = *s;
        int i;
        free( gs->Ni );
        free( gs->Ti );
        free( gs->Nii );
        free( gs->Nij );
        for( i = 0 ; i < sNo ; i++ )
        {   free( gs->Mij[i] );   }
        free( gs->Mij );

        free( gs );
        *s = NULL;
//    }
}

//void *newModelStatsG_ts( xns, gv, ivs)
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//    int sNo = gv->sNo;
//    tsGlobalStats *gs = (tsGlobalStats*)malloc( sizeof(tsGlobalStats) );
//
//    return gs;
//}

//void freeModelStatsG_ts( gs, xns, gv, ivs )
//void **gs;
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//
//    free( *gs );
//    *gs = NULL;
//}

void initializeVbHmm_ts( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    tsData *d = xn->data;
    int sNo = gv->sNo;
    tsParameters *p = gv->params;

    int i, j;

    // hyperparameter for p( pi(k) )
    p->sumUPi = 0.0;
    for( i = 0 ; i < sNo ; i++ ){
        p->uPiArr[i] = 1.0;
        p->sumUPi += p->uPiArr[i];
    }

    // hyperparameter for p( k(i,j) ) (i != j)
    for( i = 0 ; i < sNo ; i++ ){
        p->sumUKArr[i] = 0.0;
        for( j = 0 ; j < sNo ; j++ ){
            p->uKMat[i][j] = 1.0;
            if( j != i ){
                p->sumUKArr[i] += p->uKMat[i][j];
            }
        }
    }

    double meanI = (double)xn->N / d->T;

    // hyperparameter for p( I(k) )
    for( i = 0 ; i < sNo ; i++ ){
        p->aIArr[i] = 1.0;
        p->bIArr[i] = 1.0 / meanI;
    }

    initialize_indVars_ts( xn, gv, iv );

    calcStatsVars_ts( xn, gv, iv );
    maximization_ts( xn, gv, iv );
}

//void initializeVbHmmG_ts( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

void initialize_indVars_ts( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;

    int i;
    size_t n;
    double sumPar;
    for( n = 0 ; n < dLen ; n++ ){
        sumPar = 0.0;
        for( i = 0 ; i < sNo ; i++ ){
            gmMat[n][i] = enoise(1.0) + 1.0;
            sumPar += gmMat[n][i];
        }
        for( i = 0 ; i < sNo ; i++ ){
            gmMat[n][i] /= sumPar;
        }
    }
}

xnDataSet *newXnDataSet_ts( filename )
const char *filename;
{
    xnDataSet *xn = (xnDataSet*)malloc( sizeof(xnDataSet) );
    xn->name = (char*)malloc( strlen(filename) + 2 );
    strncpy( xn->name, filename, strlen(filename)+1 );
    xn->data = (tsData*)malloc( sizeof(tsData) );
    tsData *d = (tsData*)xn->data;
    d->T = 0.0;
    d->dt = NULL;
    d->time = NULL;
    return xn;
}

void freeXnDataSet_ts( xn )
xnDataSet **xn;
{
    tsData *d = (tsData*)(*xn)->data;
    free( d->dt );
    free( d->time );
    free( (*xn)->data );
    free( (*xn)->name );
    free( *xn );
    *xn = NULL;
}

double pTilde_z1_ts( i, params )
int i;
void *params;
{
    tsParameters *p = (tsParameters*)params;
    return exp( p->avgLnPi[i] );
}

double pTilde_zn_zn1_ts( i, j, params )
int i, j;
void *params;
{
    tsParameters *p = (tsParameters*)params;
    if( i == j ){
        return exp( p->avgLnI[i] - p->avgLnKI[i] );
    } else {
        return exp( p->avgLnK[i][i] + p->avgLnK[i][j] - p->avgLnKI[i] );
    }
}

double pTilde_xn_zn_ts( xnWv, n, i, params )
xnDataSet *xnWv;
size_t n;
int i;
void *params;
{
    tsParameters *p = (tsParameters*)params;
    tsData *d = (tsData*)xnWv->data;
    return exp( p->avgLnI[i] - p->avgI[i] * d->dt[n] );
}

void calcStatsVars_ts( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    tsData *d = (tsData*)xn->data;
    tsStats *s = (tsStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, ***xiMat = iv->xiMat;
    double *Ni = s->Ni, *Ti = s->Ti;
    double *Nii = s->Nii, *Nij = s->Nij, **Mij = s->Mij;
    size_t n;
    int i, j;

    for( i = 0 ; i < sNo ; i++ ){
        Ni[i] = 1e-10;
        Ti[i] = 1e-10;
        Nii[i] = 0.0;
        Nij[i] = 0.0;
        for( j = 0 ; j < sNo ; j++ ){
            Mij[i][j] = 1e-10;
        }
    }
    for( n = 0 ; n < dLen ; n++ ){
        for( i = 0 ; i < sNo ; i++ ){
            Ni[i]  += gmMat[n][i];
            Ti[i]  += gmMat[n][i] * d->dt[n];
            Nii[i] += xiMat[n][i][i];
            for( j = 0 ; j < sNo ; j++ ){
                if( j != i ){
                    Mij[i][j] += xiMat[n][i][j];
                    Nij[i]    += xiMat[n][i][j];
                }
            }
        }
    }
    for( i = 0 ; i < sNo ; i++ ){
        Nii[i] = MAX( Nii[i], 1.0 );
        Nij[i] = MAX( Nij[i], 1.0 );
    }
}

//void calcStatsVarsG_ts( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

void maximization_ts( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    tsParameters *p = (tsParameters*)gv->params;
    tsStats *s = (tsStats*)iv->stats;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgK = p->avgK, **avgLnK = p->avgLnK;
    double *avgLnKI = p->avgLnKI, *avgI = p->avgI, *avgLnI = p->avgLnI;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi, *aIArr = p->aIArr, *bIArr = p->bIArr;
    double **uKMat = p->uKMat, *sumUKArr = p->sumUKArr;
    double *Ni = s->Ni, *Ti = s->Ti;
    double *Nii = s->Nii, *Nij = s->Nij, **Mij = s->Mij;
    int i, j;

    for( i = 0 ; i < sNo ; i++ ){
        avgPi[i] = ( uPiArr[i] + gmMat[0][i] ) / ( sumUPi + 1.0 );
        avgLnPi[i] = gsl_sf_psi( uPiArr[i] + gmMat[0][i] ) - gsl_sf_psi( sumUPi + 1.0 );

        avgK[i][i] = (Ni[i] + aIArr[i]) * Nij[i] / (Nii[i] - 1.0) / (Ti[i] + bIArr[i]);
        avgLnK[i][i]  = gsl_sf_psi(Ni[i] + aIArr[i]) + gsl_sf_psi(Nij[i]);
        avgLnK[i][i] += -gsl_sf_psi(Nii[i]) - log(Ti[i] + bIArr[i]);

        avgI[i] = (Ni[i] + aIArr[i]) / (Ti[i] + bIArr[i]);
        avgLnI[i] = gsl_sf_psi(Ni[i] + aIArr[i]) - log(Ti[i] + bIArr[i]);

        avgLnKI[i]  = gsl_sf_psi(Ni[i] + aIArr[i]) + gsl_sf_psi(Nii[i] + Nij[i]);
        avgLnKI[i] += -gsl_sf_psi(Nii[i]) - log(Ti[i] + bIArr[i]);

#ifdef INTENSITY_CAP
        double meanI = (double)xnWv->N / xnWv->T;
        avgI[i]  = MIN( avgI[i], maxIntensityRatio * meanI );
        avgLnI[i]  = MIN( avgLnI[i], log(maxIntensityRatio * meanI) );
#endif
#ifdef TRANSITION_RATE_CAP
        avgK[i][i] = MIN( avgK[i][i], avgI[i]/minPhotonNumPerState );
        avgLnK[i][i]  = MIN( avgLnK[i][i], avgLnI[i] - log(minPhotonNumPerState) );
        avgLnKI[i]  = MIN( avgLnKI[i], avgLnI[i] + log(1.0 + 1.0/minPhotonNumPerState) );
#endif

        for( j = 0 ; j < sNo ; j++ ){
            if( i != j ){
                avgK[i][j] = ( uKMat[i][j] + Mij[i][j] ) / ( sumUKArr[i] + Nij[i] );
                avgLnK[i][j] = gsl_sf_psi( uKMat[i][j] + Mij[i][j] ) - gsl_sf_psi( sumUKArr[i] + Nij[i] );
            }
        }
    }
}

//void maximizationG_ts( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

double varLowerBound_ts( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    tsParameters *p = (tsParameters*)gv->params;
    tsStats *s = (tsStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, *cn = iv->cn;
    double *avgLnPi = p->avgLnPi, **avgLnK = p->avgLnK;
    double *avgLnKI = p->avgLnKI, *avgI = p->avgI, *avgLnI = p->avgLnI;
    double *uPiArr = p->uPiArr, sumUPi = p->sumUPi, *aIArr = p->aIArr, *bIArr = p->bIArr;
    double **uKMat = p->uKMat, *sumUKArr = p->sumUKArr;
    double *Ni = s->Ni, *Ti = s->Ti;
    double *Nii = s->Nii, *Nij = s->Nij, **Mij = s->Mij;
    size_t n;
    int i, j;

    double lnpPi;
    lnpPi = gsl_sf_lngamma(sumUPi);
    double Ck = 1.0, lnpKii = (double)sNo * log(Ck);
    double lnpKij = 0.0;
    double lnpI = 0.0;
    double lnqPi;
    lnqPi = gsl_sf_lngamma(sumUPi + 1.0);
    double lnqKiiI = 0.0;
    double lnqKij = 0.0;
    for( i = 0 ; i < sNo ; i++ ){
        lnpPi += (uPiArr[i]-1) * avgLnPi[i] - gsl_sf_lngamma(uPiArr[i]);

        lnpKii -= avgLnK[i][i];

        lnpI += aIArr[i] * log(bIArr[i]) - gsl_sf_lngamma(aIArr[i]);
        lnpI += (aIArr[i]-1.0)*avgLnI[i] - bIArr[i]*avgI[i];

        lnqPi += (uPiArr[i]+gmMat[0][i]-1) * (gsl_sf_psi(uPiArr[i]+gmMat[0][i]) - gsl_sf_psi(sumUPi+1));
        lnqPi -= gsl_sf_lngamma(uPiArr[i] + gmMat[0][i]);

        lnqKiiI += gsl_sf_lngamma(Nii[i] + Nij[i]) + (Ni[i] + aIArr[i]) * log(Ti[i] + bIArr[i]);
        lnqKiiI += -gsl_sf_lngamma(Ni[i] + aIArr[i]) -gsl_sf_lngamma(Nii[i]);
        lnqKiiI += -gsl_sf_lngamma(Nij[i]) + (Nij[i]-1.0)*avgLnK[i][i];
        lnqKiiI += (Ni[i] + Nii[i] + aIArr[i] - 1.0)*avgLnI[i];
        lnqKiiI += -(Nii[i] + Nij[i])*avgLnKI[i] - (Ti[i] + bIArr[i])*avgI[i];

        lnpKij += gsl_sf_lngamma(sumUKArr[i]);
        lnqKij += gsl_sf_lngamma(sumUKArr[i] + Nij[i]);
        for( j = 0 ; j < sNo ; j++ ){
            if( j != i ){
                lnpKij += (uKMat[i][j]-1)*avgLnK[i][j] - gsl_sf_lngamma(uKMat[i][j]);
                lnqKij += (uKMat[i][j]+Mij[i][j]-1)*(gsl_sf_psi(uKMat[i][j]+Mij[i][j])-gsl_sf_psi(sumUKArr[i]+Nij[i]));
                lnqKij -= gsl_sf_lngamma( uKMat[i][j] + Mij[i][j] );
            }
        }
    }

    double lnpX = 0.0;
    for( n = 0 ; n < dLen ; n++ ){
        lnpX += log( cn[n] );
    }

    double val;
    val  = lnpPi + lnpKii + lnpKij + lnpI;
    val -= lnqPi + lnqKiiI + lnqKij;
    val += lnpX;
    val += log(gsl_sf_fact(sNo));

    return val;
}

//double varLowerBoundG_ts( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

void reorderParameters_ts( xn, gv, iv )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
{
    tsParameters *p = (tsParameters*)gv->params;
    tsStats *s = (tsStats*)iv->stats;
    size_t dLen = xn->N;
    int sNo = gv->sNo;
    double **gmMat = iv->gmMat, ***xiMat = iv->xiMat;
    double *avgPi = p->avgPi, *avgLnPi = p->avgLnPi, **avgK = p->avgK;
    double **avgLnK = p->avgLnK, *avgLnKI = p->avgLnKI;
    double *avgI = p->avgI, *avgLnI = p->avgLnI;
    double *Ni = s->Ni, *Ti = s->Ti;
    size_t n;
    int i, j;

    int *index = (int*)malloc( sNo * sizeof(int) );
    double *store = (double*)malloc( sNo * sizeof(double) );
    double **s2D = (double**)malloc( sNo * sizeof(double*) );
    for( i = 0 ; i < sNo ; i++ )
    {   s2D[i] = (double*)malloc( sNo * sizeof(double) );   }

    // index indicates order of avgI values (0=biggest avgI -- sNo=smallest avgI).
    for( i = 0 ; i < sNo ; i++ ){
        index[i] = sNo - 1;
        for( j = 0 ; j < sNo ; j++ ){
            if( j != i ){
                if( avgI[i] < avgI[j] ){
                    index[i]--;
                } else if( avgI[i] == avgI[j] ){
                    if( i < j )
                    {   index[i]--;   }
                }
            }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnPi[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnPi[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgI[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgI[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnI[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnI[i] = store[i];   }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgK[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgK[i][j] = s2D[i][j];   }
    }

    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = avgLnK[i][j];   }
    }
    for( j = 0 ; j < sNo ; j++ ){
        for( i = 0 ; i < sNo ; i++ ){   avgLnK[i][j] = s2D[i][j];   }
    }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = avgLnKI[i];   }
    for( i = 0 ; i < sNo ; i++ ){   avgLnKI[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = Ni[i];   }
    for( i = 0 ; i < sNo ; i++ ){   Ni[i] = store[i];   }

    for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = Ti[i];   }
    for( i = 0 ; i < sNo ; i++ ){   Ti[i] = store[i];   }

    for( n = 0 ; n < dLen ; n++ ){
        for( i = 0 ; i < sNo ; i++ ){   store[index[i]] = gmMat[n][i];   }
        for( i = 0 ; i < sNo ; i++ ){   gmMat[n][i] = store[i];   }
    }

    for( n = 0 ; n < dLen ; n++ ){
        for( j = 0 ; j < sNo ; j++ ){
            for( i = 0 ; i < sNo ; i++ ){   s2D[index[i]][index[j]] = xiMat[n][i][j];   }
        }
        for( j = 0 ; j < sNo ; j++ ){
            for( i = 0 ; i < sNo ; i++ ){   xiMat[n][i][j] = s2D[i][j];   }
        }
    }

    for( i = 0 ; i < sNo ; i++ ){   free( s2D[i] );   }
    free( s2D );
    free( store );
    free( index );
}

//void reorderParametersG_ts( xns, gv, ivs )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//{
//}

void outputTsResults( xn, gv, iv, logFP )
xnDataSet *xn;
globalVars *gv;
indVars *iv;
FILE *logFP;
{
    tsParameters *p = (tsParameters*)gv->params;
    int sNo = gv->sNo;
    int i, j;

    fprintf(logFP, "  results: K = %d \n", sNo);

    fprintf(logFP, "   intensities: ( %g", p->avgI[0]);
    for( i = 1 ; i < sNo ; i++ )
    {   fprintf(logFP, ", %g", p->avgI[i]);   }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   pi: ( %g", p->avgPi[0]);
    for( i = 1 ; i < sNo ; i++ ){
        fprintf(logFP, ", %g", p->avgPi[i]);
    }
    fprintf(logFP, " ) \n");

    fprintf(logFP, "   k_matrix: [");
    for( i = 0 ; i < sNo ; i++ ){
        fprintf(logFP, " ( %g", p->avgK[i][0]);
        for( j = 1 ; j < sNo ; j++ )
        {   fprintf(logFP, ", %g", p->avgK[i][j]);   }
        fprintf(logFP, ")");
    }
    fprintf(logFP, " ] \n\n");

    char fn[256];
    FILE *fp;
    size_t n;

    sprintf( fn, "%s.param%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        fprintf(fp, "I, pi");
        for( i = 0 ; i < sNo ; i++ )
        {   fprintf(fp, ", K%dx", i);   }
        fprintf(fp, "\n");

        for( i = 0 ; i < sNo ; i++ ){
            fprintf(fp, "%g, %g", p->avgI[i], p->avgPi[i]);
            for( j = 0 ; j < sNo ; j++ )
            {   fprintf(fp, ", %g", p->avgK[j][i]);   }
            fprintf(fp, "\n");
        }
        fclose(fp);
    }

    sprintf( fn, "%s.Lq%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < gv->iteration ; n++ ){
            fprintf( fp, "%24.20e\n", gv->LqArr[n] );
        }
        fclose(fp);
    }

    sprintf( fn, "%s.maxS%03d", xn->name, sNo );
    if( (fp = fopen( fn, "w")) != NULL ){
        for( n = 0 ; n < xn->N ; n++ ){
            fprintf( fp, "%d\n", iv->stateTraj[n] );
        }
        fclose(fp);
    }

}

//void outputTsResultsG( xns, gv, ivs, logFP )
//xnDataBundle *xns;
//globalVars *gv;
//indVarBundle *ivs;
//FILE *logFP;
//{
//}

//

max_stars_repo_name: okamoto-kenji/varBayes-HMM

max_issues_repo_name: okamoto-kenji/varBayes-HMM

max_forks_repo_name: okamoto-kenji/varBayes-HMM


text:
#ifndef OPENMC_TALLIES_FILTER_CELL_INSTANCE_H
#define OPENMC_TALLIES_FILTER_CELL_INSTANCE_H

#include <cstdint>
#include <unordered_map>
#include <vector>

#include <gsl/gsl>

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

namespace openmc {

//==============================================================================
//! Specifies cell instances that tally events reside in.
//==============================================================================

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

  CellInstanceFilter() = default;
  CellInstanceFilter(gsl::span<CellInstance> instances);
  ~CellInstanceFilter() = default;

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

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

  void from_xml(pugi::xml_node node) override;

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

  void to_statepoint(hid_t filter_group) const override;

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

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

  const std::vector<CellInstance>& cell_instances() const { return cell_instances_; }

  void set_cell_instances(gsl::span<CellInstance> instances);

private:
  //----------------------------------------------------------------------------
  // Data members

  //! The indices of the cells binned by this filter.
  std::vector<CellInstance> cell_instances_;

  //! A map from cell/instance indices to filter bin indices.
  std::unordered_map<CellInstance, gsl::index, CellInstanceHash> map_;
};

} // namespace openmc

#endif // OPENMC_TALLIES_FILTER_CELL_INSTANCE_H

max_stars_repo_name: PullRequest-Agent/openmc

max_issues_repo_name: norberto-schmidt/openmc-surftrack

max_forks_repo_name: norberto-schmidt/openmc-surftrack


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include <algorithm>
#include <boost/iostreams/categories.hpp>
#include <gsl/gsl>
#include <iosfwd>
#include <vector>

namespace mira
{
    class memory_device
    {
    public:
        using char_type = char;
        using category = boost::iostreams::seekable_device_tag;

        memory_device(std::vector<char>& memory, std::streamoff pos)
            : m_pos{ pos }
            , m_memory{ memory }
        {}

        std::streamoff position() const
        {
            return m_pos;
        }

        void reset()
        {
            m_pos = 0;
        }

        std::streamsize read(char* s, std::streamsize n)
        {
            // Read up to n characters from the underlying data source
            // into the buffer s, returning the number of characters
            // read; return -1 to indicate EOF
            std::streamsize read = std::min<std::streamsize>(m_memory.size() - m_pos, n);
            if (read == 0)
                return -1;

            std::copy(pointer(), pointer() + gsl::narrow_cast<ptrdiff_t>(read), s);
            m_pos += read;

            return read;
        }

        std::streamsize write(const char* s, std::streamsize n)
        {
            // Write up to n characters to the underlying
            // data sink into the buffer s, returning the
            // number of characters written
            std::streamsize written = 0ll;
            if (m_pos != gsl::narrow<std::streamoff>(m_memory.size()))
            {
                written = std::min<std::streamsize>(n, m_memory.size() - m_pos);
                std::copy(s,
                          s + gsl::narrow_cast<ptrdiff_t>(written),
                          m_memory.begin() + gsl::narrow_cast<ptrdiff_t>(m_pos));
                m_pos += written;
            }

            if (written < n)
            {
                m_memory.insert(m_memory.end(), s + written, s + n);
                m_pos = gsl::narrow<std::streamoff>(m_memory.size());
            }

            return n;
        }

        std::streamoff seek(std::streamoff off, std::ios_base::seekdir way)
        {
            // Seek to position off and return the new stream
            // position. The argument way indicates how off is
            // interpretted:
            //    - std::ios_base::beg indicates an offset from the
            //      sequence beginning
            //    - std::ios_base::cur indicates an offset from the
            //      current character position
            //    - std::ios_base::end indicates an offset from the
            //      sequence end

            std::streamoff next;
            if (way == std::ios_base::beg)
            {
                next = off;
            }
            else if (way == std::ios_base::cur)
            {
                next = m_pos + off;
            }
            else if (way == std::ios_base::end)
            {
                next = m_memory.size() + off - 1;
            }
            else
            {
                throw std::ios_base::failure("bad seek direction");
            }

            // Check for errors
            if (next < 0 || next > gsl::narrow<std::streamoff>(m_memory.size()))
            {
                throw std::ios_base::failure("bad seek offset");
            }

            m_pos = next;
            return m_pos;
        }

    private:
        char* pointer()
        {
            return &m_memory[gsl::narrow_cast<ptrdiff_t>(m_pos)];
        }

        std::streamoff m_pos{0};
        std::vector<char>& m_memory;
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
//
//  Rpc.hpp
//  PretendToWork
//
//  Created by tqtifnypmb on 18/12/2017.
//  Copyright © 2017 tqtifnypmb. All rights reserved.
//

#pragma once

#include <memory>
#include <vector>

#include <gsl/gsl>
#include <uv.h>

#include "Request.h"

namespace brick
{

class Rpc {
public:
    using RpcPeer = uv_handle_t;
    using Callback = std::function<void(RpcPeer*, Request req)>;

    enum class LoopState: int {
        looping = 0,
        closing = 1,
        closed,
        ready
    };

    Rpc(const char* ip, int port, Callback req_cb);
    Rpc(const Rpc&) = delete;
    Rpc& operator=(const Rpc&) = delete;
    ~Rpc();

    void send(RpcPeer* peer, const std::string& msg);
    void loop();
    void close();
    void close(RpcPeer* peer);

private:

    void onNewConnection(uv_tcp_t* client);
    void onNewMsg(uv_stream_t* client, const std::string& msg);
    void onHandleClosed(uv_handle_t* handle);

    static void connection_cb(uv_stream_t* server, int status);
    static void read_cb(uv_stream_t* stream, ssize_t nread, const uv_buf_t* buf);
    static void close_cb(uv_handle_t*);
    static void write_cb(uv_write_t* req, int status);

    LoopState state_;

    Callback req_cb_;
    std::vector<uv_tcp_t*> clients_;
    gsl::owner<uv_loop_t*> loop_;
    gsl::owner<uv_tcp_t*> server_;
};

}   // namespace brick

max_stars_repo_name: tqtifnypmb/brick

max_issues_repo_name: tqtifnypmb/brick

max_forks_repo_name: tqtifnypmb/brick


text:
//  This random generator is a C++ wrapper for the GNU Scientific Library
//  Copyright (C) 2001 Torbjorn Vik

//  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.
#ifndef __min_fminimizer_h
#define __min_fminimizer_h

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

namespace gsl{

//! Derive this class provide a user defined function for minimisation
struct min_f
{
	//! This operator must be overridden
	virtual double operator()(const double& x)=0;

	//! This is the function gsl calls to optimize f
	static double f(double x, void *p)
	{
		return (*(min_f *)p)(x);
	}
};

//! Class for minimizing one dimensional functions.
/*!
  Usage:
       - Create with optional minimize type
	   - Set with function object and inital bounds
	   - Loop the  iterate function until convergence or maxIterations (extra facility)

	   - Recover minimum and bounds
 */
class min_fminimizer
{
 public:
	//! choose between gsl_min_fminimizer_goldensection and gsl_min_fminimizer_brent
	min_fminimizer(const gsl_min_fminimizer_type* type=gsl_min_fminimizer_brent) : s(NULL), maxIterations(100), isSet(false)
	{
		s=gsl_min_fminimizer_alloc(type);
		nIterations=0;
		if (!s)
		{
			//error
			//cout << "ERROR Couldn't allocate memory for minimizer" << endl;
			//throw ?
			exit(-1);
		}
	}
	~min_fminimizer(){if (s) gsl_min_fminimizer_free(s);}
	//! returns GSL_FAILURE if the interval does not contain a minimum
	int set(min_f& function, double minimum, double x_lower, double x_upper)
	{
		isSet=false;
		f.function = &function.f;
		f.params = &function;
		int status=	gsl_min_fminimizer_set(s, &f, minimum, x_lower, x_upper);
		if (!status)
		{
			isSet=true;
			nIterations=0;
		}
		return status;
	}
	int set_with_values(min_f& function,
						double minimum, double f_minimum,
						double x_lower,double f_lower,
						double x_upper, double f_upper)
	{
		isSet=false;
		f.function = &function.f;
		f.params = &function;
		int status=	gsl_min_fminimizer_set_with_values(s, &f, minimum, f_minimum, x_lower, f_lower, x_upper, f_upper);
		if (!status)
		{
			isSet=true;
			nIterations=0;
		}
		return status;
	}
	int iterate()
	{
		assert_set();
		int status=gsl_min_fminimizer_iterate(s);
		nIterations++;
		if (status==GSL_FAILURE)
			isConverged=true;
		return status;
	}
	double minimum(){assert_set();return gsl_min_fminimizer_minimum(s);}
	double x_upper(){assert_set();return gsl_min_fminimizer_x_upper(s);}
	double x_lower(){assert_set();return gsl_min_fminimizer_x_lower(s);}
	void SetMaxIterations(int n){maxIterations=n;}
	int GetNIterations(){return nIterations;}
	bool is_converged(){if (nIterations>=maxIterations) return true; if (isConverged) return true; return false;}
	//string name() const;

 private:
	void assert_set(){if (!isSet)exit(-1);} // Old problem of error handling: TODO

	bool isSet;
	bool isConverged;
	int nIterations;
	int maxIterations;
	gsl_min_fminimizer* s;
	gsl_function f;
};
};	 // namespace gsl

#endif //__min_fminimizer_h

max_stars_repo_name: entn-at/GlottDNN

max_issues_repo_name: entn-at/GlottDNN

max_forks_repo_name: entn-at/GlottDNN


text:
/**
 *
 * @file pdpotrf.c
 *
 *  PLASMA auxiliary routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Jakub Kurzak
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:45:10 2014
 *
 **/
#include "common.h"
#ifdef USE_MKL
#include <mkl_lapacke.h>
#else
#include <lapacke.h>
#endif
#if defined(USE_OMPEXT)
#include <omp_ext.h>
#endif

#define A(m,n) BLKADDR(A, double, m, n)

/***************************************************************************//**
 *  Parallel tile Cholesky factorization - dynamic scheduling
 **/
void plasma_pdpotrf_quark(PLASMA_enum uplo, PLASMA_desc A)
{
    int k, m, n;
    int ldak, ldam;
    int tempkm, tempmm;

    double zone  = (double) 1.0;
    double mzone = (double)-1.0;
    /*
     *  PlasmaLower
     */
    if (uplo == PlasmaLower) {
        abort();
    }
    /*
     *  PlasmaUpper
     */
    else {
        for (k = 0; k < A.nt; k++) {
            tempkm = k == A.nt-1 ? A.n-k*A.nb : A.nb;
            ldak = BLKLDD(A, k);
            double *dA = A(k, k);
#if defined(USE_OMPEXT)
omp_set_task_priority(1);
#endif
#pragma omp task depend(inout:dA[0:A.mb*A.mb])
            {
                LAPACKE_dpotrf_work(LAPACK_COL_MAJOR, lapack_const(PlasmaUpper), tempkm, dA, ldak);
            }

            for (m = k+1; m < A.nt; m++) {
                tempmm = m == A.nt-1 ? A.n-m*A.nb : A.nb;
                double *dA = A(k, k);
                double *dB = A(k, m);
#pragma omp task depend(in:dA[0:A.mb*A.mb]) depend(inout:dB[0:A.mb*A.mb])
                cblas_dtrsm(
                        CblasColMajor,
                        (CBLAS_SIDE)PlasmaLeft, (CBLAS_UPLO)PlasmaUpper,
                        (CBLAS_TRANSPOSE)PlasmaTrans, (CBLAS_DIAG)PlasmaNonUnit,
                        A.mb, tempmm,
                        zone, dA, ldak,
                        dB, ldak);
            }
            for (m = k+1; m < A.nt; m++) {
                tempmm = m == A.nt-1 ? A.n-m*A.nb : A.nb;
                ldam = BLKLDD(A, m);
                double *dA = A(k, m);
                double *dB = A(m, m);
#pragma omp task depend(in:dA[0:A.mb*A.mb]) depend(inout:dB[0:A.mb*A.mb])
                {
                    cblas_dsyrk(
                            CblasColMajor,
                            (CBLAS_UPLO)PlasmaUpper, (CBLAS_TRANSPOSE)PlasmaTrans,
                            tempmm, A.mb,
                            (-1.0), dA, ldak,
                            (1.0), dB, ldam);
                }

                for (n = k+1; n < m; n++) {
                    double *dA = A(k , n);
                    double *dB = A(k , m);
                    double *dC = A(n , m);
#pragma omp task depend(in:dA[0:A.mb*A.mb], dB[0:A.mb*A.mb]) depend(inout:dC[0:A.mb*A.mb])
                    cblas_dgemm(CblasColMajor, (CBLAS_TRANSPOSE)PlasmaTrans, (CBLAS_TRANSPOSE)PlasmaNoTrans,
                            A.mb, tempmm, A.mb,
                            mzone, dA, ldak,
                            dB, ldak,
                            zone, dC, A.mb);
                }
            }
        }
    }
}


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2018-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */

#pragma once

#include <mcbp/protocol/response.h>
#include <mcbp/protocol/status.h>
#include <nlohmann/json_fwd.hpp>
#include <platform/thread.h>
#include <chrono>
#include <gsl/gsl>
#include <mutex>
#include <queue>
#include <unordered_map>
#include <vector>

class Connection;
class AuthnAuthzServiceTask;

/**
 * The ExternalAuthManagerThread class takes care of the scheduling between
 * the authentication requests and the responses received.
 *
 * One of the problems is that we try to communicate between multiple
 * frontend threads, which could cause deadlocks. To avoid locking
 * problems we're using a dedicated thread to communicate with the
 * other threads.
 */
class ExternalAuthManagerThread : public Couchbase::Thread {
public:
    ExternalAuthManagerThread() : Couchbase::Thread("mcd:ext_auth") {
    }
    ExternalAuthManagerThread(const ExternalAuthManagerThread&) = delete;

    /**
     * The named (external) user logged on to the system
     */
    void login(const std::string& user);

    /**
     * The named (external) user logged off the system
     */
    void logoff(const std::string& user);

    /**
     * Get the list of active users defined in the system
     */
    nlohmann::json getActiveUsers() const;

    /**
     * Add the provided connection to the list of available authentication
     * providers
     */
    void add(Connection& connection);

    /**
     * Remove the connection from the list of available providers (the
     * connection may not be registered)
     */
    void remove(Connection& connection);

    /**
     * Enqueue a SASL request to the list of requests to send to the
     * RBAC provider. NOTE: This method must _NOT_ be called from the
     * worker threads as that may result in a deadlock (should be
     * called from the tasks execute() method when running on the
     * executor threads)
     *
     * @param request the task to notify when the result comes back
     */
    void enqueueRequest(AuthnAuthzServiceTask& request);

    /**
     * Received an authentication response on the provided connection
     *
     * @param connection the connection the message was received on
     * @param response
     */
    void responseReceived(const cb::mcbp::Response& response);

    /**
     * Request the authentication daemon to stop
     */
    void shutdown();

    void setPushActiveUsersInterval(std::chrono::microseconds interval) {
        activeUsersPushInterval.store(interval);
        condition_variable.notify_one();
    }

    void setRbacCacheEpoch(std::chrono::steady_clock::time_point tp);

    /// Check to see if we've got an up to date RBAC entry for the user
    bool haveRbacEntryForUser(const std::string& user) const;

protected:
    /// The main loop of the thread
    void run() override;

    /**
     * Process all of the requests added by the various SaslTasks by putting
     * them into the queue of messages to send to the Authentication Provider.
     */
    void processRequestQueue();
    /**
     * Process all of the response messages enqueued and notify the SaslTask
     */
    void processResponseQueue();

    /**
     * Iterate over all of the connections marked as closed. Generate
     * error responses for all outstanding requests; decrement the reference
     * count and tell the thread to complete its shutdown logic.
     */
    void purgePendingDeadConnections();

    /// Push the list of active users to the authentication provider
    void pushActiveUsers();

    /// Let the daemon thread run as long as this member is set to true
    bool running = true;

    /// The next value to use in the opaque field
    uint32_t next = 0;

    /// The map between the opaque field being used and the SaslAuthTask
    /// requested the operation
    /// Second should be a std pair of a connection and a task, so that
    /// when we remove a connection we can iterate over the entire
    /// map and create aborts for the one we don't have anymore (or
    /// we could redistribute them O:)
    std::unordered_map<uint32_t, std::pair<Connection*, AuthnAuthzServiceTask*>>
            requestMap;

    /// The mutex variable used to protect access to _all_ the internal
    /// members
    std::mutex mutex;
    /// The daemon thread will block and wait on this variable when there
    /// isn't any work to do
    std::condition_variable condition_variable;

    /// A list of available connections to use as authentication providers
    std::vector<Connection*> connections;

    /**
     * All of the various tasks enqueue their SASL request in this queue
     * to allow the authentication provider thread to inject them into
     * the worker thread (to avoid a possible deadlock)
     */
    std::queue<AuthnAuthzServiceTask*> incomingRequests;

    /**
     * We need to pass the response information from one of the frontend
     * thread back to the sasl task, and we'd like to use the same
     * logic when we're terminating the connections.
     */
    struct AuthResponse {
        AuthResponse(uint32_t opaque, std::string payload)
            : opaque(opaque),
              status(cb::mcbp::Status::Etmpfail),
              payload(std::move(payload)) {
        }
        AuthResponse(uint32_t opaque,
                     cb::mcbp::Status status,
                     cb::const_byte_buffer value)
            : opaque(opaque),
              status(status),
              payload(reinterpret_cast<const char*>(value.data()),
                      value.size()) {
        }
        uint32_t opaque;
        cb::mcbp::Status status;
        std::string payload;
    };

    /**
     * All of the various responses is stored within this queue to
     * be handled by the daemon thread
     */
    std::queue<std::unique_ptr<AuthResponse>> incommingResponse;

    std::vector<Connection*> pendingRemoveConnection;

    class ActiveUsers {
    public:
        void login(const std::string& user);

        void logoff(const std::string& user);

        nlohmann::json to_json() const;

    private:
        mutable std::mutex mutex;
        std::unordered_map<std::string, uint32_t> users;
    } activeUsers;

    /**
     * The interval we'll push the current active users to the external
     * authentication service.
     *
     * The list can't be considered as "absolute" as memcached is a
     * highly multithreaded application and we may have disconnects
     * happening in multiple threads at the same time and we might
     * have login messages waiting to acquire the lock for the map.
     * The intention with the list is that we'll poll the users from
     * LDAP to check if the group membership has changed since the last
     * time we checked. Given that LDAP is an external system and it
     * can't push change notifications to our system our guarantee
     * is that the change is reflected in memcached at within 2 times
     * of the interval we push these users.
     *
     * We don't need microseconds resolution "in production", but
     * we don't want to have our unit tests wait multiple seconds
     * for this message to be sent
     */
    std::atomic<std::chrono::microseconds> activeUsersPushInterval{
            std::chrono::minutes(5)};

    /**
     * The time point we last sent the list of active users to the auth
     * provider
     */
    std::chrono::steady_clock::time_point activeUsersLastSent;

    /**
     * It should be possible for the authentication provider to
     * invalidate the entire cache, but we can't simply drop the
     * current cache as that would cause all of the connected clients
     * to be disconnected. Instead the authentication provider may
     * tell us to start requesting a new RBAC entry if the cached
     * one is older than a given time point. Combined with the logic
     * that we'll push all of the connected users the authentication
     * provider can make sure that all of the connections use the
     * correct RBAC definition. Ideally we would have stored this
     * as a std::chrono::steady_clock::TimePoint, but that don't
     * play well with std::atomic, so we're storing the raw
     * number of seconds instead.
     */
    std::atomic<uint64_t> rbacCacheEpoch{{}};
};

extern std::unique_ptr<ExternalAuthManagerThread> externalAuthManager;


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <time.h>
#include <sys/time.h>
#include <gsl/gsl_math.h>

#include "cosmocalc.h"

double wtime(void)
{
  struct timeval tp;
  gettimeofday(&tp,NULL);
  return ((double) (tp.tv_sec)) + ((double) (tp.tv_usec))/1e6;
}

#define EPS 3.0e-11
void gauleg(double x1, double x2, double x[], double w[], int n)
{
  int m,j,i;
  double z1,z,xm,xl,pp,p3,p2,p1;

  m=(n+1)/2;
  xm=0.5*(x2+x1);
  xl=0.5*(x2-x1);
  for (i=1;i<=m;i++) {
    z=cos(3.141592654*(i-0.25)/(n+0.5));
    do {
      p1=1.0;
      p2=0.0;
      for (j=1;j<=n;j++) {
	p3=p2;
	p2=p1;
	p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
      }
      pp=n*(z*p1-p2)/(z*z-1.0);
      z1=z;
      z=z1-p1/pp;
    } while (fabs(z-z1) > EPS);
    x[i]=xm-xl*z;
    x[n+1-i]=xm+xl*z;
    w[i]=2.0*xl/((1.0-z*z)*pp*pp);
    w[n+1-i]=w[i];
  }
}
#undef EPS

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text:
/*
**  parse second level design file
**
** G.Lohmann
*/
#include <viaio/Vlib.h>
#include <viaio/VImage.h>

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

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

extern int VStringToken(char *, char *, int, int);

#define LEN  10000

int test_ascii(int val)
{
  if (val >= 'a' && val <= 'z') return 1;
  if (val >= 'A' && val <= 'Z') return 1;
  if (val >= '0' && val <= '9') return 1;
  if (val ==  ' ') return 1;
  if (val == '\0') return 1;
  if (val == '\n') return 1;
  if (val == '\r') return 1;
  if (val == '\t') return 1;
  if (val == '\v') return 1;
  if (val == '-') return 1;
  if (val == '+') return 1;
  if (val == '.') return 1;
  return 0;
}

int line_empty(char *buf,int len)
{
  int i;
  for (i=0; i<len; i++) {
    if (buf[i] != ' ' && buf[i] != '\n' && buf[i] != 0) return 0;
  }
  return 1;
}

int CheckBuffer(char *buf,int len)
{
  int j;

  if(strlen(buf) < 1) return 0;
  if (buf[0] == '%' || buf[0] == '#' || buf[0] == '/' || buf[0] == '\n') return 0;

  /* remove tabs */
  for (j=0; j<len; j++) {
    if (buf[0] == '\t') buf[j] = ' ';
  }
  if (line_empty(buf,len) > 0) return 0;
  return 1;
}

int VistaFormat(char *buf,int len)
{
  if (strncmp(buf,"V-data",6) == 0) return 1;
  return 0;
}

/*
** read a 2nd level design file
*/
gsl_matrix *XRead2ndLevel(VString filename)
{
  int i, j, nrows, ncols;
  double val;
  char buf[LEN], token[32];

  fprintf(stderr," Reading %s\n",filename);

  FILE *fp = VOpenInputFile (filename, TRUE);
  if (!fp) VError(" error opening %s",filename);

  nrows = ncols = 0;
  while(!feof(fp)) {
    memset(buf, 0, LEN);
    if(!fgets(buf, LEN, fp)) continue;
    if (VistaFormat(buf,LEN) > 0) VError(" Design file must be a text file");
    if (CheckBuffer(buf,LEN) < 1) continue;
    if (! test_ascii((int)buf[0])) VError(" Design file must be a text file");

    j = 0;
    while(VStringToken(buf, token, j, 30)) {
      if(!sscanf(token, "%lf", &val))
	VError("illegal text string in design file: %s", buf);
      j++;
    }

    if (ncols == 0) ncols = j;
    if (j < 1) continue;
    else if(ncols != j)
      VError(" inconsistent number of columns in row %d", nrows + 1);
    nrows++;
  }
  rewind(fp);

  /* fill design matrix */
  gsl_matrix *X = gsl_matrix_calloc(nrows,ncols);

  i = 0;
  while(!feof(fp)) {
    memset(buf, 0, LEN);
    if(!fgets(buf, LEN, fp))  continue;
    if (CheckBuffer(buf,LEN) < 1) continue;

    j = 0;
    while(VStringToken(buf, token, j, 30)) {
      sscanf(token, "%lf", &val);
      gsl_matrix_set(X,i,j,val);
      j++;
    }
    if (j < 1) continue;
    i++;
  }
  fclose(fp);

  return X;
}

/* read exchangeability file */
void XReadExchange(VString filename,int *exchange,int n)
{
  int len=1024;
  char token[32];
  char *buf = (char *)VCalloc(len,sizeof(char));
  FILE *fp = VOpenInputFile (filename, TRUE);
  if (!fp) VError(" error opening %s",filename);
  fprintf(stderr," Reading %s\n",filename);

  int nrows = 0;
  while(!feof(fp)) {
    memset(buf, 0, len);
    if(!fgets(buf, len, fp)) break;
    if (VistaFormat(buf,LEN) > 0) VError(" Exchangeability file must be a text file");
    if (CheckBuffer(buf,LEN) < 1) continue;
    if (! test_ascii((int)buf[0])) VError(" Exchangeability file must be a text file");
    nrows++;
  }
  rewind(fp);

  if (nrows != n)
    VError(" Number of rows in exchangeability file (%d) does not match number of rows in design file (%d)",nrows,n);

  int i=0,j=0,k=0;
  while(!feof(fp)) {
    memset(buf, 0, len);
    if(!fgets(buf, len, fp))  break;
    if (CheckBuffer(buf,LEN) < 1) continue;

    k = 0;
    while(VStringToken(buf, token, k, 30)) {
      if (!sscanf(token, "%d", &j))
	VError("illegal text string in exchangeability file: %s", buf);
      if (j < 0) VError(" Values in exchangeability file must be non-negative");
      if (j > n) VError(" Illegal value (%d) in exchangeability file, values must be <= %d",j,n);
      exchange[i] = j;
      k++;
    }
    if (k > 1) VError(" Exchangeability file must have exactly one column");
    i++;
  }
  VFree(buf);
  fclose(fp);
}

/*
int main()
{
  int i,j;
  gsl_matrix *X = XRead2ndLevel("bla.mat");

  fprintf(stderr,"X:\n");
  for (i=0; i<X->size1; i++) {
    for (j=0; j<X->size2; j++) {
      fprintf(stderr," %.2f",gsl_matrix_get(X,i,j));
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");
  exit(0);
}
*/

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/* Copyright (c) 2011-2012, Jérémy Fix. 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. */
/* * None of 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 AUTHOR 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 UKF_SAMPLES_H
#define UKF_SAMPLES_H

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sort_vector_double.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_randist.h>

#include <time.h>

#include <iostream>

namespace ukf
{
    namespace samples
    {
        /**
          * @short Generate 1D samples according to a uniform distribution :
          * @brief In fact, gsl_ran_choose would do the job!!
          */
        class RandomSample_1D
        {
        public:
            static void generateSample(gsl_vector * samples, unsigned int nb_samples, double min, double max, double delta)
            {
                if(samples->size != nb_samples)
                {
                    std::cerr << "[ERROR] RandomSample_1D : samples should be allocated to the size of the requested number of samples" << std::endl;
                    return;
                }

                for(unsigned int i = 0 ; i < nb_samples ; i++)
                    gsl_vector_set(samples, i, min + floor((max - min)*rand()/double(RAND_MAX)/delta)    * delta  );
            }
        };

        /**
          * @short Generate 2D samples according to a uniform distribution
          * and put them alternativaly at the odd/even positions
          */
        class RandomSample_2D
        {
        public:
            static void generateSample(gsl_vector * samples, unsigned int nb_samples, double minx, double maxx, double deltax, double miny, double maxy, double deltay)
            {
                if(samples->size != 2*nb_samples)
                {
                    std::cerr << "[ERROR] RandomSample_2D : samples should be allocated to the size of the 2 times the requested number of samples" << std::endl;
                    return;
                }
		gsl_vector_view vec_view = gsl_vector_subvector_with_stride (samples, 0, 2, nb_samples);
                RandomSample_1D::generateSample(&vec_view.vector, nb_samples, minx, maxx, deltax);
		vec_view = gsl_vector_subvector_with_stride (samples, 1, 2, nb_samples);
                RandomSample_1D::generateSample(&vec_view.vector, nb_samples, miny, maxy, deltay);
            }
        };

        /**
          * @short Generate 3D samples according to a uniform distribution
          * and put them alternativaly at the odd/even positions
          */
        class RandomSample_3D
        {
        public:
            static void generateSample(gsl_vector * samples, unsigned int nb_samples, double minx, double maxx, double deltax, double miny, double maxy, double deltay, double minz, double maxz, double deltaz)
            {
                if(samples->size != 3*nb_samples)
                {
                    std::cerr << "[ERROR] RandomSample_2D : samples should be allocated to the size of the 2 times the requested number of samples" << std::endl;
                    return;
                }
		gsl_vector_view vec_view = gsl_vector_subvector_with_stride (samples, 0, 3, nb_samples);
                RandomSample_1D::generateSample(&vec_view.vector, nb_samples, minx, maxx, deltax);
		vec_view = gsl_vector_subvector_with_stride (samples, 1, 3, nb_samples);
                RandomSample_1D::generateSample(&vec_view.vector, nb_samples, miny, maxy, deltay);
		vec_view = gsl_vector_subvector_with_stride (samples, 2, 3, nb_samples);
                RandomSample_1D::generateSample(&vec_view.vector, nb_samples, minz, maxz, deltaz);
            }
        };

        /**
          * @short Extract the indexes of the nb_samples highest values of a vector
          * Be carefull, this function is extracting the indexes as, the way it is used, it doesn't not to which sample a vector index corresponds
          */
        class MaximumIndexes_1D
        {
        public:
            static void generateSample(gsl_vector * v, gsl_vector * indexes_samples, unsigned int nb_samples)
            {
                if(indexes_samples->size != nb_samples)
                {
                    std::cerr << "[ERROR] MaximumIndexes_1D : samples should be allocated to the size of the requested number of samples" << std::endl;
                    return;
                }
                if(nb_samples > v->size)
                {
                    std::cerr << "[ERROR] MaximumIndexes_1D : the vector v must hold at least nb_samples elements" << std::endl;
                    return;
                }

                gsl_permutation * p = gsl_permutation_alloc(v->size);
                gsl_sort_vector_index(p,v);
                for(unsigned int i = 0 ; i < nb_samples ; i++)
                    gsl_vector_set(indexes_samples, i, gsl_permutation_get(p, p->size - i - 1));
                gsl_permutation_free(p);
            }
        };

        /**
          * @short Extract the indexes of the nb_samples highest values of a matrix
          * Be carefull, this function is extracting the indexes as, the way it is used, it doesn't not to which sample a vector index corresponds
          */
        class MaximumIndexes_2D
        {
        public:
            static void generateSample(gsl_matrix * m, gsl_vector * indexes_samples, unsigned int nb_samples)
            {

                if(indexes_samples->size != 2*nb_samples)
                {
                    std::cerr << "[ERROR] MaximumIndexes_2D : samples should be allocated to the size of 2 times the requested number of samples" << std::endl;
                    return;
                }
                if(nb_samples > m->size1 * m->size2 )
                {
                    std::cerr << "[ERROR] MaximumIndexes_2D : the vector v must hold at least nb_samples elements" << std::endl;
                    return;
                }

                gsl_permutation * p = gsl_permutation_alloc(m->size1 * m->size2);
                // We cast the matrix in a vector ordered in row-major and get the permutation to order the data by increasing value
		gsl_vector_view vec_view = gsl_vector_view_array(m->data,m->size1 * m->size2);
                gsl_sort_vector_index(p,&vec_view.vector);

                // The indexes stored in p are in row-major order
                // we then convert the nb_samples first indexes into matrix indexes and copy them in indexes_samples
                int k,l;
                for(unsigned int i = 0 ; i < nb_samples ; i++)
                {
                    // Line
                    k = gsl_permutation_get(p,p->size - i - 1) / m->size2;
                    // Column
                    l = gsl_permutation_get(p,p->size - i - 1) % m->size2;
                    gsl_vector_set(indexes_samples, 2 * i, k);
                    gsl_vector_set(indexes_samples, 2 * i + 1, l);
                }

                gsl_permutation_free(p);
            }
        };

        /**
          * @short Generate 1D samples according to a discrete vectorial distribution
          */
        class DistributionSample_1D
        {
        public:
            static void generateSample(gsl_vector * v, gsl_vector * indexes_samples, unsigned int nb_samples)
            {
                if(indexes_samples->size != nb_samples)
                {
                    std::cerr << "[ERROR] RandomSample_1D : samples should be allocated to the size of the requested number of samples" << std::endl;
                    return;
                }

                // Generating nb_samples samples from an array of weights for each index is easily done in the gsl with gsl_ran_discrete_* methods
                // the weights v do not need to add up to one
                gsl_rng * r = gsl_rng_alloc(gsl_rng_default);
                // Init the random seed
                gsl_rng_set(r, time(NULL));
                gsl_ran_discrete_t * ran_pre = gsl_ran_discrete_preproc(v->size, v->data);

                for(unsigned int i = 0 ; i < nb_samples ; i++)
                    gsl_vector_set(indexes_samples, i , gsl_ran_discrete(r, ran_pre));

                gsl_ran_discrete_free(ran_pre);
                gsl_rng_free(r);
            }
        };

        /**
          * @short Generate 2D samples according to a discrete matricial distribution
          */
        class DistributionSample_2D
        {
        public:
            static void generateSample(gsl_matrix * m, gsl_vector * indexes_samples, unsigned int nb_samples)
            {
                if(indexes_samples->size != 2*nb_samples)
                {
                    std::cerr << "[ERROR] RandomSample_1D : samples should be allocated to the size of the requested number of samples" << std::endl;
                    return;
                }

                // To generate 2D samples according to a matrix of weights
                // we simply cast the matrix into a vector, ask DistributionSample_1D to do the job
                // and then convert the linear row-major indexes back into 2D indexes

                gsl_vector * samples_tmp = gsl_vector_alloc(nb_samples);
		gsl_vector_view vec_view = gsl_vector_view_array(m->data,m->size1 * m->size2);
                DistributionSample_1D::generateSample(&vec_view.vector, samples_tmp, nb_samples);

                int k,l;
                for(unsigned int i = 0 ; i < nb_samples ; i++)
                {
                    k = int(gsl_vector_get(samples_tmp,i)) / m->size2;
                    l = int(gsl_vector_get(samples_tmp,i)) % m->size2;
                    gsl_vector_set(indexes_samples, 2*i, k);
                    gsl_vector_set(indexes_samples, 2*i+1, l);
                }

                gsl_vector_free(samples_tmp);
            }
        };
    }
}

#endif // UKF_SAMPLES_H