text:
/*
Simulation of a network of Izhikevich neurons, coupled by chemical synapses,
receiving uncorrelated, random input spiketrains.

The simulation was written by Ekkehard Ullner in 2013.

This is a modified version of E.U.'s original, stand-alone program,
for use e.g. via the ctypes ffi in Python.

Changelog:
- suport for networks with inhibitory synapses
- Neuron type switched from FitzHugh-Nagumo to Izhikevich
- Various simulation parameters exposed to allow arbitrary network topologies,
input patterns, neuron parameters etc.
- results are written to a provided buffer instead of text files
- the non-free, problematic[0] "Numerical Recipes" RAN1 random number generator
was replaced with a Mersenne Twister rng from the GNU scientific library.
    [0 "Random Numbers in Scientific Computing: An Introduction", Katzgraber, 2010
    http://arxiv.org/pdf/1005.4117.pdf]
- cosmetics and a minimum of documentation (i.e. the comments)

CK 2014-2016

The MIT License (MIT)

Copyright (c) 2013,2014 Ekkehard Ullner, Clemens Korndörfer

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 PARTICU
LAR 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 <string.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <time.h>
#include <assert.h>

void sagwas(){
    puts("simulation version Oct 12 2016");
}

/* function simulate
 * ------------------
 * Runs a single network simulation with the given parameters and writes the resulting
 * voltage traces into a given buffer.
 *
 * Parameters:
 *
 * N: Number of rows i in the simulated grid of neurons. (regardless of network topology,
 *    the network is a 'grid' where each node is adressed by an (i,j) coordinate.)
 * M: number of columns j in the grid
 *
 * limit: number of simulation steps to run
 *
 * KT: maximum number of coupled neighbors of any node (degree of the network)
 * connection: function pointer, where
 *             connection(i,j,l, '0' | '1') must return the first (second) coordinate of the l'th afferent to cell ij
 *             connection(i,j,l, 's') must return the strength of that connection
 *
 * recording_voltage: array of shape (M,N,limit) in which simulation results will be written in row first order
 * recording_recov: same as recording_voltage, for the neuron's recovery variable
 * recording_spikes: same as recording_voltage, for binary log of spike events
 *
 * inputc_exc: array of shape (M,N) in row-first order. modulates the rate of input pulses fed into
 *              each cell, values within [0,1]. aka "the stimulus".
 * inputc_inh: same for inhibitory input pulses
 * laminex: lambda_in_ex, the rate of excitatory input spikes to input-receiving nodes
 * lamini: lambda_in_in,  the rate of inhibitory input spikes to input-receiving nodes
 * izhi_a: parameter a of the izhikevich neuron
 * izhi_b: parameter b of the izhikevich neuron
 * izhi_reset_c: reset constant of the neuron
 * izhi_recovery_d: recovery reset constant
 * con_upstr_exc: synapse conductivity of the excitatory random inputs
 * con_upstr_inh: synapse conductivity of the inhibitory random inputs
 * # lateral synapses:
 * taus: duration of transmitter presence after spike
 * alp: rising time constant of fraction of open receptors
 * bet: decay time constant of fraction of open receptors
 * # inhibitory synapses:
 * tauni: duration of transmitter presence after spike
 * ani: rising time constant of fraction of open receptors
 * bni: decay time constant of fraction of open receptors
 * # external excitatory synapses:
 * tauna: duration of transmitter presence after spike
 * ana: rising time constant of fraction of open receptors
 * bna: decay time constant of fraction of open receptors
 * activation_noise: scale of white noise added to each neuron's activation variable
 * seed: random seed > 0
 * verbose: print a lot or not
 * delta_t: integration step width
*/

