#! /usr/bin/gawk -f
# Last edited on 2020-12-24 16:09:56 by jstolfi

BEGIN \
  {
    # Converts Nilton Kamiji's network provided to Jorge Stolfi on 2020-12-03.
    
    # User must specify with "-v" the variable {odfname},
    # the name of the table with the outdegree of each neuron.
     
    # User must also specify with "-v" the variable {dead},
    # The length of the fake refractory period (recovery time 
    # for the recharge modulator).  If {dead} is zero, there
    # will be no refractory period.
    
    # The network is based on the Potjans e Diesmann 2015 layer model of
    # a cortical column. It has 4 layers, with two populations of neurons
    # in each class, one excitatory, one inhibitory.

    abort = -1;
    
    if (odfname == "") { arg_error("must define {odfname}"); }
    if (dead == "") { arg_error("must define {dead}"); } else { dead += 0.0 }

    nnc = 8   # Number of neuron classes.
    nsc = 64  # Number of synapse classes.
    nng = 8   # Number of neuron groups.
    nsg = 64  # Number of synapse groups.

    # Define the number {nne_per_ng[ing]} of neurons per group {ing}
    # and total number of neurons {nne}.
    split("20683 5834 21915 5479 4850 1065 14395 2948", nne_per_ng)
    ashift(nng, nne_per_ng)
    nne = 0;
    for (ing = 0; ing < nng; ing++)
      { printf "group %d has %d neurons\n", ing, nne_per_ng[ing] > "/dev/stderr"
        nne += nne_per_ng[ing];
      }
    printf "\n" > "/dev/stderr"
    printf "network has  %d neurons\n", nne > "/dev/stderr"

    # Define average and dev synapse weight for each class
    # {wt_avg_per_sc[isc], wt_dev_per_sc[isc]}.  These will be scaled
    # by 0.5 becaus of Nilton's 
    
    split \
      ( \
        ( "+0.3512 +0.3512 +0.3512 +0.3512 +0.3512 +0.3512 +0.3512 +0.3512 " \
          "-1.4048 -1.4048 -1.4049 -1.4049 -1.4047 -1.4049 -1.4046 -1.4043 " \
          "+0.7024 +0.3512 +0.3512 +0.3512 +0.3512 +0.3513 +0.3512 +0.3512 " \
          "-1.4047 -1.4047 -1.4048 -1.4047 -1.4054 -1.4025 -1.4049 -1.4043 " \
          "+0.3512 +0.3512 +0.3512 +0.3512 +0.3512 +0.3512 +0.3512 +0.3513 " \
          "00.0000 00.0000 -1.4026 00.0000 -1.4048 -1.4051 -1.4050 -1.4047 " \
          "+0.3512 +0.3514 +0.3512 +0.3512 +0.3512 +0.3512 +0.3512 +0.3512 " \
          "00.0000 00.0000 00.0000 00.0000 00.0000 00.0000 -1.4047 -1.4047 " \
        ),  \
        wt_avg_per_sc \
      )
    ashift(nsc, wt_avg_per_sc)
    ahalf(nsc, wt_avg_per_sc)
    split \
      ( \
        ( "0.03530 0.03530 0.03529 0.03530 0.03529 0.03530 0.03530 0.03530 " \
          "0.14119 0.14118 0.14119 0.14120 0.14117 0.14119 0.14117 0.14115 " \
          "0.07059 0.03529 0.03529 0.03529 0.03530 0.03530 0.03529 0.03530 " \
          "0.14118 0.14117 0.14118 0.14117 0.14124 0.14093 0.14119 0.14113 " \
          "0.03529 0.03530 0.03530 0.03529 0.03530 0.03530 0.03530 0.03530 " \
          "0.00000 0.00000 0.14095 0.00000 0.14118 0.14121 0.14121 0.14118 " \
          "0.03529 0.03531 0.03529 0.03529 0.03529 0.03530 0.03530 0.03530 " \
          "0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.14118 0.14117 " \
        ),  \
        wt_dev_per_sc \
      )
    ashift(nsc, wt_dev_per_sc)
    ahalf(nsc, wt_dev_per_sc)

    # Define total synapses {nse} and properties per synapse group {isg}, namely
    # number {nse_per_sg[isg]} of synapses.
    split \
      ( \
        ( "45499805 17443694  3503670  8114254 10613575  1241436  4681225  2260836 " \
          "22323577  5018763   756561    92832  1817058   169424   556108    17207 " \
          "20253647  4105338 24482849  9933538  5507804   607667  6727570   220033 " \
          " 9670918  1690074 17413576  5223272   151900    12851  1320234     8078 " \
          " 3293578  2221213   714524    87836  2040738   319602  4112225   401638 " \
          "       0        0     7003        0  2407889   430444   305029    25218 " \
          " 2271404   353461 14624432  8810905  1438969   132414  8372649  2888426 " \
          "       0        0        0        0        0        0 10827677  1354320 " \
        ),  \
        nse_per_sg \
      )
    ashift(nsg, nse_per_sg)

    # Define the average outdegree of each synaptic bundle (number of 
    # synapses in bundle divided by the size of the pre-synaptic neuron
    # group).  This data is currently not used since we read the actual 
    # synapses from the input file. But it is useful as a consistency check.
    split \
      ( \
        ( "2199.86  843.38  169.40  392.32  513.15   60.02  226.33  109.31 " \
          "3826.46  860.26  129.68   15.91  311.46   29.04   95.32    2.95 " \
          " 924.19  187.33 1117.17  453.28  251.33   27.73  306.98   10.04 " \
          "1765.09  308.46 3178.24  953.33   27.72    2.35  240.96    1.47 " \
          " 679.09  457.98  147.32   18.11  420.77   65.90  847.88   82.81 " \
          "   0.00    0.00    6.58    0.00 2260.93  404.17  286.41   23.68 " \
          " 157.79   24.55 1015.94  612.08   99.96    9.20  581.64  200.65 " \
          "   0.00    0.00    0.00    0.00    0.00    0.00 3672.89  459.40 " \
        ),  \
        odg_per_sg \
      )
    ashift(nsg, odg_per_sg)
      
    nse = 0;
    for (isg = 0; isg < nsg; isg++)
      { printf "group %d has %d synapses\n", isg, nse_per_sg[isg] > "/dev/stderr"
        nse += nse_per_sg[isg];
      }
    printf "\n" > "/dev/stderr"
    printf "network has  %d synapses\n", nse > "/dev/stderr"
    
    # Read the table of outdegrees {nse_out_per_ne[0..nne-1]}:
    split("", nse_out_per_ne)
    read_nse_out_table(nne, nse, odfname, nse_out_per_ne)

    # Write the nmsim header:
    printf "begin nmsim_elem_net (format of 2020-12-10)\n"
    
    # Write the element counts:
    printf "  neuron_elems = %d\n", nne
    printf "  synapse_elems = %d\n", nse

    # Write the group-level network:
    write_group_net( \
        nnc,"N",dead, \
        nsc,wt_avg_per_sc,wt_dev_per_sc, \
        nng,nne_per_ng, \
        nsg,nse_per_sg \
      )

    # Write the individual neurons, and note their group in {ing_per_ne}:
    split("", ing_per_ne);
    ine = 0;
    for (ing = 0; ing < nng; ing++)
      { nne_ing = nne_per_ng[ing];
        for (ine_ing = 0; ine_ing < nne_ing; ine_ing++)
          { printf "  %d %d %d\n", ine, ing, nse_out_per_ne[ine]
            ing_per_ne[ine] = ing;
            ine++;
          }
      }

    # Start converting the synapses:
    ise = 0
    nse_read = 0;
    split("", nse_per_sg_read)
    split("", wt_tot_per_sc_read)
    split("", wtd2_tot_per_sc_read)
  }
  
