/* Last edited on 2020-12-11 00:25:56 by jstolfi */

typedef struct nmsim_elem_neuron_trace_t
  { nmsim_elem_neuron_ix_t ine;           /* Index of the neuron in the network. */
    nmsim_group_neuron_ix_t ing;          /* Index of the neuron's group in the network. */
    nmsim_class_neuron_ix_t inc;          /* Index of the neuron's class in the network. */
    nmsim_time_t tini;                    /* Discrete time of first trace entry. */
    nmsim_time_t tfin;                    /* Discrete time of last trace entry. */
    nmsim_elem_neuron_trace_entry_t *ts;  /* Full states and evolution data of the neuron. */
  } nmsim_elem_neuron_trace_t;

nmsim_elem_neuron_trace_t *nmsim_elem_neuron_trace_new
  ( nmsim_elem_neuron_ix_t ine,       /* Index of the neuron in the network. */
    nmsim_group_neuron_ix_t ing,      /* Index of the neuron's group in the network. */
    nmsim_class_neuron_ix_t inc,      /* Index of the neuron's class in the network. */
    nmsim_time_t tini,                /* Discrete time of first trace entry. */
    nmsim_time_t tfin                 /* Discrete time of last trace entry. */
  );
#define nmsim_elem_neuron_trace_read_INFO \
  "  The file begins with a line \n" \
  "      \"begin " nmsim_elem_neuron_trace_FILE_TYPE " format of " nmsim_elem_neuron_trace_VERSION "\"\n" \
  "  Then follow a set of lines in the format {NAME} = {VALUE}, where {NAME} is \"neuron_elem\", " \
  " \"neuron_group\", \"neuron_class\", \"initial_time\", and \"final time\". Then" \
  " there follows one line for each discrete time {t} in the monitored interval.  Each" \
  " line has the time {t}, followed by the state and evolution fields as described below.\n" \
  "\n" \
  "TRACE ENTRY FORMAT\n" \
  "  " nmsim_elem_neuron_trace_entry_read_INFO ""

    demand((ing >= 0) && (ing <= nmsim_group_neuron_ix_MAX), "invalid neuron group");
    demand((inc >= 0) && (inc <= nmsim_class_neuron_ix_MAX), "invalid neuron class");

    fprintf(wr, "%sneuron_group = %d\n", ind2, trne->ing);
    fprintf(wr, "%sneuron_class = %d\n", ind2, trne->inc);

    nmsim_group_neuron_ix_t ing = (nmsim_group_neuron_ix_t)
      nmsim_read_int64_param(rd, "neuron_group", 0, nmsim_group_neuron_ix_MAX);
    nmsim_class_neuron_ix_t inc = (nmsim_class_neuron_ix_t)
      nmsim_read_int64_param(rd, "neuron_class", 0, nmsim_class_neuron_ix_MAX);

/* ---------------------------------------------------------------------- */
/* nmsim_basic.h */

typedef enum {
    nmsim_vtype_VOLT,  /* Voltage value, avg, or dev in millivolts. */
    nmsim_vtype_MULT,  /* Modulator or recharge factor. */
    nmsim_vtype_SYNW,  /* Synaptic strength value, avg, or dev. */
    nmsim_vtype_INDG,  /* Average indegree. */
    nmsim_vtype_FRAC,  /* Fraction of neurons, such as firing rate. */
  } nmsim_vtype_t;
  /* Describes the nature of a value to be given to the {write_double_value}
    and {write_double_parm} procedures below. */

#define write_voltage_DIGITS 8
  /* Max number of decimal fraction digits to use when writing {double}
    values to files. */

#define write_multiplier_DIGITS 8
  /* Max number of decimal fraction digits to use when writing 
    modulator or recharge factor values. */

#define write_neuron_fraction_DIGITS 5
  /* Max number of decimal fraction digits to use when 
    writing neuron fractions, such as firing rate or 
    relative cohort size. */

#define write_indegree_DIGITS 3
  /* Max number of decimal fraction digits to use when 
    writing average indegrees, such as firing rate or 
    relative cohort size. */

void write_voltage_value(FILE *wr, double v, bool_t sign, bool_t special_0);
  /* Writes {v}, assumed to be a volatge or voltage increment in millivolts,
    rounded to three decimal fraction digits, suppressing trailing zeros. 
    
    If {sign} is true, writes an explicit "+" sign for positive values. 
  
    If {special_0} is true, writes the exact value of 0 as "0", and rounds
    any values that would print as "0" away from 0 to print as
    {+0.001} or {-0.001}, preserving their signs. Otherwise may print small
    values as "0". */

void write_multiplier_value(FILE *wr, double v);
  /* Writes {v}, assumed to be a modulator or recharge factor,
    rounded to eight decimal fraction digits, suppressing trailing zeros. 
  
    Writes the exact value of 0 as "0", and rounds
    other values that would print as "0" away from 0 to print as
    {+0.00000001} or {-0.00000001}, preserving their signs.
    
    Also writes the exact value of 1 as "1", otherwise rounds
    values that would print as "1" away from 1 to print as "0.99999999" or "1.00000001",
    depending on it being less than or greater than 1. */

void write_multiplier_value(FILE *wr, double v);
  /* Writes {v}, assumed to be a modulator or recharge factor,
    rounded to eight decimal fraction digits, suppressing trailing zeros. 
  
    Writes the exact value of 0 as "0", and rounds
    other values that would print as "0" away from 0 to print as
    {+0.00000001} or {-0.00000001}, preserving their signs.
    
    Also writes the exact value of 1 as "1", otherwise rounds
    values that would print as "1" away from 1 to print as "0.99999999" or "1.00000001",
    depending on it being less than or greater than 1. */

void write_indegree_value(FILE *wr, double v);
  /* Writes {v}, assumed to be an average indegree rounded to three decimal fraction digits, suppressing trailing zeros. 
  
    Writes the exact value of 0 as "0", and rounds
    other values that would print as "0" away from 0 to print as
    {+0.00000001} or {-0.00000001}, preserving their signs.
    
    Also writes the exact value of 1 as "1", otherwise rounds
    values that would print as "1" away from 1 to print as "0.99999999" or "1.00000001",
    depending on it being less than or greater than 1. */


/* ---------------------------------------------------------------------- */
/* nmsim_elem_net_sim.c */
    
    /* Save the trace data for time {t+1}: */
    if (etrace != NULL) { nmsim_elem_net_trace_set_V_age_H_M(etrace, t+1, nne, V, age, H, M); }