void simulate(int M,
              int N,
              int limit,
              int KT,
              double (*connection)(int,int,int,char),
              double* recording_voltage,
              double* recording_recov,
              double* recording_spikes,
              double* inputc_exc,
              double* inputc_inh,
              double laminex,
              double lamini,
              double izhi_a,
              double izhi_b,
              double izhi_reset_c,
              double izhi_recovery_d,
              double con_upstr_exc,
              double con_upstr_inh,
              double taus,
              double alp,
              double bet,
              double tauni,
              double ani,
              double bni,
              double tauna,
              double ana,
              double bna,
              double activation_noise,
              int seed,
              int verbose,
              double delta_t){

    // say hello & show the current parametrisation & input pattern.
    if (verbose){
        printf("starting simulation (Izhikevich neuron)\n");
        printf("seed: %d\nM: %d\nN: %d\nlimit: %d\nKT: %d\nlaminex: %f\nlamini: %f\nizhi_a: %f\nizhi_b: %f\nizhi_reset_c: %f\nizhi_recovery_d: %f\ncon_upstr_exc: %f\ncon_upstr_inh: %f\ntaus: %f\nalp: %f\nbet: %f\ntauni: %f\nani: %f\nbni: %f\ntauna: %f\nana: %f\nbna: %f\nactivation_noise: %f\ndelta_t: %f\n",
                seed, M, N, limit, KT, laminex, lamini, izhi_a, izhi_b, izhi_reset_c, izhi_recovery_d, con_upstr_exc, con_upstr_inh, taus, alp, bet, tauni, ani, bni, tauna, ana, bna, activation_noise, delta_t);

        int i1,j1;
        int print_exc, print_inh;
        for (j1 = 0; j1<M; j1++ ){
            for (i1 = 0; i1<N; i1++ ){
                print_exc = inputc_exc[i1 + N*j1]  >  0;
                print_inh = inputc_inh[i1 + N*j1]  >  0;
                if (print_exc && print_inh) printf(" ±");
                else if (print_exc) printf(" +");
                else if (print_inh) printf(" -");
                else printf(" .");
            }
            printf("\n");
        }
    }

    assert(seed>0);
    gsl_rng* rng  =  gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(rng, seed);

    // Initialisations & constant parameters:
    // neuron model:
    double xold[N][M],dx_h[N][M],xh[N][M],xnew[N][M]; // activation variable (previous value, first step, intermediate value, new value)
    double yold[N][M],dy_h[N][M],yh[N][M],ynew[N][M]; // recovery variable (previous value, first step, intermediate value, new value)
    double Isyn[N][M]; // collection of lateral input currents at each neuron in the current time step

    double spikes[N][M]; // recording of spike events in the current time step

    // lateral synapses:
    double rold[N][M],rh[N][M],rnew[N][M];  // fraction of open receptors (previous value, intermediate value, new value).

        // more precisely: r[i][j] is the fraction of open receptors at synapses anywhere in the network
        // that are the target of neuron (i,j). Since all synapses have the same rise and decay constants
        // alpha and beta, the effect of a spike at neuron i,j is the same at all of its target synapses.
        // We can therefore write the fractions of open receptors of all these target synapses as a
        // single variable r[i][j], associated with the source neuron i,j.

    double Tl[N][M]; // spike arrival times "T_j" (indexing follows the same logic as r[][])
    double Hs[N][M]; // transmitter concentrations "[T]_j" (...here, too)
    double Vna  = 0.0; // synaptic reversal potential (excitatory)
    double Vni  = -80; // synaptic reversal potential (inhibitory)

    // double taus = 0.01; // duration of transmitter presence after spike
    double Tm = 1.0; // maximum transmitter concentration

    // double alp = 8.0; // rising time constant of fraction of open receptors
    // double bet = 8.0; // decay time constant of fraction of open receptors

    // activatory external input synapses:
    double rnoa[N][M],rnha[N][M],rnna[N][M]; // fraction of open receptors
    double Tnla[N][M],Hna[N][M]; // spike arrival times and transmitter concentrations

    // double tauna  =  0.01; // duration of transmitter presence after spike
    double Tna  =  1.0; // maximum transmitter concentration

    // double ana  =  8.0; // rising time constant of fraction of open receptors
    // double bna  =  8.0; // decay time constant of fraction of open receptors

    // inhibitory external input synapses:
    double rnoi[N][M],rnhi[N][M],rnni[N][M]; // fraction of open receptors (previous value, intermediate value, new value)
    double Tnli[N][M],Hni[N][M]; // spike arrival times and transmitter concentrations

    // double tauni  =  0.01; // duration of transmitter presence after spike
    double Tni = 1.0; // maximum transmitter concentration

    // double ani = 8.0; // rising time constant of fraction of open receptors
    // double bni = 8.0; // decay time constant of fraction of open receptors

    // misc:
    double izhi_reset_threshold = 30; // spike detection threshold
    long step; // simulation step count
    // double delta_t = 0.001; // integration step width
    // double delta_t = 0.01; // integration step width
    double current_time = 0.0; // simulation time

    int i,j,l; // various loop variables over network nodes.

    int p1[N][M][KT],p2[N][M][KT]; // lookup tables holding the first (p1) and second (p2) coordinate of the KT'th afferent to neuron N,M
    float conductance_net[N][M][KT]; // lookup table: synapse conductance of the KT'th afferent to neuron N,M

    double in_exc[N][M]; // excitatory stimulus strength (input noise rate modulation) for neuron N,M
    double in_inh[N][M]; // inhibitory ...

    // remaining initialisations:
    for (j = 0;j<M;j++ ){
        for (i = 0;i<N;i++ ){
            in_exc[i][j] =  inputc_exc[i + N * j];
            in_inh[i][j] =  inputc_inh[i + N * j];

            xold[i][j] = 30*(gsl_rng_uniform(rng)-0.5);
            yold[i][j] = 30*(gsl_rng_uniform(rng)-0.5);
            rold[i][j] = 0.0;
            rnoa[i][j] = 0.0;
            rnoi[i][j] = 0.0;

            Tl[i][j] = -2 * taus;
            Hs[i][j] = 0.0;
            Isyn[i][j] = 0.0;
            Tnla[i][j] = -2 * tauna;
            Tnli[i][j] = -2 * tauni;
            Hna[i][j] = 0.0;
            Hni[i][j] = 0.0;
            // network connectivity lookup tables:
            for (l = 0;l<KT;l++ ){
                p1[i][j][l]  =  (int)connection(j,i,l,'0');
                p2[i][j][l]  =  (int)connection(j,i,l,'1');
                conductance_net[i][j][l] = connection(j,i,l,'s');
    }}}

    // // some more extreme verbosity for verification
    // printf("\n connectivity structure: \n");
    // for (l = 0; l<KT; l++){
    //     printf("\n\n%d. neighbour of each unit:\n\n", l);
    //     for (j = 0;j<M;j++ ){
    //         for (i = 0;i<N;i++ ){
    //             printf("(%d,%d)   ", p1[i][j][l], p2[i][j][l]);
    //         }
    //         printf("\n");
    //     }
    // }
    // for (l = 0; l<KT; l++){
    //     printf("\n\nconductance from %d. neighbour of each unit:\n\n", l);
    //     for (j = 0;j<M;j++ ){
    //         for (i = 0;i<N;i++ ){
    //             printf("%05.2f  ", conductance_net[i][j][l]);
    //         }
    //         printf("\n");
    //     }
    // }

    // begin stepwise integration by the implicit midpoint method
    for(step = 1; step<= limit;step++ ){
        if (verbose && (step % (limit/100)  ==  0)){
            printf("\r...%d%%",(int)(100 * step/(double)limit));
            fflush(stdout);
        }

        for (j = 0;j<M;j++ ){
            for (i = 0;i<N;i++ ){
                // set transmitter presence if there was a spike recently,
                // ...for synapses within the network..
                if (Tl[i][j] + taus>current_time) Hs[i][j] = 1.0 * Tm;
                else Hs[i][j] = 0.0;

                // ..for activatory external input synapses..
                if (Tnla[i][j] + tauna>current_time) Hna[i][j] = 1.0 * Tna;
                else  Hna[i][j] = 0.0;

                // ..and for inhibitory external input synapses.
                if (Tnli[i][j] + tauni>current_time) Hni[i][j] = 1.0 * Tni;
                else  Hni[i][j] = 0.0;

        }}

        // first integration step
        double x;
        double y;
        double I;
        double g;
        for (j = 0;j<M;j++ ){
            for (i = 0;i<N;i++ ){
                x = xold[i][j];
                y = yold[i][j];

                // collect lateral synaptic currents
                Isyn[i][j] = 0.0;
                for (l = 0;l<KT;l++ ){
                    if(p1[i][j][l]>=0 && p2[i][j][l]>=0){
                        // input currents to neuron (i,j) depend proportionally on
                        // r[ p1[i][j][l] ][ p2[i][j][l] ], that is, on the fraction of open
                        // receptors of all the synapses targeted by neighbour l projecting
                        // to neuron i,j. These synapses have a high fraction of open receptors
                        // if their source neuron recently fired a spike. Thus, input currents
                        // flow to neuron (i,j) if its neighbours have recently fired a spike.

                        // we use signed conductance values to encode excitatory vs inhibitory synapses.
                        // of course, conductance is in fact positive in both these synapses; what differs
                        // is the reversal potential.
                        g = conductance_net[i][j][l];
                        if (g > 0.0)
                            Isyn[i][j] -= g * rold[p1[i][j][l]][p2[i][j][l]] * (xold[i][j] - Vna);
                        else if (g < 0.0)
                            Isyn[i][j] -= -g * rold[p1[i][j][l]][p2[i][j][l]] * (xold[i][j] - Vni);

                    }
                }

                I = Isyn[i][j];

                // collect external synaptic currents
                I -= con_upstr_exc * rnoa[i][j] * (xold[i][j] - Vna);
                I -= con_upstr_inh * rnoi[i][j] * (xold[i][j] - Vni);

                // Izhikevich neuron:
                // dv/dt = 0.04v^2 + 5v + 140 - u + I
                // du/dt = a(bv - u)
                dx_h[i][j] = 0.04*x*x + 5*x + 140 - y + I;
                dy_h[i][j] = izhi_a * (izhi_b*x - y);

                // step:
                xh[i][j] = xold[i][j] + (dx_h[i][j])*delta_t;
                yh[i][j] = yold[i][j] + (dy_h[i][j])*delta_t;

                // fraction of open receptors.. dr/dt  =  alpha [T] (1 - r) -  beta r
                // ..of synapses within the network:
                rh[i][j] = rold[i][j] + (alp * Hs[i][j] * (1 - rold[i][j]) - bet * rold[i][j]) * delta_t;

                // ..of synapses receiving external noise, activatory:
                rnha[i][j] = rnoa[i][j] + (ana * Hna[i][j] * (1 - rnoa[i][j]) - bna * rnoa[i][j]) * delta_t;
                // ..of synapses receiving external noise, inhibitory:
                rnhi[i][j] = rnoi[i][j] + (ani * Hni[i][j] * (1 - rnoi[i][j]) - bni * rnoi[i][j]) * delta_t;

        }}

        // second integration step
        double dx;
        double dy;
        double Ihsyn;
        for (j = 0;j<M;j++ ){
            for (i = 0;i<N;i++ ){
                x = xh[i][j];
                y = yh[i][j];

                Ihsyn = 0.0;
                for (l = 0;l<KT;l++ ){
                    if(p1[i][j][l]>=0 && p2[i][j][l]>=0){
                        g = conductance_net[i][j][l];
                        if (g > 0)
                            Ihsyn -= g * rh[p1[i][j][l]][p2[i][j][l]] * (xh[i][j] - Vna);
                        else if (g < 0)
                            Ihsyn -= -g * rh[p1[i][j][l]][p2[i][j][l]] * (xh[i][j] - Vni);
                    }
                }

                I =  0.5 * (Isyn[i][j] + Ihsyn);

                // input currents from external synapses
                I -= 0.5 * con_upstr_exc * (rnoa[i][j] * (xold[i][j] - Vna)  +  rnha[i][j] * (xh[i][j] - Vna));
                I -= 0.5 * con_upstr_inh * (rnoi[i][j] * (xold[i][j] - Vni)  +  rnhi[i][j] * (xh[i][j] - Vni));

                // Izhikevich neuron:
                // dv/dt = 0.04v^2 + 5v + 140 - u + I
                // du/dt = a(bv -u)
                dx = 0.04*x*x + 5*x + 140 - y + I;
                dy = izhi_a * (izhi_b*x - y);

                xnew[i][j] = xold[i][j] + 0.5*(dx_h[i][j] + dx)*delta_t;
                ynew[i][j] = yold[i][j] + 0.5*(dy_h[i][j] + dy)*delta_t;

                // fraction of open receptors.. dr/dt  =  alpha [T] (1 - r) -  beta r
                // ..of synapses within the network:
                rnew[i][j] = rold[i][j] + 0.5 * (alp * Hs[i][j] * (1 - rold[i][j]) - bet * rold[i][j]  +  alp * Hs[i][j] * (1 - rh[i][j]) - bet * rh[i][j]) * delta_t;
                // ..of synapses receiving external noise, activatory:
                rnna[i][j] = rnoa[i][j] + 0.5 * (ana * Hna[i][j] * (1 - rnoa[i][j]) - bna * rnoa[i][j]  +  ana * Hna[i][j] * (1 - rnha[i][j]) - bna * rnha[i][j]) * delta_t;
                // ..of synapses receiving external noise, inhibitory:
                rnni[i][j] = rnoi[i][j] + 0.5 * (ani * Hni[i][j] * (1 - rnoi[i][j]) - bni * rnoi[i][j]  +  ani * Hni[i][j] * (1 - rnhi[i][j]) - bni * rnhi[i][j]) * delta_t;

        }}

        // identify spike times & prepare next step
        current_time = delta_t * step;
        for (j = 0;j<M;j++ ){
            for (i = 0;i<N;i++ ){
                // note spike times within the network
                if(xnew[i][j] >= izhi_reset_threshold){
                    xnew[i][j] = izhi_reset_c;
                    ynew[i][j] += izhi_recovery_d;
                    spikes[i][j] = 1;
                    Tl[i][j] = current_time;
                }
                else spikes[i][j] = 0;

                // sample random spike times for activatory external input
                if(gsl_rng_uniform(rng) <=  (laminex * (double)in_exc[i][j]) * delta_t){
                    Tnla[i][j] = current_time;
                }
                // sample random spike times for inhibitory external input
                if(gsl_rng_uniform(rng)<= (lamini * (double)in_inh[i][j]) * delta_t){
                   Tnli[i][j] = current_time;
                }

                // swap
                xold[i][j] = xnew[i][j];
                yold[i][j] = ynew[i][j];
                rold[i][j] = rnew[i][j];
                rnoa[i][j] = rnna[i][j];
                rnoi[i][j] = rnni[i][j];

                if (activation_noise > 0)
                    xold[i][j] += (gsl_rng_uniform(rng)-0.5)*activation_noise;
        }}

        // write the current network state into the provided output buffers.
        // each is a 3D array of shape (M,N,limit) in row first order,
        // so it has strides proportional to (N * limit, limit, 1).
        int linearind  =  0;
        for(int oi  =  0; oi < M; oi++ ) {
            for(int oj  =  0; oj < N; oj++ ) {
                    linearind  =  oi * limit * N  +  oj * limit  + step - 1;
                    recording_voltage[linearind]  =  xnew[oj][oi];
                    recording_recov[linearind]  = ynew[oj][oi];
                    //recording_recov[linearind] = (double) Tnla[oj][oi] == current_time; // <-uncomment to show input spikes instead
                    recording_spikes[linearind]  =  spikes[oj][oi];
            }

        }
    } // end of integration loop
    if (verbose){
        printf("\033[2K");
        fflush(stdout);
    }

    gsl_rng_free(rng);
}

