auditory: add auditory simulation file

lib/CMakeLists.txt

@@ -124,4 +124,5 @@ endfunction()
 if (NOT BUILD_SHARED_LIBS)
   irritator_add_test(test-api test/public-api.cpp)
   irritator_add_test(simulations test/simulations.cpp)
+  irritator_add_test(auditory test/auditory.cpp)
 endif ()

lib/test/auditory.cpp

+// Copyright (c) 2020 INRA Distributed under the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#include <irritator/core.hpp>
+
+#include <boost/ut.hpp>
+
+#include <fmt/format.h>
+
+#include <cstdio>
+
+#include <fstream>
+
+using namespace std;
+
+static void
+dot_graph_save(const irt::simulation& sim, std::FILE* os)
+{
+    /* With input and output port.
+
+    digraph graphname{
+        graph[rankdir = "LR"];
+        node[shape = "record"];
+        edge[];
+
+        "sum_a"[label = "sum-a | <f0> | <f1>"];
+
+        "sum_a":f0->int_a[id = 1];
+        sum_b->int_b[label = "2-10"];
+        prod->sum_b[label = "3-4"];
+        prod -> "sum_a":f0[label = "3-2"];
+        int_a->qua_a[label = "4-11"];
+        int_a->prod[label = "4-5"];
+        int_a -> "sum_a":f1[label = "4-1"];
+        int_b->qua_b[label = "5-12"];
+        int_b->prod[label = "5-6"];
+        int_b->sum_b[label = "5-3"];
+        qua_a->int_a[label = "6-7"];
+        qua_b->int_b[label = "7-9"];
+    }
+    */
+
+    !boost::ut::expect(os != nullptr);
+
+    std::fputs("digraph graphname {\n", os);
+
+    irt::output_port* output_port = nullptr;
+    while (sim.output_ports.next(output_port)) {
+        for (const irt::input_port_id dst : output_port->connections) {
+            if (auto* input_port = sim.input_ports.try_to_get(dst);
+                input_port) {
+                auto* mdl_src = sim.models.try_to_get(output_port->model);
+                auto* mdl_dst = sim.models.try_to_get(input_port->model);
+
+                if (!(mdl_src && mdl_dst))
+
+                    continue;
+                if (mdl_src->name.empty())
+                    fmt::print(os, "{} -> ", irt::get_key(output_port->model));
+                else
+                    fmt::print(os, "{} -> ", mdl_src->name.c_str());
+
+                if (mdl_dst->name.empty())
+                    fmt::print(os, "{}", irt::get_key(input_port->model));
+                else
+                    fmt::print(os, "{}", mdl_dst->name.c_str());
+
+                std::fputs(" [label=\"", os);
+
+                if (output_port->name.empty())
+                    fmt::print(
+                      os,
+                      "{}",
+                      irt::get_key(sim.output_ports.get_id(*output_port)));
+                else
+                    fmt::print(os, "{}", output_port->name.c_str());
+
+                std::fputs("-", os);
+
+                if (input_port->name.empty())
+                    fmt::print(
+                      os,
+                      "{}",
+                      irt::get_key(sim.input_ports.get_id(*input_port)));
+                else
+                    fmt::print(os, "{}", input_port->name.c_str());
+
+                std::fputs("\"];\n", os);
+            }
+        }
+    }
+}
+
+/**
+ * Reads csv file into table, exported as a vector of vector of doubles.
+ * @param inputFileName input file name (full path).
+ * @return data as vector of vector of doubles.
+ */
+vector<vector<double>> parse2DCsvFile(string inputFileName) {
+ 
+    vector<vector<double> > data;
+    ifstream inputFile(inputFileName);
+    int l = 0;
+ 
+    while (inputFile) {
+        l++;
+        string s;
+        if (!getline(inputFile, s)) break;
+        if (s[0] != '#') {
+            istringstream ss(s);
+            vector<double> record;
+ 
+            while (ss) {
+                string line;
+                if (!getline(ss, line, ','))
+                    break;
+                try {
+                    record.push_back(stof(line));
+                }
+                catch (const std::invalid_argument e) {
+                    cout << "NaN found in file " << inputFileName << " line " << l
+                         << endl;
+                    e.what();
+                }
+            }
+ 
+            data.push_back(record);
+        }
+    }
+ 
+    if (!inputFile.eof()) {
+        cerr << "Could not read file " << inputFileName << "\n";
+        __throw_invalid_argument("File not found.");
+    }
+ 
+    return data;
+}
+
+// Global data
+vector<vector<double>> sound_data = parse2DCsvFile("output_cochlea.csv");
+double samplerate = 44100.0;
+
+struct neuron {
+  irt::dynamics_id  sum;
+  irt::dynamics_id  prod;
+  irt::dynamics_id  integrator;
+  irt::dynamics_id  quantifier;
+  irt::dynamics_id  constant;
+  irt::dynamics_id  flow;
+  irt::dynamics_id  cross;
+  irt::dynamics_id  constant_cross;
+};
+
+/*typedef double Fn(double);
+
+template<long unsigned int i>
+double
+time_f(double t) noexcept
+{
+  return sound_data[i][static_cast<int>(t*samplerate)];
+}*/
+
+struct neuron 
+make_neuron(irt::simulation* sim, long unsigned int i) noexcept
+{
+  using namespace boost::ut;
+  double tau_lif = 1.5*0.001;
+  double Vr_lif = 0.0;
+  double Vt_lif = 1.0;
+  /*double ref_lif = 0.5*1e-3;
+  double sigma_lif = 0.02;*/
+
+
+
+
+  auto& sum_lif = sim->adder_2_models.alloc();
+  auto& prod_lif = sim->adder_2_models.alloc();
+  auto& integrator_lif = sim->integrator_models.alloc();
+  auto& quantifier_lif = sim->quantifier_models.alloc();
+  auto& constant_lif = sim->constant_models.alloc();
+  auto& flow_lif = sim->flow_models.alloc();
+  auto& constant_cross_lif = sim->constant_models.alloc();
+  auto& cross_lif = sim->cross_models.alloc();
+
+
+  sum_lif.default_input_coeffs[0] = -1.0;
+  sum_lif.default_input_coeffs[1] = 1.0;
+  
+  prod_lif.default_input_coeffs[0] = 1.0/tau_lif;
+  prod_lif.default_input_coeffs[1] = 0.0;
+
+  constant_lif.default_value = 1.0;
+  flow_lif.default_data = sound_data[i+1];
+  flow_lif.default_samplerate = samplerate;
+  constant_cross_lif.default_value = Vr_lif;
+  
+  integrator_lif.default_current_value = 0.0;
+
+  quantifier_lif.default_adapt_state =
+    irt::quantifier::adapt_state::possible;
+  quantifier_lif.default_zero_init_offset = true;
+  quantifier_lif.default_step_size = 0.1;
+  quantifier_lif.default_past_length = 3;
+
+  cross_lif.default_threshold = Vt_lif;
+
+  char crosslif[7];char ctecrosslif[7];char intlif[7];char flowlif[7];
+  char quantlif[7];char sumlif[7];char prodlif[7];char ctelif[7];
+
+  snprintf(crosslif, 7,"croli%ld", i);snprintf(ctecrosslif, 7,"ctcli%ld", i);snprintf(intlif, 7,"intli%ld", i);
+  snprintf(quantlif, 7,"quali%ld", i);snprintf(sumlif, 7,"sumli%ld", i);snprintf(prodlif, 7,"prdli%ld", i);
+  snprintf(ctelif, 7,"cteli%ld", i);snprintf(flowlif, 7,"flwli%ld", i);
+  
+
+  sim->alloc(sum_lif, sim->adder_2_models.get_id(sum_lif), sumlif);
+  sim->alloc(prod_lif, sim->adder_2_models.get_id(prod_lif), prodlif);
+  sim->alloc(integrator_lif, sim->integrator_models.get_id(integrator_lif), intlif);
+  sim->alloc(quantifier_lif, sim->quantifier_models.get_id(quantifier_lif), quantlif);
+  sim->alloc(constant_lif, sim->constant_models.get_id(constant_lif), ctelif);
+  sim->alloc(flow_lif, sim->flow_models.get_id(flow_lif), flowlif);
+  sim->alloc(cross_lif, sim->cross_models.get_id(cross_lif), crosslif);
+  sim->alloc(constant_cross_lif, sim->constant_models.get_id(constant_cross_lif), ctecrosslif);
+
+  struct neuron neuron_model = {sim->adder_2_models.get_id(sum_lif),
+                                sim->adder_2_models.get_id(prod_lif),
+                                sim->integrator_models.get_id(integrator_lif),
+                                sim->quantifier_models.get_id(quantifier_lif),
+                                sim->constant_models.get_id(constant_lif),
+                                sim->flow_models.get_id(flow_lif),
+                                sim->cross_models.get_id(cross_lif),                                                                                                
+                                sim->constant_models.get_id(constant_cross_lif),                                                                
+                                }; 
+
+
+  // Connections
+  expect(sim->connect(quantifier_lif.y[0], integrator_lif.x[0]) ==
+          irt::status::success);
+  expect(sim->connect(prod_lif.y[0], integrator_lif.x[1]) ==
+          irt::status::success);
+  expect(sim->connect(cross_lif.y[0], integrator_lif.x[2]) ==
+          irt::status::success);
+  expect(sim->connect(cross_lif.y[0], quantifier_lif.x[0]) ==
+          irt::status::success);
+  expect(sim->connect(cross_lif.y[0], sum_lif.x[0]) ==
+          irt::status::success);
+  expect(sim->connect(integrator_lif.y[0],cross_lif.x[0]) ==
+          irt::status::success);     
+  expect(sim->connect(integrator_lif.y[0],cross_lif.x[2]) ==
+          irt::status::success);
+  expect(sim->connect(constant_cross_lif.y[0],cross_lif.x[1]) ==
+          irt::status::success);     
+  expect(sim->connect(flow_lif.y[0], sum_lif.x[1]) ==
+          irt::status::success);  
+  expect(sim->connect(sum_lif.y[0],prod_lif.x[0]) ==
+          irt::status::success);   
+  expect(sim->connect(constant_lif.y[0],prod_lif.x[1]) ==
+          irt::status::success);
+  return neuron_model;
+}
+
+
+int
+main()
+{
+    using namespace boost::ut;
+
+   
+    "laudanski_1_simulation"_test = [] {
+        irt::simulation sim;
+
+        // Neuron constants
+        long unsigned int N = sound_data.size() - 1;
+        
+
+
+        expect(irt::is_success(sim.init(512lu, 8192lu)));
+
+
+        
+        // Neurons
+        std::vector<struct neuron> first_layer_neurons;
+        for (long unsigned int i = 0 ; i < N; i++) {
+
+          struct neuron neuron_model = make_neuron(&sim,i);
+          first_layer_neurons.emplace_back(neuron_model);
+        } 
+
+
+        dot_graph_save(sim, stdout);
+
+        irt::time t = 0.0;
+        std::FILE* os = std::fopen("output_laudanski.csv", "w");
+        !expect(os != nullptr);
+        
+        std::string s = "t,";
+        for (long unsigned int i = 0; i < N; i++)
+        {
+          s =  s + "Neuron" + std::to_string(i)
+                + ",";
+        }
+        fmt::print(os, s + "\n");
+
+        expect(irt::status::success == sim.initialize(t));
+
+        do {
+
+            irt::status st = sim.run(t);
+            expect(st == irt::status::success);
+
+            std::string s = std::to_string(t)+",";
+            for (long unsigned int i = 0; i < N; i++)
+            {
+              //s =  s + std::to_string(sim.cross_models.get(first_layer_neurons[i].cross).event)
+              s =  s + std::to_string(sim.integrator_models.get(first_layer_neurons[i].integrator).last_output_value)
+                    + ",";
+
+            }
+            fmt::print(os, s + "\n");
+
+        } while (t < sound_data[0].size()/samplerate);
+
+        std::fclose(os);
+      };
+
+   
+}
+
+