/* ---------------------------------------------------------------------- */
/* nmsim_test_042_group_net_sim.c */

nmsim_class_neuron_t *nmsim_group_net_test_parms_pick(void)
  { 
    /* Baseline at {-40} mV, shallow reset: */
    double V_B = -40.0;   
    double V_R = -30.0;
    
    /* Fixed recharge factor with no recharge modulation: */
    double c_B = 0.90;   
    double M_R = 1.00;   
    double M_mu = 1.00; 
    
    /* No output modulation: */
    double H_R = 1.00;   
    double H_mu = 1.00;  
    
    /* Firing function parameters: */
    double V_M = -20.0;
    double V_D = 10.0;
    
    fprintf(stderr, "creating firing function record...\n");
    nmsim_firing_func_t *Phi = nmsim_firing_func_new("gauss", 0, V_M, V_D);
    fprintf(stderr, "creating neuron class record...\n");
    nmsim_class_neuron_t *nclass = 
      nmsim_class_neuron_new(V_B, V_R, c_B, M_R, M_mu, G_R, G_mu, H_R, H_mu, Phi);
    nmsim_test_plot_firing_func(nclass);
    return nclass;
  }

nmsim_group_net_t *nmsim_group_net_test_new(uint64_t nng, nmsim_elem_neuron_count_t N)
  { 
    /* Choose the neuron parameters: */
    nmsim_class_neuron_t *nclass = nmsim_group_net_test_nclass_pick();
    
    /* Create the pops: */
    fprintf(stderr, "allocating the meanfield net model...\n");
    nmsim_group_net_t *net = nmsim_group_net_new(nng);
    assert(net->N == N);
    fprintf(stderr, "setting the parameters of each neuron class...\n");
    for (nmsim_pop_ix_t i = 0; i < N; i++) { net->nclass[i] = nclass; }
    
    /* Create random synapses: */
    double pr = 1.0/sqrt((double)N); /* Average in- and out-degree. */
    fprintf(stderr, "creating synapses with pr = %6.4f...\n", pr);
    double swt_tot_avg = 5.0;
    double swt_tot_dev = 0.1*swt_tot_avg;
    fprintf(stderr, "tot input weight avg = %+8.3f dev = %7.3f\n", swt_tot_avg, swt_tot_dev);
    nmsim_group_net_add_random_synapses(net, 0, N-1, 0, N-1, pr, swt_tot_avg, swt_tot_dev);
    fprintf(stderr, "network created.\n");
    return net;
  }
  
void nmsim_test_net_init(nmsim_group_net_t *net)
  {
    fprintf(stderr, "initializing neuron group states...\n");
    nmsim_elem_neuron_count_t N = net->N;
    for (nmsim_pop_ix_t i = 0; i < N; i++) 
      { nmsim_pop_state_t *state = &(net->state[i]); 
        nmsim_class_neuron_t *nclass = net->nclass[i];
        double rho = 0.005; /* For now, assume firing every 100 steps. */
        nmsim_pop_state_throw(state, nclass, rho);
        if ((i < 10) || (i >= N-10))
          { fprintf(stderr, "  ngrp[%lu]: V = %+6.2f age = %lu\n", i, state->V, state->age); }
      }
  }
  
void nmsim_test_net_simulate(nmsim_elem_net_t *net, uint64_t T)
  {
    fprintf(stderr, "simulating for T = %lu steps...\n", T);
    nmsim_elem_neuron_count_t N = net->N;
    
    /* Open trace file: */
    char *fname = NULL;
    asprintf(&fname, "out/net-%09lu-trace.txt", N);
    FILE *wr = open_write(fname, TRUE);
    
    /* Decide how many neurons {NW} to trace: */
    nmsim_elem_neuron_count_t NW_max = 20;
    nmsim_elem_neuron_count_t NW = (N < NW_max ? N : NW_max);
    
    /* Allocate work vectors: */
    bool_t *X = notnull(malloc(N*sizeof(bool_t)), "no mem"); /* Firing indicators. */ 
    double *dV = notnull(malloc(N*sizeof(double)), "no mem"); /* Input potential increments. */ 
    
    /* Run simulation: */
    for (uint64_t t = 0; t < T; t++)
      { 
        bool_t debug = (N <= 10) && ((t < 10) || (t >= T-10));
        
        if (debug) { fprintf(stderr, "iteration %lu -> %lu\n", t, t+1); }

        if (debug) { nmsim_test_debug_ages("ai", net); }
        if (debug) { nmsim_test_debug_potentials("Vi", net); }
        /* First pass: compute the firing indicators {X[0..N-1]} and activity {rho}: */
        nmsim_elem_neuron_count_t NX = 0; /* Count of neurons that fired during step. */
        for (nmsim_elem_neuron_ix_t i = 0; i < N; i++) 
          { nmsim_neuron_state_t *state = &(net->state[i]); 
            nmsim_class_neuron_t *nclass = net->nclass[i];
            double pr;
            nmsim_firing_func_eval(nclass->Phi, state->V, &pr, NULL);
            X[i] = (drandom() < pr);
            NX += (int)(X[i]);
          }
        if (debug) { nmsim_test_debug_firings("X ", net, X); }
        double rho = ((double)NX)/((double)N);

        /* Second pass: Compute the potential increments {dv[0..N-1]}: */
        nmsim_elem_net_tot_inputs(net, X, NULL, dV);
        
        /* Third pass: reset or decay the previous potentials,: */
        for (nmsim_elem_neuron_ix_t i = 0; i < net->N; i++) 
          { nmsim_class_neuron_t *nclass = net->nclass[i];
            nmsim_neuron_state_t *state = &(net->state[i]);
            if (X[i]) 
              { state->V = nclass->V_R; }
            else
              { double L = nmsim_neuron_state_potential_decay_modulator(nclass, state->age);
                double c = L*nclass->c_B;
                state->V = c * state->V;
              }
          }
        if (debug) { nmsim_test_debug_potentials("Vm", net); }

        /* Third pass: add input increments and update state: */
        for (nmsim_elem_neuron_ix_t i = 0; i < net->N; i++) 
          { nmsim_neuron_state_t *state = &(net->state[i]);
            state->V += dV[i];
            state->age = (X[i] ? 0 : state->age + 1);  
         }
        if (debug) { nmsim_test_debug_potential_increments("dV", net, dV); }
        if (debug) { nmsim_test_debug_potentials("Vf", net); }
        if (debug) { nmsim_test_debug_ages("af", net); }

        if (debug) { fprintf(stderr, "\n"); }
        
        /* Write trace: */
        nmsim_test_write_trace(wr, NW, N, t, rho, X, NULL, dV, net->state);
        
      }
  }