// helper function to play with the random number generator.
void rantest(long seed,int N,double *  out){
    gsl_rng *  rng  =  gsl_rng_alloc(gsl_rng_mt19937); // initialize a mersenne twister rng
    gsl_rng_set(rng, seed); // seed it
    double number  =  - 1;

    for(int i = 0; i<N; i++ ){
        number  =  gsl_rng_uniform(rng); // sample a number
        printf("%ld  %f\n",seed,number);
        out[i]  =  number;
    }
    puts(" - - - ");

    gsl_rng_free(rng);
}

max_stars_repo_name: cknd/synchrony

max_issues_repo_name: cknd/synchrony

max_forks_repo_name: cknd/synchrony


text:
#include <bindings.cmacros.h>
#include <gsl/gsl_multimin.h>

BC_INLINE2(GSL_MULTIMIN_FN_EVAL,gsl_multimin_function*,gsl_vector*,double)
BC_INLINE2(GSL_MULTIMIN_FN_EVAL_F,gsl_multimin_function_fdf*,gsl_vector*,double)
BC_INLINE3VOID(GSL_MULTIMIN_FN_EVAL_DF,gsl_multimin_function_fdf*,gsl_vector*,gsl_vector*)
BC_INLINE4VOID(GSL_MULTIMIN_FN_EVAL_F_DF,gsl_multimin_function_fdf*,gsl_vector*,double*,gsl_vector*)