(abort >= 0) { exit(abort); }

/N_layers/ { next; }

/^ *[0-9]+ +[0-9]+ *$/ { next; }

/^ *$/ { next; }

/^ *[0-9]+[,] *[0-9]+[,] *[-.0-9]+[,] *[.0-9]+ *$/ \
  {
    gsub(/^ +/, "", $0);
    gsub(/ +$/, "", $0);
    gsub(/ *[,] */, " ", $0);
    ine_pre = $1 - 1;
    ine_pos = $2 - 1;
    wt_ise = 0.500*($3 + 0.0)/250.0; # Weight (mV).
    dt_ise = $4 + 0.0; # Delay (ms)

    # Determne the synapse group, count:
    ing_pre = ing_per_ne[ine_pre];
    ing_pos = ing_per_ne[ine_pos];
    isg = nng*ing_pre + ing_pos
    isc = isg
    nse_per_sg_read[isg] ++;
    wt_tot_per_sc_read[isc] += wt_ise;
    wtd_ise = wt_ise - wt_avg_per_sc[isc];
    wtd2_tot_per_sc_read[isc] += wtd_ise*wtd_ise;
    nse_read++;

    # Ignore delay for now:
    printf "  %d %d %d %d  %.4f\n", ise, isg, ine_pre, ine_pos, wt_ise
    ise++;
    next
  }