void nmsim_test_write_trace
  ( FILE *wr, 
    nmsim_elem_neuron_count_t NW, 
    nmsim_elem_neuron_count_t N, 
    uint64_t t, 
    double rho, 
    bool_t X[], 
    double J[],
    double dV[], 
    nmsim_neuron_state_t state[]
  )
  {
    fprintf(wr, "%10lu %10.8f", t, rho);
    for (nmsim_elem_neuron_ix_t kW = 0; kW < NW; kW++) 
      { uint64_t m = (uint64_t)N-1;
        uint64_t mW = (uint64_t)NW-1;
        nmsim_elem_neuron_ix_t i = (NW == 1 ? 0 : (nmsim_elem_neuron_ix_t)((kW*m)/mW));
        fprintf(wr, "  ");
        fprintf(wr, " %d", (X == NULL ? 0 : (int)(X[i])));
        fprintf(wr, " %+6.3f", (J == NULL ? 0.0 : J[i]));
        fprintf(wr, " %+6.3f", (dV == NULL ? 0.0 : dV[i]));
        fprintf(wr, " %+6.3f", (state == NULL ? 0.0 : state[i].V));
        fprintf(wr, " %lu", (state == NULL ? 0 : state[i].age));
      }
    fprintf(wr, "\n");
    fflush(wr);
  }

void nmsim_test_debug_firings(char *lab, nmsim_elem_net_t *net, bool_t X[])
  { fprintf(stderr, "  %s = ", lab);        
    for (nmsim_elem_neuron_ix_t i = 0; i < net->N; i++) 
      { fprintf(stderr, " %6d", (int)(X[i])); }
    fprintf(stderr, "\n");
  }
          
void nmsim_test_debug_potentials(char *lab, nmsim_elem_net_t *net)
  { fprintf(stderr, "  %s = ", lab);        
    for (nmsim_elem_neuron_ix_t i = 0; i < net->N; i++) 
      { nmsim_neuron_state_t *state = &(net->state[i]);
        fprintf(stderr, " %6.2f", state->V);
      }
    fprintf(stderr, "\n");
  }
           
void nmsim_test_debug_ages(char *lab, nmsim_elem_net_t *net)
  { fprintf(stderr, "  %s = ", lab);        
    for (nmsim_elem_neuron_ix_t i = 0; i < net->N; i++) 
      { nmsim_neuron_state_t *state = &(net->state[i]);
        fprintf(stderr, " %6lu", state->age);
      }
    fprintf(stderr, "\n");
  }
           
void nmsim_test_debug_potential_increments(char *lab, nmsim_elem_net_t *net, double dV[])
  { fprintf(stderr, "  %s = ", lab);        
    for (nmsim_elem_neuron_ix_t i = 0; i < net->N; i++) 
      { fprintf(stderr, " %6.2f", dV[i]); }
    fprintf(stderr, "\n");
  }

/* ---------------------------------------------------------------------- */
/* nmsim_group_net_trace.{h,c} */

/* TRACE SAVING PROCEDURES 
  
  These procedures receive vectors {V,X,age,H,M,I,J} indexed by
  neuron group number {ing} in {0..nng-1}, for time {t}
  or the time step from {t} to {t+1}.  They store those values
  in the corresponding fields of the traces of monitored neuron groups.
  
  More precisely, for each non-null neuron group trace {tr = gtrace.tr[i]}
  with {t} in {tr.tini .. tr.tfin}, the procedures store in 
  state record {st = tr->st[t - tr.tini]} the value of {V[tr.ing]}
  in the {.V} field, {X[tr.ing]} in the {.X} field, etc.. */

void nmsim_group_neuron_trace_set_t_V_X
  ( nmsim_group_net_trace_t *gtrace,
    nmsim_time_t t,                /* Time at start of step. */
    nmsim_group_neuron_count_t nng, /* Number of neuron groups in net. */
    double V[],                    /* Membrane potential of each neuron group at {t}. */
    bool_t X[]                     /* Firing indicators between {t} and {t+1}. */
  );

void nmsim_group_neuron_trace_set_age_H_I_J
  ( nmsim_group_net_trace_t *gtrace,
    nmsim_time_t t,                /* Time at start of step. */
    nmsim_group_neuron_count_t nng, /* Number of neuron groups in net. */
    nmsim_step_count_t age[],      /* Firing age of each neuron group at time {t}. */
    double H[],                    /* Output modulator of each neuron group at time {t}. */
    double I[],                    /* External input of each neuron group from {t} to {t+1}. */
    double J[]                     /* Total input of each neuron group from {t} to {t+1}. */
  );

void nmsim_group_neuron_trace_set_M
  ( nmsim_group_net_trace_t *gtrace,
    nmsim_time_t t,                /* Time at start of step. */
    nmsim_group_neuron_count_t nng, /* Number of neuron groups in net. */
    double M[]                     /* Recharge modulator of each neuron group at time {t}. */
  );

void nmsim_group_neuron_trace_set_t_V_X
  ( nmsim_group_net_trace_t *etrace,
    nmsim_time_t t,                 /* Time at start of step. */
    nmsim_group_neuron_count_t nng, /* Number of neuron groups in net. */
    double V[],                     /* Membrane potential of each neuron group at {t}. */
    bool_t X[]                      /* Firing indicators between {t} and {t+1}. */
  )
  { if ((etrace != NULL) && (t >= etrace->tini) && (t <= etrace->tfin))
      { for (nmsim_group_neuron_ix_t i = 0; i < etrace->nng; i++)
          { nmsim_group_neuron_trace_t *tr = etrace->tr[i];
            if ((tr != NULL) && (t >= tr->tini) && (t <= tr->tfin))
              { nmsim_group_neuron_ix_t ing = tr->ing; /* Index of monitored neuron group. */
                assert((ing >= 0) && (ing < nng));
                nmsim_group_neuron_state_t *st = &(tr->st[t - tr->tini]); /* State rec for time {t}. */
                st->t = t;
                st->V = V[ing];
                st->X = X[ing];
              }
          }
      }
  }