ext: c


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

#ifndef WSS_BOARD_PREMIUMS_H
#define WSS_BOARD_PREMIUMS_H

#include <gsl/string_span>

#include "board.h"

enum class premium {
    normal,
    dl,
    tl,
    dw,
    tw
};

constexpr std::array<int, 5> letter_multipliers {
        1,
        2,
        3,
        1,
        1
};

constexpr std::array<int, 5> word_multipliers {
        1,
        1,
        1,
        2,
        3
};

auto load_board_premiums(gsl::cstring_span<> filename)
-> std::optional<board<premium>>;

#endif //WSS_BOARD_PREMIUMS_H

ext: h

max_stars_repo_name: johnmcfarlane/wss

max_issues_repo_name: johnmcfarlane/wss

max_forks_repo_name: johnmcfarlane/wss


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

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_filter.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_statistics.h>

/* compute filtered data by explicitely constructing window, sorting it and finding median */
int
slow_movmedian(const gsl_filter_end_t etype, const gsl_vector * x, gsl_vector * y, const int K)
{
  const size_t H = K / 2;
  const size_t n = x->size;
  double *window = malloc(K * sizeof(double));
  size_t i;

  for (i = 0; i < n; ++i)
    {
      size_t wsize = gsl_movstat_fill(etype, x, i, H, H, window);
      double yi = gsl_stats_median(window, 1, wsize);

      gsl_vector_set(y, i, yi);
    }

  free(window);

  return GSL_SUCCESS;
}

