text:
#pragma once

#include <utility>
#include <iostream>
#include <initializer_list>
#include <exception>

#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>

#include "utils/fcmp.h"

namespace gsl_wrapper
{

  class Vector
  {
  public:
    // Constructors and destructor
    Vector(size_t vec_size);
    Vector(gsl_vector *gsl_vec_ptr);

    Vector(const Vector &copy_from);
    Vector(Vector &&move_from);

    Vector(std::initializer_list<double> args);

    ~Vector();

    // Member functions
    auto get_gsl_vector() const -> gsl_vector *;
    auto size() const -> size_t;
    auto begin() const -> double *;
    auto end() const -> double *;

    // Operators
    auto operator=(const Vector &copy_from) -> Vector &;
    auto operator=(Vector &&move_from) -> Vector &;

    auto operator==(const Vector &comparasion_vector) -> bool;
    auto operator!=(const Vector &comparasion_vector) -> bool;

    auto operator[](const size_t index) -> double &;
    auto operator[](const size_t index) const -> const double &;

    auto operator+(const Vector &add_to) const -> Vector;
    auto operator-(const Vector &sub) const -> Vector;
    auto operator*(const double number) const -> Vector;

    // Friend declarations
    friend auto operator<<(std::ostream &stream, const Vector &to_print) -> std::ostream &;
    friend auto operator*(const double number, const Vector &vec) -> Vector;

  private:
    gsl_vector *m_vector_ptr;
    size_t m_vector_size;
  };

  inline Vector::Vector(size_t vec_size)
      : m_vector_ptr{gsl_vector_calloc(vec_size)},
        m_vector_size{vec_size}
  {
  }

  inline Vector::Vector(gsl_vector *gsl_vec_ptr)
      : m_vector_ptr{gsl_vec_ptr},
        m_vector_size{gsl_vec_ptr->size}
  {
  }

  inline Vector::Vector(const Vector &copy_from)
      : m_vector_ptr{gsl_vector_calloc(copy_from.m_vector_size)},
        m_vector_size{copy_from.m_vector_size}
  {
    gsl_vector_memcpy(m_vector_ptr, copy_from.m_vector_ptr);
  }

  inline Vector::Vector(Vector &&move_from)
      : m_vector_ptr{std::exchange(move_from.m_vector_ptr, nullptr)},
        m_vector_size{std::exchange(move_from.m_vector_size, 0)}
  {
  }

  inline Vector::Vector(std::initializer_list<double> args)
      : m_vector_ptr{gsl_vector_calloc(args.size())},
        m_vector_size{args.size()}
  {
    size_t i = 0;
    for (auto &&el : args)
    {
      gsl_vector_set(m_vector_ptr, i++, el);
    }
  }

  inline Vector::~Vector()
  {
    gsl_vector_free(m_vector_ptr);
  }

  inline auto Vector::get_gsl_vector() const -> gsl_vector *
  {
    return m_vector_ptr;
  }

  inline auto Vector::size() const -> size_t
  {
    return m_vector_size;
  }

  inline auto Vector::begin() const -> double *
  {
    return m_vector_ptr->data;
  }

  inline auto Vector::end() const -> double *
  {
    return m_vector_ptr->data + m_vector_size;
  }

  inline auto Vector::operator=(const Vector &copy_from) -> Vector &
  {
    // Prevent self copy
    if (m_vector_ptr == copy_from.m_vector_ptr)
      return *this;

    gsl_vector_free(m_vector_ptr);
    m_vector_ptr = gsl_vector_calloc(copy_from.m_vector_size);
    gsl_vector_memcpy(m_vector_ptr, copy_from.m_vector_ptr);
    m_vector_size = copy_from.m_vector_size;

    return *this;
  }

  inline auto Vector::operator=(Vector &&move_from) -> Vector &
  {
    // Prevent self move
    if (m_vector_ptr == move_from.m_vector_ptr)
      return *this;

    gsl_vector_free(m_vector_ptr);
    m_vector_ptr = std::exchange(move_from.m_vector_ptr, nullptr);
    m_vector_size = std::exchange(move_from.m_vector_size, 0);

    return *this;
  }

  inline auto Vector::operator==(const Vector &comparasion_vector) -> bool
  {
    if (m_vector_size != comparasion_vector.m_vector_size)
      return false;

    for (size_t i = 0; i < m_vector_size; i++)
    {
      bool test = ::gsl_wrapper::utils::equal((*this)[i], comparasion_vector[i]);
      if (!test)
        return false;
    }
    return true;
  }

  inline auto Vector::operator!=(const Vector &comparasion_vector) -> bool
  {
    return !(*this == comparasion_vector);
  }

  inline auto Vector::operator[](const size_t index) -> double &
  {
    if (index >= m_vector_size)
      throw std::range_error{"Accesing vector elements out of bounds"};
    return *gsl_vector_ptr(m_vector_ptr, index);
  }

  inline auto Vector::operator[](const size_t index) const -> const double &
  {
    return *gsl_vector_const_ptr(m_vector_ptr, index);
  }

  inline auto operator<<(std::ostream &stream, const Vector &to_print) -> std::ostream &
  {
    for (size_t i = 0; i < to_print.m_vector_size - 1; i++)
    {
      stream << to_print[i] << " ";
    }
    stream << to_print[to_print.m_vector_size - 1];

    return stream;
  }

  inline auto Vector::operator+(const Vector &add_to) const -> Vector
  {
    if (m_vector_size != add_to.m_vector_size)
      throw std::range_error{"Adding vector of diffrent sizes"};

    Vector result(m_vector_size);

    for (size_t i = 0; i < m_vector_size; i++)
    {
      result[i] = (*this)[i] + add_to[i];
    }

    return result;
  }

  inline auto Vector::operator-(const Vector &sub) const -> Vector
  {
    if (m_vector_size != sub.m_vector_size)
      throw std::range_error{"Subtracting vector of diffrent sizes"};

    Vector result(m_vector_size);

    for (size_t i = 0; i < m_vector_size; i++)
    {
      result[i] = (*this)[i] - sub[i];
    }

    return result;
  }

  inline auto Vector::operator*(const double number) const -> Vector
  {
    Vector result = *this;
    for (auto &&el : result)
    {
      el *= number;
    }

    return result;
  }

  inline auto operator*(const double number, const Vector &vec) -> Vector
  {
    Vector result = vec;
    for (auto &&el : result)
    {
      el *= number;
    }

    return result;
  }

}


text:
//Normalizes RMS of each vector in X according to dim.
//See C++ command-line code for more info.

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <lapacke.h>

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

int rms_scale_s (float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const float sr, const float tau, const float target_dB_SPL, const float max_dB_SPL);
int rms_scale_d (double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const double sr, const double tau, const double target_dB_SPL, const double max_dB_SPL);

int rms_scale_s (float *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const float sr, const float tau, const float target_dB_SPL, const float max_dB_SPL)
{
    if (dim>3u) { fprintf(stderr,"error in rms_scale_s: dim must be in [0 3]\n"); return 1; }
    if (sr<=0.0f) { fprintf(stderr,"error in rms_scale_s: sr (sample rate) must be positive\n"); return 1; }
    if (tau<=0.0f) { fprintf(stderr,"error in rms_scale_s: tau (time constant) must be positive\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t L = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;

    if (N<2u || L<2u) {}
    else
    {
        const float rms_prctile = 0.9f;
        const float target_RMS = 20e-6f * powf(10.0f, target_dB_SPL/20.0f);
        const float max_RMS = 20e-6f * powf(10.0f, max_dB_SPL/20.0f);
        const float a1 = expf(-1.0f/(sr*tau));
        const float b0 = 1.0f - a1;
        const size_t tranz = (L>5u) ? 5u : 0u;  //nsamps to zero in case transient
        const size_t rms_p = (size_t)roundf((1.0f-rms_prctile)*(float)L);  //index for rms_prctile
        float rms, scale, scalemx;

        float *Y;
        if (!(Y=(float *)malloc(L*sizeof(float)))) { fprintf(stderr,"error in rms_scale_s: problem with malloc. "); perror("malloc"); return 1; }

        if (L==N)
        {
            *Y = b0 * *X * *X;
            ++X; ++Y;
            for (size_t l=L-1u; l>0u; --l, ++X, ++Y) { *Y = b0**X**X + a1**(Y-1); }
            Y -= L;
            for (size_t l=tranz; l>0u; --l, ++Y) { *Y = 0.0f; }  //in case transient
            Y -= tranz;
            if (LAPACKE_slasrt_work('D',(int)L,Y)) { fprintf(stderr,"error in rms_scale_s: problem with LAPACKE function\n"); }
            rms = sqrtf(*Y);
            scalemx = (rms>FLT_EPSILON) ? max_RMS/rms : 1.0f;
            rms = sqrtf(*(Y+rms_p));
            scale = (rms>FLT_EPSILON) ? target_RMS/rms : 1.0f;
            scale = (scale<scalemx) ? scale : scalemx;
            for (size_t l=L; l>0u; --l) { *--X *= scale; }
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/L, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=V; v>0u; --v, X+=L)
                {
                    *Y = b0 * *X * *X;
                    ++X; ++Y;
                    for (size_t l=L-1u; l>0u; --l, ++X, ++Y) { *Y = b0**X**X + a1**(Y-1); }
                    Y -= L;
                    for (size_t l=tranz; l>0u; --l, ++Y) { *Y = 0.0f; }  //in case transient
                    Y -= tranz;
                    if (LAPACKE_slasrt_work('D',(int)L,Y)) { fprintf(stderr,"error in rms_scale_s: problem with LAPACKE function\n"); }
                    rms = sqrtf(*Y);
                    scalemx = (rms>FLT_EPSILON) ? max_RMS/rms : 1.0f;
                    rms = sqrtf(*(Y+rms_p));
                    scale = (rms>FLT_EPSILON) ? target_RMS/rms : 1.0f;
                    scale = (scale<scalemx) ? scale : scalemx;
                    for (size_t l=L; l>0u; --l) { *--X *= scale; }
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B*(L-1u))
                {
                    for (size_t b=B; b>0u; --b, ++X)
                    {
                        *Y = b0 * *X * *X;
                        X+=K; ++Y;
                        for (size_t l=L-1u; l>0u; --l, X+=K, ++Y) { *Y = b0**X**X + a1**(Y-1); }
                        Y -= L;
                        for (size_t l=tranz; l>0u; --l, ++Y) { *Y = 0.0f; }  //in case transient
                        Y -= tranz;
                        if (LAPACKE_slasrt_work('D',(int)L,Y)) { fprintf(stderr,"error in rms_scale_s: problem with LAPACKE function\n"); }
                        rms = sqrtf(*Y);
                        scalemx = (rms>FLT_EPSILON) ? max_RMS/rms : 1.0f;
                        rms = sqrtf(*(Y+rms_p));
                        scale = (rms>FLT_EPSILON) ? target_RMS/rms : 1.0f;
                        scale = (scale<scalemx) ? scale : scalemx;
                        for (size_t l=L; l>0u; --l) { X-=K; *X *= scale; }
                    }
                }
            }
        }

        free(Y);
    }

    return 0;
}