void nmsim_group_neuron_trace_set_age_H_I_J
  ( nmsim_group_net_trace_t *etrace,
    nmsim_time_t t,                /* Time at start of step. */
    nmsim_group_neuron_count_t nng, /* Number of neuron groups in net. */
    nmsim_step_count_t age[],      /* Firing age of each neuron group at time {t}. */
    double H[],                    /* Output modulator of each neuron group at time {t}. */
    double I[],                    /* External input of each neuron group from {t} to {t+1}. */
    double J[]                     /* Total input of each neuron group from {t} to {t+1}. */
  )
  { if ((etrace != NULL) && (t >= etrace->tini) && (t <= etrace->tfin))
      { for (nmsim_group_neuron_ix_t i = 0; i < etrace->nng; i++)
          { nmsim_group_neuron_trace_t *tr = etrace->tr[i];
            if ((tr != NULL) && (t >= tr->tini) && (t <= tr->tfin))
              { nmsim_group_neuron_ix_t ing = tr->ing; /* Index of monitored neuron group. */
                assert((ing >= 0) && (ing < nng));
                nmsim_group_neuron_state_t *st = &(tr->st[t - tr->tini]); /* State rec for time {t}. */
                assert(st->t == t);
                st->age = age[ing];
                st->H = H[ing];
                st->I = I[ing];
                st->J = J[ing];
              }
          }
      }
  }

void nmsim_group_neuron_trace_set_M
  ( nmsim_group_net_trace_t *etrace,
    nmsim_time_t t,                /* Time at start of step. */
    nmsim_group_neuron_count_t nng, /* Number of neuron groups in net. */
    double M[]                     /* Potential decay modulator of each neuron group at time {t}. */
  )
  { if ((etrace != NULL) && (t >= etrace->tini) && (t <= etrace->tfin))
      { for (nmsim_group_neuron_ix_t i = 0; i < etrace->nng; i++)
          { nmsim_group_neuron_trace_t *tr = etrace->tr[i];
            if ((tr != NULL) && (t >= tr->tini) && (t <= tr->tfin))
              { nmsim_group_neuron_ix_t ing = tr->ing; /* Index of monitored neuron group. */
                assert((ing >= 0) && (ing < nng));
                nmsim_group_neuron_state_t *st = &(tr->st[t - tr->tini]); /* State rec for time {t}. */
                assert(st->t == t);
                st->M = M[ing];
              }
          }
      }
  }

/* ---------------------------------------------------------------------- */
/* nmsim_elem_net.{h,c} */    
/*     
    On input, the variable {pnet->nse} must contain the number {nse}
    of previously defined synapses. On output, it will be incremented
    with the number of synapses created by the function call. The
    vector {pnet->bparms} must be allocated before the call (even if
    with zero size) and will be extended by the function if and when needed.
    
    The parameter {pr} is basically the probability of adding a synapse
    between any neuron {ia} in {ine0..ine1} and any neuron {ib} in {jne0..jne1}.
    In particular, if {pr} is zero, no synapses are added and the call is
    a no-op.  
    
    If {pr} is positive, usually at most one synapse is added, with
    probability {pr}, for each of those pairs {ia,ib}. Then the number
    of synapses added by the call is a random variable with binomial
    distribution, max value {na*nb} and mean value {pr*na*nb} where
    {na = ine1-ine0+1} and {nb = jne1-jne0+1}. On average, the call will
    increment the output degree of each neuron in {ine0..ine1} by
    {pr*nb} synapses, and the input degree of each neuron in
    {jne0..jne1} by {pr*na} synapses.
    
    Exceptionally, if this process gives zero new synapses to a neuron
    {ib}, but {pr} is not zero, one synapse is added anyway, from a
    neuron {ia} randomly chosen among {ine0..ine1}. Thus, if {pr} is not
    zero, the procedure always adds at least {nb} synapses. This
    exceptional case is unlikely to occur if {pr*na} is substantially
    greater than 1.
    
    The weights of the added synapses will have a log-normal
    distribution. Specifically, the log {zwt} of the absolute weight
    will be a random variable with normal distribution.  The
    sign of the weight will be the same as that of {swt_tot_avg}.
    
    The mean {zavg} and deviation {zdev} of the log-weight {zwt} will
    be chosen so that the SUM of the weights of the synapses entering
    each neuron {ib} in {jne0..jne1} is a random variable with mean
    {swt_tot_avg} and deviation {swt_tot_dev}. In particular, if
    {swt_tot_dev} is zero, every new synapse entering {ib} will have
    the same weight {swt = swt_tot_avg/din} where {din} is the number
    of new input synapses added to {ib}. */

void nmsim_elem_net_add_random_synapses
  ( nmsim_elem_net_t *net,
    nmsim_neuron_elem_ix_t ine0,
    nmsim_neuron_elem_ix_t ine1,
    nmsim_neuron_elem_ix_t jne0,
    nmsim_neuron_elem_ix_t jne1,
    double pr,
    double swt_tot_avg,
    double swt_tot_dev
  );
  /* Appends to the vector {net->syn[0..nse-1]} a number of random synapses
    from neurons {ine0..ine1} to neurons {jne0..jne1}.
    
    On input, the variable {net->nse} must contain the number {nse}
    of previously defined synapses. On output, it will be incremented
    with the number of synapses created by the function call. The
    vector {net->syn} must be allocated before the call (even if
    with zero size) and will be extended by the function if and when needed.
    
    The parameter {pr} is basically the probability of adding a synapse
    between any neuron {ia} in {ine0..ine1} and any neuron {ib} in {jne0..jne1}.
    In particular, if {pr} is zero, no synapses are added and the call is
    a no-op.  
    
    If {pr} is positive, usually at most one synapse is added, with
    probability {pr}, for each of those pairs {ia,ib}. Then the number
    of synapses added by the call is a random variable with binomial
    distribution, max value {na*nb} and mean value {pr*na*nb} where
    {na = ine1-ine0+1} and {nb = jne1-jne0+1}. On average, the call will
    increment the output degree of each neuron in {ine0..ine1} by
    {pr*nb} synapses, and the input degree of each neuron in
    {jne0..jne1} by {pr*na} synapses.
    
    Exceptionally, if this process gives zero new synapses to a neuron
    {ib}, but {pr} is not zero, one synapse is added anyway, from a
    neuron {ia} randomly chosen among {ine0..ine1}. Thus, if {pr} is not
    zero, the procedure always adds at least {nb} synapses. This
    exceptional case is unlikely to occur if {pr*na} is substantially
    greater than 1.
    
    The weights of the added synapses will have a log-normal
    distribution. Specifically, the log {zwt} of the absolute weight
    will be a random variable with normal distribution.  The
    sign of the weight will be the same as that of {swt_tot_avg}.
    
    The mean {zavg} and deviation {zdev} of the log-weight {zwt} will
    be chosen so that the SUM of the weights of the synapses entering
    each neuron {ib} in {jne0..jne1} is a random variable with mean
    {swt_tot_avg} and deviation {swt_tot_dev}. In particular, if
    {swt_tot_dev} is zero, every new synapse entering {ib} will have
    the same weight {swt = swt_tot_avg/din} where {din} is the number
    of new input synapses added to {ib}. */