static void
test_median_proc(const double tol, const size_t n, const size_t K,
                 const gsl_filter_end_t etype, gsl_rng *rng_p)
{
  gsl_filter_median_workspace *w = gsl_filter_median_alloc(K);
  gsl_vector *x = gsl_vector_alloc(n);
  gsl_vector *y = gsl_vector_alloc(n);
  gsl_vector *z = gsl_vector_alloc(n);
  char buf[2048];

  /* test moving median with random input */
  random_vector(x, rng_p);

  /* y = median(x) with slow brute force algorithm */
  slow_movmedian(etype, x, y, K);

  /* z = median(x) */
  gsl_filter_median(etype, x, z, w);

  /* test y = z */
  sprintf(buf, "n=%zu K=%zu endtype=%u median random", n, K, etype);
  compare_vectors(tol, z, y, buf);

  /* z = median(x) in-place */
  gsl_vector_memcpy(z, x);
  gsl_filter_median(etype, z, z, w);

  sprintf(buf, "n=%zu K=%zu endtype=%u median random in-place", n, K, etype);
  compare_vectors(tol, z, y, buf);

  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(z);
  gsl_filter_median_free(w);
}

static void
test_median(gsl_rng * rng_p)
{
  test_median_proc(GSL_DBL_EPSILON, 1000, 1, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 1000, 3, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 1000, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 100, 301, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 17, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 9, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 901, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 201, GSL_MOVSTAT_END_PADZERO, rng_p);

  test_median_proc(GSL_DBL_EPSILON, 1000, 1, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 1000, 3, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 1000, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 100, 301, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 17, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 9, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 901, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 201, GSL_MOVSTAT_END_PADVALUE, rng_p);

  test_median_proc(GSL_DBL_EPSILON, 1000, 1, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 1000, 3, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 1000, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 100, 301, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 17, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 5000, 9, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 901, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_median_proc(GSL_DBL_EPSILON, 50, 201, GSL_MOVSTAT_END_TRUNCATE, rng_p);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/*
    Licensed to the OpenAirInterface (OAI) Software Alliance under one or more
    contributor license agreements.  See the NOTICE file distributed with
    this work for additional information regarding copyright ownership.
    The OpenAirInterface Software Alliance licenses this file to You under
    the OAI Public License, Version 1.1  (the "License"); you may not use this file
    except in compliance with the License.
    You may obtain a copy of the License at

        http://www.openairinterface.org/?page_id=698

    Unless required by applicable law or agreed to in writing, software
    distributed under the License is distributed on an "AS IS" BASIS,
    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    See the License for the specific language governing permissions and
    limitations under the License.
    -------------------------------------------------------------------------------
    For more information about the OpenAirInterface (OAI) Software Alliance:
        contact@openairinterface.org
*/

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

#include "SIMULATION/TOOLS/sim.h"
#include "PHY/TOOLS/tools_defs.h"
#include "assertions.h"

// NEW code with lookup table for sin/cos based on delay profile (TO BE TESTED)

double **cos_lut = NULL, * *sin_lut = NULL;

//#if 1

int init_freq_channel(channel_desc_t *desc, uint16_t nb_rb, int16_t n_samples)
{
    double delta_f, freq; // 90 kHz spacing
    double delay;
    int16_t f;
    uint8_t l;

    if((n_samples & 1) == 0)
    {
        fprintf(stderr, "freq_channel_init: n_samples has to be odd\n");
        return(-1);
    }

    cos_lut = (double **)malloc(n_samples * sizeof(double *));
    sin_lut = (double **)malloc(n_samples * sizeof(double *));
    delta_f = nb_rb * 180000 / (n_samples - 1);

    for(f = -(n_samples >> 1); f <= (n_samples >> 1); f++)
    {
        freq = delta_f * (double)f * 1e-6; // due to the fact that delays is in mus
        cos_lut[f + (n_samples >> 1)] = (double *)malloc((int)desc->nb_taps * sizeof(double));
        sin_lut[f + (n_samples >> 1)] = (double *)malloc((int)desc->nb_taps * sizeof(double));

        for(l = 0; l < (int)desc->nb_taps; l++)
        {
            if(desc->nb_taps == 1)
            {
                delay = desc->delays[l];
            }
            else
            {
                delay = desc->delays[l] + NB_SAMPLES_CHANNEL_OFFSET / desc->sampling_rate;
            }

            cos_lut[f + (n_samples >> 1)][l] = cos(2 * M_PI * freq * delay);
            sin_lut[f + (n_samples >> 1)][l] = sin(2 * M_PI * freq * delay);
            //printf("values cos:%d, sin:%d\n", cos_lut[f][l], sin_lut[f][l]);
        }
    }

    return(0);
}

int freq_channel(channel_desc_t *desc, uint16_t nb_rb, int16_t n_samples)
{
    int16_t f, f2, d;
    uint8_t aarx, aatx, l;
    double *clut, *slut;
    static int freq_channel_init = 0;
    static int n_samples_max = 0;

    // do some error checking
    // n_samples has to be a odd number because we assume the spectrum is symmetric around the DC and includes the DC
    if((n_samples & 1) == 0)
    {
        fprintf(stderr, "freq_channel: n_samples has to be odd\n");
        return(-1);
    }

    // printf("no of taps:%d,",(int)desc->nb_taps);

    if(freq_channel_init == 0)
    {
        // we are initializing the lut for the largets possible n_samples=12*nb_rb+1
        // if called with n_samples<12*nb_rb+1, we decimate the lut
        n_samples_max = 12 * nb_rb + 1;

        if(init_freq_channel(desc, nb_rb, n_samples_max) == 0)
        {
            freq_channel_init = 1;
        }
        else
        {
            return(-1);
        }
    }

    d = (n_samples_max - 1) / (n_samples - 1);
    //printf("no_samples=%d, n_samples_max=%d, d=%d\n",n_samples,n_samples_max,d);
    start_meas(&desc->interp_freq);

    for(f = -n_samples_max / 2, f2 = -n_samples / 2; f < n_samples_max / 2; f += d, f2++)
    {
        clut = cos_lut[n_samples_max / 2 + f];
        slut = sin_lut[n_samples_max / 2 + f];

        for(aarx = 0; aarx < desc->nb_rx; aarx++)
        {
            for(aatx = 0; aatx < desc->nb_tx; aatx++)
            {
                desc->chF[aarx + (aatx * desc->nb_rx)][n_samples / 2 + f2].x = 0.0;
                desc->chF[aarx + (aatx * desc->nb_rx)][n_samples / 2 + f2].y = 0.0;

                for(l = 0; l < (int)desc->nb_taps; l++)
                {
                    desc->chF[aarx + (aatx * desc->nb_rx)][n_samples / 2 + f2].x += (desc->a[l][aarx + (aatx * desc->nb_rx)].x * clut[l] +
                            desc->a[l][aarx + (aatx * desc->nb_rx)].y * slut[l]);
                    desc->chF[aarx + (aatx * desc->nb_rx)][n_samples / 2 + f2].y += (-desc->a[l][aarx + (aatx * desc->nb_rx)].x * slut[l] +
                            desc->a[l][aarx + (aatx * desc->nb_rx)].y * clut[l]);
                }
            }
        }
    }

    stop_meas(&desc->interp_freq);
    return(0);
}

double compute_pbch_sinr(channel_desc_t *desc,
                         channel_desc_t *desc_i1,
                         channel_desc_t *desc_i2,
                         double snr_dB, double snr_i1_dB,
                         double snr_i2_dB,
                         uint16_t nb_rb)
{
    double avg_sinr, snr = pow(10.0, .1 * snr_dB), snr_i1 = pow(10.0, .1 * snr_i1_dB), snr_i2 = pow(10.0, .1 * snr_i2_dB);
    uint16_t f;
    uint8_t aarx, aatx;
    double S;
    struct complex S_i1;
    struct complex S_i2;
    avg_sinr = 0.0;

    //  printf("nb_rb %d\n",nb_rb);
    for(f = (nb_rb - 6); f < (nb_rb + 6); f++)
    {
        S = 0.0;
        S_i1.x = 0.0;
        S_i1.y = 0.0;
        S_i2.x = 0.0;
        S_i2.y = 0.0;

        for(aarx = 0; aarx < desc->nb_rx; aarx++)
        {
            for(aatx = 0; aatx < desc->nb_tx; aatx++)
            {
                S    += (desc->chF[aarx + (aatx * desc->nb_rx)][f].x * desc->chF[aarx + (aatx * desc->nb_rx)][f].x +
                         desc->chF[aarx + (aatx * desc->nb_rx)][f].y * desc->chF[aarx + (aatx * desc->nb_rx)][f].y);
                //  printf("%d %d chF[%d] => (%f,%f)\n",aarx,aatx,f,desc->chF[aarx+(aatx*desc->nb_rx)][f].x,desc->chF[aarx+(aatx*desc->nb_rx)][f].y);

                if(desc_i1)
                {
                    S_i1.x += (desc->chF[aarx + (aatx * desc->nb_rx)][f].x * desc_i1->chF[aarx + (aatx * desc->nb_rx)][f].x +
                               desc->chF[aarx + (aatx * desc->nb_rx)][f].y * desc_i1->chF[aarx + (aatx * desc->nb_rx)][f].y);
                    S_i1.y += (desc->chF[aarx + (aatx * desc->nb_rx)][f].x * desc_i1->chF[aarx + (aatx * desc->nb_rx)][f].y -
                               desc->chF[aarx + (aatx * desc->nb_rx)][f].y * desc_i1->chF[aarx + (aatx * desc->nb_rx)][f].x);
                }

                if(desc_i2)
                {
                    S_i2.x += (desc->chF[aarx + (aatx * desc->nb_rx)][f].x * desc_i2->chF[aarx + (aatx * desc->nb_rx)][f].x +
                               desc->chF[aarx + (aatx * desc->nb_rx)][f].y * desc_i2->chF[aarx + (aatx * desc->nb_rx)][f].y);
                    S_i2.y += (desc->chF[aarx + (aatx * desc->nb_rx)][f].x * desc_i2->chF[aarx + (aatx * desc->nb_rx)][f].y -
                               desc->chF[aarx + (aatx * desc->nb_rx)][f].y * desc_i2->chF[aarx + (aatx * desc->nb_rx)][f].x);
                }
            }
        }

        //    printf("snr %f f %d : S %f, S_i1 %f, S_i2 %f\n",snr,f-nb_rb,S,snr_i1*sqrt(S_i1.x*S_i1.x + S_i1.y*S_i1.y),snr_i2*sqrt(S_i2.x*S_i2.x + S_i2.y*S_i2.y));
        avg_sinr += (snr * S / (desc->nb_tx + snr_i1 * sqrt(S_i1.x * S_i1.x + S_i1.y * S_i1.y) + snr_i2 * sqrt(S_i2.x * S_i2.x + S_i2.y * S_i2.y)));
    }

    //  printf("avg_sinr %f (%f,%f,%f)\n",avg_sinr/12.0,snr,snr_i1,snr_i2);
    return(10 * log10(avg_sinr / 12.0));
}

double compute_sinr(channel_desc_t *desc,
                    channel_desc_t *desc_i1,
                    channel_desc_t *desc_i2,
                    double snr_dB, double snr_i1_dB,
                    double snr_i2_dB,
                    uint16_t nb_rb)
{
    double avg_sinr, snr = pow(10.0, .1 * snr_dB), snr_i1 = pow(10.0, .1 * snr_i1_dB), snr_i2 = pow(10.0, .1 * snr_i2_dB);
    uint16_t f;
    uint8_t aarx, aatx;
    double S;
    struct complex S_i1;
    struct complex S_i2;
    DevAssert(nb_rb > 0);
    avg_sinr = 0.0;

    //  printf("nb_rb %d\n",nb_rb);
    for(f = 0; f < 2 * nb_rb; f++)
    {
        S = 0.0;
        S_i1.x = 0.0;
        S_i1.y = 0.0;
        S_i2.x = 0.0;
        S_i2.y = 0.0;

        for(aarx = 0; aarx < desc->nb_rx; aarx++)
        {
            for(aatx = 0; aatx < desc->nb_tx; aatx++)
            {
                S    += (desc->chF[aarx + (aatx * desc->nb_rx)][f].x * desc->chF[aarx + (aatx * desc->nb_rx)][f].x +
                         desc->chF[aarx + (aatx * desc->nb_rx)][f].y * desc->chF[aarx + (aatx * desc->nb_rx)][f].y);

                if(desc_i1)
                {
                    S_i1.x += (desc->chF[aarx + (aatx * desc->nb_rx)][f].x * desc_i1->chF[aarx + (aatx * desc->nb_rx)][f].x +
                               desc->chF[aarx + (aatx * desc->nb_rx)][f].y * desc_i1->chF[aarx + (aatx * desc->nb_rx)][f].y);
                    S_i1.y += (desc->chF[aarx + (aatx * desc->nb_rx)][f].x * desc_i1->chF[aarx + (aatx * desc->nb_rx)][f].y -
                               desc->chF[aarx + (aatx * desc->nb_rx)][f].y * desc_i1->chF[aarx + (aatx * desc->nb_rx)][f].x);
                }

                if(desc_i2)
                {
                    S_i2.x += (desc->chF[aarx + (aatx * desc->nb_rx)][f].x * desc_i2->chF[aarx + (aatx * desc->nb_rx)][f].x +
                               desc->chF[aarx + (aatx * desc->nb_rx)][f].y * desc_i2->chF[aarx + (aatx * desc->nb_rx)][f].y);
                    S_i2.y += (desc->chF[aarx + (aatx * desc->nb_rx)][f].x * desc_i2->chF[aarx + (aatx * desc->nb_rx)][f].y -
                               desc->chF[aarx + (aatx * desc->nb_rx)][f].y * desc_i2->chF[aarx + (aatx * desc->nb_rx)][f].x);
                }
            }
        }

        //        printf("f %d : S %f, S_i1 %f, S_i2 %f\n",f-nb_rb,snr*S,snr_i1*sqrt(S_i1.x*S_i1.x + S_i1.y*S_i1.y),snr_i2*sqrt(S_i2.x*S_i2.x + S_i2.y*S_i2.y));
        avg_sinr += (snr * S / (desc->nb_tx + snr_i1 * sqrt(S_i1.x * S_i1.x + S_i1.y * S_i1.y) + snr_i2 * sqrt(S_i2.x * S_i2.x + S_i2.y * S_i2.y)));
    }

    //  printf("avg_sinr %f (%f,%f,%f)\n",avg_sinr/12.0,snr,snr_i1,snr_i2);
    return(10 * log10(avg_sinr / (nb_rb * 2)));
}

int pbch_polynomial_degree = 6;
double pbch_awgn_polynomial[7] = {-7.2926e-05, -2.8749e-03, -4.5064e-02, -3.5301e-01, -1.4655e+00, -3.6282e+00, -6.6907e+00};

void load_pbch_desc(FILE *pbch_file_fd)
{
    int i, ret;
    char dummy[25];
    ret = fscanf(pbch_file_fd, "%d", &pbch_polynomial_degree);

    if(ret < 0)
    {
        printf("fscanf failed: %s\n", strerror(errno));
        exit(EXIT_FAILURE);
    }

    if(pbch_polynomial_degree > 6)
    {
        printf("Illegal degree for pbch interpolation polynomial %d\n", pbch_polynomial_degree);
        exit(-1);
    }

    printf("PBCH polynomial : ");

    for(i = 0; i <= pbch_polynomial_degree; i++)
    {
        ret = fscanf(pbch_file_fd, "%24s", dummy);

        if(ret < 0)
        {
            printf("fscanf failed: %s\n", strerror(errno));
            exit(EXIT_FAILURE);
        }

        pbch_awgn_polynomial[i] = strtod(dummy, NULL);
        printf("%f ", pbch_awgn_polynomial[i]);
    }

    printf("\n");
}

double pbch_bler(double sinr)
{
    int i;
    double log10_bler = pbch_awgn_polynomial[pbch_polynomial_degree];
    double sinrpow = sinr;
    double bler = 0.0;

    //  printf("log10bler %f\n",log10_bler);
    if(sinr < -10.0)
    {
        bler = 1.0;
    }
    else if(sinr >= 0.0)
    {
        bler = 0.0;
    }
    else
    {
        for(i = 1; i <= pbch_polynomial_degree; i++)
        {
            //    printf("sinrpow %f\n",sinrpow);
            log10_bler += (pbch_awgn_polynomial[pbch_polynomial_degree - i] * sinrpow);
            sinrpow *= sinr;
            //    printf("log10bler %f\n",log10_bler);
        }

        bler = pow(10.0, log10_bler);
    }

    //printf ("sinr %f bler %f\n",sinr,bler);
    return(bler);
}