// { printf "** line %d: bad format\n  [[%s]]\n", FNR, $0 > "/dev/stderr"; exit(1); }

END \
  {
    if (abort >= 0) { exit(abort); }
    
    printf "end nmsim_elem_net\n"

    # Print synapse counts per group and total:
    for (isg = 0; isg < nsg; isg++)
      { print_compare_count(("synapse group " isg), nse_per_sg[isg], nse_per_sg_read[isg])
        printf "\n" > "/dev/stderr"
      }
    printf "\n" > "/dev/stderr"
    print_compare_count("total synapses", nse, nse_read)
    printf "\n" > "/dev/stderr"
    printf "\n" > "/dev/stderr"

    # Print synapse fanout, weight average and deviation per group:
    for (isc = 0; isc < nsc; isc++)
      { printf "synapse group %d:", isc > "/dev/stderr"
        isg = isc
        nse_isg = nse_per_sg_read[isg]
        # Observed fanout factor:
        ing_pos = isg % nng;
        ing_pre = int((isg - ing_pos)/nng)
        odg_isg = nse_isg/nne_per_ng[ing_pre]
        print_compare_float(" K avg", odg_per_sg[isg], odg_isg, 0.01);
        # Observed average weight:
        wt_tot_isc_read = wt_tot_per_sc_read[isc];
        wt_avg_isc_read = (nse_isg <= 0 ? 0.0 : wt_tot_isc_read/nse_isg)
        print_compare_float(" wt avg", wt_avg_per_sc[isc], wt_avg_isc_read, 0.0002);
        # Observed weight deviation:
        wtd2_tot_isc_read = wtd2_tot_per_sc_read[isc];
        wt_dev_isc_read = (nse_isg <= 1 ? 0.0 : sqrt(wtd2_tot_isc_read/(nse_isg - 1)))
        print_compare_float(" wt dev", wt_dev_per_sc[isc], wt_dev_isc_read, 0.02);
        printf "\n" > "/dev/stderr"
      }
    printf "\n" > "/dev/stderr"
  }

function ashift(n,arr, i) 
  { # Assumes {arr} is an array with elements indexed {1..n}.
    # Shifts them so that they are indexed {0..n-1}.
    for (i = 0; i < n; i++) { arr[i] = arr[i+1]; }
    delete arr[n];
  }

function ahalf(n,arr, i) 
  { # Assumes {arr} is an array with elements indexed {0..n-1}.
    # Multiplies every element by 0.5.
    for (i = 0; i < n; i++) { arr[i] = 0.500 * arr[i]; }
  }
  
function print_compare_count(label,n_exp,n_read)
  {
    printf "%s: ", label > "/dev/stderr"
    printf "expected %d  read %d", n_exp, n_read > "/dev/stderr"
    if (n_exp != n_read) { printf " ** error" > "/dev/stderr" }
  }

function print_compare_float(label,x_exp,x_read,tol,  d)
  {
    printf "%s: ", label > "/dev/stderr"
    printf "expected %10.6f  read %10.6f", x_exp, x_read > "/dev/stderr"
    d = x_exp - x_read;
    if (d < 0) { d = -d; }
    if (d > tol) { printf " ** error" > "/dev/stderr" }
  }

function write_group_net \
  ( nnc,phiclass,dead, \
    nsc,wt_avg_per_sc,wt_dev_per_sc, \
    nng,nne_per_ng, \
    nsg,nse_per_sg, \
    ing,inc_ing,nne_ing, \
    isg,ing_pre,ing_pos \
  )
  {
    printf "  begin nmsim_group_net (format of 2020-12-10)\n"

    printf "    neuron_groups = %d\n", nng
    printf "    synapse_groups = %d\n", nsg

    write_class_net(nnc,phiclass,dead,nsc,wt_avg_per_sc,wt_dev_per_sc);

    # Write neuron groups:
    for (ing = 0; ing < nng; ing++)
      { inc_ing = ing; # Each neuron group is one class.
        nne_ing = nne_per_ng[ing];
        nsg_out = nng; # Bundles to all neurons (but some bundles have zero synapses).
        printf "      %d %d %d %d\n", ing, inc_ing, nne_ing, nsg_out
      }

    if (nsg > 0)
      { isg = 0;
        for (ing_pre = 0; ing_pre < nng; ing_pre++)
          { for (ing_pos = 0; ing_pos < nng; ing_pos++)
              { isc_isg = isg; # Each synapse group is one synapse class.
                nse_g = nse_per_sg[isg]
                printf "      %d %d  %d %d  %d\n", isg, isc_isg, ing_pre, ing_pos, nse_g
                isg++
              }
          }
        if (isg != nsg) { printf "** bug {nsg}\n" > "/dev/stderr"; exit(1); }
      }
    printf "  end nmsim_group_net\n"
  }

function write_class_net(nnc,phiclass,dead,nsc,wt_avg_per_sc,wt_dev_per_sc,  inc,isc)
  {
    # Write the class-level network:
    printf "    begin nmsim_class_net (format of 2020-12-10)\n"

    printf "      neuron_classes = %d\n", nnc
    printf "      synapse_classes = %d\n", nsc

    for (inc = 0; inc < nnc; inc++)
      { write_neuron_class(inc,dead,phiclass) }

    for (isc = 0; isc < nsc; isc++)
      { write_synapse_class(isc,wt_avg_per_sc[isc],wt_dev_per_sc[isc]) }

     printf "    end nmsim_class_net\n"
  }

function  write_neuron_class(inc,dead,phiclass)
  {
    printf "      neuron_class = %d\n", inc
    printf "      begin nmsim_class_neuron (format of 2019-01-08)\n"
    printf "        V_R = 0\n"
    printf "        V_B = 0\n"
    printf "        V_tau = 10\n"
    if (dead <= 0)
      { # No refractory period: */
        printf "        M_R = 1\n"
        printf "        M_tau = 0\n"
      }
    else
      { # With fake refractory period: */
        printf "        M_R = 0\n"
        printf "        M_tau = %.8f\n", dead
      }
    printf "        H_R = 1\n"
    printf "        H_tau = 0\n"
    printf "        Phi_class = %s\n", phiclass
    printf "        V_M = +20\n"
    printf "        V_D = +5\n"
    printf "      end nmsim_class_neuron\n"
  }

function  write_synapse_class(isc,wt_avg_isc,wt_dev_isc)
  {
    printf "      synapse_class = %d\n", isc
    printf "      begin nmsim_class_synapse (format of 2019-01-10)\n"
    printf "        W_avg = %+7.4f\n", wt_avg_isc
    printf "        W_dev = %4.2f\n", wt_dev_isc
    printf "      end nmsim_class_synapse\n"
  }

function read_nse_out_table(nne,nse,fname,tbl,   ine,nlin,lin,fld,nfld) 
  {
    ine=0; # Number of table entries read.
    nlin=0; # Number of lines read (including blanks and comments).
    while((getline lin < fname) > 0) { 
      nlin++;
      if (! match(lin, /^[ \011]*([#]|$)/))
        { nfld = split(lin, fld, " ");
          if ((nfld >= 3) && (fld[3] ~ /^[#]/)) { nfld = 2; }
          if (nfld != 2) { tbl_error(fname, nlin, ("bad table entry = \"" lin "\"")); }
          if ((fld[1] + 0) != ine) { tbl_error(fname, nlin, ("invalid neuron index = \"" lin "\"")); }
          tbl[ine] = fld[2] + 0;
          ine++;
        }
    }
    if ((ERRNO != "0") && (ERRNO != "")) { tbl_error(fname, nlin, ERRNO); }
    close (fname);
    if (nlin == 0) { arg_error(("file \"" fname "\" empty or missing")); }
    # printf "loaded %6d map pairs\n", ine > "/dev/stderr"
  }

function tbl_error(f,n,msg)
  { printf "%s:%d: %s\n", f, n, msg > "/dev/stderr";
    abort = 1;
    exit(1)
  }

function data_error(msg,lin)
  { printf "%s:%d: ** %s\n", FILENAME, FNR, msg > "/dev/stderr";
    if (lin != "") { printf "  [[%s]]\n", lin > "/dev/stderr"; }
    abort = 1;
    exit(1)
  }

function arg_error(msg)
  { printf "** %s\n", msg > "/dev/stderr";
    abort = 1;
    exit(1)
  }