nmsim_elem_net_t *nmsim_elem_net_throw
  ( nmsim_group_net_t *gnet,
    nmsim_neuron_elem_count_t nne, 
    nmsim_synapse_elem_count_t nse
  )
  {
    demand((nne >= 0) && (nne <= nmsim_neuron_elem_count_MAX), "invalid neuron count");
    demand((nse >= 0) && (nse <= nmsim_synapse_elem_count_MAX), "invalid synapse count");

    nmsim_elem_net_t *net = nmsim_elem_net_new(gnet);

    /* Add neurons: */
    for (nmsim_neuron_elem_ix_t ine = 0; ine < nne; ine++)
      { nmsim_neuron_elem_t neu = nmsim_neuron_elem_throw(gnet->nng - 1);
        nmsim_neuron_elem_ix_t jne = nmsim_elem_net_add_neuron(net, &neu);
        assert(jne == ine);
      }
    
    /* Add synapses: */
    demand((nse == 0) || (nne > 0), "no neurons for synapses");
    for (nmsim_synapse_elem_ix_t ise = 0; ise < nse; ise++)
      { /* Pick a random synapse group {isg}, and get its parameters {sclass}: */
        nmsim_synapse_group_ix_t isg = (nmsim_synapse_group_ix_t)int64_abrandom(0,gnet->nsg-1);
        nmsim_synapse_class_ix_t isc = gnet->sgrp.e[isg].isc;
        nmsim_synapse_class_t *sclass = gnet->cnet->sclass.e[isc];
        /* Throw a synapse of that group: */
        double W_avg = sclass->W_avg;
        double W_dev = sclass->W_dev;
        nmsim_synapse_elem_t syn = nmsim_synapse_elem_throw(gnet->nsg - 1, net->nne - 1, W_avg, W_dev);
        syn.isg = isg; /* !!! UGLY !!! */
        nmsim_synapse_elem_ix_t jse = nmsim_elem_net_add_synapse(net, &syn);
        assert(jse == ise);
      }
  }
  
/* ---------------------------------------------------------------------- */
/* nmsim_neuron_class.h */

double nmsim_neuron_elem_state_leakage_modulator(nmsim_neuron_class_t *nclass, uint64_t age);
double nmsim_neuron_elem_state_input_modulator(nmsim_neuron_class_t *nclass, uint64_t age);
double nmsim_neuron_elem_state_output_modulator(nmsim_neuron_class_t *nclass, uint64_t age);  /*
 Return the leakage modulator {L}, input modulator {G}, and output modulator {H} of a 
    neuron with parameters {nclass} and given {age}. */ 

void nmsim_neuron_elem_state_throw(nmsim_neuron_elem_state_t *state, nmsim_neuron_class_t *nclass, double rho);
  /* Sets {*state} to a random state compatible with {nclass}.
    Assumes that the average activity of the neuron in the past was {rho},
    and therefore sets the firing age to about {1/rho}. */ 
    
double nmsim_neuron_elem_state_exp_dyn(double R, double mu, uint64_t age);
  /* The value of a modulator with reset-decay parameters {R,mu}
    for a neuron with specified firing {age}. */

double nmsim_neuron_elem_state_leakage_modulator(nmsim_neuron_class_t *nclass, uint64_t age)
  { return nmsim_neuron_elem_state_exp_dyn(nclass->M_R, nclass->M_mu, age); }
  
double nmsim_neuron_elem_state_input_modulator(nmsim_neuron_class_t *nclass, uint64_t age)
  { return nmsim_neuron_elem_state_exp_dyn(nclass->G_R, nclass->G_mu, age); }

double nmsim_neuron_elem_state_output_modulator(nmsim_neuron_class_t *nclass, uint64_t age)
  { return nmsim_neuron_elem_state_exp_dyn(nclass->H_R, nclass->H_mu, age); }
  
double nmsim_neuron_elem_state_exp_dyn(double R, double mu, uint64_t age)
  {
    /* Should use a table and have a better cutoff. */
    if ((mu = 1.0) || (R == 1.0) || (age > 1000000)) { return 1.0; }
    double P = 1.0 - (1.0 - R)*pow(mu, (double)age);
    return P;
  }

void nmsim_neuron_elem_state_throw(nmsim_neuron_elem_state_t *state, nmsim_neuron_class_t *nclass, double rho) 
  { 
    demand((rho > 1.0e-8) && (rho <= 1.0), "invalid mean activity");
    /* Random potential: */
    double V_sat = nmsim_firing_func_eval_inv(nclass->Phi, 0.999);
    double V_min = fmin(nclass->V_B, nclass->V_R);
    double V_max = fmax(fmax(nclass->V_B, nclass->V_R), V_sat);
    state->V = dabrandom(V_min, V_max);

    /* Random age: */
    uint64_t max_age = ((uint64_t)floor(2.0/rho + 0.5)) - 1;
    state-> age = uint64_abrandom(0,max_age);
  }

/* ---------------------------------------------------------------------- */
/* nmsim_firing_func.{h,c} */

typedef struct nmsim_firing_func_t nmsim_firing_func_t;
  /* A descriptor of of a firing function {Phi}. */

struct nmsim_firing_func_t *nmsim_firing_func_new
  ( char *class, 
    double deg, 
    double V_M, 
    double V_D
  );
  /* Returns a descriptor for a firing function of the given {class},
    with midpoint potential {V_M} and slope parameter
    {V_D} (in millivoltds), and degree {deg}. The current classes are
    described in {nmsim_firing_func_INFO}. */
    
  /* (which should have been multiplied by the firing modulator) */

struct nmsim_firing_func_t *nmsim_firing_func_new
  ( char *class, 
    double deg, 
    double V_M,
    double V_D
  )
  {
    nmsim_firing_func_t *Phi = NULL;
    if ((strcmp(class, "Gauss") == 0) || (strcmp(class, "gauss") == 0))
      { demand(deg == 0, "invalid degree for Gauss Phi");
        Phi = nmsim_firing_func_gauss_new(V_M, V_D);
      }
    else
      { demand(FALSE, "invalid class"); }
    return Phi;
  }

void nmsim_firing_func_eval
  ( nmsim_firing_func_t Phi,
    double V, 
    double *prP,
    double *dprP
  )
  {
    Phi->eval(V, Phi->deg, Phi->V_M, Phi->V_D, prP, dprP);
  }
 
double nmsim_firing_func_eval_inv(struct nmsim_firing_func_t *Phi, double pr)
  {
    demand((pr >= 0) && (pr <= 1), "invalid probability");
    return Phi->eval_inv(pr, Phi->deg, Phi->V_M, Phi->V_D);
  }

nmsim_firing_func_t *nmsim_firing_func_new_gen
  ( nmsim_firing_func_eval_proc_t *eval, 
    nmsim_firing_func_eval_inv_proc_t *eval_inv, 
    double deg, 
    double V_M, 
    double V_D,
    char *descr
  )
  {
    nmsim_firing_func_t *Phi = notnull(malloc(sizeof(nmsim_firing_func_t)), "no mem");
    (*Phi) = (nmsim_firing_func_t){ 
      .eval = eval, .eval_inv = eval_inv,
      .deg = deg, .V_M = V_M, .V_D = V_D,
      .descr = descr
    };
    return Phi;
  }

nmsim_firing_func_t *nmsim_firing_func_gauss_new(double V_M, double V_D)
  {
    char *descr = NULL;
    asprintf(&descr, "PhiGauss[V_M=%.3f,V_D=%.3f]", V_M, V_D);
    nmsim_firing_func_t *Phi = nmsim_firing_func_new_gen
      ( &nmsim_firing_func_gauss_eval, 
        &nmsim_firing_func_gauss_eval_inv,
        0, V_M, V_D,
        descr
      );
    return Phi;
  }
    
void nmsim_firing_func_free(nmsim_firing_func_t *Phi);
  /* Releases the storage of the record {*Phi}. Does nothing if {Phi==NULL}. */

nmsim_firing_func_t *nmsim_firing_func_new
  ( nmsim_firing_func_class_t class,
    double V_M,
    double V_D
  )
  { nmsim_firing_func_t *Phi = notnull(malloc(sizeof(nmsim_firing_func_t)), "no mem");
    (*Phi) = (nmsim_firing_func_t) {
      .class = class, .V_M = V_M, .V_D = V_D
    };
    return Phi;
  }

void nmsim_firing_func_free(nmsim_firing_func_t *Phi)
  { 
    if (Phi != NULL) { free(Phi); }
  }


/* ---------------------------------------------------------------------- */
/* nmsim_bundle_parms.{h,c} */

    /* INTERNAL IMPLEMENTATIONS */

    double read_param(char *name, double vmin, double vmax)
      { double v = nget_double(rd, name);
        if ((v < vmin) || (v > vmax))
          { fprintf
              ( stderr, "** parameter {%s} = %24.16e is out of range [ %24.16e _ %24.16e ]\n", 
                name, v, vmin, vmax
              ); 
            demand(FALSE, "aborted");
          }
        fget_eol(rd);
        return v;
      }

/* ---------------------------------------------------------------------- */
/* nmsim_pop_parms.{h,c} */

    double read_parm(char *name, double vmin, double vmax)
      { double v = nget_double(rd, name);
        if ((v < vmin) || (v > vmax))
          { fprintf
              ( stderr, "** parameter {%s} = %24.16e is out of range [ %24.16e _ %24.16e ]\n", 
                name, v, vmin, vmax
              ); 
            demand(FALSE, "aborted");
          }
        fget_eol(rd);
        return v;
      }
    
/* ---------------------------------------------------------------------- */
/* nmsim_elem_net.{h,c} */

    /* These arrays are indexed with {i} in {0..N-1} and synapse index in {0..deg_in[i]-1}: */
    nmsim_neuron_elem_ix_t **src_in; /* Indices of source neurons of input synapses. */
    double **swt_in;    /* Rest weights of input synapses. */
    /* These parameters are indexed with {i} in {0..N-1} and synapse index in {0..deg_ot[i]-1}: */
    nmsim_neuron_elem_ix_t **dst_ot; /* Indices of destination neurons of output synapses. */

void nmsim_elem_net_set_synapses
  ( nmsim_elem_net_t *net,
    nmsim_synapse_elem_count_t M,
    nmsim_neuron_elem_ix_vec_t *pre,
    nmsim_neuron_elem_ix_vec_t *pos,
    double_vec_t *swt
  );
  /* Sets the synapses of network {net}to be the {M} synapses 
    in the vectors {pre,pos,swt}.  Any previous synapses of {net} are
    discarded.
    
    Specifically, for each {k} in {0..M-1}, adds to {net}
    a synapse from neyron {pre[k]} to neuron {pos[k]}, with 
    resting weight {swt[k]}.  Alloes synapses from a neuron to itself and
    multiple synapses between the same pair of neurons.
    
    The synapse lists {net.syn_in} and {net.syn_ot} are left unsorted. */

    net->src_in = notnull(malloc(N*sizeof(nmsim_neuron_elem_ix_t *)), "no mem");
    net->swt_in = notnull(malloc(N*sizeof(double*)), "no mem");
    net->dst_ot = notnull(malloc(N*sizeof(nmsim_neuron_elem_ix_t *)), "no mem");

        net->src_in[i] = NULL;
        net->swt_in[i] = NULL;
        net->dst_ot[i] = NULL;