int rms_scale_d (double *X, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t dim, const double sr, const double tau, const double target_dB_SPL, const double max_dB_SPL)
{
    if (dim>3u) { fprintf(stderr,"error in rms_scale_d: dim must be in [0 3]\n"); return 1; }
    if (sr<=0.0) { fprintf(stderr,"error in rms_scale_d: sr (sample rate) must be positive\n"); return 1; }
    if (tau<=0.0) { fprintf(stderr,"error in rms_scale_d: tau (time constant) must be positive\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t L = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;

    if (N<2u || L<2u) {}
    else
    {
        const double rms_prctile = 0.9;
        const double target_RMS = 20e-6 * pow(10.0, target_dB_SPL/20.0);
        const double max_RMS = 20e-6 * pow(10.0, max_dB_SPL/20.0);
        const double a1 = exp(-1.0/(sr*tau));
        const double b0 = 1.0 - a1;
        const size_t tranz = (L>5u) ? 5u : 0u;  //nsamps to zero in case transient
        const size_t rms_p = (size_t)round((1.0-rms_prctile)*(double)L);  //index for rms_prctile
        double rms, scale, scalemx;

        double *Y;
        if (!(Y=(double *)malloc(L*sizeof(double)))) { fprintf(stderr,"error in rms_scale_d: problem with malloc. "); perror("malloc"); return 1; }

        if (L==N)
        {
            *Y = b0 * *X * *X;
            ++X; ++Y;
            for (size_t l=L-1u; l>0u; --l, ++X, ++Y) { *Y = b0**X**X + a1**(Y-1); }
            Y -= L;
            for (size_t l=tranz; l>0u; --l, ++Y) { *Y = 0.0; }  //in case transient
            Y -= tranz;
            if (LAPACKE_dlasrt_work('D',(int)L,Y)) { fprintf(stderr,"error in rms_scale_d: problem with LAPACKE function\n"); }
            rms = sqrt(*Y);
            scalemx = (rms>DBL_EPSILON) ? max_RMS/rms : 1.0;
            rms = sqrt(*(Y+rms_p));
            scale = (rms>DBL_EPSILON) ? target_RMS/rms : 1.0;
            scale = (scale<scalemx) ? scale : scalemx;
            for (size_t l=L; l>0u; --l) { *--X *= scale; }
        }
        else
        {
            const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
            const size_t B = (iscolmajor && dim==0u) ? C*S*H : K;
            const size_t V = N/L, G = V/B;

            if (K==1u && (G==1u || B==1u))
            {
                for (size_t v=V; v>0u; --v, X+=L)
                {
                    *Y = b0 * *X * *X;
                    ++X; ++Y;
                    for (size_t l=L-1u; l>0u; --l, ++X, ++Y) { *Y = b0**X**X + a1**(Y-1); }
                    Y -= L;
                    for (size_t l=tranz; l>0u; --l, ++Y) { *Y = 0.0; }  //in case transient
                    Y -= tranz;
                    if (LAPACKE_dlasrt_work('D',(int)L,Y)) { fprintf(stderr,"error in rms_scale_d: problem with LAPACKE function\n"); }
                    rms = sqrt(*Y);
                    scalemx = (rms>DBL_EPSILON) ? max_RMS/rms : 1.0;
                    rms = sqrt(*(Y+rms_p));
                    scale = (rms>DBL_EPSILON) ? target_RMS/rms : 1.0;
                    scale = (scale<scalemx) ? scale : scalemx;
                    for (size_t l=L; l>0u; --l) { *--X *= scale; }
                }
            }
            else
            {
                for (size_t g=G; g>0u; --g, X+=B*(L-1u))
                {
                    for (size_t b=B; b>0u; --b, ++X)
                    {
                        *Y = b0 * *X * *X;
                        X+=K; ++Y;
                        for (size_t l=L-1u; l>0u; --l, X+=K, ++Y) { *Y = b0**X**X + a1**(Y-1); }
                        Y -= L;
                        for (size_t l=tranz; l>0u; --l, ++Y) { *Y = 0.0; }  //in case transient
                        Y -= tranz;
                        if (LAPACKE_dlasrt_work('D',(int)L,Y)) { fprintf(stderr,"error in rms_scale_d: problem with LAPACKE function\n"); }
                        rms = sqrt(*Y);
                        scalemx = (rms>DBL_EPSILON) ? max_RMS/rms : 1.0;
                        rms = sqrt(*(Y+rms_p));
                        scale = (rms>DBL_EPSILON) ? target_RMS/rms : 1.0;
                        scale = (scale<scalemx) ? scale : scalemx;
                        for (size_t l=L; l>0u; --l) { X-=K; *X *= scale; }
                    }
                }
            }
        }

        free(Y);
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif

ext: c


text:
// THE VIRTUAL BRAIN -- C
// A fast implementation of TVB-style brain network models based on
// DYNAMIC MEAN FIELD MODEL Deco et al. 2014 Journal of Neuroscience
//
//  m.schirner@fu-berlin.de
//  michael.schirner@charite.de
//
// MIT LICENSE
// Copyright 2020 Michael Schirner
//
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

#include <stdio.h>
#include <xmmintrin.h>
#include <emmintrin.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <math.h>
#include <pthread.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
//#include <semaphore.h>

/*

 Definition of pthread barrier -- since OSX doesn't seem to have it

 */

 typedef int pthread_barrierattr_t;
 typedef struct
 {
 pthread_mutex_t mutex;
 pthread_cond_t cond;
 int count;
 int tripCount;
 } pthread_barrier_t;

 int pthread_barrier_init(pthread_barrier_t *barrier, const pthread_barrierattr_t *attr, unsigned int count)
 {
 if(count == 0)
 {
 //errno = EINVAL;
 return -1;
 }
 if(pthread_mutex_init(&barrier->mutex, 0) < 0)
 {
 return -1;
 }
 if(pthread_cond_init(&barrier->cond, 0) < 0)
 {
 pthread_mutex_destroy(&barrier->mutex);
 return -1;
 }
 barrier->tripCount = count;
 barrier->count = 0;

 return 0;
 }

 int pthread_barrier_destroy(pthread_barrier_t *barrier)
 {
 pthread_cond_destroy(&barrier->cond);
 pthread_mutex_destroy(&barrier->mutex);
 return 0;
 }

 int pthread_barrier_wait(pthread_barrier_t *barrier)
 {
 pthread_mutex_lock(&barrier->mutex);
 ++(barrier->count);
 if(barrier->count >= barrier->tripCount)
 {
 barrier->count = 0;
 pthread_cond_broadcast(&barrier->cond);
 pthread_mutex_unlock(&barrier->mutex);
 return 1;
 }
 else
 {
 pthread_cond_wait(&barrier->cond, &(barrier->mutex));
 pthread_mutex_unlock(&barrier->mutex);
 return 0;
 }
 }

/*

 End: Definition of pthread barrier

 */

pthread_mutex_t mutex_thrcount;

struct Xi_p{
    float **Xi_elems;
};

struct SC_capS{
    float *cap;
};

struct SC_inpregS{
    int *inpreg;
};

pthread_barrier_t mybarrier_base, mybarrier1, mybarrier2, mybarrier3;

#define REAL float
//#define REAL double

int divRoundClosest(const int n, const int d)
{
    return ((n < 0) ^ (d < 0)) ? ((n - d/2)/d) : ((n + d/2)/d);
}

int divRoundUp(const int x, const int y)
{
    return (x + y - 1) / y;
}

/* Compute Pearson's correlation coefficient */
float corr(float *x, float *y, int n){
    int i;
    float mx=0, my=0;

    /* Calculate the mean of the two series x[], y[] */
    for (i=0; i<n; i++) {
        mx += x[i];
        my += y[i];
    }
    mx /= n;
    my /= n;

    /* Calculate the correlation */
    float sxy = 0, sxsq = 0, sysq = 0, tmpx, tmpy;
    for (i=0; i<n; i++) {
        tmpx = x[i] - mx;
        tmpy = y[i] - my;
        sxy += tmpx*tmpy;
        sxsq += tmpx*tmpx;
        sysq += tmpy*tmpy;
    }

    return (sxy / (sqrt(sxsq)*sqrt(sysq)));
}

//float gaussrand_ret()
//{
//    static double V1=0.0, V2=0.0, S=0.0, U1=0.0, U2=0.0;
//    S=0.0;
//    do {
//        U1 = (double)rand() / RAND_MAX;
//        U2 = (double)rand() / RAND_MAX;
//        V1 = 2 * U1 - 1;
//        V2 = 2 * U2 - 1;
//        S = V1 * V1 + V2 * V2;
//    } while(S >= 1 || S == 0);
//
//    return (float)(V1 * sqrt(-2 * log(S) / S));
//}
//
//static inline void gaussrand(float *randnum)
//{
//    static double V1=0.0, V2=0.0, S=0.0, U1=0.0, U2=0.0;
//
//    S=0.0;
//    do {
//        U1 = (double)rand() / RAND_MAX;
//        U2 = (double)rand() / RAND_MAX;
//        V1 = 2 * U1 - 1;
//        V2 = 2 * U2 - 1;
//        S = V1 * V1 + V2 * V2;
//    } while(S >= 1 || S == 0);
//
//    randnum[0] = (float)(V1 * sqrt(-2 * log(S) / S));
//    randnum[1] = (float)(V2 * sqrt(-2 * log(S) / S));
//
//    S=0.0;
//    do {
//        U1 = (double)rand() / RAND_MAX;
//        U2 = (double)rand() / RAND_MAX;
//        V1 = 2 * U1 - 1;
//        V2 = 2 * U2 - 1;
//        S = V1 * V1 + V2 * V2;
//    } while(S >= 1 || S == 0);
//
//    randnum[2] = (float)(V1 * sqrt(-2 * log(S) / S));
//    randnum[3] = (float)(V2 * sqrt(-2 * log(S) / S));
//}

int importGlobalConnectivity(char *SC_cap_filename, char *SC_dist_filename, char *SC_inputreg_filename, int regions, float **region_activity, struct Xi_p **reg_globinp_p, float global_trans_v, int **n_conn_table, float G_J_NMDA, struct SC_capS **SC_cap, float **SC_rowsums, struct SC_inpregS **SC_inpreg)
{

    int i,j,k, maxdelay=0, tmpint;
    float *region_activity_p;
    double tmp, tmp2;
    struct Xi_p *reg_globinp_pp;
    struct SC_capS      *SC_capp;
    struct SC_inpregS   *SC_inpregp;

    // Open SC files
    FILE *file_cap, *file_dist;
    file_cap=fopen(SC_cap_filename, "r");
    file_dist=fopen(SC_dist_filename, "r");
    if (file_cap==NULL || file_dist==NULL)
    {
        printf( "\nERROR: Could not open SC files. Terminating... \n\n");
        exit(0);
    }

    // Allocate a counter that stores number of region input connections for each region and the SCcap array
    *SC_rowsums = (float *)_mm_malloc(regions*sizeof(float),16);
    *n_conn_table = (int *)_mm_malloc(regions*sizeof(int),16);
    *SC_cap = (struct SC_capS *)_mm_malloc(regions*sizeof(struct SC_capS),16);
    SC_capp = *SC_cap;
    *SC_inpreg = (struct SC_inpregS *)_mm_malloc(regions*sizeof(struct SC_inpregS),16);
    SC_inpregp = *SC_inpreg;
    if(*n_conn_table==NULL || SC_capp==NULL || SC_rowsums==NULL || SC_inpregp==NULL){
        printf("Running out of memory. Terminating.\n");fclose(file_dist);fclose(file_cap);exit(2);
    }

    // Read out the maximal fiber length and the degree of each node and rewind SC files
    int n_entries = 0, curr_col = 0, curr_row = 0, curr_row_nonzero = 0;
    double tmp_max = -9999;
    while(fscanf(file_dist,"%lf",&tmp) != EOF && fscanf(file_cap,"%lf",&tmp2) != EOF){
        if (tmp_max < tmp) tmp_max = tmp;
        n_entries++;

        if (tmp2 > 0.0) curr_row_nonzero++;

        curr_col++;
        if (curr_col == regions) {
            curr_col = 0;
            (*n_conn_table)[curr_row] = curr_row_nonzero;
            curr_row_nonzero = 0;
            curr_row++;
        }
    }

    if (n_entries < regions * regions) {
        printf( "ERROR: Unexpected end-of-file in file. File contains less input than expected. Terminating... \n");
        exit(0);
    }

    rewind(file_dist);
    rewind(file_cap);

    maxdelay = (int)(((tmp_max/global_trans_v)*10)+0.5); // *10 for getting from m/s to 10kHz sampling, +0.5 for rounding by casting
    if (maxdelay < 1) maxdelay = 1; // Case: no time delays

    // Allocate ringbuffer that contains region activity for each past time-step until maxdelay and another ringbuffer that contains pointers to the first ringbuffer
    *region_activity = (float *)_mm_malloc(maxdelay*regions*sizeof(float),16);
    region_activity_p = *region_activity;
    *reg_globinp_p = (struct Xi_p *)_mm_malloc(maxdelay*regions*sizeof(struct Xi_p),16);
    reg_globinp_pp = *reg_globinp_p;
    if(region_activity_p==NULL || reg_globinp_p==NULL){
        printf("Running out of memory. Terminating.\n");fclose(file_dist);exit(2);
    }
    for (j=0; j<maxdelay*regions; j++) {
        region_activity_p[j]=0.001;
    }

    // Read SC files and set pointers for each input region and correspoding delay for each ringbuffer time-step
    int ring_buff_position;
    for (i=0; i<regions; i++) {

        if ((*n_conn_table)[i] > 0) {
            // SC strength and inp region numbers
            SC_capp[i].cap          = (float *)_mm_malloc(((*n_conn_table)[i])*sizeof(float),16);
            SC_inpregp[i].inpreg    = (int *)_mm_malloc(((*n_conn_table)[i])*sizeof(int),16);
            if(SC_capp[i].cap==NULL || SC_inpregp[i].inpreg==NULL){
                printf("Running out of memory. Terminating.\n");exit(2);
            }

            // Allocate memory for input-region-pointer arrays for each time-step in ringbuffer
            for (j=0; j<maxdelay; j++){
                reg_globinp_pp[i+j*regions].Xi_elems=(float **)_mm_malloc(((*n_conn_table)[i])*sizeof(float *),16);
                if(reg_globinp_pp[i+j*regions].Xi_elems==NULL){
                    printf("Running out of memory. Terminating.\n");exit(2);
                }
            }

            float sum_caps=0.0;
            // Read incoming connections and set pointers
            curr_row_nonzero = 0;
            for (j=0; j<regions; j++) {
                if(fscanf(file_cap,"%lf",&tmp) != EOF && fscanf(file_dist,"%lf",&tmp2) != EOF){
                    if (tmp > 0.0) {
                        tmpint = (int)(((tmp2/global_trans_v)*10)+0.5); //  *10 for getting from m/s or mm/ms to 10kHz sampling, +0.5 for rounding by casting
                        if (tmpint < 0 || tmpint > maxdelay){
                            printf( "\nERROR: Negative or too high (larger than maximum specified number) connection length/delay %d -> %d. Terminating... \n\n",i,j);exit(0);
                        }
                        if (tmpint <= 0) tmpint = 1; // If time delay is smaller than integration step size, than set time delay to one integration step

                        SC_capp[i].cap[curr_row_nonzero] = (float)tmp * G_J_NMDA;
                        //sum_caps                += (float)tmp;
                        sum_caps                += SC_capp[i].cap[curr_row_nonzero];
                        SC_inpregp[i].inpreg[curr_row_nonzero]  =  j;

                        ring_buff_position=maxdelay*regions - tmpint*regions + j;
                        for (k=0; k<maxdelay; k++) {
                            reg_globinp_pp[i+k*regions].Xi_elems[curr_row_nonzero]=&region_activity_p[ring_buff_position];
                            ring_buff_position += regions;
                            if (ring_buff_position > (maxdelay*regions-1)) ring_buff_position -= maxdelay*regions;
                        }

                        curr_row_nonzero++;
                    }
                } else{
                    printf( "\nERROR: Unexpected end-of-file in file %s or %s. File contains less input than expected. Terminating... \n\n", SC_cap_filename, SC_dist_filename);
                    exit(0);
                }

            }
            if (sum_caps <= 0) {
                printf( "\nERROR: Sum of connection strenghts is negative or zero. sum-caps node %d = %f. Terminating... \n\n",i,sum_caps);exit(0);
            }
            (*SC_rowsums)[i] = sum_caps;
        }
    }

    fclose(file_dist);fclose(file_cap);
    return maxdelay;
}

/*
 Initialize thread barriers
 */
void initialize_thread_barriers(int n_threads) {
    pthread_barrier_init(&mybarrier1, NULL, n_threads);
    pthread_barrier_init(&mybarrier2, NULL, n_threads);
    pthread_barrier_init(&mybarrier3, NULL, n_threads);
    return;
}

/* create thread argument struct for thr_func() */
typedef struct _thread_data_t {
    int                 tid, rand_num_seed;
    int                 nodes, nodes_vec, fake_nodes, n_threads, vectorization_grade, time_steps, BOLD_TR, BOLD_ts_len;
    float               *J_i;
    int                 reg_act_size;
    float               *region_activity;
    float               model_dt;
    __m128              _gamma;
    __m128              _one;
    __m128              _imintau_E;
    __m128              _dt;
    __m128              _sigma_sqrt_dt;
    __m128              _sigma;
    __m128              _gamma_I;
    __m128              _imintau_I;
    __m128              _min_d_I;
    __m128              _b_I;
    __m128              _J_NMDA;
    __m128              _w_I__I_0;
    __m128              _a_I;
    __m128              _min_d_E;
    __m128              _b_E;
    __m128              _a_E;
    __m128              _w_plus_J_NMDA;
    __m128              _w_E__I_0;
    struct SC_capS      *SC_cap;
    struct SC_inpregS   *SC_inpreg;
    struct Xi_p         *reg_globinp_p;
    int                 *n_conn_table;
    float               *BOLD_ex;
    char                *output_file;
} thread_data_t;

/*

 Function for multi-threaded FIC tuning

 */
void *run_simulation(void *arg)
{

    /* Local function parameters (no concurrency) */
    int j, i_node_vec, i_node_vec_local, k, int_i, ts;
    float tmpglobinput, tmpglobinput_FFI;
    __m128 _tmp_H_E, _tmp_H_I, _tmp_I_I, _tmp_I_E;
    float tmp_exp_E[4]          __attribute__((aligned(16)));
    float tmp_exp_I[4]          __attribute__((aligned(16)));
    __m128          *_tmp_exp_E         = (__m128*)tmp_exp_E;
    __m128          *_tmp_exp_I         = (__m128*)tmp_exp_I;
    float rand_number[4]        __attribute__((aligned(16)));
    __m128          *_rand_number       = (__m128*)rand_number;
    int     ring_buf_pos                = 0;

    /* Global function parameters (some with concurrency) */
    thread_data_t *thr_data = (thread_data_t *)arg;
    //printf("Hello from thread %d, process %d\n", thr_data->tid, thr_data->rank);
    int     t_id                    = thr_data->tid;
    int     rand_num_seed           = thr_data->rand_num_seed;
    int     nodes                   = thr_data->nodes;
    int     fake_nodes              = thr_data->fake_nodes;
    const   int nodes_vec           = thr_data->nodes_vec;
    int     n_threads               = thr_data->n_threads;
    int     vectorization_grade     = thr_data->vectorization_grade;
    int     reg_act_size            = thr_data->reg_act_size;
    int     *n_conn_table           = thr_data->n_conn_table;
    float   *J_i                    = thr_data->J_i;
    float   *region_activity        = thr_data->region_activity;
    float   *BOLD_ex                = thr_data->BOLD_ex;
    const   __m128 _gamma           = thr_data->_gamma;
    const   __m128 _one             = thr_data->_one;
    const   __m128 _imintau_E       = thr_data->_imintau_E;
    const   __m128 _dt              = thr_data->_dt;
    const   __m128 _sigma_sqrt_dt   = thr_data->_sigma_sqrt_dt;
    const   __m128 _gamma_I         = thr_data->_gamma_I;
    const   __m128 _imintau_I       = thr_data->_imintau_I;
    const   __m128 _min_d_I         = thr_data->_min_d_I;
    const   __m128 _b_I             = thr_data->_b_I;
    const   __m128 _J_NMDA          = thr_data->_J_NMDA;
    const   __m128 _w_I__I_0        = thr_data->_w_I__I_0;
    const   __m128 _a_I             = thr_data->_a_I;
    const   __m128 _min_d_E         = thr_data->_min_d_E;
    const   __m128 _b_E             = thr_data->_b_E;
    const   __m128 _a_E             = thr_data->_a_E;
    const   __m128 _w_plus_J_NMDA   = thr_data->_w_plus_J_NMDA;
    const   __m128 _w_E__I_0        = thr_data->_w_E__I_0;
    struct  SC_capS *SC_cap         = thr_data->SC_cap;
    struct  Xi_p *reg_globinp_p     = thr_data->reg_globinp_p;
    int     time_steps              = thr_data->time_steps;
    int     BOLD_TR                 = thr_data->BOLD_TR;
    float   model_dt                = thr_data->model_dt;
    int     BOLD_ts_len             = thr_data->BOLD_ts_len;

    /* Parallelism: divide work among threads */
    int nodes_vec_mt        = divRoundUp(nodes_vec, n_threads);
    int nodes_mt            = nodes_vec_mt  * vectorization_grade;

    int start_nodes_vec_mt  = t_id         * nodes_vec_mt;
    int start_nodes_mt      = t_id         * nodes_mt;
    int end_nodes_vec_mt    = (t_id + 1)   * nodes_vec_mt;
    int end_nodes_mt        = (t_id + 1)   * nodes_mt;
    int end_nodes_mt_glob   = end_nodes_mt; // end_nodes_mt_glob may differ from end_nodes_mt in the last thread (fake-nodes vs nodes)

    /* Correct settings for last thread */
    if (end_nodes_mt > nodes) {
        end_nodes_vec_mt    = nodes_vec;
        end_nodes_mt        = fake_nodes;
        end_nodes_mt_glob   = nodes;
        nodes_mt            = end_nodes_mt - start_nodes_mt;
        nodes_vec_mt        = end_nodes_vec_mt - start_nodes_vec_mt;
    }

    printf("thread %d: start: %d end: %d size: %d\n", thr_data->tid, start_nodes_mt, end_nodes_mt, nodes_mt);

    /* Check whether splitting makes sense (i.e. each thread has something to do). Terminate if not. */
    if (nodes_vec_mt <= 1) {
        printf("The splitting is ineffective (e.g. fewer nodes than threads or some threads have nothing to do). Use different numbers of threads. Terminating.\n");
        exit(0);
    }

    /* Set up barriers */
    if (t_id == 0) {
        initialize_thread_barriers(n_threads);
    }
    pthread_barrier_wait(&mybarrier_base);

    /* Initialize random number generator */
    //const gsl_rng_type * T;
    gsl_rng *r = gsl_rng_alloc (gsl_rng_mt19937);
    gsl_rng_set (r, (t_id+rand_num_seed)); // Random number seed -- otherwise every node gets the same random nums
    srand((unsigned)(t_id+rand_num_seed));

    // All threads
    float   *meanFR                 = (float *)_mm_malloc(nodes_mt * sizeof(float),16);  // summation array for mean firing rate over all threads
    float   *meanFR_INH             = (float *)_mm_malloc(nodes_mt * sizeof(float),16);  // summation array for mean firing rate over all threads
    float   *global_input           = (float *)_mm_malloc(nodes_mt * sizeof(float),16);
    float   *global_input_FFI       = (float *)_mm_malloc(nodes_mt * sizeof(float),16);
    float   *S_i_E                  = (float *)_mm_malloc(nodes_mt * sizeof(float),16);
    float   *S_i_I                  = (float *)_mm_malloc(nodes_mt * sizeof(float),16);
    //float   *r_i_E                    = (float *)_mm_malloc(nodes_mt * sizeof(float),16);
    //float   *r_i_I                    = (float *)_mm_malloc(nodes_mt * sizeof(float),16);
    float   *J_i_local              = (float *)_mm_malloc(nodes_mt * sizeof(float),16);

    if (meanFR == NULL || meanFR_INH == NULL || global_input == NULL || global_input_FFI == NULL || S_i_E==NULL || S_i_I==NULL || J_i_local==NULL) {
        printf( "ERROR: Running out of memory. Aborting... \n");
        exit(0);
    }

    // Cast to vector/SSE arrays
    __m128          *_meanFR         = (__m128*)meanFR;
    __m128          *_meanFR_INH     = (__m128*)meanFR_INH;
    __m128          *_global_input   = (__m128*)global_input;
    __m128          *_global_input_FFI = (__m128*)global_input_FFI;
    __m128          *_S_i_E          = (__m128*)S_i_E;
    __m128          *_S_i_I          = (__m128*)S_i_I;
    //__m128          *_r_i_E             = (__m128*)r_i_E;
    //__m128          *_r_i_I             = (__m128*)r_i_I;
    __m128          *_J_i_local      = (__m128*)J_i_local;

    //Balloon-Windkessel model parameters / arrays
    float rho = 0.34, alpha = 0.32, tau = 0.98, y = 1.0/0.41, kappa = 1.0/0.65;
    float V_0 = 0.02, k1 = 7 * rho, k2 = 2.0, k3 = 2 * rho - 0.2, ialpha = 1.0/alpha, itau = 1.0/tau, oneminrho = (1.0 - rho);
    float f_tmp;
    float *BOLD       = (float *)_mm_malloc(nodes_mt * BOLD_ts_len * sizeof(float),16);  // resulting BOLD data
    float *bw_x_ex    = (float *)_mm_malloc(nodes_mt * sizeof(float),16);  // State-variable 1 of BW-model (exc. pop.)
    float *bw_f_ex    = (float *)_mm_malloc(nodes_mt * sizeof(float),16);  // State-variable 2 of BW-model (exc. pop.)
    float *bw_nu_ex   = (float *)_mm_malloc(nodes_mt * sizeof(float),16);  // State-variable 3 of BW-model (exc. pop.)
    float *bw_q_ex    = (float *)_mm_malloc(nodes_mt * sizeof(float),16);  // State-variable 4 of BW-model (exc. pop.)
    //float *bw_x_in    = (float *)_mm_malloc(nodes_mt * sizeof(float),16);  // State-variable 1 of BW-model (inh. pop.)
    //float *bw_f_in    = (float *)_mm_malloc(nodes_mt * sizeof(float),16);  // State-variable 2 of BW-model (inh. pop.)
    //float *bw_nu_in   = (float *)_mm_malloc(nodes_mt * sizeof(float),16);  // State-variable 3 of BW-model (inh. pop.)
    //float *bw_q_in    = (float *)_mm_malloc(nodes_mt * sizeof(float),16);  // State-variable 4 of BW-model (inh. pop.)
    //if (bw_x_ex == NULL || bw_f_ex == NULL || bw_nu_ex==NULL || bw_q_ex==NULL || bw_x_in==NULL || bw_f_in==NULL || bw_nu_in==NULL || bw_q_in==NULL) {}
    if (bw_x_ex == NULL || bw_f_ex == NULL || bw_nu_ex==NULL || bw_q_ex==NULL) {
        printf( "ERROR: Running out of memory. Aborting... \n");
        exit(0);
    }

    // Reset arrays and variables
    ring_buf_pos     = 0;
    for (j = 0; j < nodes_mt; j++) {
        meanFR[j]           = 0.0;
        meanFR_INH[j]       = 0.0;
        global_input[j]     = 0.00;
        global_input_FFI[j] = 0.00;
        S_i_E[j]            = 0.00;
        S_i_I[j]            = 0.00;
        J_i_local[j]        = J_i[j+start_nodes_mt];
    }
    if (t_id == 0) {
        for (j=0; j<reg_act_size; j++) {
            region_activity[j]=0.00;
        }
    }

    // Reset Balloon-Windkessel model parameters and arrays
    for (j = 0; j < nodes_mt; j++) {
        bw_x_ex[j] = 0.0;
        bw_f_ex[j] = 1.0;
        bw_nu_ex[j] = 1.0;
        bw_q_ex[j] = 1.0;
        //bw_x_in[j] = 0.0;
        //bw_f_in[j] = 1.0;
        //bw_nu_in[j] = 1.0;
        //bw_q_in[j] = 1.0;
    }

    /* Barrier: only start simulating when arrays were cleared */
    pthread_barrier_wait(&mybarrier1);

    /*  Simulation loop */
    int ts_bold_i = 0;
    for (ts = 0; ts < time_steps; ts++) {

        if (t_id == 0){
            printf("%.1f %% \r", ((float)ts / (float)time_steps) * 100.0f );
        }

        // dt = 0.1 ms => 10 integration iterations per ms
        for (int_i = 0; int_i < 10; int_i++) {
            /* Barrier: only start integrating input when all other threads are finished copying their results to buffer */
            pthread_barrier_wait(&mybarrier2);

            /* Compute global coupling */
            i_node_vec_local = 0;
            for(j=start_nodes_mt; j<end_nodes_mt_glob; j++){
                tmpglobinput     = 0;
                tmpglobinput_FFI = 0;
                for (k=0; k<n_conn_table[j]; k++) {
                    tmpglobinput     += *reg_globinp_p[j+ring_buf_pos].Xi_elems[k] * SC_cap[j].cap[k];
                    //tmpglobinput_FFI += *reg_globinp_p[j+ring_buf_pos].Xi_elems[k] * SC_cap_FFI[j].cap[k];
                }
                global_input[i_node_vec_local]     = tmpglobinput;
                //global_input_FFI[i_node_vec_local] = tmpglobinput_FFI;

                i_node_vec_local++;
            }

            i_node_vec_local = 0;
            for (i_node_vec = start_nodes_vec_mt; i_node_vec < end_nodes_vec_mt; i_node_vec++) {

                // Excitatory population firing rate
                _tmp_I_E    = _mm_sub_ps(_mm_mul_ps(_a_E,_mm_add_ps(_mm_add_ps(_w_E__I_0,_mm_mul_ps(_w_plus_J_NMDA, _S_i_E[i_node_vec_local])),_mm_sub_ps(_global_input[i_node_vec_local],_mm_mul_ps(_J_i_local[i_node_vec_local], _S_i_I[i_node_vec_local])))),_b_E);

                *_tmp_exp_E     = _mm_mul_ps(_min_d_E, _tmp_I_E);
                tmp_exp_E[0]    = tmp_exp_E[0] != 0 ? expf(tmp_exp_E[0]) : 0.9;
                tmp_exp_E[1]    = tmp_exp_E[1] != 0 ? expf(tmp_exp_E[1]) : 0.9;
                tmp_exp_E[2]    = tmp_exp_E[2] != 0 ? expf(tmp_exp_E[2]) : 0.9;
                tmp_exp_E[3]    = tmp_exp_E[3] != 0 ? expf(tmp_exp_E[3]) : 0.9;
                _tmp_H_E        = _mm_div_ps(_tmp_I_E, _mm_sub_ps(_one, *_tmp_exp_E));
                _meanFR[i_node_vec_local] = _mm_add_ps(_meanFR[i_node_vec_local],_tmp_H_E);
                //_r_i_E[i_node_vec_local]  = _tmp_H_E;

                // Inhibitory population firing rate
                _tmp_I_I = _mm_sub_ps(_mm_mul_ps(_a_I,_mm_sub_ps(_mm_add_ps(_mm_add_ps(_w_I__I_0,_global_input_FFI[i_node_vec_local]),_mm_mul_ps(_J_NMDA, _S_i_E[i_node_vec_local])), _S_i_I[i_node_vec_local])),_b_I);
                *_tmp_exp_I   = _mm_mul_ps(_min_d_I, _tmp_I_I);
                tmp_exp_I[0]  = tmp_exp_I[0] != 0 ? expf(tmp_exp_I[0]) : 0.9;
                tmp_exp_I[1]  = tmp_exp_I[1] != 0 ? expf(tmp_exp_I[1]) : 0.9;
                tmp_exp_I[2]  = tmp_exp_I[2] != 0 ? expf(tmp_exp_I[2]) : 0.9;
                tmp_exp_I[3]  = tmp_exp_I[3] != 0 ? expf(tmp_exp_I[3]) : 0.9;
                _tmp_H_I  = _mm_div_ps(_tmp_I_I, _mm_sub_ps(_one, *_tmp_exp_I));
                _meanFR_INH[i_node_vec_local] = _mm_add_ps(_meanFR_INH[i_node_vec_local],_tmp_H_I);
                //_r_i_I[i_node_vec_local] = _tmp_H_I;

                //gaussrand(rand_number);
                rand_number[0] = (float)gsl_ran_gaussian(r, 1.0);
                rand_number[1] = (float)gsl_ran_gaussian(r, 1.0);
                rand_number[2] = (float)gsl_ran_gaussian(r, 1.0);
                rand_number[3] = (float)gsl_ran_gaussian(r, 1.0);

                _S_i_I[i_node_vec_local] = _mm_add_ps(_mm_add_ps(_mm_mul_ps(_sigma_sqrt_dt, *_rand_number),_S_i_I[i_node_vec_local]),_mm_mul_ps(_dt,_mm_add_ps(_mm_mul_ps(_imintau_I, _S_i_I[i_node_vec_local]),_mm_mul_ps(_tmp_H_I,_gamma_I))));

                //gaussrand(rand_number);
                rand_number[0] = (float)gsl_ran_gaussian(r, 1.0);
                rand_number[1] = (float)gsl_ran_gaussian(r, 1.0);
                rand_number[2] = (float)gsl_ran_gaussian(r, 1.0);
                rand_number[3] = (float)gsl_ran_gaussian(r, 1.0);

                _S_i_E[i_node_vec_local] = _mm_add_ps(_mm_add_ps(_mm_mul_ps(_sigma_sqrt_dt, *_rand_number),_S_i_E[i_node_vec_local]),_mm_mul_ps(_dt, _mm_add_ps(_mm_mul_ps(_imintau_E, _S_i_E[i_node_vec_local]),_mm_mul_ps(_mm_mul_ps(_mm_sub_ps(_one, _S_i_E[i_node_vec_local]),_gamma),_tmp_H_E))));

                i_node_vec_local++;
            }

            /* ensure that synaptic activity is within its boundaries */
            for(j=0; j<nodes_mt; j++){
                S_i_E[j] = S_i_E[j] >= 0.0 ? S_i_E[j] : 0.0;
                S_i_E[j] = S_i_E[j] <= 1.0 ? S_i_E[j] : 1.0;
                S_i_I[j] = S_i_I[j] >= 0.0 ? S_i_I[j] : 0.0;
                S_i_I[j] = S_i_I[j] <= 1.0 ? S_i_I[j] : 1.0;
            }

            /* Copy result of this thread into global region_activity buffer */
            memcpy(&region_activity[ring_buf_pos+start_nodes_mt], S_i_E, nodes_mt*sizeof( float ));

            /* Shift region_activity ring-buffer start */
            ring_buf_pos = ring_buf_pos<(reg_act_size-nodes) ? (ring_buf_pos+nodes) : 0;

        }

        /*
         Compute BOLD for that time-step (subsampled to 1 ms)
         */
        for (j = 0; j < nodes_mt; j++) {
            bw_x_ex[j]  = bw_x_ex[j]  +  model_dt * (S_i_E[j] - kappa * bw_x_ex[j] - y * (bw_f_ex[j] - 1.0));
            f_tmp       = bw_f_ex[j]  +  model_dt * bw_x_ex[j];
            bw_nu_ex[j] = bw_nu_ex[j] +  model_dt * itau * (bw_f_ex[j] - powf(bw_nu_ex[j], ialpha));
            bw_q_ex[j]  = bw_q_ex[j]  +  model_dt * itau * (bw_f_ex[j] * (1.0 - powf(oneminrho,(1.0/bw_f_ex[j]))) / rho  - powf(bw_nu_ex[j],ialpha) * bw_q_ex[j] / bw_nu_ex[j]);
            bw_f_ex[j]  = f_tmp;

            /*
             bw_x_in[j]  = bw_x_in[j]  +  model_dt * (S_i_I[j] - kappa * bw_x_in[j] - y * (bw_f_in[j] - 1.0));
             f_tmp       = bw_f_in[j]  +  model_dt * bw_x_in[j];
             bw_nu_in[j] = bw_nu_in[j] +  model_dt * itau * (bw_f_in[j] - powf(bw_nu_in[j], ialpha));
             bw_q_in[j]  = bw_q_in[j]  +  model_dt * itau * (bw_f_in[j] * (1.0 - powf(oneminrho,(1.0/bw_f_in[j]))) / rho  - powf(bw_nu_in[j],ialpha) * bw_q_in[j] / bw_nu_in[j]);
             bw_f_in[j]  = f_tmp;
             */
        }

        if (ts % BOLD_TR == 0) {
            for (j = 0; j < nodes_mt; j++) {
                BOLD[ts_bold_i + j * BOLD_ts_len] = 100 / rho * V_0 * (k1 * (1 - bw_q_ex[j]) + k2 * (1 - bw_q_ex[j]/bw_nu_ex[j]) + k3 * (1 - bw_nu_ex[j]));
            }
            ts_bold_i++;
        }
    }

    /* Copy results of this thread into shared memory buffer */
    memcpy(&BOLD_ex[start_nodes_mt*BOLD_ts_len], BOLD, nodes_mt * BOLD_ts_len * sizeof( float ));

    /* Wait until other threads finished before writing out result */
    pthread_barrier_wait(&mybarrier3);

    /* Write out fMRI time series */
    if (t_id == 0){
        FILE *FCout = fopen(thr_data->output_file, "w");
        for (j = 0; j < nodes; j++) {
            for (k = 0; k < ts_bold_i; k++) {
                fprintf(FCout, "%.5f ",BOLD_ex[j*BOLD_ts_len + k]);
            }
            fprintf(FCout, "\n");
        }
        fclose(FCout);
    }

    gsl_rng_free (r);
    pthread_exit(NULL);
}

/*
 Usage: tvbii <paramfile> <subject_id>
 */

int main(int argc, char *argv[])
{
    /*
     Check whether right number of arguments was passed to program call
     */
    if (argc != 4 || atoi(argv[3]) <= 0) {
        printf("\nERROR: Invalid arguments.\n\nUsage: tvbii <paramset_file> <sub_id> <#threads>\n\nTerminating... \n\n");
        printf("\nProvided arguments:\n");
        int i;
        for (i=0; i<argc; i++) {
            printf("%s\n", argv[i]);
        }
        exit(0);
    }

    /*
     Get current time and do some initializations
     */
    time_t  start = time(NULL);
    int     i, j;
    int     n_threads = atoi(argv[3]);

    /*
     Global model and integration parameters
     */
    const float dt                  = 0.1;      // Integration step length dt = 0.1 ms
    const float sqrt_dt             = sqrtf(dt);// Noise in the diffusion part scaled by sqrt of dt
    const float model_dt            = 0.001;    // Time-step of model (sampling-rate=1000 Hz)
    const int   vectorization_grade = 4;        // How many operations can be done simultaneously. Depends on CPU Architecture and available intrinsics.
    int         time_steps          = 667*1.94*1000;    // Simulation length
    int         nodes               = 84;    // Number of surface vertices; must be a multiple of vectorization grade
    int         fake_nodes          = 84;
    float       global_trans_v      = 1.0;     // Global transmission velocity (m/s); Local time-delays can be ommited since smaller than integration time-step
    float       G                   = 0.5;        // Global coupling strength
    int         BOLD_TR             = 1940;     // TR of BOLD data
    int         rand_num_seed       = 1403;

    /*
     Local model: DMF-Parameters from Deco et al. JNeuro 2014
     */
    float w_plus  = 1.4;          // local excitatory recurrence
    //float I_ext   = 0;            // External stimulation
    float J_NMDA  = 0.15;         // (nA) excitatory synaptic coupling
    //float J_i     = 1.0;          // 1 for no-FIC, !=1 for Feedback Inhibition Control
    const float a_E     = 310;          // (n/C)
    const float b_E     = 125;          // (Hz)
    const float d_E     = 0.16;         // (s)
    const float a_I     = 615;          // (n/C)
    const float b_I     = 177;          // (Hz)
    const float d_I     = 0.087;        // (s)
    const float gamma   = 0.641/1000.0; // factor 1000 for expressing everything in ms
    const float tau_E   = 100;          // (ms) Time constant of NMDA (excitatory)
    const float tau_I   = 10;           // (ms) Time constant of GABA (inhibitory)
    float       sigma   = 0.00316228;   // (nA) Noise amplitude
    const float I_0     = 0.382;        // (nA) overall effective external input
    const float w_E     = 1.0;          // scaling of external input for excitatory pool
    const float w_I     = 0.7;          // scaling of external input for inhibitory pool
    const float gamma_I = 1.0/1000.0;   // for expressing inhib. pop. in ms
    float       tmpJi   = 1.0;          // Feedback inhibition J_i

    /* Read parameters from input file. Input file is a simple text file that contains one line with parameters and white spaces in between. */
    FILE *file;
    char param_file[300];memset(param_file, 0, 300*sizeof(char));
    snprintf(param_file, sizeof(param_file), "/input/%s",argv[1]);
    file=fopen(param_file, "r");
    if (file==NULL){
        printf( "\nERROR: Could not open file %s. Terminating... \n\n", param_file);
        exit(0);
    }
    if(fscanf(file,"%d",&nodes) != EOF && fscanf(file,"%f",&G) != EOF && fscanf(file,"%f",&J_NMDA) != EOF && fscanf(file,"%f",&w_plus) != EOF && fscanf(file,"%f",&tmpJi) != EOF && fscanf(file,"%f",&sigma) != EOF && fscanf(file,"%d",&time_steps) != EOF && fscanf(file,"%d",&BOLD_TR) != EOF && fscanf(file,"%f",&global_trans_v) != EOF && fscanf(file,"%d",&rand_num_seed) != EOF){
    } else{
        printf( "\nERROR: Unexpected end-of-file in file %s. File contains less input than expected. Terminating... \n\n", param_file);
        exit(0);
    }
    fclose(file);

    char output_file[300];memset(output_file, 0, 300*sizeof(char));
    snprintf(output_file, sizeof(output_file), "/output/%s_%s_fMRI.txt",argv[2],argv[1]);

    /* Add fake regions to make region count a multiple of vectorization grade */
    if (nodes % vectorization_grade != 0){
        printf( "\nWarning: Specified number of nodes (%d) is not a multiple of vectorization degree (%d). Will add some fake nodes... \n\n", nodes, vectorization_grade);

        int remainder   = nodes%vectorization_grade;
        if (remainder > 0) {
            fake_nodes = nodes + (vectorization_grade - remainder);
        }
    } else{
        fake_nodes = nodes;
    }

    /* Initialize random number generator */
    srand((unsigned)rand_num_seed);

    /* Derived parameters */
    const float sigma_sqrt_dt = sqrt_dt * sigma;
    const int   nodes_vec     = fake_nodes/vectorization_grade;
    const float min_d_E       = -1.0 * d_E;
    const float min_d_I       = -1.0 * d_I;
    const float imintau_E     = -1.0 / tau_E;
    const float imintau_I     = -1.0 / tau_I;
    const float w_E__I_0      = w_E * I_0;
    const float w_I__I_0      = w_I * I_0;
    const float one           = 1.0;
    const float w_plus__J_NMDA= w_plus * J_NMDA;
    const float G_J_NMDA      = G * J_NMDA;
          float TR            = (float)BOLD_TR / 1000;   // (s) TR of fMRI data
          int   BOLD_ts_len   = time_steps / (TR / model_dt) + 1;  // Length of BOLD time-series written to HDD

    /*
     Import and setup global and local connectivity
     */
    int         *n_conn_table;
    float       *region_activity, *SC_rowsums;
    struct Xi_p *reg_globinp_p;
    struct SC_capS      *SC_cap;
    struct SC_inpregS   *SC_inpreg;
    char cap_file[300];memset(cap_file, 0, 300*sizeof(char));
    snprintf(cap_file, sizeof(cap_file), "/input/%s_SC_weights.txt",argv[2]);
    char dist_file[300];memset(dist_file, 0, 300*sizeof(char));
    snprintf(dist_file, sizeof(dist_file), "/input/%s_SC_distances.txt",argv[2]);
    char reg_file[300];memset(reg_file, 0, 300*sizeof(char));
    snprintf(reg_file, sizeof(reg_file), "/input/%s_SC_regionids.txt",argv[2]);

    int         maxdelay = importGlobalConnectivity(cap_file, dist_file, reg_file, nodes, &region_activity, &reg_globinp_p, global_trans_v, &n_conn_table, G_J_NMDA, &SC_cap, &SC_rowsums, &SC_inpreg);

    int         reg_act_size = nodes * maxdelay;

    /*
     Initialize and/or cast to vector-intrinsics types for variables & parameters
     */
    float *J_i      = (float *)_mm_malloc(fake_nodes * sizeof(float),16); // (nA) inhibitory synaptic coupling
    float *BOLD_ex  = (float *)_mm_malloc(nodes * BOLD_ts_len * sizeof(float),16);

    if (J_i == NULL || BOLD_ex == NULL) {
        printf( "ERROR: Running out of memory. Aborting... \n");
        _mm_free(J_i);_mm_free(BOLD_ex);
        return 1;
    }
    // Initialize state variables / parameters
    for (j = 0; j < fake_nodes; j++) {
        J_i[j]              = tmpJi;
    }
    const __m128    _dt                 = _mm_load1_ps(&dt);
    const __m128    _sigma_sqrt_dt      = _mm_load1_ps(&sigma_sqrt_dt);
    const __m128    _w_plus_J_NMDA      = _mm_load1_ps(&w_plus__J_NMDA);
    const __m128    _a_E                = _mm_load1_ps(&a_E);
    const __m128    _b_E                = _mm_load1_ps(&b_E);
    const __m128    _min_d_E            = _mm_load1_ps(&min_d_E);
    const __m128    _a_I                = _mm_load1_ps(&a_I);
    const __m128    _b_I                = _mm_load1_ps(&b_I);
    const __m128    _min_d_I            = _mm_load1_ps(&min_d_I);
    const __m128    _gamma              = _mm_load1_ps(&gamma);
    const __m128    _gamma_I            = _mm_load1_ps(&gamma_I);
    const __m128    _imintau_E          = _mm_load1_ps(&imintau_E);
    const __m128    _imintau_I          = _mm_load1_ps(&imintau_I);
    const __m128    _w_E__I_0           = _mm_load1_ps(&w_E__I_0);
    const __m128    _w_I__I_0           = _mm_load1_ps(&w_I__I_0);
    float           tmp_sigma           = sigma*dt;// pre-compute dt*sigma for the integration of sigma*randnumber in equations (9) and (10) of Deco2014
    const __m128    _sigma              = _mm_load1_ps(&tmp_sigma);
    const __m128    _one                = _mm_load1_ps(&one);
    const __m128    _J_NMDA             = _mm_load1_ps(&J_NMDA);

    /*
     Start multithread simulation
     */

    /* Create threads */
    pthread_barrier_init(&mybarrier_base, NULL, n_threads);
    pthread_t       thr[n_threads];
    thread_data_t   thr_data[n_threads];
    int rc;
    for (i = 0; i < n_threads; ++i) {
        thr_data[i].tid             =	i               ;
        thr_data[i].n_threads       =	n_threads       ;
        thr_data[i].vectorization_grade = vectorization_grade;
        thr_data[i].nodes			=	nodes			;
        thr_data[i].fake_nodes      =   fake_nodes      ;
        thr_data[i].J_i             =	J_i            	;
        thr_data[i].reg_act_size	=	reg_act_size	;
        thr_data[i].region_activity	=	region_activity	;
        thr_data[i]._gamma			=	_gamma			;
        thr_data[i]._one			=	_one			;
        thr_data[i]._imintau_E		=	_imintau_E		;
        thr_data[i]._dt             =	_dt             ;
        thr_data[i]._sigma_sqrt_dt  =   _sigma_sqrt_dt  ;
        thr_data[i]._sigma			=	_sigma			;
        thr_data[i]._gamma_I		=	_gamma_I		;
        thr_data[i]._imintau_I		=	_imintau_I		;
        thr_data[i]._min_d_I		=	_min_d_I		;
        thr_data[i]._b_I			=	_b_I			;
        thr_data[i]._J_NMDA         =	_J_NMDA         ;
        thr_data[i]._w_I__I_0		=	_w_I__I_0		;
        thr_data[i]._a_I			=	_a_I			;
        thr_data[i]._min_d_E		=	_min_d_E		;
        thr_data[i]._b_E			=	_b_E			;
        thr_data[i]._a_E			=	_a_E			;
        thr_data[i]._w_plus_J_NMDA	=	_w_plus_J_NMDA	;
        thr_data[i]._w_E__I_0		=	_w_E__I_0		;
        thr_data[i].nodes_vec		=	nodes_vec		;
        thr_data[i].SC_cap			=	SC_cap			;
        thr_data[i].SC_inpreg       =	SC_inpreg       ;
        thr_data[i].reg_globinp_p	=	reg_globinp_p	;
        thr_data[i].n_conn_table	=	n_conn_table	;
        thr_data[i].time_steps      =   time_steps      ;
        thr_data[i].BOLD_TR         =   BOLD_TR         ;
        thr_data[i].model_dt        =   model_dt        ;
        thr_data[i].BOLD_ex         =   BOLD_ex         ;
        thr_data[i].rand_num_seed   =   rand_num_seed   ;
        thr_data[i].BOLD_ts_len     =   BOLD_ts_len     ;
        thr_data[i].output_file     =   output_file     ;

        if ((rc = pthread_create(&thr[i], NULL, run_simulation, &thr_data[i]))) {
            fprintf(stderr, "error: pthread_create, rc: %d\n", rc);
            return EXIT_FAILURE;
        }
    }

    /* block until all threads complete */
    for (i = 0; i < n_threads; ++i) {
        pthread_join(thr[i], NULL);
    }
    printf("Threads finished. Back to main thread.\n");

    printf("MT-TVBii finished. Execution took %.2f s for %d nodes. Goodbye!\n", (float)(time(NULL) - start), nodes);

    return 0;
}


text:
#include <stdio.h>
#include <gsl/gsl_ieee_utils.h>

int main() {
    //gsl_ieee_printf_double
    double a = 1.1;
    while (a > 0){
        gsl_ieee_printf_double(&a);
        printf("\r\n");
        a /= 2.0;
    }

    float b = 1.1;
    while (b > 0){
        gsl_ieee_printf_float(&b);
        printf("\r\n");
        b /= 2.0;
    }

    return 0;
}


text:
#ifndef alphas_H_INCLUDED
#define alphas_H_INCLUDED

//for fixed alpha
//#define g2 (4.6)

#include <gsl/gsl_spline.h>

//http://pdg.lbl.gov/2010/reviews/rpp2010-rev-qcd.pdf

#define Mz 91.1876
#define LamQCD 0.217
//fixednf = 0 for nf(mu)
#define fixednf 3
//extrapolates below mu0*LamQCD to value of (mu0-1.0)*LamQCD at mu=0
#define mu0 4.
#define pi 3.14159265359
#define Nc (3.)
#define Cf ((Nc*Nc-1.)/(2.*Nc))

void TabulateAlpha();
double alpha(double mu);

double nf(double mu);
double alpha_b0(double nf);
double alpha_b1(double nf);
double alpha_b2(double nf);
double alpha_s0(double mu);
double alpha_s1(double mu);
double alpha_s2(double mu);

gsl_interp_accel *accAlpha;
gsl_spline *Alpha;

#endif

ext: h

max_stars_repo_name: kdusling/mpc

max_issues_repo_name: kdusling/mpc

max_forks_repo_name: kdusling/mpc


text:
#if !defined (RAND_H)
#define RAND_H

#include <cstdlib>
#include <iostream>
#include <cmath>
#include <vector>
#include <list>
#include <set>
#include <limits>
#include <algorithm>

#include <gsl/gsl_sf_zeta.h>

using namespace std;

const double TOLERANCE = numeric_limits<double>::epsilon();
const int    INTMAX    = numeric_limits<int>::max();

const int BIGPRIME = 899999963;
const double PI    = 3.1415926535897932384626433832795028841971693993751;

typedef double (*FUNC)(const double, const double, const double);
inline double fx(const double htilde, const double y, const double u)
{
  return (pow(htilde,y+1) - (y+1)*htilde + fabs(2*u-1)*y);
}

inline int Getnode(list<int>& X, int index)
{
  int pos=0;
  for(list<int>::iterator p = X.begin(); p!= X.end(); p++, pos++) {
    if(pos==index)
      return (*p);
  }

  return -1;
}

class Rand {
  double c,cd,cm,u[97];
  int i97,j97 ;
  bool outputReady;
  double output;

  int idum;

 public:
  Rand(int seed) {
    ranmarin(fabs(seed%BIGPRIME));
    outputReady=false;
    idum=-1*seed;

    flag_setFvector = false;
    flag_setPvector = false;
  }

  Rand() {
    ranmarin(1);
    outputReady=false;
    idum=-1;

    flag_setFvector = false;
    flag_setPvector = false;
  }

  void seed(int seed) {
    ranmarin(fabs(seed%BIGPRIME));
    outputReady=false;
    idum=-1*seed;

    flag_setFvector = false;
    flag_setPvector = false;
  }

  double ran() { return rand()/(double)RAND_MAX;}

  double uniform()
  {
    double uni;

    uni=u[i97]-u[j97];
    if (uni<0.0) uni+=1.0;
    u[i97]=uni;
    if (--i97<0) i97=96;
    if (--j97<0) j97=96;
    c-=cd;
    if (c<0.0) c+=cm;
    uni-=c;
    if (uni<0.0) uni+=1.0;
    return(uni);
  }

  double ran1()
  {
    const int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32;
    const int NDIV=(1+(IM-1)/NTAB);
    const double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS);
    static int iy=0;
    static vector<int> iv(NTAB);
    int j,k;
    double temp;

    if (idum <= 0 || !iy) {
      if (-idum < 1) idum=1;
      else idum = -idum;
      for (j=NTAB+7;j>=0;j--) {
	k=idum/IQ;
	idum=IA*(idum-k*IQ)-IR*k;
	if (idum < 0) idum += IM;
	if (j < NTAB) iv[j] = idum;
      }
      iy=iv[0];
    }
    k=idum/IQ;
    idum=IA*(idum-k*IQ)-IR*k;
    if (idum < 0) idum += IM;
    j=iy/NDIV;
    iy=iv[j];
    iv[j] = idum;
    if ((temp=AM*iy) > RNMX) return RNMX;
    else return temp;
  }

  int discrete(int min, int max)
  {
    int ret= (int)(floor(ran1()*(max+1-min)+min));
    return ret;
  }

  void discrete(int min, int max, int& a, int& b)
  {
    a = (int)(floor(ran1()*(max+1-min)+min));

    do {
      b = (int)(floor(ran1()*(max+1-min)+min));
    } while (a==b);

  }

  void discrete(int min, int max, int& a, int& b, int& c)
  {
    a = (int)(floor(ran1()*(max+1-min)+min));

    do {
      b = (int)(floor(ran1()*(max+1-min)+min));
    } while (b==a);

    do {
      c = (int)(floor(ran1()*(max+1-min)+min));
    } while (c==b || c==a);

  }

  void discrete(int min, int max, int& a, int& b, int& c, int& d)
  {
    a = (int)(floor(ran1()*(max+1-min)+min));

    do {
      b = (int)(floor(ran1()*(max+1-min)+min));
    } while (b==a);

    do {
      c = (int)(floor(ran1()*(max+1-min)+min));
    } while (c==b || c==a);

    do {
      d = (int)(floor(ran1()*(max+1-min)+min));
    } while (d==c || d==b || d==a);

  }

  void discrete(int min, int max, int n, set<int>& Set)
  {
    Set.clear();
    for(int i=0; i<n; ) {
      int a = (int)(floor(ran1()*(max+1-min)+min));
      if(Set.find(a)==Set.end()) {
	Set.insert(a);
	i++;
      }
    }
  }

  void discrete(int min, int max, int n, vector<int>& Vec)
  {
    set<int> Set;

    Vec.clear();
    for(int i=0; i<n; ) {
      int a = (int)(floor(ran1()*(max+1-min)+min));
      if(Set.find(a)==Set.end()) {
	Set.insert(a);
	Vec.push_back(a);
	i++;
      }
    }
  }

  int discrete_logbin(int xmin, int xmax)
  {
    double ratio = 1.15;
    if (xmin==0)
      xmin=1;

    int Nbins = (int)(floor(log(xmax/xmin)/log(ratio))) + 1;
    int i = discrete(0, Nbins-1);
    int a = xmin * pow(ratio, (double)i);
    int b = a*ratio;
    if(b>xmax)
      b = xmax;

    int x = discrete(a,b);
    return x;
  }

  void randomsubset(int N, int Nc, vector<int>& S)
  {
    list<int> Lbig;
    for(int i=0; i<N; i++) {
      Lbig.push_back(i);
    }

    int N0= N;
    int M = (2*Nc>N0) ? (N0-Nc) : Nc;

    list<int> Lsmall;
    for(int i=0; i<M; i++) {
      int index = discrete(0, N-1);

      int j = Getnode(Lbig, index);
      Lsmall.push_back(j);
      Lbig.remove(j);
      N--;
    }

    S.clear();
    S.resize(Nc);

    if(2*Nc>N0)
      copy(Lbig.begin(), Lbig.end(), S.begin());
    else
      copy(Lsmall.begin(), Lsmall.end(), S.begin());

  }

  void randomsubset_knuth(int n, int m, vector<int>& S)
  {
    for (int i = 0; i < n; i++)
      if ((discrete(1,n) % (n-i)) < m) {
	S.push_back(i);
	m--;
      }
  }

  void randomsubset_shuf(int n, int m, vector<int>& S)
  {
    int i,j,t;
    vector<int> x(n);
    for (i = 0; i < n; i++)
      x[i] = i;
    for (i = 0; i < m; i++) {
      j = discrete(i, n-1);
      t = x[i];
      x[i] = x[j];
      x[j] = t;
    }
    sort(x.begin(), x.begin()+m);
    for (i = 0; i< m; i++)
      S.push_back(x[i]);
  }

  void randomsubset_shuffle(int n, int m, vector<int>& S)
  {
    vector<int> x(n);
    for(int i=0; i<n; i++)
      x[i] = i;

    random_shuffle(x.begin(), x.end());

    for (int i=0; i<m; i++)
      S.push_back(x[i]);
  }

  double gaussian(double sd)
  {
    double v1,v2,Rsq,z;
    if(outputReady)
      {
	outputReady = false;
	return output;
      }

    do
      {

	v1=(uniform()*2.0)-1.0;
	v2=(uniform()*2.0)-1.0;

	Rsq=v1*v1+v2*v2;
      } while (Rsq>=1.0);

    z=sqrt(-2.0*log(Rsq)/Rsq);
    output = v1*z*sd;
    outputReady=true;
    return v2*z*sd;
  }

  double lorentzian(double sd)
  {
    double v1,v2,z;
    do {
      v2 = uniform();
      v1 = 2.0 * uniform() - 1.0;
    } while ((v1*v1+v2*v2)>1.0);
    z = v2/v1;
    return sd*z/2.0;
  }

  double powerlaw_continuous(double alpha, double xmin)
  {
    return ( xmin*pow(1-ran1(), 1./(1-alpha) ) );
  }

  double powerlaw_continuous_interval(double a, double beta)
  {
    return pow( (1- (1-ran1())*(1-pow(a, 1-beta))), 1/(1-beta) );
  }

  int powerlaw_discrete_approximate(double alpha, double xmin)
  {
    double r = ran1();
    int k = (int)( (xmin-0.5)*pow(1-r, 1./(1-alpha) ) + 0.5 );
    if(k<0) {
      cout << 1-r << ',' << k << endl;
      k =  INTMAX;
    }
    return k;

  }

  int powerlaw_discrete(double alpha, int xmin)
  {
    if(!flag_setPvector)
      SetPvector(alpha,xmin);

    double r = ran1();
    int j = Findj(P, 1-r);
    return j;

  }

  int powerlaw_exponentialcutoff_discrete(double gamma, double kappa)
  {
    if(!flag_setPvector)
      SetPvector_SFEC(gamma, kappa);

    double r = ran1();

    int j = Findj(P, 1-r);
    return j;

  }

  int arbitrary_discrete(vector<double>& CDF)
  {
    double r = ran1();
    int j = Findj(CDF, r);
    return j+1;
  }

  int arbitrary_discrete(vector<long double>& CDF)
  {
    double r = ran1();
    int j = Findj(CDF, r);
    return j+1;
  }

  double lognormal(double sd)
  {
    double v1,v2,Rsq,z;
    if(outputReady)
      {
	outputReady = false;
	return output;
      }

    do
      {
	v1=(ran1()*2.0)-1.0;
	v2=(ran1()*2.0)-1.0;
	Rsq=v1*v1+v2*v2;
      } while (Rsq>=1.0);

    z=sqrt(-2.0*log(Rsq)/Rsq);
    output = exp(v1*z*sd);
    outputReady=true;
    return exp(v2*z*sd);
  }

  double bimodal(double R)
  {

    return (ran1()>0.5) ? (R) : (-R);
  }

  double rectangular(double R)
  {

    return ((ran1()*2-1)*R);
  }

  double parabolic(double R)
  {

    double x = ran1();

    int sign =  (x>0.5) ? (-1) : (+1);

    double R2 = R*R;
    double R3 = R2*R;

    double q     = 4*R3*(x-0.5);
    double p     = -3*R2;
    double Disc  = q*q/4. + p*p*p/27.;
    double r     = R3;
    double theta = atan(sqrt(-Disc)/(-q/2.));
    double ronethird = pow(r,1/3.);

    double h1 = sign * 2 * ronethird * cos(theta/3.);
    double h2 = sign * 2 * ronethird * cos((theta+2*PI)/3.);
    double h3 = sign * 2 * ronethird * cos((theta-2*PI)/3.);

    double y = 0;
    if(h1<R && h1>(-R))
      y = h1;
    else if(h2<R && h2>(-R))
      y = h2;
    else if(h3<R && h3>(-R))
      y = h3;
    else
      cout << "range is wrong!";

    return y;
  }

  double antiparabolic(double R)
  {

    double x = ran1();

    int sign =  (x>0.5) ? (+1) : (-1);

    double h = sign * pow(sign*(2*x-1),1./3.) * R;
    return h;
  }

  double DMI_lookuptable(double y, double R)

  {
    if(!flag_setFvector)
      SetFvector(y,R);

    double u = ran1();
    int j = Findj(F, u);
    double h = -R + j*dh;

    return h;
  }

  double DMI(double y, double R)
  {
    double u = ran1();
    int sign =  (u>0.5) ? (+1) : (-1);

    const double hacc = 1e-12;
    double htilde = rtbis(&fx, y, u, 0, 1, hacc);
    double h = sign*htilde*R;
    return h;
  }

  int exponential(double kappa)
  {
    double u = ran1();

    int k = -kappa * log(1-u);
    return k;
  }

 private:

  bool flag_setFvector;
  vector<double> F;
  double dh;

  double P_DMI(double h, double y, double R)
  {
    if(h<-R || h>R) return 0;
    else            return (y+1)/(2*y*R)*(1-pow(fabs(h)/R,y));
  }

  void SetFvector(double y, double R)
  {
    int N = 10000001;
    dh = R/(double)(N-1);
    F.resize(2*N-1,0);

    F[0]     = 0;
    F[2*N-2] = 1 - F[0];
    for(int i=1;i<N;i++)
      {
	double h1 = -R + i*dh;
	double h2 = h1 + dh;

	F[i]       = F[i-1] + 0.5* (P_DMI(h1,y,R)+P_DMI(h2,y,R))*dh;
	F[2*N-2-i] = 1- F[i];
      }

    flag_setFvector = true;
  }

  int Findj(vector<double>& xx, const double x)
  {
    int j,jl,jm,ju;
    int n=xx.size();

    jl=-1; ju=n;

    bool ascnd = (xx[n-1] >= xx[0]);
    while (ju-jl > 1)
      {
	jm=(ju+jl) >> 1;
	if ((x >= xx[jm]) == ascnd) jl=jm;
	else    	          ju=jm;
      }

    if (x == xx[0]) j=0;
    else if (x == xx[n-1]) j=n-2;
    else j=jl;

    return j;
  }

  int Findj(vector<long double>& xx, const double x)
  {
    int j,jl,jm,ju;
    int n=xx.size();

    jl=-1; ju=n;

    bool ascnd = (xx[n-1] >= xx[0]);
    while (ju-jl > 1)
      {
	jm=(ju+jl) >> 1;
	if ((x >= xx[jm]) == ascnd) jl=jm;
	else    	          ju=jm;
      }

    if (x == xx[0]) j=0;
    else if (x == xx[n-1]) j=n-2;
    else j=jl;

    return j;
  }

  double rtbis(FUNC func, const double y, const double u, const double x1, const double x2, const double xacc)
  {
    const int JMAX=40;
    int j;
    double dx,f,fmid,xmid,rtb;

    f=func(x1,y,u);
    fmid=func(x2,y,u);
    if (f*fmid >= 0.0)
      cout << "Root must be bracketed for bisection in rtbis\n";

    rtb = f < 0.0 ? (dx=x2-x1,x1) : (dx=x1-x2,x2);
    for (j=0;j<JMAX;j++) {
      fmid=func(xmid=rtb+(dx *= 0.5),y,u);
      if (fmid <= 0.0) rtb=xmid;
      if (fabs(dx) < xacc || fmid == 0.0) return rtb;
    }
    cout << "Too many bisections in rtbis\n";
    return 0.0;
  }

  void ranmarin(int ijkl)
  {
    int i,ii,j,jj,k,l,m ;
    double s,t ;

    int ij=ijkl/30082;
    int kl=ijkl-30082*ij;

    i=((ij/177)%177)+2 ;
    j=(ij%177)+2 ;
    k=((kl/169)%178)+1 ;
    l=kl%169 ;
    for (ii=0;ii<97;ii++) {
      s=0.0 ;
      t=0.5 ;
      for (jj=0;jj<24;jj++) {
	m=(((i*j)%179)*k)%179 ;
	i=j;
	j=k;
	k=m;
	l=(53*l+1)%169;
	if (((l*m)%64)>=32) s+=t;
	t*=0.5;
      }
      u[ii]=s;
    }
    c=362436.0/16777216.0;
    cd=7654321.0/16777216.0;
    cm=16777213.0/16777216.0;
    i97=96;
    j97=32;
  }

  bool flag_setPvector;
  vector<double> P;

  double Hzeta(double s, double q)
  {
    double sum = 0;
    double nterm = 0;
    int n;
    for(n=0; ;n++)
      {
	nterm = pow((n+q),-s);
	sum += nterm;

	if(fabs(nterm/sum) < TOLERANCE)
	  break;
      }
    cout << "nmax= " << n << endl;
    return sum;
  }

  void SetPvector(double alpha, int xmin)
  {

    double a = gsl_sf_hzeta(alpha, xmin);

    int N = 10000000+xmin;

    P.resize(N,2);

    P[xmin] = 1;
    for(int i=xmin;i<N-1;i++)
      {

	P[i+1] = P[i] - 1/(pow((double)i,alpha)*a);
      }

    flag_setPvector = true;
  }

double polylog(double s, double z)
{
  double sum = 0;
  double zn = z;
  int nmax = 1;
  for(int n=1; ; n++) {
    double term = zn/pow((double)n,s);
    sum += term;
    zn *= z;
    if(term < TOLERANCE) {
      nmax = n;
      break;
    }
  }

  return sum;
}

  void SetPvector_SFEC(double gamma, double kappa)
  {
    double enkappa= exp(-1/kappa);
    double C = polylog(gamma, enkappa);

    int N = 10000000;

    P.resize(N,2);

    P[1] = 1;
    double a = enkappa;
    for(int i=1;i<=N;i++) {
      P[i+1] = P[i] - pow((double)i,-gamma) * a/ C ;
      a *= enkappa;
    }

    flag_setPvector = true;
  }

};

#endif

max_stars_repo_name: spxuw/RFIM

max_issues_repo_name: spxuw/RFIM

max_forks_repo_name: spxuw/RFIM


text:
#ifndef GNUPLOT_H
#define GNUPLOT_H

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

#include "constants.h"

void writePlotData(gsl_vector *f, const std::string& filename);

#endif // GNUPLOT_H

max_stars_repo_name: ichi-rika/glottal-inverse

max_issues_repo_name: ichi-rika/glottal-inverse

max_forks_repo_name: ichi-rika/glottal-inverse


text:
/**
 * @file lapacke.h
 * @brief Header file to automatically handle Intel MKL or LAPACKE includes.
 *
 * `npypacke` can be linked with a reference LAPACKE implementation, OpenBLAS,
 * or Intel MKL, but Intel MKL uses a different header file, `mkl.h`, while
 * OpenBLAS uses the LAPACKE header `lapacke.h`. This file includes the correct
 * header file given an appropriate preprocessor macro is defined and defines
 * some types appropriately so that the user can simply write in terms of the
 * Intel MKL interface, whether or not Intel MKL is actually linked.
 */

#ifndef NPY_LPK_LAPACKE_H
#define NPY_LPK_LAPACKE_H
// if linking with Netlib LAPACKE
#if defined(LAPACKE_INCLUDE)
#include <lapacke.h>
// define MKL_INT used in extension modules as in mkl.h
#ifndef MKL_INT
#define MKL_INT int
#endif /* MKL_INT */
// else if linking with OpenBLAS CBLAS
#elif defined(OPENBLAS_INCLUDE)
#include <complex.h>
#include <lapacke.h>
// OpenBLAS has blasint typedef, so define MKL_INT using blasint. for LAPACKE
// routines, should cast to lapack_int (same size as blasint in OpenBLAS)
#ifndef MKL_INT
#define MKL_INT blasint
#endif /* MKL_INT */
// else if linking with Intel MKL LAPACKE implementation
#elif defined(MKL_INCLUDE)
#include <mkl.h>
// else error, no LAPACKE includes specified
#else
// silence error squiggles in VS Code (__INTELLISENSE__ always defined)
#ifndef __INTELLISENSE__
#error "no LAPACKE includes specified. try -D(LAPACKE|OPENBLAS|MKL)_INCLUDE"
#endif /* __INTELLISENSE__ */
#endif

// silence error squiggles in VS Code. use mkl.h since it also define MKL types
#ifdef __INTELLISENSE__
#include <mkl.h>
#endif /* __INTELLISENSE__ */

#endif /* NPY_LPK_LAPACKE_H */

max_issues_repo_name: phetdam/numpy-lapacke-demo

max_forks_repo_name: phetdam/numpy-lapacke-demo


text:
#pragma once

#include "SimpleEvaluator.h"
#include "BooleanDAG.h"
#include "Interface.h"
#include "SuggestionGenerator.h"
#include <map>
#include <set>
#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#else
#include "FakeGSL.h"
#endif

class SmartSuggestionGenerator1: public SuggestionGenerator {
    Interface* interf;
    BooleanDAG* dag;
    map<string, int>& ctrls;
    SimpleEvaluator* seval;
    const vector<vector<int>>& dependentInputs;
    const vector<vector<int>>& dependentCtrls;
    vector<int> inputsToAsserts; // TODO: this is pretty local hack
    int counter;
public:
    SmartSuggestionGenerator1(BooleanDAG* _dag, Interface* _interface, map<string, int>& _ctrls, const vector<vector<int>>& _dependentInputs, const vector<vector<int>>& _dependentCtrls): dag(_dag), interf(_interface), ctrls(_ctrls), dependentInputs(_dependentInputs), dependentCtrls(_dependentCtrls) {
        seval = new SimpleEvaluator(*dag, ctrls);
        seval->setInputs(interf);
        counter = 0;
        inputsToAsserts.resize(dag->size(), -1);
        for (int i = 0; i < dag->size(); i++) {
            bool_node* n = (*dag)[i];
            if (n->type == bool_node::ASSERT) {
                const auto& inputs = dependentInputs[i];
				for (auto it = inputs.begin(); it != inputs.end(); ++it) {
                    inputsToAsserts[*it] = i;
                }
            }
        }
    }

    pair<int, int> getSuggestion(const gsl_vector* state) {
        return make_pair(-1, 0);
    }

    pair<int, int> getUnsatSuggestion(const gsl_vector* state) {
        return make_pair(-1, 0);
    }

    set<bool_node*> getResponsibleNodes(bool_node* node) {
        Assert(node->type == bool_node::ASSERT, "dsfqhuopwqe");
        set<bool_node*> responsibleNodes;
        vector<bool_node*> nodes; // queue
        nodes.push_back(node);

        while(nodes.size() > 0) {
            bool_node* n = nodes.back();
            nodes.pop_back();
            Assert(n->type == bool_node::LT || n->type == bool_node::ASSERT || n->type == bool_node::NOT || n->type == bool_node::AND || n->type == bool_node::OR, "Unexpected node type");
            if (n->type == bool_node::LT) {
                responsibleNodes.insert(n);
            }
            if (n->type == bool_node::ASSERT || n->type == bool_node::NOT) {
                nodes.push_back(n->mother());
            }
            if (n->type == bool_node::AND || n->type == bool_node::OR) {
                double dist = seval->d(n);
                double mdist = seval->d(n->mother());
                double fdist = seval->d(n->father());
                if (dist * mdist > 0.0) {
                    nodes.push_back(n->mother());
                }
                if (dist * fdist > 0.0) {
                    nodes.push_back(n->father());
                }
            }
        }
        return responsibleNodes;
    }

    virtual IClause* getConflictClause(int level, LocalState * state) {
        return NULL;
    }

    virtual SClause* getUnsatClause(const gsl_vector* state) {
        return NULL;
    }

    virtual SClause* getUnsatClause(const gsl_vector* state, const gsl_vector* initState) {
        return NULL;
    }

    vector<tuple<int, int, int>> getUnsatSuggestions(const gsl_vector* state) {
        vector<tuple<int, int, int>> suggestions;
        cout << "state: " << Util::print(state) << endl;
        seval->run(state);

        double minNegDist = GradUtil::MAXVAL;
        int minNegDistId = -1;
        for (int i = 0; i < dag->size(); i++) {
            bool_node* n = (*dag)[i];
            if (n->type == bool_node::ASSERT) {
                double dist = seval->d(n);
                if (dist < 0.0) {
                    if (dist < minNegDist) {
                        minNegDist = dist;
                        minNegDistId = n->id;
                    }
                }
            }
        }

        cout << "MinNegNode: " << (*dag)[minNegDistId]->lprint() << " with error " << minNegDist << endl;

        /*cout << "python test.py " << counter++ << " \"" << Util::print(state) << "\" data.txt";
        string  assertMsg = ((ASSERT_node*)(*dag)[minNegDistId])->getMsg();
        size_t start = assertMsg.find("(");
        if (start != string::npos) {
            size_t end = assertMsg.find(")");
            cout << " \"" << assertMsg.substr(start, end - start + 1) << "\"";
        }

        const set<int>& inputConstraints = interface->getInputConstraints();
        cout << " \"";
        for (auto it = inputConstraints.begin(); it != inputConstraints.end(); it++) {
            string line = to_string(dependentCtrls[*it][0]);
            string assertMsg = ((ASSERT_node*)(*dag)[inputsToAsserts[*it]])->getMsg();
            string point = "";
            size_t start = assertMsg.find("(");
            if (start != string::npos) {
                size_t end = assertMsg.find(")");
                point = assertMsg.substr(start, end - start + 1);
            }
            cout << line << ":" << point << ":" << interface->getValue(*it) << ";" ;
        }
        cout << "\"" << endl;*/

        map<int, int> ctrlsToScore;
        for (int i = 0; i < dag->size(); i++) {
            bool_node* n = (*dag)[i];
            if (n->type == bool_node::ASSERT) {
                bool satisfied = seval->d(n) >= 0.0;
                const set<bool_node*>& responsibleNodes = getResponsibleNodes(n);
                for (auto it = responsibleNodes.begin(); it != responsibleNodes.end(); it++) {
                    const vector<int>& ctrls = dependentCtrls[(*it)->id];
                    int cost = satisfied ? 1 : -1;
                    for (int j = 0; j < ctrls.size(); j++) {
                        if (ctrlsToScore.find(ctrls[j]) == ctrlsToScore.end()) {
                            ctrlsToScore[ctrls[j]] = cost;
                        } else {
                            ctrlsToScore[ctrls[j]] += cost;
                        }
                    }
                }
            }
        }

        for (auto it = ctrlsToScore.begin(); it != ctrlsToScore.end(); it++) {
            cout << "(" <<  it->first << "," << it->second << "); ";
        }
        cout << endl;

        const set<bool_node*>& influentialInputs = getResponsibleNodes((*dag)[minNegDistId]); // this should be a subset of dependentInputs of minNegDistId

        vector<tuple<int, int, int>> s;
        for (auto it = influentialInputs.begin(); it != influentialInputs.end(); it++) {
            int inputId = (*it)->id;
            double inputDist = seval->d(*it);
            int val = inputDist >= 0.0 ? 0 : 1;
            const vector<int>& ctrls = dependentCtrls[inputId];
            int cost = 0;
            for (int j = 0; j < ctrls.size(); j++) {
                if (ctrlsToScore.find(ctrls[j]) != ctrlsToScore.end()) {
                    cost += ctrlsToScore[ctrls[j]];
                }
            }
            s.push_back(make_tuple(cost, inputId, val));
            cout << (*dag)[inputId]->lprint() << " " << seval->d((*dag)[inputId]) << "  [" << Util::print(dependentCtrls[inputId])  << "] " << cost << endl;
        }
        sort(s.begin(), s.end());
        reverse(s.begin(), s.end());
        for (int k = 0; k < s.size(); k++) {
            int nodeid = get<1>(s[k]);
            if (!interf->hasValue(nodeid)) {
                suggestions.push_back(make_tuple(0, nodeid, get<2>(s[k])));
            }
        }
        int lastIdx = suggestions.size() - 1;
        if (lastIdx > 0) {
            cout << "Suggesting " << ((*dag)[get<1>(suggestions[lastIdx])])->lprint() << " " << get<2>(suggestions[lastIdx]) << endl;
        }

        return suggestions;

    }
};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend