text:
/*
Simulation of a network of FitzHugh-Nagumo neurons, coupled by excitatory
chemical synapses, receiving decorrelated, random input spiketrains with
both excitatory and inhibitory components.

The simulation was written by Ekkehard Ullner in 2013.

This is a modified version of E.U.'s original, stand-alone program,
for use as a library (e.g. via the ctypes ffi in Python).
Changelog:
- Various simulation parameters are 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 and problematic[0] "Numerical Recipes" RAN1 random number generator
was replaced with a Mersenne Twister rng from the GNU scientific library.
- cosmetics and a minimum of documentation (i.e. the comments)

C Korndoerfer 2014

[0 "Random Numbers in Scientific Computing: An Introduction", Katzgraber, 2010
http://arxiv.org/pdf/1005.4117.pdf]

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: max. simulation time
 * rec_resolution: to store only every ~'th simulation step.
 *
 * 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
 *
 * outputbuffer: array of shape (M,N,limit) in which simulation results will be written in row first order
 *
 * inputc: 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".
 * laminex: firing rate "lambda_in,ex", the base rate of excitatory input spikes modulated by inputc
 * lamini: firing rate "lambda_in,in", the base rate of inhibitory input spikes modulated by inputc
 *
 * double fhn_a: parameter a of the FHN neuron
 * double fhn_eps: parameter epsilon of the FHN neuron
 *
 * double con_upstr_exc: synapse conductance "g^{up,ex}" of the excitatory random external inputs.
 * double con_upstr_inh: synapse conductance "g^{up,in}" of the inhibitory random external inputs.
 *
 * seed: a random seed > 0
 *
 * delta_t : integration step width
 *
 * Return: None (simulation results are passed back as a side effect, by filling the provided output buffer)
*/

void simulate(int M,
              int N,
              int limit,
              int rec_resolution,
              int KT,
              double (*connection)(int,int,int,char),
              double* outputbuffer,
              double* debug_outbuf,
              double* inputc_exc,
              double* inputc_inh,
              double laminex,
              double lamini,
              double fhn_a,
              double fhn_eps,
              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,
              int seed,
              int verbose,
              double delta_t){

    // say hello & show the current input pattern.
    if (verbose){
        printf(" simulation with seed %d\n",seed);
        printf("laminex: %f\nlamini: %f\nfhn_a: %f\nfhn_eps: %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\n", laminex, lamini, fhn_a, fhn_eps, con_upstr_exc,
               con_upstr_inh, taus, alp, bet, tauni, ani, bni, tauna, ana, bna);
        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],xh[N][M],xnew[N][M]; // activation variable (previous value, intermediate value, new value)
    double yold[N][M],yh[N][M],ynew[N][M]; // recovery variable (previous value, intermediate value, new value)
    double Isyn[N][M], Ihsyn[N][M]; // lateral input current (Ihsyn: intermediate value)

    // 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 Vsyn = 0.0; // synaptic reversal potential

    // 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 Vna  =  0.0; // synaptic reversal potential

    // 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 Vni  =  - 2.1; // synaptic reversal potential

    // 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 xth = 1.0; // 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, in some irrelevant real - numbered unit

    long i,j,l; // various loop variables over network nodes. don't ask me why those were chosen as long ints.

    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] = -1.05 + 0.2 * (gsl_rng_uniform(rng) - 0.5);
            yold[i][j] = -0.66 + 0.2 * (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;
            Ihsyn[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');
    }}}

    // begin stepwise integration.
    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;

                // collect lateral synaptic currents (for first integration step)
                Isyn[i][j] = 0.0;
                for (l = 0;l<KT;l++ ){
                    if(p1[i][j][l]<N && p2[i][j][l]<M)
                        Isyn[i][j]  +=   conductance_net[i][j][l]  *  rold[p1[i][j][l]][p2[i][j][l]] * (xold[i][j] - Vsyn);
                        // Here, we add input currents to neuron (i,j) that depend mostly on
                        // r[ p1[i][j][l] ][ p2[i][j][l] ], that is, on the fraction of open
                        // receptors of those synapses targeted by neighbours 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.
                }

        }}

        // first integration step
        for (j = 0;j<M;j++ ){
            for (i = 0;i<N;i++ ){
                //FHN neuron: eps * (dx/dt)  =  x -  x^3 / 3 -  y
                xh[i][j]  =  xold[i][j]  +  ( (1./fhn_eps) * (xold[i][j] -  xold[i][j] * xold[i][j] * xold[i][j]/3. -  yold[i][j]) )  *  delta_t;
                xh[i][j] +=  - (1./fhn_eps) * Isyn[i][j] * delta_t;
                //FHN neuron: (dy/dt)  =  x + a
                yh[i][j] = yold[i][j] + (xold[i][j] + fhn_a) * 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;

                // directly add input currents from the external noise synapses to the neuron's activation variable:
                xh[i][j] -=  (1./fhn_eps) * con_upstr_exc * rnoa[i][j] * (xold[i][j] - Vna) * delta_t;
                xh[i][j] -=  (1./fhn_eps) * con_upstr_inh * rnoi[i][j] * (xold[i][j] - Vni) * delta_t;
        }}

        // collect updated lateral synaptic currents
        for (j = 0;j<M;j++ ){
            for (i = 0;i<N;i++ ){
                Ihsyn[i][j] = 0.0;
                for (l = 0;l<KT;l++ ){
                    if(p1[i][j][l]<N && p2[i][j][l]<M)
                        Ihsyn[i][j] +=  conductance_net[i][j][l] * rh[p1[i][j][l]][p2[i][j][l]] * (xh[i][j] - Vsyn);
                }
        }}

        // second integration step
        for (j = 0;j<M;j++ ){
            for (i = 0;i<N;i++ ){
                //FHN neuron: eps * (dx/dt)  =  x -  x^3 / 3 -  y
                xnew[i][j] = xold[i][j] + 0.5 * ((1./fhn_eps) * (xold[i][j] - xold[i][j] * xold[i][j] * xold[i][j]/3. - yold[i][j]  +  xh[i][j] - xh[i][j] * xh[i][j] * xh[i][j]/3. - yh[i][j])) * delta_t;
                xnew[i][j] +=  - (1./fhn_eps) * 0.5 * (Isyn[i][j]  +  Ihsyn[i][j]) * delta_t;
                //FHN neuron: (dy/dt)  =  x + a
                ynew[i][j] = yold[i][j] + 0.5 * (xold[i][j] + fhn_a  +  xh[i][j] + fhn_a) * 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;

                // directly add input currents from the external noise synapses to the neuron's activation variable:
                xnew[i][j] -=  0.5 * (1./fhn_eps) * con_upstr_exc * (rnoa[i][j] * (xold[i][j] - Vna)  +  rnha[i][j] * (xh[i][j] - Vna)) * delta_t;
                xnew[i][j] -=  0.5 * (1./fhn_eps) * con_upstr_inh * (rnoi[i][j] * (xold[i][j] - Vni)  +  rnhi[i][j] * (xh[i][j] - Vni)) * 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(xold[i][j]<xth && xnew[i][j] >= xth){
                    Tl[i][j] = current_time;
                }

                // 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];
        }}

        // write the current network state into the provided output buffer.
        // This is a 3D array of shape (M,N,limit) in row first order,
        // so it has strides proportional to (N * limit, limit, 1).
        // ..and we write only every rec_resolution'th step.
        int limit_rec  =  limit/rec_resolution;
        if (step % rec_resolution  ==  0) {
            int linearind  =  0;
            for(int oi  =  0; oi < M; oi++ ) {
                for(int oj  =  0; oj < N; oj++ ) {
                        linearind  =  oi * limit_rec * N  +  oj * limit_rec  +  (step/rec_resolution) - 1;
                        outputbuffer[linearind]  =  xnew[oj][oi];
                        debug_outbuf[linearind]  =  ynew[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);
}

ext: c

max_stars_repo_name: cknd/synchrony

max_issues_repo_name: cknd/synchrony

max_forks_repo_name: cknd/synchrony


text:
/* Starting from version 7.8, MATLAB BLAS expects ptrdiff_t arguments for integers */
#if MATLAB_VERSION >= 0x0708
#include <stddef.h>
#include <stdlib.h>
#endif
#include <string.h>

/* Define MX_HAS_INTERLEAVED_COMPLEX for version <9.4 */
#ifndef MX_HAS_INTERLEAVED_COMPLEX
#define MX_HAS_INTERLEAVED_COMPLEX 0
#endif

/* Starting from version 7.6, MATLAB BLAS is seperated */
#if MATLAB_VERSION >= 0x0705
#include <blas.h>
#else
#define dcabs1 FORTRAN_WRAPPER(dcabs1)
extern doublereal dcabs1(
        doublereal *z
        );
#endif
#include <lapack.h>
#include "f2c.h"

#define cgeqpb FORTRAN_WRAPPER(cgeqpb)
int cgeqpb(integer *, integer *, integer *, integer *,
        complex *, integer *, complex *, integer *, integer *,
        real *, real *, integer *, real *, complex *, integer *,
        real *, integer *);

#define cgeqpc FORTRAN_WRAPPER(cgeqpc)
int cgeqpc(integer *, integer *, integer *, integer *, complex *,
        integer *, complex *, integer *, integer *, integer *, real *,
        integer *, integer *, complex *, complex *, real *, real *,
        complex *, integer *, real *);

#define cgeqpw FORTRAN_WRAPPER(cgeqpw)
int cgeqpw(integer *, integer *, integer *, integer *, integer *,
        complex *, integer *, integer *, real *, complex *, real *, real *,
        complex *, complex *, real *);

#define cgeqpx FORTRAN_WRAPPER(cgeqpx)
int cgeqpx(integer*, integer*, integer*, integer*, complex*,
        integer*, complex*, integer*, integer*, real*, real*, integer*,
        real*, complex*, integer*, real*, integer*);

#define cgeqpy FORTRAN_WRAPPER(cgeqpy)
int cgeqpy(integer*, integer*, integer*, integer*, complex*,
        integer*, complex*, integer*, integer*, real*, real*, integer*,
        real*, complex*, integer*, real*, integer*);

#define clasmx FORTRAN_WRAPPER(clasmx)
doublereal clasmx(integer *);

#define clauc1 FORTRAN_WRAPPER(clauc1)
logical clauc1(integer *, complex *, real *, complex *, complex *, real *);

#define ctrqpx FORTRAN_WRAPPER(ctrqpx)
int ctrqpx(integer *, integer *, integer *, integer *, complex *,
        integer *, complex *, integer *, integer *, real *, real *,
        integer *, real *, complex *, real *, integer *);

#define ctrqpy FORTRAN_WRAPPER(ctrqpy)
int ctrqpy(integer *, integer *, integer *, integer *,
        complex *, integer *, complex *c, integer *, integer *,
        real *, real *, integer *, real *, complex *, real *, integer *);

#define ctrqxc FORTRAN_WRAPPER(ctrqxc)
int ctrqxc(integer *, integer *, integer *, integer *, complex *,
        integer *, complex *, integer *, integer *, integer *, real *,
        real *, real *, complex *, real *, integer *);

#define cglbif FORTRAN_WRAPPER(cglbif)
int cglbif(integer *, integer *, integer *, integer *,
        complex *, integer *, complex *, integer *, integer *,
        real *, integer *, logical *, complex *, real *, integer *);

#define ccniif FORTRAN_WRAPPER(ccniif)
int ccniif(integer *, integer *, integer *, integer *,
        complex *, integer *, complex *, integer *, integer *,
        real *, real *, integer *, logical *, complex *, real *, integer *);

#define cgret FORTRAN_WRAPPER(cgret)
int cgret(integer *, integer *, integer *, integer *,
       complex *, integer *, complex *, integer *, complex *,
       real *, integer *);

#define shess FORTRAN_WRAPPER(shess)
int chess(integer *, integer *, integer *, integer *,
        complex *, integer *, complex *, integer *,
       complex *, real *, integer *);

#define ctrqyc FORTRAN_WRAPPER(ctrqyc)
int ctrqyc(integer *, integer *, integer *, integer *,
        complex *, integer *, complex *, integer *, integer *,
        integer *, real *, real *, real *, complex *, real *, integer *);

#define ctrrnk FORTRAN_WRAPPER(ctrrnk)
int ctrrnk(integer *, complex *, integer *, real *, integer *, complex *, integer *);

#define dgeqpb FORTRAN_WRAPPER(dgeqpb)
int dgeqpb(integer *, integer *, integer *,
        integer *, doublereal *, integer *, doublereal *, integer *,
        integer *, doublereal *, doublereal *, integer *, doublereal *,
        doublereal *, integer *, integer *);

#define dgeqpc FORTRAN_WRAPPER(dgeqpc)
int dgeqpc(integer *, integer *, integer *,
        integer *, doublereal *, integer *, doublereal *, integer *,
        integer *, integer *, doublereal *, integer *, integer *,
        doublereal *, doublereal *, doublereal *, doublereal *,
        doublereal *, integer *);

#define dgeqpw FORTRAN_WRAPPER(dgeqpw)
int dgeqpw(integer *, integer *, integer *,
        integer *, integer *, doublereal *, integer *, integer *,
        doublereal *, doublereal *, doublereal *, doublereal *,
        doublereal *, doublereal *);

#define dgeqpx FORTRAN_WRAPPER(dgeqpx)
int dgeqpx(integer*, integer*, integer*, integer*, doublereal*,
        integer*, doublereal*, integer*, integer*, doublereal*,
        doublereal*, integer*, doublereal*, doublereal*, integer*,
        integer*);

#define dgeqpy FORTRAN_WRAPPER(dgeqpy)
int dgeqpy(integer*, integer*, integer*, integer*, doublereal*,
        integer*, doublereal*, integer*, integer*, doublereal*,
        doublereal*, integer*, doublereal*, doublereal*, integer*,
        integer*);

#define dlasmx FORTRAN_WRAPPER(dlasmx)
doublereal dlasmx(integer *);

#define dlauc1 FORTRAN_WRAPPER(dlauc1)
logical dlauc1(integer *, doublereal *, doublereal *, doublereal *,
        doublereal *, doublereal *);

#define dtrqpx FORTRAN_WRAPPER(dtrqpx)
int dtrqpx(integer *, integer *, integer *, integer *,
        doublereal *, integer *, doublereal *, integer *, integer *,
        doublereal *, doublereal *, integer *, doublereal *, doublereal *,
        integer *, integer *);

#define dtrqpy FORTRAN_WRAPPER(dtrqpy)
int dtrqpy(integer *, integer *, integer *, integer *,
        doublereal *, integer *, doublereal *, integer *, integer *,
        doublereal *, doublereal *, integer *, doublereal *, doublereal *,
        integer *, integer *);

#define dtrqxc FORTRAN_WRAPPER(dtrqxc)
int dtrqxc(integer *, integer *, integer *, integer *, doublereal *,
        integer*, doublereal *, integer *, integer *, integer *,
        doublereal *, doublereal *, doublereal *, doublereal *, integer *);

#define dcniif FORTRAN_WRAPPER(dcniif)
int dcniif(integer *, integer *, integer *, integer *,
        doublereal *, integer *, doublereal *, integer *, integer *,
        doublereal *, doublereal *, integer *, logical *, doublereal *,
        integer *);

#define dglbif FORTRAN_WRAPPER(dglbif)
int dglbif(integer *, integer *, integer *, integer *,
        doublereal *, integer *, doublereal *, integer *, integer *,
        doublereal *, integer *, logical *, doublereal *, integer *);

#define dgret FORTRAN_WRAPPER(dgret)
int dgret(integer *, integer *, integer *, integer *,
        doublereal *, integer *, doublereal *, integer *,
        doublereal *, integer *);

#define dhess FORTRAN_WRAPPER(dhess)
int dhess(integer *, integer *, integer *, integer *,
       doublereal *, integer *, doublereal *, integer *,
       doublereal *, integer *);

#define dtrqyc FORTRAN_WRAPPER(dtrqyc)
int dtrqyc(integer *, integer *, integer *, integer *, doublereal *,
        integer*, doublereal *, integer *, integer *, integer *,
        doublereal *, doublereal *, doublereal *, doublereal *, integer *);

#define dtrrnk FORTRAN_WRAPPER(dtrrnk)
int dtrrnk(integer *, doublereal *, integer *,
        doublereal *, integer *, doublereal *, integer *);

#define sgeqpb FORTRAN_WRAPPER(sgeqpb)
int sgeqpb(integer *, integer *, integer *, integer *, real *,
        integer *, real *, integer *, integer *, real *, real *, integer *,
        real *, real *, integer *, integer *);

#define sgeqpc FORTRAN_WRAPPER(sgeqpc)
int sgeqpc(integer *, integer *, integer *, integer *, real *,
        integer *, real *, integer *, integer *, integer *, real *,
        integer *, integer *, real *, real *, real *, real *, real *, integer *);

#define sgeqpw FORTRAN_WRAPPER(sgeqpw)
int sgeqpw(integer *, integer *, integer *, integer *,
        integer *, real *, integer *, integer *, real *, real *, real *,
        real *, real *, real *);

#define sgeqpx FORTRAN_WRAPPER(sgeqpx)
int sgeqpx(integer*, integer*, integer*, integer*, real*, integer*,
        real*, integer*, integer*, real*, real*, integer*, real*, real*,
        integer*, integer*);

#define sgeqpy FORTRAN_WRAPPER(sgeqpy)
int sgeqpy(integer*, integer*, integer*, integer*, real*, integer*,
        real*, integer*, integer*, real*, real*, integer*, real*, real*,
        integer*, integer*);

#define slasmx FORTRAN_WRAPPER(slasmx)
doublereal slasmx(integer *);

#define slauc1 FORTRAN_WRAPPER(slauc1)
logical slauc1(integer *, real *, real *, real *, real *, real *);

#define strqpx FORTRAN_WRAPPER(strqpx)
int strqpx(integer *, integer *, integer *, integer *, real *,
        integer *, real *, integer *, integer *, real *, real *, integer *,
        real *, real *, integer *, integer *);

#define strqpy FORTRAN_WRAPPER(strqpy)
int strqpy(integer *, integer *, integer *, integer *, real *,
        integer *, real *, integer *, integer *, real *, real *, integer *,
        real *, real *, integer *, integer *);

#define strqxc FORTRAN_WRAPPER(strqxc)
int strqxc(integer *, integer *, integer *, integer *, real *,
        integer *, real *, integer *, integer *, integer *, real *,
        real *, real *, real *, integer *);

#define sglbif FORTRAN_WRAPPER(sglbif)
int sglbif(integer *, integer *, integer *, integer *,
        real *, integer *, real *, integer *, integer *,
        real *, integer *, logical *, real *, integer *);

#define scniif FORTRAN_WRAPPER(scniif)
int scniif(integer *, integer *, integer *, integer *,
        real *, integer *, real *, integer *, integer *,
        real *, real *, integer *, logical *, real *, integer *);

#define sgret FORTRAN_WRAPPER(sgret)
int sgret(integer *, integer *, integer *, integer *,
        real *, integer *, real *, integer *, real *, integer *);

#define shess FORTRAN_WRAPPER(shess)
int shess(integer *, integer *, integer *, integer *,
        real *, integer *, real *, integer *, real *, integer *);

#define strqyc FORTRAN_WRAPPER(strqyc)
int strqyc(integer *, integer *, integer *, integer *, real *,
        integer *, real *, integer *, integer *, integer *, real *,
        real *, real *, real *, integer *);

#define strrnk FORTRAN_WRAPPER(strrnk)
int strrnk(integer *, real *, integer *, real *,
        integer *, real *, integer *);

#define zgeqpb FORTRAN_WRAPPER(zgeqpb)
int zgeqpb(integer *, integer *, integer *, integer *, doublecomplex *,
        integer *, doublecomplex *, integer *, integer *, doublereal *,
        doublereal *, integer *, doublereal *, doublecomplex *, integer *,
        doublereal *, integer *);

#define zgeqpc FORTRAN_WRAPPER(zgeqpc)
int zgeqpc(integer *, integer *, integer *, integer *, doublecomplex *,
        integer *, doublecomplex *, integer *, integer *, integer *, doublereal *,
        integer *, integer *, doublecomplex *, doublecomplex *, doublereal *, doublereal *,
        doublecomplex *, integer *, doublereal *);

#define zgeqpw FORTRAN_WRAPPER(zgeqpw)
int zgeqpw(integer *, integer *, integer *, integer *, integer *,
        doublecomplex *, integer *, integer *, doublereal *, doublecomplex *,
        doublereal *, doublereal *, doublecomplex *, doublecomplex *, doublereal *);

#define zgeqpx FORTRAN_WRAPPER(zgeqpx)
int zgeqpx(integer*, integer*, integer*, integer*, doublecomplex*,
        integer*, doublecomplex*, integer*, integer*, doublereal*,
        doublereal*, integer*, doublereal*, doublecomplex*, integer*,
        doublereal*, integer*);

#define zgeqpy FORTRAN_WRAPPER(zgeqpy)
int zgeqpy(integer*, integer*, integer*, integer*, doublecomplex*,
        integer*, doublecomplex*, integer*, integer*, doublereal*,
        doublereal*, integer*, doublereal*, doublecomplex*, integer*,
        doublereal*, integer*);

#define zlasmx FORTRAN_WRAPPER(zlasmx)
doublereal zlasmx(integer *i);

#define zlauc1 FORTRAN_WRAPPER(zlauc1)
logical zlauc1(integer *, doublecomplex *, doublereal *,
        doublecomplex *, doublecomplex *, doublereal *);

#define ztrqpx FORTRAN_WRAPPER(ztrqpx)
int ztrqpx(integer *, integer *, integer *, integer *,
        doublecomplex *, integer *, doublecomplex *, integer*,
        integer *, doublereal *, doublereal *, integer *, doublereal *,
        doublecomplex *, doublereal *, integer *);

#define ztrqpy FORTRAN_WRAPPER(ztrqpy)
int ztrqpy(integer *, integer *, integer *, integer *,
        doublecomplex *, integer *, doublecomplex *, integer*,
        integer *, doublereal *, doublereal *, integer *, doublereal *,
        doublecomplex *, doublereal *, integer *);

#define ztrqxc FORTRAN_WRAPPER(ztrqxc)
int ztrqxc(integer *, integer *, integer *, integer *, doublecomplex *,
        integer *, doublecomplex *, integer *, integer *, integer *,
        doublereal *, doublereal *, doublereal *, doublecomplex *,
        doublereal *, integer *);

#define zcniif FORTRAN_WRAPPER(zcniif)
int zcniif(integer *, integer *, integer *, integer *,
        doublecomplex *, integer *, doublecomplex *, integer *, integer *,
        doublereal *, doublereal *, integer *, logical *, doublecomplex *,
        doublereal *, integer *);

#define zglbif FORTRAN_WRAPPER(zglbif)
int zglbif(integer *, integer *, integer *, integer *,
        doublecomplex *, integer *, doublecomplex *, integer *, integer *,
        doublereal *, integer *, logical *, doublecomplex *, doublereal *,
        integer *);

#define zgret FORTRAN_WRAPPER(zgret)
int zgret(integer *, integer *, integer *, integer *,
        doublecomplex *, integer *,doublecomplex *, integer *,
        doublecomplex *, doublereal *, integer *);

#define dzhess FORTRAN_WRAPPER(dhess)
int zhess(integer *, integer *, integer *, integer *,
        doublecomplex *, integer *, doublecomplex *, integer *,
        doublecomplex *, doublereal *, integer *);

#define ztrqyc FORTRAN_WRAPPER(ztrqyc)
int ztrqyc(integer *, integer *, integer *, integer *, doublecomplex *,
        integer *, doublecomplex *, integer *, integer *, integer *,
        doublereal *, doublereal *, doublereal *, doublecomplex *,
        doublereal *, integer *);

#define ztrrnk FORTRAN_WRAPPER(ztrrnk)
int ztrrnk(integer *, doublecomplex *, integer *,
        doublereal *, integer *, doublecomplex *, integer *);

ext: h

max_stars_repo_name: iwoodsawyer/rrqr

max_issues_repo_name: iwoodsawyer/rrqr

max_forks_repo_name: iwoodsawyer/rrqr


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header for the conversion between time scales.
 *
 */

#ifndef _TIMECONVERSION_H
#define _TIMECONVERSION_H

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>

#include "constants.h"

/*************************/
/****** Prototypes ******/

/* Function computing the gmst angle from the gps time */
double gmst_angle_from_gpstime(const double gpstime); /* gpstime in seconds */

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _TIMECONVERSION_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


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

/* Some functions and tests based on test.c by G. Jungman.
*/

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_odeiv.h>
#include "odeiv_util.h"

/* Maximum number of ODE equations */
#define MAXEQ 4

/* RHS for f=2. Solution y = 2 * t + t0 */

int
rhs_linear (double t, const double y[], double f[], void *params)
{
  f[0] = 2.0;

  return GSL_SUCCESS;
}

int
jac_linear (double t, const double y[], double *dfdy, double dfdt[],
            void *params)
{
  dfdy[0] = 0.0;
  dfdt[0] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_lin = {
  rhs_linear,
  jac_linear,
  1,
  0
};

/* RHS for f=y. Equals y=exp(t) with initial value y(0)=1.0 */

int
rhs_exp (double t, const double y[], double f[], void *params)
{
  f[0] = y[0];

  return GSL_SUCCESS;
}

int
jac_exp (double t, const double y[], double *dfdy, double dfdt[],
         void *params)
{
  dfdy[0] = y[0];
  dfdt[0] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_exp = {
  rhs_exp,
  jac_exp,
  1,
  0
};

/* RHS for f0 = -y1, f1 = y0
   equals y = [cos(t), sin(t)] with initial values [1, 0]
*/

int
rhs_sin (double t, const double y[], double f[], void *params)
{
  f[0] = -y[1];
  f[1] = y[0];

  return GSL_SUCCESS;
}

int
jac_sin (double t, const double y[], double *dfdy, double dfdt[],
         void *params)
{
  dfdy[0] = 0.0;
  dfdy[1] = -1.0;
  dfdy[2] = 1.0;
  dfdy[3] = 0.0;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_sin = {
  rhs_sin,
  jac_sin,
  2,
  0
};

/*  Sine/cosine with random failures */

int
rhs_xsin (double t, const double y[], double f[], void *params)
{
  static int n = 0;

  n++;

  if (n > 40 && n < 65) {
    f[0] = GSL_NAN;
    f[1] = GSL_NAN;
    return GSL_EFAILED;
  }

  f[0] = -y[1];
  f[1] = y[0];

  return GSL_SUCCESS;
}

int
jac_xsin (double t, const double y[], double *dfdy, double dfdt[],
         void *params)
{
  static int n = 0;

  n++;

  if (n > 50 && n < 55) {
    dfdy[0] = GSL_NAN;
    dfdy[1] = GSL_NAN;
    dfdy[2] = GSL_NAN;
    dfdy[3] = GSL_NAN;

    dfdt[0] = GSL_NAN;
    dfdt[1] = GSL_NAN;
    return GSL_EFAILED;
  }

  dfdy[0] = 0.0;
  dfdy[1] = -1.0;
  dfdy[2] = 1.0;
  dfdy[3] = 0.0;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_xsin = {
  rhs_xsin,
  jac_xsin,
  2,
  0
};

/* RHS for classic stiff example
   dy0 / dt =  998 * y0 + 1998 * y1    y0(0) = 1.0
   dy1 / dt = -999 * y0 - 1999 * y1    y1(0) = 0.0

   solution is
   y0 = 2 * exp(-t) - exp(-1000 * t)
   y1 = - exp(-t) + exp(-1000 * t)
*/

int
rhs_stiff (double t, const double y[], double f[], void *params)
{
  f[0] = 998.0 * y[0] + 1998.0 * y[1];
  f[1] = -999.0 * y[0] - 1999.0 * y[1];

  return GSL_SUCCESS;
}

int
jac_stiff (double t, const double y[], double *dfdy, double dfdt[],
           void *params)
{
  dfdy[0] = 998.0;
  dfdy[1] = 1998.0;
  dfdy[2] = -999.0;
  dfdy[3] = -1999.0;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_stiff = {
  rhs_stiff,
  jac_stiff,
  2,
  0
};

/* van Der Pol oscillator:
   f0 = y1                           y0(0) = 1.0
   f1 = -y0 + mu * y1 * (1 - y0^2)   y1(0) = 0.0
*/

int
rhs_vanderpol (double t, const double y[], double f[], void *params)
{
  const double mu = 10.0;

  f[0] = y[1];
  f[1] = -y[0] + mu * y[1] * (1.0 - y[0]*y[0]);

  return GSL_SUCCESS;
}

int
jac_vanderpol (double t, const double y[], double *dfdy, double dfdt[],
	       void *params)
{
  const double mu = 10.0;

  dfdy[0] = 0.0;
  dfdy[1] = 1.0;
  dfdy[2] = -2.0 * mu * y[0] * y[1] - 1.0;
  dfdy[3] = mu * (1.0 - y[0] * y[0]);

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_vanderpol = {
  rhs_vanderpol,
  jac_vanderpol,
  2,
  0
};

/* The Oregonator - chemical Belusov-Zhabotinskii reaction
   y0(0) = 1.0, y1(0) = 2.0, y2(0) = 3.0
*/

int
rhs_oregonator (double t, const double y[], double f[], void *params)
{
  const double c1=77.27;
  const double c2=8.375e-6;
  const double c3=0.161;

  f[0] = c1 * (y[1] + y[0] * (1 - c2 * y[0] - y[1]));
  f[1] = 1/c1 * (y[2] - y[1] * (1 + y[0]));
  f[2] = c3 * (y[0] - y[2]);

  return GSL_SUCCESS;
}

int
jac_oregonator (double t, const double y[], double *dfdy, double dfdt[],
		void *params)
{
  const double c1=77.27;
  const double c2=8.375e-6;
  const double c3=0.161;

  dfdy[0] = c1 * (1 - 2 * c2 * y[0] - y[1]);
  dfdy[1] = c1 * (1 - y[0]);
  dfdy[2] = 0.0;

  dfdy[3] = 1/c1 * (-y[1]);
  dfdy[4] = 1/c1 * (-1 - y[0]);
  dfdy[5] = 1/c1;

  dfdy[6] = c3;
  dfdy[7] = 0.0;
  dfdy[8] = -c3;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  dfdt[2] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_oregonator = {
  rhs_oregonator,
  jac_oregonator,
  3,
  0
};

/* Volterra-Lotka predator-prey model

   f0 = (a - b * y1) * y0     y0(0) = 3.0
   f1 = (-c + d * y0) * y1    y1(0) = 1.0
 */

int
rhs_vl (double t, const double y[], double f[], void *params)
{
  const double a = 1.0;
  const double b = 1.0;
  const double c = 1.0;
  const double d = 1.0;

  f[0] = (a - b * y[1]) * y[0];
  f[1] = (-c + d * y[0]) * y[1];

  return GSL_SUCCESS;
}

int
jac_vl (double t, const double y[], double *dfdy, double dfdt[],
            void *params)
{
  const double a = 1.0;
  const double b = 1.0;
  const double c = 1.0;
  const double d = 1.0;

  dfdy[0] = a - b * y[1];
  dfdy[1] = -b * y[0];
  dfdy[2] = d * y[1];
  dfdy[3] = -c + d * y[0];

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_vl = {
  rhs_vl,
  jac_vl,
  2,
  0
};

/* Stiff trigonometric example

   f0 = -50 * (y0 - cos(t))    y0(0) = 0.0
 */

int
rhs_stifftrig (double t, const double y[], double f[], void *params)
{
  f[0] = -50 * (y[0] - cos(t));

  return GSL_SUCCESS;
}

int
jac_stifftrig (double t, const double y[], double *dfdy, double dfdt[],
            void *params)
{
  dfdy[0] = -50;

  dfdt[0] = -50 * sin(t);

  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_stifftrig = {
  rhs_stifftrig,
  jac_stifftrig,
  1,
  0
};

/* E5 - a stiff badly scaled chemical problem by Enright, Hull &
   Lindberg (1975): Comparing numerical methods for stiff systems of
   ODEs. BIT, vol. 15, pp. 10-48.

   f0 = -a * y0 - b * y0 * y2                            y0(0) = 1.76e-3
   f1 = a * y0 - m * c * y1 * y2                         y1(0) = 0.0
   f2 = a * y0 - b * y0 * y2 - m * c * y1 * y2 + c * y3  y2(0) = 0.0
   f3 = b * y0 * y2 - c * y3                             y3(0) = 0.0
 */

int
rhs_e5 (double t, const double y[], double f[], void *params)
{
  const double a = 7.89e-10;
  const double b = 1.1e7;
  const double c = 1.13e3;
  const double m = 1.0e6;

  f[0] = -a * y[0] - b * y[0] * y[2];
  f[1] = a * y[0] - m * c * y[1] * y[2];
  f[3] = b * y[0] * y[2] - c * y[3];
  f[2] = f[1] - f[3];

  return GSL_SUCCESS;
}

int
jac_e5 (double t, const double y[], double *dfdy, double dfdt[],
            void *params)
{
  const double a = 7.89e-10;
  const double b = 1.1e7;
  const double c = 1.13e3;
  const double m = 1.0e6;

  dfdy[0] = -a - b * y[2];
  dfdy[1] = 0.0;
  dfdy[2] = -b * y[0];
  dfdy[3] = 0.0;

  dfdy[4] = a;
  dfdy[5] = -m * c * y[2];
  dfdy[6] = -m * c * y[1];
  dfdy[7] = 0.0;

  dfdy[8] = a - b * y[2];
  dfdy[9] = -m * c * y[2];
  dfdy[10] = -b * y[0] - m * c * y[1];
  dfdy[11] = c;

  dfdy[12] = b * y[2];
  dfdy[13] = 0.0;
  dfdy[14] = b * y[0];
  dfdy[15] = -c;

  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  dfdt[2] = 0.0;
  dfdt[3] = 0.0;

  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_e5 = {
  rhs_e5,
  jac_e5,
  4,
  0
};

void
test_odeiv_stepper (const gsl_odeiv_step_type *T, const gsl_odeiv_system *sys,
		    const double h, const double t, const char desc[],
		    const double ystart[], const double yfin[],
		    const double relerr)
{
  /* tests stepper T with one fixed length step advance of system sys
     and compares with given values yfin
  */

  double y[MAXEQ] = {0.0};
  double yerr[MAXEQ] = {0.0};
  size_t ne = sys->dimension;
  size_t i;

  gsl_odeiv_step *step = gsl_odeiv_step_alloc (T, ne);

  DBL_MEMCPY (y, ystart, MAXEQ);

  {
    int s = gsl_odeiv_step_apply (step, t, h, y, yerr, 0, 0, sys);
    if (s != GSL_SUCCESS)
      {
	gsl_test(s, "test_odeiv_stepper: %s step_apply returned %d", desc, s);
      }
  }

  for (i = 0; i < ne; i++)
    {
      gsl_test_rel (y[i], yfin[i], relerr,
		    "%s %s step(%d)",
		    gsl_odeiv_step_name (step), desc,i);
    }

  gsl_odeiv_step_free (step);
}

void
test_stepper (const gsl_odeiv_step_type *T)
{
  /* Tests stepper T with a step of selected systems */

  double y[MAXEQ] = {0.0};
  double yfin[MAXEQ] = {0.0};

  /* Step length */
  double h;

  /* Required tolerance */
  double err_target;

  /* linear */
  h = 1e-1;
  err_target = 1e-10;
  y[0] = 0.58;
  yfin[0] = y[0] + 2 * h;
  test_odeiv_stepper (T, &rhs_func_lin, h, 0.0, "linear",
		      y, yfin, err_target);

  /* exponential */
  h = 1e-4;
  err_target = 1e-8;
  y[0] = exp(2.7);
  yfin[0] = exp(2.7 + h);
  test_odeiv_stepper (T, &rhs_func_exp, h, 2.7, "exponential",
		      y, yfin, err_target);
  /* cosine-sine */
  h = 1e-3;
  err_target = 1e-6;
  y[0] = cos(1.2);
  y[1] = sin(1.2);
  yfin[0] = cos(1.2 + h);
  yfin[1] = sin(1.2 + h);
  test_odeiv_stepper (T, &rhs_func_sin, h, 1.2, "cosine-sine",
		      y, yfin, err_target);

  /* classic stiff */
  h = 1e-7;
  err_target = 1e-4;
  y[0] = 1.0;
  y[1] = 0.0;

  {
    const double e1 = exp (-h);
    const double e2 = exp (-1000.0 * h);
    yfin[0] = 2.0 * e1 - e2;
    yfin[1] = -e1 + e2;
  }

  test_odeiv_stepper (T, &rhs_func_stiff, h, 0.0, "classic_stiff",
		      y, yfin, err_target);
}

void
test_evolve_system (const gsl_odeiv_step_type * T,
                    const gsl_odeiv_system * sys,
                    double t0, double t1, double hstart,
                    double y[], double yfin[],
                    double err_target, const char *desc)
{
  /* Tests system sys with stepper T. Step length is controlled by
     error estimation from the stepper.
  */

  int steps = 0;
  size_t i;

  double t = t0;
  double h = hstart;

  /* Tolerance factor in testing errors */
  const double factor = 10;

  gsl_odeiv_step * step = gsl_odeiv_step_alloc (T, sys->dimension);

  gsl_odeiv_control *c =
    gsl_odeiv_control_standard_new (err_target, err_target, 1.0, 0.0);

  gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension);

  while (t < t1)
    {
      int s = gsl_odeiv_evolve_apply (e, c, step, sys, &t, t1, &h, y);

      if (s != GSL_SUCCESS && sys != &rhs_func_xsin)
	{
	  gsl_test(s, "%s evolve_apply returned %d",
		   gsl_odeiv_step_name (step), s);
	  break;
	}

      if (steps > 100000)
	{
	  gsl_test(GSL_EFAILED,
		   "%s evolve_apply reached maxiter",
		   gsl_odeiv_step_name (step));
	  break;
	}

      steps++;
    }

  /* err_target is target error of one step. Test if stepper has made
     larger error than (tolerance factor times) the number of steps
     times the err_target */

  for (i = 0; i < sys->dimension; i++)
    {
      gsl_test_abs (y[i], yfin[i], factor * e->count * err_target,
		    "%s %s evolve(%d)",
		    gsl_odeiv_step_name (step), desc, i);
    }

  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (step);
}

int
sys_driver (const gsl_odeiv_step_type * T,
	    const gsl_odeiv_system * sys,
	    double t0, double t1, double hstart,
	    double y[], double epsabs, double epsrel,
	    const char desc[])
{
  /* This function evolves a system sys with stepper T from t0 to t1.
     Step length is varied via error control with possibly different
     absolute and relative error tolerances.
  */

  int s = 0;
  int steps = 0;

  double t = t0;
  double h = hstart;

  gsl_odeiv_step * step = gsl_odeiv_step_alloc (T, sys->dimension);

  gsl_odeiv_control *c =
    gsl_odeiv_control_standard_new (epsabs, epsrel, 1.0, 0.0);
  gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension);

  while (t < t1)
    {
      s = gsl_odeiv_evolve_apply (e, c, step, sys, &t, t1, &h, y);

      if (s != GSL_SUCCESS)
	{
	  gsl_test(s, "sys_driver: %s evolve_apply returned %d",
		   gsl_odeiv_step_name (step), s);
	  break;
	}

      if (steps > 1e7)
	{
	  gsl_test(GSL_EMAXITER,
		   "sys_driver: %s evolve_apply reached maxiter at t=%g",
		   gsl_odeiv_step_name (step), t);
	  s = GSL_EMAXITER;
	  break;
	}

      steps++;
    }

  gsl_test(s, "%s %s [%g,%g], %d steps completed",
	   gsl_odeiv_step_name (step), desc, t0, t1, steps);

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

  return s;
}

void
test_compare_vanderpol (void)
{
  /* Compares output of van Der Pol oscillator with several steppers */

  /* system dimension */
  const size_t sd = 2;

  const gsl_odeiv_step_type *steppers[20];
  const gsl_odeiv_step_type **T;

  /* Required error tolerance for each stepper. */
  double err_target[20];

  /* number of ODE solvers */
  const size_t ns = 11;

  /* initial values for each ode-solver */
  double y[11][2];
  double *yp = &y[0][0];

  size_t i, j, k;
  int status = 0;

  /* Parameters for the problem and stepper  */
  const double start = 0.0;
  const double end = 100.0;
  const double epsabs = 1e-8;
  const double epsrel = 1e-8;
  const double initstepsize = 1e-5;

  /* Initialize */

  steppers[0] = gsl_odeiv_step_rk2;
  err_target[0] = 1e-6;
  steppers[1] = gsl_odeiv_step_rk4;
  err_target[1] = 1e-6;
  steppers[2] = gsl_odeiv_step_rkf45;
  err_target[2] = 1e-6;
  steppers[3] = gsl_odeiv_step_rkck;
  err_target[3] = 1e-6;
  steppers[4] = gsl_odeiv_step_rk8pd;
  err_target[4] = 1e-6;
  steppers[5] = gsl_odeiv_step_rk2imp;
  err_target[5] = 1e-5;
  steppers[6] = gsl_odeiv_step_rk2simp;
  err_target[6] = 1e-5;
  steppers[7] = gsl_odeiv_step_rk4imp;
  err_target[7] = 1e-6;
  steppers[8] = gsl_odeiv_step_bsimp;
  err_target[8] = 1e-7;
  steppers[9] = gsl_odeiv_step_gear1;
  err_target[9] = 1e-2;
  steppers[10] = gsl_odeiv_step_gear2;
  err_target[10] = 1e-6;
  steppers[11] = 0;

  T = steppers;

  for (i = 0; i < ns; i++)
    {
      y[i][0] = 1.0;
      y[i][1] = 0.0;
    }

  /* Call each solver for the problem */

  i = 0;
  while (*T != 0)
    {
      {
	int s = sys_driver (*T, &rhs_func_vanderpol,
			    start, end, initstepsize, &yp[i],
			    epsabs, epsrel, "vanderpol");
	if (s != GSL_SUCCESS)
	  {
	    status++;
	  }
      }

      T++;
      i += sd;
    }

  if (status != GSL_SUCCESS)
    {
      return;
    }

  /* Compare results */

  T = steppers;

  for (i = 0; i < ns; i++)
    for (j = i+1; j < ns; j++)
      for (k = 0; k < sd; k++)
	{
	  const double val1 = yp[sd * i + k];
	  const double val2 = yp[sd * j + k];
	  gsl_test_abs (val1, val2,
			( GSL_MAX(err_target[i], err_target[j]) ),
			"%s/%s vanderpol",
			T[i]->name, T[j]->name);
	}

}

void
test_compare_oregonator (void)
{
  /* Compares output of the Oregonator with several steppers */

  /* system dimension */
  const size_t sd = 3;

  const gsl_odeiv_step_type *steppers[20];
  const gsl_odeiv_step_type **T;

  /* Required error tolerance for each stepper. */
  double err_target[20];

  /* number of ODE solvers */
  const size_t ns = 2;

  /* initial values for each ode-solver */
  double y[2][3];
  double *yp = &y[0][0];

  size_t i, j, k;
  int status = 0;

  /* Parameters for the problem and stepper  */
  const double start = 0.0;
  const double end = 360.0;
  const double epsabs = 1e-8;
  const double epsrel = 1e-8;
  const double initstepsize = 1e-5;

  /* Initialize */

  steppers[0] = gsl_odeiv_step_rk2simp;
  err_target[0] = 1e-6;
  steppers[1] = gsl_odeiv_step_bsimp;
  err_target[1] = 1e-6;
  steppers[2] = 0;

  T = steppers;

  for (i = 0; i < ns; i++)
    {
      y[i][0] = 1.0;
      y[i][1] = 2.0;
      y[i][2] = 3.0;
    }

  /* Call each solver for the problem */

  i = 0;
  while (*T != 0)
    {
      {
	int s = sys_driver (*T, &rhs_func_oregonator,
			    start, end, initstepsize, &yp[i],
			    epsabs, epsrel, "oregonator");

	if (s != GSL_SUCCESS)
	  {
	    status++;
	  }
      }

      T++;
      i += sd;
    }

  if (status != GSL_SUCCESS)
    {
      return;
    }

  /* Compare results */

  T = steppers;

  for (i = 0; i < ns; i++)
    for (j = i+1; j < ns; j++)
      for (k = 0; k < sd; k++)
	{
	  const double val1 = yp[sd * i + k];
	  const double val2 = yp[sd * j + k];
	  gsl_test_rel (val1, val2,
			( GSL_MAX(err_target[i], err_target[j]) ),
			"%s/%s oregonator",
			T[i]->name, T[j]->name);
	}

}

void
test_evolve_linear (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[1];
  double yfin[1];

  y[0] = 1.0;
  yfin[0] = 9.0;
  test_evolve_system (T, &rhs_func_lin, 0.0, 4.0, h, y, yfin, err,
		      "linear[0,4]");
}

void
test_evolve_exp (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[1];
  double yfin[1];

  y[0] = 1.0;
  yfin[0] = exp (2.0);
  test_evolve_system (T, &rhs_func_exp, 0.0, 2.0, h, y, yfin, err,
		      "exp[0,2]");
}

void
test_evolve_sin (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  yfin[0] = cos (2.0);
  yfin[1] = sin (2.0);
  test_evolve_system (T, &rhs_func_sin, 0.0, 2.0, h, y, yfin, err,
		      "sine[0,2]");
}

void
test_evolve_xsin (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  yfin[0] = cos (2.0);
  yfin[1] = sin (2.0);
  test_evolve_system (T, &rhs_func_xsin, 0.0, 2.0, h, y, yfin, err,
		      "sine[0,2] w/errors");
}

void
test_evolve_stiff1 (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  {
    double arg = 1.0;
    double e1 = exp (-arg);
    double e2 = exp (-1000.0 * arg);
    yfin[0] = 2.0 * e1 - e2;
    yfin[1] = -e1 + e2;
  }
  test_evolve_system (T, &rhs_func_stiff, 0.0, 1.0, h, y, yfin, err,
                      "stiff[0,1]");
}

void
test_evolve_stiff5 (const gsl_odeiv_step_type * T, double h, double err)
{
  double y[2];
  double yfin[2];

  y[0] = 1.0;
  y[1] = 0.0;
  {
    double arg = 5.0;
    double e1 = exp (-arg);
    double e2 = exp (-1000.0 * arg);
    yfin[0] = 2.0 * e1 - e2;
    yfin[1] = -e1 + e2;
  }
  test_evolve_system (T, &rhs_func_stiff, 0.0, 5.0, h, y, yfin, err,
                      "stiff[0,5]");
}

int
main (void)
{
  int i;

  struct ptype
  {
    const gsl_odeiv_step_type *type;
    double h;
  }
  p[20];

  p[0].type = gsl_odeiv_step_rk2;
  p[0].h = 1.0e-3;
  p[1].type = gsl_odeiv_step_rk4;
  p[1].h = 1.0e-3;
  p[2].type = gsl_odeiv_step_rkf45;
  p[2].h = 1.0e-3;
  p[3].type = gsl_odeiv_step_rkck;
  p[3].h = 1.0e-3;
  p[4].type = gsl_odeiv_step_rk8pd;
  p[4].h = 1.0e-3;
  p[5].type = gsl_odeiv_step_rk2imp;
  p[5].h = 1.0e-3;
  p[6].type = gsl_odeiv_step_rk2simp;
  p[6].h = 1.0e-3;
  p[7].type = gsl_odeiv_step_rk4imp;
  p[7].h = 1.0e-3;
  p[8].type = gsl_odeiv_step_bsimp;
  p[8].h = 1.0e-3;
  p[9].type = gsl_odeiv_step_gear1;
  p[9].h = 1.0e-3;
  p[10].type = gsl_odeiv_step_gear2;
  p[10].h = 1.0e-3;
  p[11].type = 0;

  gsl_ieee_env_setup ();

  for (i = 0; p[i].type != 0; i++)
    {
      test_stepper(p[i].type);
    }

  for (i = 0; p[i].type != 0; i++)
    {
      test_evolve_linear (p[i].type, p[i].h, 1e-10);
      test_evolve_exp (p[i].type, p[i].h, 1e-6);
      test_evolve_sin (p[i].type, p[i].h, 1e-8);
      test_evolve_xsin (p[i].type, p[i].h, 1e-8);
      test_evolve_stiff1 (p[i].type, p[i].h, 1e-7);
      test_evolve_stiff5 (p[i].type, p[i].h, 1e-7);
    }

  test_compare_vanderpol();
  test_compare_oregonator();

  exit (gsl_test_summary ());
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * 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.
 */

#pragma once

#ifndef widget_c3b57c97_cea3_4794_b683_762a7b11af49_h
#define widget_c3b57c97_cea3_4794_b683_762a7b11af49_h

#include <gslib/std.h>
#include <gslib/dvt.h>
#include <gslib/uuid.h>
#include <ariel/sysop.h>
#include <ariel/painter.h>

__ariel_begin__

class wsys_manager;

class widget:
    public dvt_holder
{
    friend class wsys_manager;

public:
    widget(wsys_manager* m);
    virtual ~widget();
    virtual void* query_interface(const uuid& uid);
    virtual bool create(widget* ptr, const gchar* name, const rect& rc, uint style);
    virtual void close();
    virtual void show(bool b);
    virtual void enable(bool b);
    virtual void move(const rect& rc);
    virtual void refresh(const rect& rc, bool imm = false);
    virtual void draw(painter* paint) {}
    virtual widget* capture(bool b);
    virtual widget* focus();
    virtual bool hit_test(const point& pt) { return (_style & sm_hitable) != 0; }

public:
    enum laytag
    {
        lay_before,
        lay_after,
        lay_first,
        lay_last,
    };
    virtual void lay(widget* ptr, laytag t);

public:
    virtual void on_press(uint um, unikey uk, const point& pt);
    virtual void on_click(uint um, unikey uk, const point& pt);
    virtual void on_hover(uint um, const point& pt);
    virtual void on_leave(uint um, const point& pt);
    virtual void on_keydown(uint um, unikey uk) {}
    virtual void on_keyup(uint um, unikey uk) {}
    virtual void on_char(uint um, uint ch) {}
    virtual void on_caret() {}
    virtual void on_capture(bool b) {}
    virtual void on_focus(bool b) {}
    virtual void on_accelerator(unikey key, uint mask) {}

protected:
    wsys_manager*   _manager;
    string          _name;
    uint            _style;
    rect            _pos;
    bool            _visible;
    bool            _enabled;

public:
    const string& get_name() const { return _name; }
    const rect& get_rect() const { return _pos; }
    rectf get_rectf() const { return to_rectf(get_rect()); }
    bool is_visible() const { return _visible && (_style & sm_visible); }
    bool is_enabled() const { return _enabled; }
    bool is_focused() const;
    void hide() { show(false); }
    void disable() { enable(false); }
    point& to_global(point& pt) const;
    rect& to_global(rect& rc) const;
    point& to_local(point& pt) const;
    rect& to_local(rect& rc) const;
    void move(const point& pt);
    void resize(int w, int h);
    void refresh(bool imm);
    int get_width() const { return _pos.width(); }
    int get_height() const { return _pos.height(); }
    wsys_manager* get_manager() const { return _manager; }
    bool register_accelerator(unikey key, uint mask);
    widget* unregister_accelerator(unikey key, uint mask);

protected:
    widget*         _prev;
    widget*         _next;
    widget*         _child;
    widget*         _last_child;
    widget*         _parent;

public:
    widget* get_parent() const { return _parent; }
    widget* get_prev() const { return _prev; }
    widget* get_next() const { return _next; }
    widget* get_child() const { return _child; }
    widget* get_last_child() const { return _last_child; }

public:
    template<class _lamb>
    static int traverse_widget(widget* w, _lamb trav_to_stop)
    {
        assert(w);
        for(auto* p = w; p; p = p->get_next()) {
            if(auto r = trav_to_stop(p))
                return r;
        }
        return 0;
    }
    template<class _lamb>
    static int traverse_widget_reversed(widget* w, _lamb trav_to_stop)
    {
        assert(w);
        for(auto* p = w; p; p = p->get_prev()) {
            if(auto r = trav_to_stop(p))
                return r;
        }
        return 0;
    }
    template<class _lamb>
    int traverse_child_widget(_lamb trav_to_stop) { return traverse_widget(_child, trav_to_stop); }
    template<class _lamb>
    int traverse_child_widget_reversed(_lamb trav_to_stop) { return traverse_widget_reversed(_last_child, trav_to_stop); }
};

class button:
    public widget
{
public:
    typedef widget superref;
    typedef button self;

public:
    button(wsys_manager* m);
    virtual ~button();
    virtual void draw(painter* paint) override;
    virtual void enable(bool b) override;
    virtual void on_press(uint um, unikey uk, const point& pt) override;
    virtual void on_click(uint um, unikey uk, const point& pt) override;
    virtual void on_hover(uint um, const point& pt) override;
    virtual void on_leave(uint um, const point& pt) override;

protected:
    virtual void set_press();
    virtual void set_normal();
    virtual void set_hover();
    virtual void set_gray();

protected:
    texture2d*      _source;
    texture2d*      _bkground;
    bool            _4states;

public:
    void set_image(texture2d* img, bool fs = false);

    enum btnstate
    {
        bs_none     = 0,
        bs_press,
        bs_hover,
        bs_normal,
    };
    void set_btnstate(btnstate bs) { _btnstate = bs; }
    btnstate get_btnstate() const { return _btnstate; }

protected:
    btnstate        _btnstate;
};

class edit_line:
    public widget
{
public:
    typedef widget superref;
    typedef edit_line self;

public:
    edit_line(wsys_manager* m);
    virtual void draw(painter* paint) override;
    virtual void on_press(uint um, unikey uk, const point& pt) override;
    virtual void on_click(uint um, unikey uk, const point& pt) override;
    virtual void on_hover(uint um, const point& pt) override;
    virtual void on_char(uint um, uint ch) override;
    virtual void on_keydown(uint um, unikey uk) override;
    virtual void on_caret() override;
    virtual void on_focus(bool b) override;

protected:
    bool            _caret_on;
    int             _caretpos;
    string          _textbuf;
    int             _sel_start;
    int             _sel_end;

public:
    virtual void set_text(const gchar* str);
    virtual void set_caret(int n);
    virtual void set_select(int start, int end);
    virtual void replace_select(const gchar* str);

protected:
    virtual void draw_background(painter* paint);
    virtual void draw_normal_text(painter* paint);
    virtual void draw_select_text(painter* paint);
    virtual void draw_caret(painter* paint);

public:
    void set_font(const font& ft) { _font = ft; }
    void set_bkground(texture2d* ptr) { _bkground = ptr; }
    void set_text_color(color cr) { _txtcolor = cr; }
    void set_select_color(color cr) { _selcolor = cr; }
    void set_caret_color(color cr) { _crtcolor = cr; }
    const gchar* get_text() const { return _textbuf.c_str(); }
    void del_select();
    int hit_char(point pt);
    int prev_char(int pos);
    int next_char(int pos);
    int prev_logic_char(int pos);
    int next_logic_char(int pos);
    bool no_select() const { return _sel_start == -1 || _sel_start == _sel_end; }

protected:
    texture2d*      _bkground;
    color           _txtcolor;
    color           _selcolor;
    color           _crtcolor;
    font            _font;

protected:
    int trim_if_overrun(bool alarm);
};

typedef system_driver wsys_driver;
typedef system_notify wsys_notify;

class timer:
    public dvt_holder
{
public:
    timer(wsys_manager* mgr);
    void start(int elapse);
    void start_single(int elapse);
    void set_user_data(uint usd) { _userdata = usd; }
    uint get_user_data() const { return _userdata; }

public:
    virtual ~timer();
    virtual void on_timer(uint tid) {}
    virtual void on_notified();             /* CANNOT be used by connect notify */

protected:
    wsys_driver*    _driver;
    uint            _userdata;
    int             _elapse;
    bool            _single;

private:
    void initialize(wsys_driver* drv);
};

struct accel_key
{
    unikey          key;
    uint            mask;

public:
    accel_key(): key((unikey)-1), mask(0) {}        /* invalid */
    accel_key(unikey k, uint m): key(k), mask(m) {}
    bool is_valid() const { return key != (unikey)-1; }
    bool operator == (const accel_key& that) const { return key == that.key && mask == that.mask; }
    bool from_string(const string& str);
    const string& to_string(string& str) const;
};

class wsys_manager:
    public wsys_notify
{
public:
    wsys_manager();
    void set_wsysdrv(wsys_driver* drv);
    void set_painter(painter* paint);
    painter* get_painter() const { return _painter; }
    wsys_driver* get_driver() const { return _driver; }
    void initialize(const rect& rc);

protected:
    wsys_driver*    _driver;
    painter*        _painter;
    dirty_list      _dirty;
    int             _width;
    int             _height;
    timer*          _caret;

public:
    void set_dimension(int w, int h);
    int get_width() const { return _width; }
    int get_height() const { return _height; }
    void refresh() { _dirty.set_whole(); }
    void refresh(const rect& rc, bool imm = false);
    void update();
    void update(widget* w);
    widget* hit_test(point pt) { return hit_test(_root, pt); }
    widget* hit_test(widget* ptr, point pt);
    widget* hit_proc(const point& pt, point& pt1);
    widget* set_capture(widget* ptr, bool b);
    widget* set_focus(widget* ptr);
    void on_caret(uint);
    widget* get_capture() const { return _capture; }
    widget* get_focus() const { return _focus; }
    void set_cursor(cursor_type curty);
    void reset_cursor() { set_cursor(cur_arrow); }

public:
    virtual ~wsys_manager();
    virtual void on_show(bool b) override;
    virtual void on_create(wsys_driver* ptr, const rect& rc) override;
    virtual void on_close() override;
    virtual void on_resize(const rect& rc) override;
    virtual void on_paint(const rect& rc) override;
    virtual void on_halt() override;
    virtual void on_resume() override;
    virtual bool on_mouse_down(uint um, unikey uk, const point& pt) override;
    virtual bool on_mouse_up(uint um, unikey uk, const point& pt) override;
    virtual bool on_mouse_move(uint um, const point& pt) override;
    virtual bool on_key_down(uint um, unikey uk) override;
    virtual bool on_key_up(uint um, unikey uk) override;
    virtual bool on_char(uint um, uint ch) override;
    virtual void on_timer(uint tid) override;

protected:
    widget*         _hover;
    widget*         _capture;
    widget*         _focus;
    widget*         _root;

private:
    typedef unordered_map<string, widget*> widget_map;
    widget_map      _widget_map;

public:
    template<class _ctor>
    _ctor* add_widget(widget* ptr, const gchar* name, const rect& rc, uint style)
    {
        if(name && name[0] && _widget_map.find(name) != _widget_map.end())
            return nullptr;
        if(!ptr && _root)
            ptr = _root;
        _ctor* p = new _ctor(this);
        assert(p);
        if(!p->create(ptr, name ? name : _t(""), rc, style)) {
            delete p;
            return nullptr;
        }
        if(name && name[0]) {
            _widget_map.insert(
                std::make_pair(name, p)
                );
        }
        if(!ptr && !_root)
            _root = p;
        return p;
    }
    widget* get_root() const { return _root; }
    widget* find_widget(const string& name);
    bool remove_widget(widget* ptr);
    bool remove_widget(const string& name);

protected:
    bool remove_widget_internal(widget* ptr);

public:
    void set_ime(widget* ptr, point pt, const font& ft);
    void set_clipboard(clipfmt fmt, const void* ptr, int size);
    int get_clipboard(clipfmt fmt, const void*& ptr);
    int get_clipboard(clipboard_list& cl, int c);

protected:
    struct accel_key_hasher { size_t operator()(const accel_key& kval) const { return hash_bytes((const byte*)&kval, sizeof(kval)); } };
    typedef unordered_map<accel_key, widget*, accel_key_hasher> accel_map;
    friend bool try_proceed_accelerator(const accel_map&, unikey, uint);
    accel_map       _accel_map;

public:
    bool register_accelerator(widget* w, unikey key, uint mask);
    widget* unregister_accelerator(unikey key, uint mask);
};

extern const uuid uuid_widget;

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#pragma once
#include "Type/peColor.h"
#include "Type\peBitmask.h"
#include "peCoreDefs.h"

#include <glm/detail/type_vec2.hpp>
#undef max
#include <gsl.h>
#include <random>

#pragma warning(push)
#pragma warning(disable : 4251)

/*
 * Everything in here is taking more or less one-to-one from 'Physically based
 * rendering' (Pharr, Humphreys, 2010)
 */

namespace pe {
struct Sample;
class peCoordSys;

//! \brief Represents a color spectrum. We will use RGB with floating-point
//! precision here
using Spectrum_t = RGB_32BitFloat;

//! \brief Types of BRDFs and BTDFs
enum class BxDFType {
  Reflection = 1 << 0,
  Transmission = 1 << 1,
  Diffuse = 1 << 2,
  Glossy = 1 << 3,
  Specular = 1 << 4,
  AllTypes = Diffuse | Glossy | Specular,
  AllReflection = Reflection | AllTypes,
  AllTransmission = Transmission | AllTypes,
  All = AllReflection | AllTransmission
};

template <> struct EnableEnumBitmask<BxDFType> : std::true_type {};

//! \brief BRDF or BTDF base class
class PE_CORE_API peBxDF {
public:
  virtual ~peBxDF() = default;
  explicit peBxDF(BxDFType type);

  bool HasFlags(BxDFType flags) const;

  auto Type() const { return _type; }

  //! \brief Evaluate this BxDF for the given outgoing and incoming vectors
  //! \param wo Outgoing vector
  //! \param wi Incoming vector
  //! \returns Value of the distribution function for the two vectors
  virtual Spectrum_t Eval(const glm::vec3 &wo, const glm::vec3 &wi) const = 0;

  //! \brief Sampling function for BxDF distributions that utilize delta
  //! distributions. Here, the incident direction wi cannot
  //!        be known to the user, hence the BxDF computes it
  //! \param wo Outgoing vector
  //! \param wi Incoming vector, will be computed by the BxDF
  //! \param rnd1 First uniform random variable
  //! \param rnd2 Second uniform random variable
  //! \param pdf Probability density function for this BxDF and the given pair
  //! of vectors \returns Value of the distribution function for the two vectors
  virtual Spectrum_t Sample_f(const glm::vec3 &wo, glm::vec3 &wi,
                              const float rnd1, const float rnd2,
                              float &pdf) const;

  //! \brief Computes the hemishperical-directional reflectance, which is the
  //! total reflection in the given direction due to constant illumination over
  //! the hemisphere. Not every BxDF will be able to compute this in closed
  //! form, so they will use something like a monte carlo algorithm to compute
  //! it, hence the <paramref name="samples"/> member
  //! \param wo Outgoing vector
  //! \param samples Samples
  //! \returns Hemispherical reflection around <paramref name="wo"/>
  virtual Spectrum_t rho(const glm::vec3 &wo,
                         const gsl::span<glm::vec2> &samples) const;

  //! \brief Computes the hemispherical-hemispherical reflectance, which is the
  //! fraction of incident light reflected by the surface when the incident
  //! light is the same from all directions. Not every BxDF will be able to
  //! compute this in closed form, so they will use something like a monte carlo
  //! algorithm to compute it
  //! \param samples1 First set of samples for monte carlo method
  //! \param samples2 Second set of samples for monte carlo method
  //! \returns Direction-independent reflectance
  virtual Spectrum_t rho(const gsl::span<glm::vec2> &samples1,
                         const gsl::span<glm::vec2> &samples2) const;

  //! \brief Returns the probability density function for the given pairs of
  //! vectors \param wo Outgoing vector in shading space \param wi Incoming
  //! vector in shading space \returns PDF
  virtual float Pdf(const glm::vec3 &wo, const glm::vec3 &wi) const;

private:
  const BxDFType _type;
};

#pragma region Fresnel

//! \brief Helper class to encapsulate fresnel reflectance
struct PE_CORE_API peFresnel {
  virtual ~peFresnel() {}
  virtual Spectrum_t Eval(float cosi) const = 0;
};

//! \brief Fresnel reflectance for a conductor
class PE_CORE_API peFresnelConductor : public peFresnel {
public:
  //! \brief Initializes this FresnelConductor structure with the given index of
  //! refraction and absorption
  //! \param eta Index of refraction of conductor
  //! \param k Absorption index of conductor
  peFresnelConductor(const Spectrum_t &eta, const Spectrum_t &k);
  Spectrum_t Eval(float cosi) const override;

private:
  const Spectrum_t _eta;
  const Spectrum_t _k;
};

//! \brief Fresnel reflectance for dielectric material
class PE_CORE_API peFresnelDielectric : public peFresnel {
public:
  //! \brief Initializes this FresnelDielectric structure with the given
  //! incident and transmitted medium indices of refraction
  //! \param etaIncident Index of refraction for incident material
  //! \param etaTransmitted Index of refraction for transmitted material
  peFresnelDielectric(float etaIncident, float etaTransmitted);
  Spectrum_t Eval(float cosi) const override;

private:
  const float _etaIndicent, _etaTransmitted;
};

#pragma endregion

//! \brief BRDF for specular reflection
class PE_CORE_API peSpecularReflection : public peBxDF {
public:
  peSpecularReflection(const Spectrum_t &spectrum, const peFresnel &fresnel);

  Spectrum_t Eval(const glm::vec3 &wo, const glm::vec3 &wi) const override;
  Spectrum_t Sample_f(const glm::vec3 &wo, glm::vec3 &wi, const float rnd1,
                      const float rnd2, float &pdf) const override;

private:
  const Spectrum_t _color;
  const peFresnel &_fresnel;
};

//! \brief Lambertian diffuse reflection
class PE_CORE_API peLambert : public peBxDF {
public:
  explicit peLambert(const Spectrum_t &color);
  Spectrum_t Eval(const glm::vec3 &wo, const glm::vec3 &wi) const override;

  Spectrum_t rho(const glm::vec3 &wo,
                 const gsl::span<glm::vec2> &samples) const override;
  Spectrum_t rho(const gsl::span<glm::vec2> &samples1,
                 const gsl::span<glm::vec2> &samples2) const override;

private:
  const Spectrum_t _color;
};

//! \brief Random values for BSDF sampling
struct PE_CORE_API BSDFSample {
  BSDFSample() = default;
  BSDFSample(const glm::vec2 &dir, float component);
  template <typename Rnd> explicit BSDFSample(Rnd &rng) {
    std::uniform_real_distribution<float> dist{0.f, 1.f};
    dir = {dist(rng), dist(rng)};
    component = dist(rng);
  }

  glm::vec2 dir;
  float component;
};

//! \brief Bidirectional scattering distribution function. Determines surface
//! properties of an object
class PE_CORE_API BSDF {
public:
  BSDF();

  BSDF(const BSDF &other);
  BSDF(BSDF &&other) noexcept;

  BSDF &operator=(const BSDF &);
  BSDF &operator=(BSDF &&) noexcept;

  void Add(peBxDF const *bxdf);

  uint32_t NumBxDFs() const;
  uint32_t NumBxDFsWithFlags(BxDFType flags) const;

  //! \brief Evaluate the BSDF for the given set of incoming and outgoing
  //! vectors \param outgoingWorld Outgoing vector in world space \param
  //! incomingWorld Incoming vector in world space \param shadingCoordSys
  //! Shading coordinate system \param flags Types of BxDFs to sample \param
  //! geometricNormal The normal of the actual geometry for the evaluated point
  //! \returns Evaluated spectrum
  Spectrum_t Eval(const glm::vec3 &outgoingWorld,
                  const glm::vec3 &incomingWorld,
                  const peCoordSys &shadingCoordSys,
                  const glm::vec3 &geometricNormal, BxDFType flags) const;

  //! \brief Sample this BSDF by evaluating a random BxDF
  Spectrum_t Sample_f(const glm::vec3 &wo, glm::vec3 &wi,
                      const peCoordSys &shadingCoordSys,
                      const glm::vec3 &geometryNormal, const BSDFSample &sample,
                      float &pdf, BxDFType flags, BxDFType &sampledType) const;

  //! \brief Sums up the hemispherical-hemispherical reflectance values of
  //! all assigned BxDFs \param rnd Random number generator for monte carlo
  //! method \param flags Types of BxDFs to sample \param sqrtSamples
  //! Sampling parameter for monte carlo \returns Summed reflectance
  Spectrum_t rho(std::default_random_engine &rnd,
                 BxDFType flags = BxDFType::All,
                 uint32_t sqrtSamples = 6) const;

  //! \brief Sums up the hemispherical-directional reflectance values of all
  //! assigned BxDFs \param wo Outgoing direction \param rnd Random number
  //! generator for monte carlo method \param flags Types of BxDFs to sample
  //! \param sqrtSamples Sampling parameter for monte carlo
  //! \returns Summed reflectance
  Spectrum_t rho(const glm::vec3 &wo, std::default_random_engine &rnd,
                 BxDFType flags = BxDFType::All,
                 uint32_t sqrtSamples = 6) const;

  //! \brief Probability density function of this BSDF for the given pair of
  //! vectors \param wo Outgoing vector in world space \param wi Incoming vector
  //! in world space
  float Pdf(const glm::vec3 &wo, const glm::vec3 &wi,
            const peCoordSys &shadingCoordSys,
            BxDFType flags = BxDFType::All) const;

private:
  constexpr static size_t MaxBxDF = 8;
  uint32_t _numBxdfs;
  std::array<peBxDF const *, MaxBxDF> _bxdfs;
};

#pragma region HelperFunctions

//! \brief Computes fresnel reflectance for dielectric materials
//! \param cosIncident Cosine of angle of incident direction
//! \param cosTransmitted Cosine of angle of transmitted direction
//! \param etaIncident Index of refraction for incident medium
//! \param etaTransmitted Index of refraction for transmitted medium
//! \returns Fresnel reflectance
Spectrum_t PE_CORE_API FresnelDielectric(float cosIncident,
                                         float cosTransmitted,
                                         const Spectrum_t &etaIncident,
                                         const Spectrum_t &etaTransmitted);

//! \brief Computes the fresnel reflectance for a conducting material
//! \param cosIndicent Cosine of the angle of incident direction
//! \param eta Index of refraction of the conductor
//! \param k Absorption coefficient
//! \returns Fresnel reflectance
Spectrum_t PE_CORE_API FresnelConductor(float cosIndicent,
                                        const Spectrum_t &eta,
                                        const Spectrum_t &k);

//! \brief Helper to evaluate fresnel shading for dieletric materials
//! \param cosIncident Cosine of angle of incident direction
//! \param etaIncident Index of refraction of incident medium
//! \param etaTransmitted Index of refraction of transmitted medium
//! \returns Fresnel reflectance
Spectrum_t PE_CORE_API EvalFresnelDielectric(float cosIncident,
                                             float etaIncident,
                                             float etaTransmitted);

#pragma endregion

} // namespace pe

#pragma warning(pop)

max_stars_repo_name: Mortano/Prismatic

max_issues_repo_name: Mortano/Prismatic

max_forks_repo_name: Mortano/Prismatic


text:
//
// Created by karl_ on 2020-5-3.
//

#include <stdlib.h>
#include <math.h>
#include <cblas.h>
#include <core/color-space.h>
#include <imgproc/histogram-equalization.h>
#include "color-convert.h"

// ref: https://www.math.uci.edu/icamp/courses/math77c/demos/hist_eq.pdf
static void
histeq(double const *src, unsigned src_stride, double *dst, unsigned dst_stride, size_t area, unsigned range) {
    unsigned *cdf;
    size_t i;
    cdf = calloc(range, sizeof(unsigned));
    for (i = 0; i < area; ++i) {
        cdf[(int) (src[i * src_stride] * range)] += 1;
    }
    for (i = 1; i < range; ++i) {
        cdf[i] += cdf[i - 1];
    }
    for (i = 0; i < area; ++i) {
        dst[i * dst_stride] = floor((double) (range - 1) * cdf[(int) (src[i * src_stride] * range)] / area) / range;
    }
    free(cdf);
}

static gboolean histeq_grayscale(CoreImage *src, CoreImage **dst) {
    gdouble *src_data_cast = NULL;
    gdouble *dst_data_cast = NULL;
    gpointer src_data, dst_data;
    CoreSize *size = NULL;
    guint8 channel;
    gsize i, area;
    CoreColorSpace color_space;
    CorePixelType pixel_type;
    gdouble range, inv_range;
    g_return_val_if_fail(src != NULL, FALSE);
    g_return_val_if_fail(dst != NULL, FALSE);

    channel = core_image_get_channel(src);
    g_return_val_if_fail(channel == 1, FALSE);

    src = g_object_ref(src);
    src_data = core_image_get_data(src);
    size = core_image_get_size(src);
    area = core_size_get_area(size);
    color_space = core_image_get_color_space(src);
    pixel_type = core_image_get_pixel_type(src);
    range = core_pixel_get_range(pixel_type);
    inv_range = 1.0 / range;

    /* cast to double and normalize pixel data to [0, 1] */
    if (core_pixel_is_double(pixel_type)) {
        src_data_cast = src_data;
    } else if (core_pixel_is_uint8(pixel_type)) {
        src_data_cast = g_malloc(sizeof(gdouble) * area);
        for (i = 0; i < area; ++i) {
            src_data_cast[i] = ((guint8 *) src_data)[i] * inv_range;
        }
    } else {
        g_return_val_if_fail(FALSE, FALSE);
    }
    dst_data_cast = g_malloc(sizeof(gdouble) * area);

    /* run the algorithm */
    histeq(src_data_cast, 1, dst_data_cast, 1, area, 256);

    /* cast back to uchar */
    if (core_pixel_is_double(pixel_type)) {
        /* exactly do nothing */
    } else if (core_pixel_is_uint8(pixel_type)) {
        /* reuse allocated memory */
        for (i = 0; i < area; ++i) {
            ((guint8 *) dst_data_cast)[i] = dst_data_cast[i] * 255.0;
        }
    } else {
        g_return_val_if_fail(FALSE, FALSE);
    }
    dst_data = dst_data_cast;
    if (core_pixel_is_uint8(pixel_type)) {
        /* release useless memory, can be removed for speed */
        dst_data = g_realloc(dst_data, sizeof(guint8) * area);
    }

    if (*dst == NULL) {
        *dst = core_image_new_with_data(dst_data, color_space, pixel_type, size, FALSE);
    } else {
        core_image_assign_data(*dst, dst_data, color_space, pixel_type, size, FALSE);
    }

    /* free allocated additional memory */
    if (core_pixel_is_uint8(pixel_type)) {
        g_free(src_data_cast);
    }
    g_object_unref(size);
    g_object_unref(src);
    return TRUE;
}

static gboolean histeq_hsl(CoreImage *src, CoreImage **dst) {
    gdouble *src_data, *dst_data;
    CoreSize *size;
    gsize area, block_size;

    src_data = core_image_get_data(src);
    size = core_image_get_size(src);
    area = core_size_get_area(size);
    block_size = area * 3 * sizeof(gdouble);
    dst_data = g_malloc(block_size);
    memcpy(dst_data, src_data, block_size);
    histeq(src_data + 2, 3, dst_data + 2, 3, area, 255);

    if (*dst == NULL) {
        *dst = core_image_new_with_data(dst_data, CORE_COLOR_SPACE_HSL, CORE_PIXEL_D3, size, FALSE);
    } else {
        core_image_assign_data(*dst, dst_data, CORE_COLOR_SPACE_HSL, CORE_PIXEL_D3, size, FALSE);
    }

    g_object_unref(size);
}

gboolean imgproc_histogram_equalization(CoreImage *src, CoreImage **dst) {
    CoreColorSpace src_color_space = core_image_get_color_space(src);
    CoreImage *dummy = NULL;
    gboolean ret;

    if (src_color_space == CORE_COLOR_SPACE_GRAY_SCALE) {
        ret = histeq_grayscale(src, dst);
    } else {
        if (src_color_space == CORE_COLOR_SPACE_RGB) {
            imgproc_to_HSL(src, &dummy);
        } else {
            dummy = g_object_ref(src);
        }
        ret = histeq_hsl(dummy, dst);
        g_object_unref(dummy);
        if (src_color_space == CORE_COLOR_SPACE_RGB) {
            imgproc_to_RGB(*dst, dst);
        }
    }
    return ret;
}