void nmsim_elem_net_create_synapses
  ( nmsim_elem_net_t *net,
    nmsim_synapse_elem_count_t M,
    double swt_avg,
    double swt_dev,
    double pr_neg
  )
  {
    bool_t debug = TRUE;
    nmsim_neuron_elem_count_t N = net->N;
    /* Pass 0: clear degree counts, pick inhibitors, make sure that synapse lists are {NULL}: */
    bool_t *inhib = notnull(malloc(N*sizeof(bool_t)), "no mem"); /* Inhibitory flag. */
    nmsim_neuron_elem_count_t N_inhib = 0; /* Count of inhibitory neurons. */
    for (nmsim_neuron_elem_ix_t i = 0; i < N; i++)
      { net->deg_in[i] = 0;
        net->deg_ot[i] = 0;
        if (net->src_in[i] != NULL) { free(net->src_in[i]); net->src_in[i] = NULL; }
        if (net->swt_in[i] != NULL) { free(net->swt_in[i]); net->swt_in[i] = NULL; }
        if (net->dst_ot[i] != NULL) { free(net->dst_ot[i]); net->dst_ot[i] = NULL; }
        inhib[i] = (drandom() <= pr_neg);
        if (inhib[i]) { N_inhib++; }
      }
    if (debug) { fprintf(stderr, "selected %lu inhibitors out of %lu neurons\n", N_inhib, N); }
    
    /* Pass 1: generate the synapses as pairs {pre[k],pos[k]} for {k} in {0..M-1}: */
    /* Also accumulate degrees: */
    if (debug) { fprintf(stderr, "choosing %lu pre/post neuron pairs\n", M); }
    nmsim_neuron_elem_ix_t *pre = notnull(malloc(M*sizeof(nmsim_neuron_elem_ix_t)), "no mem");
    nmsim_neuron_elem_ix_t *pos = notnull(malloc(M*sizeof(nmsim_neuron_elem_ix_t)), "no mem");
    for (nmsim_synapse_elem_ix_t k = 0; k < M; k++)
      { nmsim_neuron_elem_ix_t j = uint64_abrandom(0, N-1); /* Pre-synaptic (source) neuron. */
        nmsim_neuron_elem_ix_t i = uint64_abrandom(0, N-1); /* Post-synaptic (dest) neuron. */
        assert((i < N) && (j < N));
        pre[k] = j; net->deg_ot[j]++;
        pos[k] = i; net->deg_in[i]++;
      }

    /* Pass 2: allocate input and output synapse tables, reset degrees: */
    if (debug) { fprintf(stderr, "allocating input/output tables\n"); }
    nmsim_synapse_elem_count_t M_in = 0, M_ot = 0; /* Total in- and out-degrees. */
    for (nmsim_neuron_elem_ix_t i = 0; i < N; i++) 
      { 
        /* Allocate input source and weight table: */
        nmsim_neuron_elem_count_t din = net->deg_in[i];
        assert(din <= M);
        nmsim_neuron_elem_ix_t *src = notnull(malloc(din*sizeof(nmsim_neuron_elem_ix_t)), "no mem");
        double *swt = notnull(malloc(din*sizeof(double)), "no mem");
        for(nmsim_synapse_elem_ix_t k = 0; k < din; k++) { src[k] = i; swt[k] = NAN; /* For now. */ }
        M_in += din;
        net->src_in[i] = src;
        net->swt_in[i] = swt;
        net->deg_in[i] = 0;

        /* Allocate output destination table: */
        nmsim_neuron_elem_count_t dot = net->deg_ot[i];
        assert(dot <= M);
        nmsim_neuron_elem_ix_t *dst = notnull(malloc(dot*sizeof(nmsim_neuron_elem_ix_t)), "no mem");
        for(nmsim_synapse_elem_ix_t k = 0; k < dot; k++) { dst[k] = i; /* For now. */ }
        M_ot += dot;
        net->dst_ot[i] = dst;
        net->deg_ot[i] = 0;
        
        if (debug && ((i < 10) || (i >= N-10))) 
          { fprintf(stderr, "  neuron[%lu]: deg_in = %lu deg_ot = %lu\n", i, din, dot); }
      }
    assert((M_in == M) && (M_ot == M)); /* Paranoia. */

    /* Pass 3: collect synapses, set weights, increment degrees. */
    if (debug) { fprintf(stderr, "collecting the synapses...\n"); }
    double swt_avg_log = log(swt_avg);
    double swt_dev_log = log((swt_avg + swt_dev)/swt_avg);
    for (nmsim_synapse_elem_ix_t k = 0; k < M; k++)
      { /* Grab neuron indices: */
        nmsim_neuron_elem_ix_t j = pre[k]; /* Pre-synaptic (source) neuron. */
        nmsim_neuron_elem_ix_t i = pos[k]; /* Post-synaptic (dest) neuron. */
        
        /* Compute the synapse's rest weight {wtk}, store as the input weight of {j}: */
        double wtk = exp(swt_avg_log + swt_dev_log*dgaussrand());
        if (inhib[i]) { wtk = - wtk; }

        /* Store indices and weight in the synapse tables: */
        nmsim_neuron_elem_ix_t ki = net->deg_in[i]; 
        net->src_in[i][ki] = j;
        net->swt_in[i][ki] = wtk;
        net->deg_in[i]++;
        
        nmsim_neuron_elem_ix_t kj = net->deg_ot[j];
        net->dst_ot[j][kj] = i;
        net->deg_ot[j]++;
      }
      
    if (debug) { fprintf(stderr, "done creating synapses.\n"); }

    /* Free auxiliary storage: */
    free(inhib);
    free(pre);
    free(pos);
  }

/* ---------------------------------------------------------------------- */

void nmsim_elem_net_define_neuron
  ( nmsim_elem_net_t *enet, 
    nmsim_elem_neuron_ix_t ine,
    nmsim_elem_neuron_t *neu
  );
  /* Stores into the network {enet} the attributes {neu}
    of the neuron with index {ine}.
    
    The neuron is stored into {enet->neu[ine]}; this vector
    must have been allocated with {enet.nne} elements.
    The index and the fields {ing,nse_out,ise_ini} of {neu} are 
    checked only for range, which must be {0..enet.nne-1},
    {0..gnet.nng-1}, {0..enet.nse-1}, and {0..enet.nse},
    respectively. */
    
void nmsim_elem_net_define_synapse
  ( nmsim_elem_net_t *enet,
    nmsim_elem_synapse_ix_t ise,
    nmsim_elem_synapse_t *syn
  );
  /* Stores into the network {enet} the attributes
    {*syn} of the synapse with index {ise}.
    
    The synapse is stored into {enet->syn[ise]}; this vector 
    must have been allocated with {enet.nse} elements.
    The caller must ensure that, once all synapses
    have been defined, they are in non-decreasing order of 
    the pre-synaptic index {.ine_pre}. */

void nmsim_elem_net_define_neuron
  ( nmsim_elem_net_t *enet, 
    nmsim_elem_neuron_ix_t ine,
    nmsim_elem_neuron_t *neu
  )
  { demand((ine >= 0) && (ine < enet->nne), "invalid neuron index");
    demand((neu->ing >= 0) && (neu->ing < enet->gnet->nng), "invalid neuron group");
    demand((neu->nse_out >= 0) && (neu->nse_out < enet->nse), "invalid output synapse count");
    demand((neu->ise_ini >= 0) && (neu->ise_ini <= enet->nse), "invalid output synapse index");
    enet->neu[ine] = (*neu);
  }
  
void nmsim_elem_net_define_synapse
  ( nmsim_elem_net_t *enet,
    nmsim_elem_synapse_ix_t ise,
    nmsim_elem_synapse_t *syn
  )
  { demand((ise >= 0) && (ise < eset->nse), "invalid synapse index");
    demand((syn->isg >= 0) && (syn->isg < enet->gnet->nsg), "invalid synapse group index");
    demand((syn->ine_pre >= 0) && (syn->ine_pre < enet->nne), "invalid pre-synaptic neuron index");
    demand((syn->ine_pos >= 0) && (syn->ine_pos < enet->nne), "invalid post_synaptic neuron index");
    enet->syn[ise] = (*syn);
  }

/* ---------------------------------------------------------------------- */
nmsim_group_neuron_ix_t nmsim_group_net_add_neuron_group
  ( nmsim_group_net_t *gnet, 
    nmsim_group_neuron_t *ngrp
  );
  /* Adds to the network {gnet} a new neuron group with attributes
    {*ngrp}. The population is stored in the {gnet->ngrp} vector, which
    is expanded if necessary, and incrementng {gnet->nng}. Returns the
    assigned population index {ing}. */
    
nmsim_group_synapse_ix_t nmsim_group_net_add_synapse_group
  ( nmsim_group_net_t *gnet,
    nmsim_group_synapse_t *sgrp
  );
  /* Adds to the network {gnet} a new synaptic bundle with attributes
    {*sgrp}. The bundle is stored in the {gnet->sgrp} vector, which is
    expanded if necessary, and incrementng {gnet->nsg}. Returns the
    assigned bundle index {isg}. */

nmsim_group_neuron_ix_t nmsim_group_net_add_neuron_group
  ( nmsim_group_net_t *gnet, 
    nmsim_group_neuron_t *ngrp
  )
  { demand((ngrp->inc >= 0) && (ngrp->inc < gnet->cnet->nnc), "invalid neuron class");
    demand((ngrp->nne >= 0) && (ngrp->nne <= nmsim_elem_neuron_count_MAX), "invalid neuron count");
    nmsim_group_neuron_ix_t ing = gnet->nng;
    demand(ing <= nmsim_group_neuron_ix_MAX, "too many neuron  groups");
    nmsim_group_neuron_vec_expand(&(gnet->ngrp), ing);
    gnet->ngrp.e[ing] = (*ngrp);
    (gnet->nng)++;
    return ing;
  }
  
nmsim_group_synapse_ix_t nmsim_group_net_add_synapse_group
  ( nmsim_group_net_t *gnet,
    nmsim_group_synapse_t *sgrp
  )
  { demand((sgrp->isc >= 0) && (sgrp->isc < gnet->cnet->nsc), "invalid synapse class index");
    demand((sgrp->ing_pre >= 0) && (sgrp->ing_pre < gnet->nng), "invalid pre-synaptic neuron group index");
    demand((sgrp->ing_pos >= 0) && (sgrp->ing_pos < gnet->nng), "invalid post_synaptic neuron group index");
    demand(sgrp->K >= 0.0, "invalid in-degree");
    nmsim_group_synapse_ix_t isg = gnet->nsg;
    demand(isg <= nmsim_group_synapse_ix_MAX, "too many synapse groups");
    nmsim_group_synapse_vec_expand(&(gnet->sgrp), isg);
    gnet->sgrp.e[isg] = (*sgrp);
    (gnet->nsg)++;
    return isg;
  }
nmsim_class_neuron_ix_t nmsim_class_net_add_neuron_class
  ( nmsim_class_net_t *cnet, 
    nmsim_class_neuron_t *nclass
  );
  /* Adds to the network {cnet} a new neuron class, described by the
    given neuron parameter record {*nclass}. Assigns to those
    parameters the next neuron class index {inc}and stores {parm} in
    the {cnet->nclass} vector, expanding it if needed, and incrementng
    {cnet->nnc}. Returns the class index {inc}. */

nmsim_class_synapse_ix_t nmsim_class_net_add_synapse_class
  ( nmsim_class_net_t *cnet, 
    nmsim_class_synapse_t *sclass
  );
  /* Adds to the network {cnet} a new synapse class, described by the
    given synapse parameter record {*sclass}. Assigns to those
    parameters the next synapse class index {isc} and stores {parm} in
    the {cnet->sclass} vector, expanding it if needed, and incrementng
    {cnet->nsc}. Returns the index {isc}. */

nmsim_class_neuron_ix_t nmsim_class_net_add_neuron_class
  ( nmsim_class_net_t *cnet, 
    nmsim_class_neuron_t *nclass
  )
  { nmsim_class_neuron_ix_t inc = cnet->nnc;
    demand(inc <= nmsim_class_neuron_ix_MAX, "too many neuron classes");
    nmsim_class_neuron_ref_vec_expand(&(cnet->nclass), inc);
    cnet->nclass[inc] = nclass;
    (cnet->nnc)++;
    return inc;
  }
  
nmsim_class_synapse_ix_t nmsim_class_net_add_synapse_class
  ( nmsim_class_net_t *cnet, 
    nmsim_class_synapse_t *sclass
  )
  { nmsim_class_synapse_ix_t isc = cnet->nsc;
    demand(isc <= nmsim_class_synapse_ix_MAX, "too many synapse classes");
    nmsim_class_synapse_ref_vec_expand(&(cnet->sclass), isc);
    cnet->sclass[isc] = sclass;
    (cnet->nsc)++;
    return isc;
  }