From 95b3733225d9273890f98316bfebf8b493a075da Mon Sep 17 00:00:00 2001
From: Dario Izzo <dario.izzo@gmail.com>
Date: Tue, 31 Oct 2023 13:55:47 +0100
Subject: [PATCH 01/14] ffnn first impl

---
 include/heyoka/model/ffnn.hpp |  36 ++++++++++
 src/model/ffnn.cpp            | 124 ++++++++++++++++++++++++++++++++++
 2 files changed, 160 insertions(+)
 create mode 100644 include/heyoka/model/ffnn.hpp
 create mode 100644 src/model/ffnn.cpp

diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
new file mode 100644
index 000000000..bd29a21b6
--- /dev/null
+++ b/include/heyoka/model/ffnn.hpp
@@ -0,0 +1,36 @@
+// Copyright 2020, 2021, 2022, 2023 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com)
+//
+// This file is part of the heyoka library.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef HEYOKA_MODEL_FFNN_HPP
+#define HEYOKA_MODEL_FFNN_HPP
+
+#include <tuple>
+#include <utility>
+#include <vector>
+
+#include <heyoka/config.hpp>
+#include <heyoka/detail/visibility.hpp>
+#include <heyoka/expression.hpp>
+
+HEYOKA_BEGIN_NAMESPACE
+
+namespace model
+{
+// This c++ function returns the symbolic expressions of the `n_out` output neurons in a feed forward neural network,
+// as a function of the `n_in` input expressions.
+//
+// The expression will contain the weights and biases of the neural network flattened into `pars` with the following conventions:
+//
+// from the left to right layer of parameters: [flattened weights1, flattened weights2,  ... , biases1, bises2, ...]
+HEYOKA_DLL_PUBLIC std::vector<expression>
+ffnn_impl(const std::vector<expression> & in, unsigned n_out, const std::vector<unsigned> &n_neurons_per_hidden_layer,
+          const std::vector<std::function<expression(const expression &)>> &activations,
+          const std::vector<expression> &pars);
+} // namespace model
+
+HEYOKA_END_NAMESPACE
diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
new file mode 100644
index 000000000..9905e1d09
--- /dev/null
+++ b/src/model/ffnn.cpp
@@ -0,0 +1,124 @@
+// Copyright 2020, 2021, 2022, 2023 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com)
+//
+// This file is part of the heyoka library.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <cmath>
+#include <stdexcept>
+#include <utility>
+#include <variant>
+#include <vector>
+
+#include <fmt/core.h>
+
+#include <heyoka/config.hpp>
+#include <heyoka/expression.hpp>
+#include <heyoka/math/pow.hpp>
+#include <heyoka/math/sqrt.hpp>
+#include <heyoka/model/cr3bp.hpp>
+#include <heyoka/number.hpp>
+
+HEYOKA_BEGIN_NAMESPACE
+
+namespace model
+{
+namespace detail
+{
+// From the i-th neuron of the layer_id and the incoming j-th connection, returns the corresponding weight position in
+// the flattened structure. NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
+std::vector<expression>::size_type flattenw(std::uint32_t i, std::uint32_t j,
+                                            const std::vector<std::uint32_t> &n_neurons, std::uint32_t layer_id)
+{
+    assert(layer_id > 0);
+    // The weight for the jth-neuron of the ith layer will be placed after all previous layers.
+    // We start counting how many in flattened.
+    std::uint32_t counter = 0;
+    for (std::uint32_t k = 1; k < layer_id; ++k) {
+        counter += n_neurons[k] * n_neurons[k - 1];
+    }
+    // We then add the weights used for the previous neurons on the same layer.
+    counter += i * n_neurons[layer_id - 1];
+    // And return the index corresponding to the j-th weight.
+    return counter + j;
+}
+
+// From the i-th neuron of the layer_id, returns the corresponding bias position in the
+// flattened structure.
+// NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
+std::vector<expression>::size_type flattenb(std::uint32_t i, const std::vector<std::uint32_t> &n_neurons,
+                                            // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
+                                            std::uint32_t layer_id, std::uint32_t n_net_w)
+{
+    assert(layer_id > 0);
+    // The weight for the jth-neuron of the ith layer will be placed after all previous layers.
+    // We start counting how many in flattened.
+    std::uint32_t counter = 0;
+    for (std::uint32_t k = 1; k < layer_id; ++k) {
+        counter += n_neurons[k];
+    }
+    // And return the index corresponding to the i-th bias.
+    return counter + i + n_net_w;
+}
+
+std::vector<expression> compute_layer(std::uint32_t layer_id, const std::vector<expression> &inputs,
+                                      const std::vector<std::uint32_t> &n_neurons,
+                                      const std::function<expression(const expression &)> &activation,
+                                      const std::vector<expression> &net_wb, std::uint32_t n_net_w)
+{
+    assert(layer_id > 0);
+    auto n_neurons_prev_layer = boost::numeric_cast<std::uint32_t>(inputs.size());
+    auto n_neurons_curr_layer = n_neurons[layer_id];
+
+    std::vector<expression> retval(n_neurons_curr_layer, 0_dbl);
+    for (std::uint32_t i = 0u; i < n_neurons_curr_layer; ++i) {
+        for (std::uint32_t j = 0u; j < n_neurons_prev_layer; ++j) {
+            retval[i] += net_wb[flattenw(i, j, n_neurons, layer_id)] * inputs[j];
+        }
+        retval[i] += net_wb[flattenb(i, n_neurons, layer_id, n_net_w)];
+        retval[i] = activation(retval[i]);
+    }
+    return retval;
+}
+} // namespace detail
+
+HEYOKA_DLL_PUBLIC std::vector<expression> ffnn_impl(
+    // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
+    const std::vector<expression> &in, std::uint32_t n_out,
+    const std::vector<std::uint32_t> &n_neurons_per_hidden_layer,
+    const std::vector<std::function<expression(const expression &)>> &activations,
+    const std::vector<expression> &net_wb)
+{
+    // Sanity check (should be a throw check?)
+    assert(n_neurons_per_hidden_layer.size() + 1 == activations.size());
+
+    // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons)
+    auto n_hidden_layers = boost::numeric_cast<std::uint32_t>(n_neurons_per_hidden_layer.size());
+    // Number of neuronal layers (counting input and output)
+    auto n_layers = n_hidden_layers + 2;
+    // Number o
+    auto n_in = boost::numeric_cast<std::uint32_t>(in.size());
+    // Number of neurons per neuronal layer
+    std::vector<std::uint32_t> n_neurons = n_neurons_per_hidden_layer;
+    n_neurons.insert(n_neurons.begin(), n_in);
+    n_neurons.insert(n_neurons.end(), n_out);
+    // Number of network parameters
+    std::uint32_t n_net_wb = 0u;
+    std::uint32_t n_net_w = 0u;
+    for (std::uint32_t i = 1u; i < n_layers; ++i) {
+        n_net_wb += n_neurons[i - 1] * n_neurons[i];
+        n_net_w += n_neurons[i - 1] * n_neurons[i];
+        n_net_wb += n_neurons[i];
+    }
+    // Sanity check (should be a throw check?)
+    assert(net_wb.size() == n_net_wb);
+    std::vector<expression> retval{};
+    for (std::uint32_t i = 1u; i < n_layers; ++i) {
+        retval = detail::compute_layer(i, retval, n_neurons, activations[i], net_wb, n_net_w);
+    }
+    return retval;
+}
+} // namespace model
+HEYOKA_END_NAMESPACE
\ No newline at end of file

From b5858c38f224317776fefd5215d9b608ff10928b Mon Sep 17 00:00:00 2001
From: Dario Izzo <dario.izzo@gmail.com>
Date: Tue, 31 Oct 2023 14:04:38 +0100
Subject: [PATCH 02/14] bug fixes

---
 include/heyoka/model/ffnn.hpp |  2 ++
 src/model/ffnn.cpp            | 10 +---------
 2 files changed, 3 insertions(+), 9 deletions(-)

diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
index bd29a21b6..f1c9c616a 100644
--- a/include/heyoka/model/ffnn.hpp
+++ b/include/heyoka/model/ffnn.hpp
@@ -34,3 +34,5 @@ ffnn_impl(const std::vector<expression> & in, unsigned n_out, const std::vector<
 } // namespace model
 
 HEYOKA_END_NAMESPACE
+
+#endif
diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
index 9905e1d09..261cf08e6 100644
--- a/src/model/ffnn.cpp
+++ b/src/model/ffnn.cpp
@@ -7,19 +7,11 @@
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #include <cmath>
-#include <stdexcept>
-#include <utility>
-#include <variant>
 #include <vector>
 
-#include <fmt/core.h>
-
 #include <heyoka/config.hpp>
 #include <heyoka/expression.hpp>
-#include <heyoka/math/pow.hpp>
-#include <heyoka/math/sqrt.hpp>
-#include <heyoka/model/cr3bp.hpp>
-#include <heyoka/number.hpp>
+#include <heyoka/model/ffnn.hpp>
 
 HEYOKA_BEGIN_NAMESPACE
 

From 1d7f8516f0883a567a95d6c3d3ac39db96b59499 Mon Sep 17 00:00:00 2001
From: Dario Izzo <dario.izzo@gmail.com>
Date: Tue, 31 Oct 2023 14:06:10 +0100
Subject: [PATCH 03/14] fix of types and hpps

---
 include/heyoka/model/ffnn.hpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
index f1c9c616a..0c2cc38da 100644
--- a/include/heyoka/model/ffnn.hpp
+++ b/include/heyoka/model/ffnn.hpp
@@ -9,6 +9,7 @@
 #ifndef HEYOKA_MODEL_FFNN_HPP
 #define HEYOKA_MODEL_FFNN_HPP
 
+#include <atomic>
 #include <tuple>
 #include <utility>
 #include <vector>
@@ -28,9 +29,9 @@ namespace model
 //
 // from the left to right layer of parameters: [flattened weights1, flattened weights2,  ... , biases1, bises2, ...]
 HEYOKA_DLL_PUBLIC std::vector<expression>
-ffnn_impl(const std::vector<expression> & in, unsigned n_out, const std::vector<unsigned> &n_neurons_per_hidden_layer,
-          const std::vector<std::function<expression(const expression &)>> &activations,
-          const std::vector<expression> &pars);
+ffnn_impl(const std::vector<expression> &, std::uint32_t, const std::vector<std::uint32_t> &,
+          const std::vector<std::function<expression(const expression &)>> &,
+          const std::vector<expression> &);
 } // namespace model
 
 HEYOKA_END_NAMESPACE

From 0effd03303ca1fee96cb153856883206720d7af7 Mon Sep 17 00:00:00 2001
From: Dario Izzo <dario.izzo@gmail.com>
Date: Tue, 31 Oct 2023 16:28:05 +0100
Subject: [PATCH 04/14] WIP

---
 CMakeLists.txt                |  1 +
 include/heyoka/model/ffnn.hpp |  1 +
 src/model/ffnn.cpp            | 55 ++++++++++++++++++++++++++++++-----
 test/CMakeLists.txt           |  1 +
 4 files changed, 50 insertions(+), 8 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 56e649501..352d0a478 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -261,6 +261,7 @@ set(HEYOKA_SRC_FILES
     "${CMAKE_CURRENT_SOURCE_DIR}/src/model/mascon.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/model/vsop2013.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/model/cr3bp.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/model/ffnn.cpp"
     # Math functions.
     "${CMAKE_CURRENT_SOURCE_DIR}/src/math/kepE.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/math/kepF.cpp"
diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
index 0c2cc38da..e1c0b456a 100644
--- a/include/heyoka/model/ffnn.hpp
+++ b/include/heyoka/model/ffnn.hpp
@@ -28,6 +28,7 @@ namespace model
 // The expression will contain the weights and biases of the neural network flattened into `pars` with the following conventions:
 //
 // from the left to right layer of parameters: [flattened weights1, flattened weights2,  ... , biases1, bises2, ...]
+//
 HEYOKA_DLL_PUBLIC std::vector<expression>
 ffnn_impl(const std::vector<expression> &, std::uint32_t, const std::vector<std::uint32_t> &,
           const std::vector<std::function<expression(const expression &)>> &,
diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
index 261cf08e6..4ef1c60da 100644
--- a/src/model/ffnn.cpp
+++ b/src/model/ffnn.cpp
@@ -6,9 +6,14 @@
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
+#include <algorithm>
 #include <cmath>
+#include <iostream>
 #include <vector>
 
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+
 #include <heyoka/config.hpp>
 #include <heyoka/expression.hpp>
 #include <heyoka/model/ffnn.hpp>
@@ -65,8 +70,13 @@ std::vector<expression> compute_layer(std::uint32_t layer_id, const std::vector<
     auto n_neurons_curr_layer = n_neurons[layer_id];
 
     std::vector<expression> retval(n_neurons_curr_layer, 0_dbl);
+    fmt::print("nneurons: {}", n_neurons_prev_layer);
+    std::cout << std::endl;
+
     for (std::uint32_t i = 0u; i < n_neurons_curr_layer; ++i) {
         for (std::uint32_t j = 0u; j < n_neurons_prev_layer; ++j) {
+            fmt::print("layer, i, j: {}, {}, {}", layer_id, i, j);
+            std::cout << std::endl;
             retval[i] += net_wb[flattenw(i, j, n_neurons, layer_id)] * inputs[j];
         }
         retval[i] += net_wb[flattenb(i, n_neurons, layer_id, n_net_w)];
@@ -76,27 +86,46 @@ std::vector<expression> compute_layer(std::uint32_t layer_id, const std::vector<
 }
 } // namespace detail
 
-HEYOKA_DLL_PUBLIC std::vector<expression> ffnn_impl(
+std::vector<expression> ffnn_impl(
     // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
     const std::vector<expression> &in, std::uint32_t n_out,
     const std::vector<std::uint32_t> &n_neurons_per_hidden_layer,
     const std::vector<std::function<expression(const expression &)>> &activations,
     const std::vector<expression> &net_wb)
 {
-    // Sanity check (should be a throw check?)
-    assert(n_neurons_per_hidden_layer.size() + 1 == activations.size());
+    // Sanity checks
+    if (n_neurons_per_hidden_layer.size() + 1 != activations.size()) {
+        throw std::invalid_argument(fmt::format(
+            "The number of hidden layers, as detected from the inputs, was {}, while"
+            "the number of activation function supplied was {}. A FFNN needs exactly one more activation function "
+            "than the number of hidden layers.",
+            n_neurons_per_hidden_layer.size(), activations.size()));
+    }
+    if (in.empty()) {
+        throw std::invalid_argument("The inputs provided to the ffnn seem to be an empty vector.");
+    }
+    if (n_out == 0) {
+        throw std::invalid_argument("The number of network outputs cannot be zero.");
+    }
+    if (!std::all_of(n_neurons_per_hidden_layer.begin(), n_neurons_per_hidden_layer.end(),
+                     [](std::uint32_t item) { return item > 0; })) {
+        throw std::invalid_argument("The number of neurons for each hidden layer must be greater than zero!");
+    }
+    if (n_neurons_per_hidden_layer.empty()) { // TODO(darioizzo): maybe this is actually a wanted corner case, remove?
+        throw std::invalid_argument("The number of hidden layers cannot be zero.");
+    }
 
     // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons)
     auto n_hidden_layers = boost::numeric_cast<std::uint32_t>(n_neurons_per_hidden_layer.size());
     // Number of neuronal layers (counting input and output)
     auto n_layers = n_hidden_layers + 2;
-    // Number o
+    // Number of inputs
     auto n_in = boost::numeric_cast<std::uint32_t>(in.size());
     // Number of neurons per neuronal layer
     std::vector<std::uint32_t> n_neurons = n_neurons_per_hidden_layer;
     n_neurons.insert(n_neurons.begin(), n_in);
     n_neurons.insert(n_neurons.end(), n_out);
-    // Number of network parameters
+    // Number of network parameters (wb: weights and biases, w: only weights)
     std::uint32_t n_net_wb = 0u;
     std::uint32_t n_net_w = 0u;
     for (std::uint32_t i = 1u; i < n_layers; ++i) {
@@ -104,10 +133,20 @@ HEYOKA_DLL_PUBLIC std::vector<expression> ffnn_impl(
         n_net_w += n_neurons[i - 1] * n_neurons[i];
         n_net_wb += n_neurons[i];
     }
-    // Sanity check (should be a throw check?)
-    assert(net_wb.size() == n_net_wb);
-    std::vector<expression> retval{};
+    // Sanity check
+    if (net_wb.size() != n_net_wb) {
+        throw std::invalid_argument(fmt::format(
+            "The number of network parameters, detected from its structure to be {}, does not match the size of"
+            "the corresponding expressions {} ",
+            n_net_wb, net_wb.size()));
+    }
+
+    // Now we build the expressions recursively going from layer to layer (L = f(Wx+b)))
+
+    std::vector<expression> retval = in;
     for (std::uint32_t i = 1u; i < n_layers; ++i) {
+        fmt::print("{},{}", i, n_neurons[i]);
+        std::cout << std::endl;
         retval = detail::compute_layer(i, retval, n_neurons, activations[i], net_wb, n_net_w);
     }
     return retval;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 7463f3ab1..2a3874e3e 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -140,6 +140,7 @@ ADD_HEYOKA_TESTCASE(model_fixed_centres)
 ADD_HEYOKA_TESTCASE(model_rotating)
 ADD_HEYOKA_TESTCASE(model_mascon)
 ADD_HEYOKA_TESTCASE(model_cr3bp)
+ADD_HEYOKA_TESTCASE(model_ffnn)
 ADD_HEYOKA_TESTCASE(step_callback)
 ADD_HEYOKA_TESTCASE(llvm_state_mem_cache)
 

From 67b8b448e662b964aa4db2e0b29a32735d8808ea Mon Sep 17 00:00:00 2001
From: Dario Izzo <dario.izzo@gmail.com>
Date: Tue, 31 Oct 2023 16:32:30 +0100
Subject: [PATCH 05/14] WIP

---
 test/model_ffnn.cpp | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)
 create mode 100644 test/model_ffnn.cpp

diff --git a/test/model_ffnn.cpp b/test/model_ffnn.cpp
new file mode 100644
index 000000000..b09c6a113
--- /dev/null
+++ b/test/model_ffnn.cpp
@@ -0,0 +1,27 @@
+// Copyright 2020, 2021, 2022, 2023 Francesco Biscani (bluescarni@gmail.com), Dario Izzo (dario.izzo@gmail.com)
+//
+// This file is part of the heyoka library.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+
+#include <heyoka/config.hpp>
+#include <heyoka/expression.hpp>
+#include <heyoka/math.hpp>
+#include <heyoka/model/ffnn.hpp>
+
+#include "catch.hpp"
+
+using namespace heyoka;
+
+TEST_CASE("impl")
+{
+    auto [x] = make_vars("x");
+    auto my_net = model::ffnn_impl({x}, 2, {2, 2}, {heyoka::tanh, heyoka::tanh, [](auto ret) { return ret; }},
+                                   {1_dbl, 2_dbl, 3_dbl, 4_dbl, 5_dbl, 6_dbl, 7_dbl, 8_dbl, 9_dbl, 0_dbl, 1_dbl, 2_dbl,
+                                    3_dbl, 4_dbl, 5_dbl, 6_dbl});
+}
\ No newline at end of file

From 79ac7f977a07683e6e850f756eb279e862f68090 Mon Sep 17 00:00:00 2001
From: Dario Izzo <darioizzo@gmail.com>
Date: Tue, 31 Oct 2023 19:12:27 +0100
Subject: [PATCH 06/14] WIP

---
 src/model/ffnn.cpp  | 20 ++++++++++++++++----
 test/model_ffnn.cpp |  3 ++-
 2 files changed, 18 insertions(+), 5 deletions(-)

diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
index 4ef1c60da..a0e2c6420 100644
--- a/src/model/ffnn.cpp
+++ b/src/model/ffnn.cpp
@@ -70,18 +70,30 @@ std::vector<expression> compute_layer(std::uint32_t layer_id, const std::vector<
     auto n_neurons_curr_layer = n_neurons[layer_id];
 
     std::vector<expression> retval(n_neurons_curr_layer, 0_dbl);
-    fmt::print("nneurons: {}", n_neurons_prev_layer);
+    fmt::print("net_wb: {}\n", net_wb.size());
     std::cout << std::endl;
 
     for (std::uint32_t i = 0u; i < n_neurons_curr_layer; ++i) {
         for (std::uint32_t j = 0u; j < n_neurons_prev_layer; ++j) {
-            fmt::print("layer, i, j: {}, {}, {}", layer_id, i, j);
+            fmt::print("layer, i, j, idx: {}, {}, {}\n", layer_id, i, j, flattenw(i, j, n_neurons, layer_id));
             std::cout << std::endl;
-            retval[i] += net_wb[flattenw(i, j, n_neurons, layer_id)] * inputs[j];
+            retval[i] += 1_dbl;//net_wb[flattenw(i, j, n_neurons, layer_id)] * inputs[j];
         }
-        retval[i] += net_wb[flattenb(i, n_neurons, layer_id, n_net_w)];
+        fmt::print("idxb {}\n", flattenb(i, n_neurons, layer_id, n_net_w));
+        std::cout << std::endl;
+
+        retval[i]+= 1_dbl; //net_wb[flattenb(i, n_neurons, layer_id, n_net_w)];
+
+        fmt::print("\n{}\n", retval[i]);
+        fmt::print("Here1");
+
+        std::cout << std::endl;
         retval[i] = activation(retval[i]);
+        fmt::print("Here2");
+        std::cout << std::endl;
     }
+    fmt::print("Here3");
+
     return retval;
 }
 } // namespace detail
diff --git a/test/model_ffnn.cpp b/test/model_ffnn.cpp
index b09c6a113..5ac59511e 100644
--- a/test/model_ffnn.cpp
+++ b/test/model_ffnn.cpp
@@ -20,8 +20,9 @@ using namespace heyoka;
 
 TEST_CASE("impl")
 {
+    auto linear = [](expression ret) -> expression { return ret; };
     auto [x] = make_vars("x");
-    auto my_net = model::ffnn_impl({x}, 2, {2, 2}, {heyoka::tanh, heyoka::tanh, [](auto ret) { return ret; }},
+    auto my_net = model::ffnn_impl({x}, 2, {2, 2}, {heyoka::tanh, heyoka::tanh, heyoka::tanh},
                                    {1_dbl, 2_dbl, 3_dbl, 4_dbl, 5_dbl, 6_dbl, 7_dbl, 8_dbl, 9_dbl, 0_dbl, 1_dbl, 2_dbl,
                                     3_dbl, 4_dbl, 5_dbl, 6_dbl});
 }
\ No newline at end of file

From 7af74f3c7592c836b363cc4b61771a25806dcd64 Mon Sep 17 00:00:00 2001
From: Dario Izzo <darioizzo@gmail.com>
Date: Tue, 31 Oct 2023 23:04:13 +0100
Subject: [PATCH 07/14] bug fix and simplification

---
 include/heyoka/model/ffnn.hpp | 14 ++++----
 src/model/ffnn.cpp            | 68 ++++++++---------------------------
 test/model_ffnn.cpp           |  2 +-
 3 files changed, 23 insertions(+), 61 deletions(-)

diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
index e1c0b456a..fde4e126b 100644
--- a/include/heyoka/model/ffnn.hpp
+++ b/include/heyoka/model/ffnn.hpp
@@ -25,14 +25,16 @@ namespace model
 // This c++ function returns the symbolic expressions of the `n_out` output neurons in a feed forward neural network,
 // as a function of the `n_in` input expressions.
 //
-// The expression will contain the weights and biases of the neural network flattened into `pars` with the following conventions:
+// The expression will contain the weights and biases of the neural network flattened into `pars` with the following
+// conventions:
 //
-// from the left to right layer of parameters: [flattened weights1, flattened weights2,  ... , biases1, bises2, ...]
+// from the left to right layer of parameters: [W01, W12,W23, ..., B1,B2,B3,....] where the weight matrices Wij are
+// to be considered as flattened (row first) and so are the bias vectors.
 //
-HEYOKA_DLL_PUBLIC std::vector<expression>
-ffnn_impl(const std::vector<expression> &, std::uint32_t, const std::vector<std::uint32_t> &,
-          const std::vector<std::function<expression(const expression &)>> &,
-          const std::vector<expression> &);
+HEYOKA_DLL_PUBLIC std::vector<expression> ffnn_impl(const std::vector<expression> &, std::uint32_t,
+                                                    const std::vector<std::uint32_t> &,
+                                                    const std::vector<std::function<expression(const expression &)>> &,
+                                                    const std::vector<expression> &);
 } // namespace model
 
 HEYOKA_END_NAMESPACE
diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
index a0e2c6420..2b6a9570b 100644
--- a/src/model/ffnn.cpp
+++ b/src/model/ffnn.cpp
@@ -24,46 +24,11 @@ namespace model
 {
 namespace detail
 {
-// From the i-th neuron of the layer_id and the incoming j-th connection, returns the corresponding weight position in
-// the flattened structure. NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
-std::vector<expression>::size_type flattenw(std::uint32_t i, std::uint32_t j,
-                                            const std::vector<std::uint32_t> &n_neurons, std::uint32_t layer_id)
-{
-    assert(layer_id > 0);
-    // The weight for the jth-neuron of the ith layer will be placed after all previous layers.
-    // We start counting how many in flattened.
-    std::uint32_t counter = 0;
-    for (std::uint32_t k = 1; k < layer_id; ++k) {
-        counter += n_neurons[k] * n_neurons[k - 1];
-    }
-    // We then add the weights used for the previous neurons on the same layer.
-    counter += i * n_neurons[layer_id - 1];
-    // And return the index corresponding to the j-th weight.
-    return counter + j;
-}
-
-// From the i-th neuron of the layer_id, returns the corresponding bias position in the
-// flattened structure.
-// NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
-std::vector<expression>::size_type flattenb(std::uint32_t i, const std::vector<std::uint32_t> &n_neurons,
-                                            // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
-                                            std::uint32_t layer_id, std::uint32_t n_net_w)
-{
-    assert(layer_id > 0);
-    // The weight for the jth-neuron of the ith layer will be placed after all previous layers.
-    // We start counting how many in flattened.
-    std::uint32_t counter = 0;
-    for (std::uint32_t k = 1; k < layer_id; ++k) {
-        counter += n_neurons[k];
-    }
-    // And return the index corresponding to the i-th bias.
-    return counter + i + n_net_w;
-}
-
 std::vector<expression> compute_layer(std::uint32_t layer_id, const std::vector<expression> &inputs,
                                       const std::vector<std::uint32_t> &n_neurons,
                                       const std::function<expression(const expression &)> &activation,
-                                      const std::vector<expression> &net_wb, std::uint32_t n_net_w)
+                                      const std::vector<expression> &net_wb, std::uint32_t n_net_w,
+                                      std::uint32_t &wcounter, std::uint32_t &bcounter)
 {
     assert(layer_id > 0);
     auto n_neurons_prev_layer = boost::numeric_cast<std::uint32_t>(inputs.size());
@@ -75,25 +40,20 @@ std::vector<expression> compute_layer(std::uint32_t layer_id, const std::vector<
 
     for (std::uint32_t i = 0u; i < n_neurons_curr_layer; ++i) {
         for (std::uint32_t j = 0u; j < n_neurons_prev_layer; ++j) {
-            fmt::print("layer, i, j, idx: {}, {}, {}\n", layer_id, i, j, flattenw(i, j, n_neurons, layer_id));
+            fmt::print("layer, i, j, idx: {}, {}, {}, {}\n", layer_id, i, j, wcounter);
             std::cout << std::endl;
-            retval[i] += 1_dbl;//net_wb[flattenw(i, j, n_neurons, layer_id)] * inputs[j];
+            // Add the weight and update the weight counter
+            retval[i] += net_wb[wcounter] * inputs[j];
+            ++wcounter;
         }
-        fmt::print("idxb {}\n", flattenb(i, n_neurons, layer_id, n_net_w));
-        std::cout << std::endl;
-
-        retval[i]+= 1_dbl; //net_wb[flattenb(i, n_neurons, layer_id, n_net_w)];
-
-        fmt::print("\n{}\n", retval[i]);
-        fmt::print("Here1");
-
+        fmt::print("idxb {}\n", bcounter + n_net_w);
         std::cout << std::endl;
+        // Add the bias and update the counter
+        retval[i] += net_wb[bcounter + n_net_w];
+        ++bcounter;
+        // Activation function
         retval[i] = activation(retval[i]);
-        fmt::print("Here2");
-        std::cout << std::endl;
     }
-    fmt::print("Here3");
-
     return retval;
 }
 } // namespace detail
@@ -156,10 +116,10 @@ std::vector<expression> ffnn_impl(
     // Now we build the expressions recursively going from layer to layer (L = f(Wx+b)))
 
     std::vector<expression> retval = in;
+    std::uint32_t wcounter = 0;
+    std::uint32_t bcounter = 0;
     for (std::uint32_t i = 1u; i < n_layers; ++i) {
-        fmt::print("{},{}", i, n_neurons[i]);
-        std::cout << std::endl;
-        retval = detail::compute_layer(i, retval, n_neurons, activations[i], net_wb, n_net_w);
+        retval = detail::compute_layer(i, retval, n_neurons, activations[i - 1], net_wb, n_net_w, wcounter, bcounter);
     }
     return retval;
 }
diff --git a/test/model_ffnn.cpp b/test/model_ffnn.cpp
index 5ac59511e..98e658290 100644
--- a/test/model_ffnn.cpp
+++ b/test/model_ffnn.cpp
@@ -22,7 +22,7 @@ TEST_CASE("impl")
 {
     auto linear = [](expression ret) -> expression { return ret; };
     auto [x] = make_vars("x");
-    auto my_net = model::ffnn_impl({x}, 2, {2, 2}, {heyoka::tanh, heyoka::tanh, heyoka::tanh},
+    auto my_net = model::ffnn_impl({x}, 2, {2, 2}, {heyoka::tanh, heyoka::tanh, linear},
                                    {1_dbl, 2_dbl, 3_dbl, 4_dbl, 5_dbl, 6_dbl, 7_dbl, 8_dbl, 9_dbl, 0_dbl, 1_dbl, 2_dbl,
                                     3_dbl, 4_dbl, 5_dbl, 6_dbl});
 }
\ No newline at end of file

From 253450fdae74e981f9772dfc5b7458a5f52dd23d Mon Sep 17 00:00:00 2001
From: Dario Izzo <darioizzo@gmail.com>
Date: Wed, 1 Nov 2023 07:50:00 +0100
Subject: [PATCH 08/14] [SKIP CI] WIP

---
 src/model/ffnn.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
index 2b6a9570b..aa5dcca31 100644
--- a/src/model/ffnn.cpp
+++ b/src/model/ffnn.cpp
@@ -113,7 +113,7 @@ std::vector<expression> ffnn_impl(
             n_net_wb, net_wb.size()));
     }
 
-    // Now we build the expressions recursively going from layer to layer (L = f(Wx+b)))
+    // Now we build the expressions recursively transvering from layer to layer (L = f(Wx+b)))
 
     std::vector<expression> retval = in;
     std::uint32_t wcounter = 0;

From c3ecf67a3b6ac9684382d985d6489a6afb2ff689 Mon Sep 17 00:00:00 2001
From: Dario Izzo <dario.izzo@gmail.com>
Date: Wed, 1 Nov 2023 17:49:00 +0100
Subject: [PATCH 09/14] [CI SKIP] WIP

---
 include/heyoka/kw.hpp         |  7 +++
 include/heyoka/model/ffnn.hpp | 93 ++++++++++++++++++++++++++++++++++-
 src/model/ffnn.cpp            | 59 +++++++++-------------
 test/model_ffnn.cpp           | 72 +++++++++++++++++++++++++--
 4 files changed, 188 insertions(+), 43 deletions(-)

diff --git a/include/heyoka/kw.hpp b/include/heyoka/kw.hpp
index c28481287..6cbadce24 100644
--- a/include/heyoka/kw.hpp
+++ b/include/heyoka/kw.hpp
@@ -30,6 +30,13 @@ IGOR_MAKE_NAMED_ARGUMENT(parallel_mode);
 IGOR_MAKE_NAMED_ARGUMENT(prec);
 IGOR_MAKE_NAMED_ARGUMENT(mu);
 
+// kwargs for the ffnn
+IGOR_MAKE_NAMED_ARGUMENT(inputs);
+IGOR_MAKE_NAMED_ARGUMENT(nn_hidden);
+IGOR_MAKE_NAMED_ARGUMENT(n_out);
+IGOR_MAKE_NAMED_ARGUMENT(activations);
+IGOR_MAKE_NAMED_ARGUMENT(nn_wb);
+
 } // namespace kw
 
 HEYOKA_END_NAMESPACE
diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
index fde4e126b..0acd1ddf3 100644
--- a/include/heyoka/model/ffnn.hpp
+++ b/include/heyoka/model/ffnn.hpp
@@ -15,6 +15,7 @@
 #include <vector>
 
 #include <heyoka/config.hpp>
+#include <heyoka/detail/igor.hpp>
 #include <heyoka/detail/visibility.hpp>
 #include <heyoka/expression.hpp>
 
@@ -22,6 +23,88 @@ HEYOKA_BEGIN_NAMESPACE
 
 namespace model
 {
+namespace detail
+{
+template <typename... KwArgs>
+auto ffnn_common_opts(KwArgs &&...kw_args)
+{
+    igor::parser p{kw_args...};
+
+    static_assert(!p.has_unnamed_arguments(), "This function accepts only named arguments");
+
+    // Network inputs. Mandatory
+    auto inputs = [&p]() {
+        if constexpr (p.has(kw::inputs)) {
+            return std::vector<expression>{p(kw::inputs)};
+        } else {
+            static_assert(::heyoka::detail::always_false_v<KwArgs...>,
+                          "The 'inputs' keyword argument is necessary but it was not provided");
+        }
+    }();
+
+    // Number of hidden neurons per hidden layer. Mandatory
+    auto nn_hidden = [&p]() {
+        if constexpr (p.has(kw::nn_hidden)) {
+            return std::vector<std::uint32_t>{p(kw::nn_hidden)};
+        } else {
+            static_assert(::heyoka::detail::always_false_v<KwArgs...>,
+                          "The 'nn_hidden' keyword argument is necessary but it was not provided");
+        }
+    }();
+
+    // Number of network outputs. Mandatory
+    auto n_out = [&p]() {
+        if constexpr (p.has(kw::n_out)) {
+            return std::uint32_t{p(kw::n_out)};
+        } else {
+            static_assert(::heyoka::detail::always_false_v<KwArgs...>,
+                          "The 'n_out' keyword argument is necessary but it was not provided");
+        }
+    }();
+
+    // Network activation functions. Mandatory
+    auto activations = [&p]() {
+        if constexpr (p.has(kw::activations)) {
+            return std::vector<std::function<expression(const expression &)>>{p(kw::activations)};
+        } else {
+            static_assert(::heyoka::detail::always_false_v<KwArgs...>,
+                          "The 'activations' keyword argument is necessary but it was not provided");
+        }
+    }();
+
+    // Network weights and biases. Defaults to heyoka parameters.
+    auto nn_wb = [&p, &nn_hidden, &inputs, n_out]() {
+        if constexpr (p.has(kw::nn_wb)) {
+            return std::vector<expression> {p(kw::nn_wb)};
+        } else {
+            // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons)
+            auto n_hidden_layers = boost::numeric_cast<std::uint32_t>(nn_hidden.size());
+            // Number of neuronal layers (counting input and output)
+            auto n_layers = n_hidden_layers + 2;
+            // Number of inputs
+            auto n_in = boost::numeric_cast<std::uint32_t>(inputs.size());
+            // Number of neurons per neuronal layer
+            std::vector<std::uint32_t> n_neurons = nn_hidden;
+            n_neurons.insert(n_neurons.begin(), n_in);
+            n_neurons.insert(n_neurons.end(), n_out);
+            // Number of network parameters (wb: weights and biases, w: only weights)
+            std::uint32_t n_wb = 0u;
+            for (std::uint32_t i = 1u; i < n_layers; ++i) {
+                n_wb += n_neurons[i - 1] * n_neurons[i];
+                n_wb += n_neurons[i];
+            }
+            std::vector<expression> retval(n_wb);
+            for (decltype(retval.size()) i = 0; i < retval.size(); ++i) {
+                retval[i] = heyoka::par[i];
+            }
+            return retval;
+        }
+    }();
+
+    return std::tuple{std::move(inputs), std::move(nn_hidden), std::move(n_out), std::move(activations),
+                      std::move(nn_wb)};
+}
+
 // This c++ function returns the symbolic expressions of the `n_out` output neurons in a feed forward neural network,
 // as a function of the `n_in` input expressions.
 //
@@ -31,10 +114,16 @@ namespace model
 // from the left to right layer of parameters: [W01, W12,W23, ..., B1,B2,B3,....] where the weight matrices Wij are
 // to be considered as flattened (row first) and so are the bias vectors.
 //
-HEYOKA_DLL_PUBLIC std::vector<expression> ffnn_impl(const std::vector<expression> &, std::uint32_t,
-                                                    const std::vector<std::uint32_t> &,
+HEYOKA_DLL_PUBLIC std::vector<expression> ffnn_impl(const std::vector<expression> &, const std::vector<std::uint32_t> &,
+                                                    std::uint32_t,
                                                     const std::vector<std::function<expression(const expression &)>> &,
                                                     const std::vector<expression> &);
+} // namespace detail
+
+inline constexpr auto ffnn = [](const auto &...kw_args) -> std::vector<expression> {
+    return std::apply(detail::ffnn_impl, detail::ffnn_common_opts(kw_args...));
+};
+
 } // namespace model
 
 HEYOKA_END_NAMESPACE
diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
index aa5dcca31..ad1bb674e 100644
--- a/src/model/ffnn.cpp
+++ b/src/model/ffnn.cpp
@@ -20,14 +20,12 @@
 
 HEYOKA_BEGIN_NAMESPACE
 
-namespace model
-{
-namespace detail
+namespace model::detail
 {
 std::vector<expression> compute_layer(std::uint32_t layer_id, const std::vector<expression> &inputs,
                                       const std::vector<std::uint32_t> &n_neurons,
                                       const std::function<expression(const expression &)> &activation,
-                                      const std::vector<expression> &net_wb, std::uint32_t n_net_w,
+                                      const std::vector<expression> &nn_wb, std::uint32_t n_net_w,
                                       std::uint32_t &wcounter, std::uint32_t &bcounter)
 {
     assert(layer_id > 0);
@@ -35,66 +33,54 @@ std::vector<expression> compute_layer(std::uint32_t layer_id, const std::vector<
     auto n_neurons_curr_layer = n_neurons[layer_id];
 
     std::vector<expression> retval(n_neurons_curr_layer, 0_dbl);
-    fmt::print("net_wb: {}\n", net_wb.size());
-    std::cout << std::endl;
-
     for (std::uint32_t i = 0u; i < n_neurons_curr_layer; ++i) {
         for (std::uint32_t j = 0u; j < n_neurons_prev_layer; ++j) {
-            fmt::print("layer, i, j, idx: {}, {}, {}, {}\n", layer_id, i, j, wcounter);
-            std::cout << std::endl;
+
             // Add the weight and update the weight counter
-            retval[i] += net_wb[wcounter] * inputs[j];
+            retval[i] += nn_wb[wcounter] * inputs[j];
             ++wcounter;
         }
-        fmt::print("idxb {}\n", bcounter + n_net_w);
-        std::cout << std::endl;
+
         // Add the bias and update the counter
-        retval[i] += net_wb[bcounter + n_net_w];
+        retval[i] += nn_wb[bcounter + n_net_w];
         ++bcounter;
         // Activation function
         retval[i] = activation(retval[i]);
     }
     return retval;
 }
-} // namespace detail
 
-std::vector<expression> ffnn_impl(
-    // NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
-    const std::vector<expression> &in, std::uint32_t n_out,
-    const std::vector<std::uint32_t> &n_neurons_per_hidden_layer,
-    const std::vector<std::function<expression(const expression &)>> &activations,
-    const std::vector<expression> &net_wb)
+std::vector<expression> ffnn_impl(const std::vector<expression> &in, const std::vector<std::uint32_t> &nn_hidden,
+                                  std::uint32_t n_out,
+                                  const std::vector<std::function<expression(const expression &)>> &activations,
+                                  const std::vector<expression> &nn_wb)
 {
     // Sanity checks
-    if (n_neurons_per_hidden_layer.size() + 1 != activations.size()) {
+    if (nn_hidden.size() + 1 != activations.size()) {
         throw std::invalid_argument(fmt::format(
-            "The number of hidden layers, as detected from the inputs, was {}, while"
+            "The number of hidden layers, as detected from the inputs, was {}, while "
             "the number of activation function supplied was {}. A FFNN needs exactly one more activation function "
             "than the number of hidden layers.",
-            n_neurons_per_hidden_layer.size(), activations.size()));
+            nn_hidden.size(), activations.size()));
     }
     if (in.empty()) {
-        throw std::invalid_argument("The inputs provided to the ffnn seem to be an empty vector.");
+        throw std::invalid_argument("The inputs provided to the FFNN is an empty vector.");
     }
     if (n_out == 0) {
         throw std::invalid_argument("The number of network outputs cannot be zero.");
     }
-    if (!std::all_of(n_neurons_per_hidden_layer.begin(), n_neurons_per_hidden_layer.end(),
-                     [](std::uint32_t item) { return item > 0; })) {
+    if (!std::all_of(nn_hidden.begin(), nn_hidden.end(), [](std::uint32_t item) { return item > 0; })) {
         throw std::invalid_argument("The number of neurons for each hidden layer must be greater than zero!");
     }
-    if (n_neurons_per_hidden_layer.empty()) { // TODO(darioizzo): maybe this is actually a wanted corner case, remove?
-        throw std::invalid_argument("The number of hidden layers cannot be zero.");
-    }
 
     // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons)
-    auto n_hidden_layers = boost::numeric_cast<std::uint32_t>(n_neurons_per_hidden_layer.size());
+    auto n_hidden_layers = boost::numeric_cast<std::uint32_t>(nn_hidden.size());
     // Number of neuronal layers (counting input and output)
     auto n_layers = n_hidden_layers + 2;
     // Number of inputs
     auto n_in = boost::numeric_cast<std::uint32_t>(in.size());
     // Number of neurons per neuronal layer
-    std::vector<std::uint32_t> n_neurons = n_neurons_per_hidden_layer;
+    std::vector<std::uint32_t> n_neurons = nn_hidden;
     n_neurons.insert(n_neurons.begin(), n_in);
     n_neurons.insert(n_neurons.end(), n_out);
     // Number of network parameters (wb: weights and biases, w: only weights)
@@ -106,20 +92,19 @@ std::vector<expression> ffnn_impl(
         n_net_wb += n_neurons[i];
     }
     // Sanity check
-    if (net_wb.size() != n_net_wb) {
+    if (nn_wb.size() != n_net_wb) {
         throw std::invalid_argument(fmt::format(
-            "The number of network parameters, detected from its structure to be {}, does not match the size of"
-            "the corresponding expressions {} ",
-            n_net_wb, net_wb.size()));
+            "The number of network parameters, detected from its structure to be {}, does not match the size of "
+            "the corresponding expressions: {}.",
+            n_net_wb, nn_wb.size()));
     }
 
     // Now we build the expressions recursively transvering from layer to layer (L = f(Wx+b)))
-
     std::vector<expression> retval = in;
     std::uint32_t wcounter = 0;
     std::uint32_t bcounter = 0;
     for (std::uint32_t i = 1u; i < n_layers; ++i) {
-        retval = detail::compute_layer(i, retval, n_neurons, activations[i - 1], net_wb, n_net_w, wcounter, bcounter);
+        retval = detail::compute_layer(i, retval, n_neurons, activations[i - 1], nn_wb, n_net_w, wcounter, bcounter);
     }
     return retval;
 }
diff --git a/test/model_ffnn.cpp b/test/model_ffnn.cpp
index 98e658290..8dc2172e6 100644
--- a/test/model_ffnn.cpp
+++ b/test/model_ffnn.cpp
@@ -15,14 +15,78 @@
 #include <heyoka/model/ffnn.hpp>
 
 #include "catch.hpp"
+#include "heyoka/kw.hpp"
 
 using namespace heyoka;
 
 TEST_CASE("impl")
 {
+    // A linear layer, just because
     auto linear = [](expression ret) -> expression { return ret; };
-    auto [x] = make_vars("x");
-    auto my_net = model::ffnn_impl({x}, 2, {2, 2}, {heyoka::tanh, heyoka::tanh, linear},
-                                   {1_dbl, 2_dbl, 3_dbl, 4_dbl, 5_dbl, 6_dbl, 7_dbl, 8_dbl, 9_dbl, 0_dbl, 1_dbl, 2_dbl,
-                                    3_dbl, 4_dbl, 5_dbl, 6_dbl});
+    // We also define a few symbols
+    auto [x, y, z] = make_vars("x", "y", "z");
+
+    // First, we test malformed cases and their throws.
+    // 1 - number of activations function is wrong
+    REQUIRE_THROWS_AS(model::detail::ffnn_impl({x}, {1}, 2, {heyoka::tanh, heyoka::tanh, linear},
+                                               {1_dbl, 2_dbl, 3_dbl, 4_dbl, 5_dbl, 6_dbl}),
+                      std::invalid_argument);
+    // 2 - number of inputs is zero
+    REQUIRE_THROWS_AS(
+        model::detail::ffnn_impl({}, {1}, 2, {heyoka::tanh, heyoka::tanh}, {1_dbl, 2_dbl, 3_dbl, 4_dbl, 5_dbl, 6_dbl}),
+        std::invalid_argument);
+    // 3 - number of outputs is zero
+    REQUIRE_THROWS_AS(model::detail::ffnn_impl({x}, {1}, 0, {heyoka::tanh, heyoka::tanh}, {1_dbl, 2_dbl, 3_dbl, 4_dbl}),
+                      std::invalid_argument);
+    // 4 - One of the hidden layers has zero neurons
+    REQUIRE_THROWS_AS(model::detail::ffnn_impl({x}, {1, 0}, 2, {heyoka::tanh, heyoka::tanh, linear}, {1_dbl, 2_dbl}),
+                      std::invalid_argument);
+    // 5 - Wrong number of weights/biases
+    REQUIRE_THROWS_AS(
+        model::detail::ffnn_impl({x}, {1}, 1, {heyoka::tanh, heyoka::tanh}, {1_dbl, 2_dbl, 3_dbl, 5_dbl, 6_dbl}),
+        std::invalid_argument);
+
+    // We now check some hand coded networks
+    {
+        auto my_net = model::detail::ffnn_impl({x}, {}, 1, {linear}, {1_dbl, 2_dbl});
+        REQUIRE(my_net[0] == expression(2_dbl + x));
+    }
+    {
+        auto my_net = model::detail::ffnn_impl({x}, {}, 1, {heyoka::tanh}, {1_dbl, 2_dbl});
+        REQUIRE(my_net[0] == expression(heyoka::tanh(2_dbl + x)));
+    }
+    {
+        auto my_net = model::detail::ffnn_impl({x}, {1}, 1, {heyoka::tanh, linear}, {1_dbl, 2_dbl, 3_dbl, 4_dbl});
+        REQUIRE(my_net[0] == expression(4_dbl + (2_dbl * heyoka::tanh(3_dbl + x))));
+    }
+    {
+        auto my_net = model::detail::ffnn_impl({x}, {1}, 1, {heyoka::tanh, heyoka::sin}, {1_dbl, 2_dbl, 3_dbl, 4_dbl});
+        REQUIRE(my_net[0] == expression(heyoka::sin(4_dbl + (2_dbl * heyoka::tanh(3_dbl + x)))));
+    }
+    {
+        auto my_net = model::detail::ffnn_impl({x, y}, {2}, 1, {heyoka::sin, heyoka::cos},
+                                               {1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl});
+        REQUIRE(my_net[0] == expression(heyoka::cos(1_dbl + 2_dbl * heyoka::sin(1_dbl + x + y))));
+    }
+}
+
+TEST_CASE("igor_iface")
+{
+    auto [x, y, z] = make_vars("x", "y", "z");
+    {
+        auto igor_v = model::ffnn(kw::inputs = {x, y}, kw::nn_hidden = std::vector<std::uint32_t>{2u}, kw::n_out = 1u, 
+                 kw::activations = std::vector<std::function<expression(const expression &)>>{heyoka::sin, heyoka::cos},
+                 kw::nn_wb = {1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl});
+        auto vanilla_v = model::detail::ffnn_impl({x, y}, {2u}, 1u, {heyoka::sin, heyoka::cos},
+                                                  {1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl});
+        REQUIRE(igor_v == vanilla_v);
+    }
+    // We test the expected setting for the default weights+biases expressions to par[i].
+    {
+        auto igor_v = model::ffnn(kw::inputs = {x, y}, kw::nn_hidden = std::vector<std::uint32_t>{2u}, kw::n_out = 1u, 
+                 kw::activations = std::vector<std::function<expression(const expression &)>>{heyoka::sin, heyoka::cos});
+        auto vanilla_v = model::detail::ffnn_impl({x, y}, {2u}, 1u, {heyoka::sin, heyoka::cos},
+                                                  {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8]});
+        REQUIRE(igor_v == vanilla_v);
+    }
 }
\ No newline at end of file

From a4a5ae9f98d7ecb746e4a83afa2de02452158b12 Mon Sep 17 00:00:00 2001
From: Dario Izzo <darioizzo@gmail.com>
Date: Wed, 1 Nov 2023 18:26:57 +0100
Subject: [PATCH 10/14] [skip ci] - trying to make sense of igor

---
 include/heyoka/model/ffnn.hpp |  7 ++++---
 test/model_ffnn.cpp           | 17 ++++++++++-------
 2 files changed, 14 insertions(+), 10 deletions(-)

diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
index 0acd1ddf3..b9831a156 100644
--- a/include/heyoka/model/ffnn.hpp
+++ b/include/heyoka/model/ffnn.hpp
@@ -74,8 +74,9 @@ auto ffnn_common_opts(KwArgs &&...kw_args)
 
     // Network weights and biases. Defaults to heyoka parameters.
     auto nn_wb = [&p, &nn_hidden, &inputs, n_out]() {
+        std::vector<expression> retval{};
         if constexpr (p.has(kw::nn_wb)) {
-            return std::vector<expression> {p(kw::nn_wb)};
+            retval = p(kw::nn_wb);
         } else {
             // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons)
             auto n_hidden_layers = boost::numeric_cast<std::uint32_t>(nn_hidden.size());
@@ -93,12 +94,12 @@ auto ffnn_common_opts(KwArgs &&...kw_args)
                 n_wb += n_neurons[i - 1] * n_neurons[i];
                 n_wb += n_neurons[i];
             }
-            std::vector<expression> retval(n_wb);
+            retval.resize(n_wb);
             for (decltype(retval.size()) i = 0; i < retval.size(); ++i) {
                 retval[i] = heyoka::par[i];
             }
-            return retval;
         }
+        return retval;
     }();
 
     return std::tuple{std::move(inputs), std::move(nn_hidden), std::move(n_out), std::move(activations),
diff --git a/test/model_ffnn.cpp b/test/model_ffnn.cpp
index 8dc2172e6..9cd2fe141 100644
--- a/test/model_ffnn.cpp
+++ b/test/model_ffnn.cpp
@@ -74,19 +74,22 @@ TEST_CASE("igor_iface")
 {
     auto [x, y, z] = make_vars("x", "y", "z");
     {
-        auto igor_v = model::ffnn(kw::inputs = {x, y}, kw::nn_hidden = std::vector<std::uint32_t>{2u}, kw::n_out = 1u, 
-                 kw::activations = std::vector<std::function<expression(const expression &)>>{heyoka::sin, heyoka::cos},
-                 kw::nn_wb = {1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl});
+        auto igor_v = model::ffnn(
+            kw::inputs = {x, y}, kw::nn_hidden = std::vector<std::uint32_t>{2u}, kw::n_out = 1u,
+            kw::activations = std::vector<std::function<expression(const expression &)>>{heyoka::sin, heyoka::cos},
+            kw::nn_wb = {1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl});
         auto vanilla_v = model::detail::ffnn_impl({x, y}, {2u}, 1u, {heyoka::sin, heyoka::cos},
                                                   {1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl});
         REQUIRE(igor_v == vanilla_v);
     }
     // We test the expected setting for the default weights+biases expressions to par[i].
     {
-        auto igor_v = model::ffnn(kw::inputs = {x, y}, kw::nn_hidden = std::vector<std::uint32_t>{2u}, kw::n_out = 1u, 
-                 kw::activations = std::vector<std::function<expression(const expression &)>>{heyoka::sin, heyoka::cos});
-        auto vanilla_v = model::detail::ffnn_impl({x, y}, {2u}, 1u, {heyoka::sin, heyoka::cos},
-                                                  {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8]});
+        auto igor_v = model::ffnn(
+            kw::inputs = {x, y}, kw::nn_hidden = std::vector<std::uint32_t>{2u}, kw::n_out = 1u,
+            kw::activations = std::vector<std::function<expression(const expression &)>>{heyoka::sin, heyoka::cos});
+        auto vanilla_v
+            = model::detail::ffnn_impl({x, y}, {2u}, 1u, {heyoka::sin, heyoka::cos},
+                                       {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8]});
         REQUIRE(igor_v == vanilla_v);
     }
 }
\ No newline at end of file

From 3b4e8c3f34df6a83795f00611167b0e55ddab762 Mon Sep 17 00:00:00 2001
From: Dario Izzo <dario.izzo@gmail.com>
Date: Wed, 1 Nov 2023 19:10:06 +0100
Subject: [PATCH 11/14] final changes to make igor happy

---
 include/heyoka/model/ffnn.hpp | 24 ++++++++++++------------
 src/model/ffnn.cpp            |  2 +-
 2 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
index b9831a156..a358e5f4a 100644
--- a/include/heyoka/model/ffnn.hpp
+++ b/include/heyoka/model/ffnn.hpp
@@ -74,20 +74,20 @@ auto ffnn_common_opts(KwArgs &&...kw_args)
 
     // Network weights and biases. Defaults to heyoka parameters.
     auto nn_wb = [&p, &nn_hidden, &inputs, n_out]() {
+        // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons)
+        auto n_hidden_layers = boost::numeric_cast<std::uint32_t>(nn_hidden.size());
+        // Number of neuronal layers (counting input and output)
+        auto n_layers = n_hidden_layers + 2;
+        // Number of inputs
+        auto n_in = boost::numeric_cast<std::uint32_t>(inputs.size());
+        // Number of neurons per neuronal layer
+        std::vector<std::uint32_t> n_neurons = nn_hidden;
+        n_neurons.insert(n_neurons.begin(), n_in);
+        n_neurons.insert(n_neurons.end(), n_out);
         std::vector<expression> retval{};
         if constexpr (p.has(kw::nn_wb)) {
             retval = p(kw::nn_wb);
         } else {
-            // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons)
-            auto n_hidden_layers = boost::numeric_cast<std::uint32_t>(nn_hidden.size());
-            // Number of neuronal layers (counting input and output)
-            auto n_layers = n_hidden_layers + 2;
-            // Number of inputs
-            auto n_in = boost::numeric_cast<std::uint32_t>(inputs.size());
-            // Number of neurons per neuronal layer
-            std::vector<std::uint32_t> n_neurons = nn_hidden;
-            n_neurons.insert(n_neurons.begin(), n_in);
-            n_neurons.insert(n_neurons.end(), n_out);
             // Number of network parameters (wb: weights and biases, w: only weights)
             std::uint32_t n_wb = 0u;
             for (std::uint32_t i = 1u; i < n_layers; ++i) {
@@ -95,8 +95,8 @@ auto ffnn_common_opts(KwArgs &&...kw_args)
                 n_wb += n_neurons[i];
             }
             retval.resize(n_wb);
-            for (decltype(retval.size()) i = 0; i < retval.size(); ++i) {
-                retval[i] = heyoka::par[i];
+            for (decltype(retval.size()) i = 0u; i < retval.size(); ++i) {
+                retval[i] = heyoka::par[boost::numeric_cast<std::uint32_t>(i)];
             }
         }
         return retval;
diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
index ad1bb674e..1a5120765 100644
--- a/src/model/ffnn.cpp
+++ b/src/model/ffnn.cpp
@@ -108,5 +108,5 @@ std::vector<expression> ffnn_impl(const std::vector<expression> &in, const std::
     }
     return retval;
 }
-} // namespace model
+} // namespace model::detail
 HEYOKA_END_NAMESPACE
\ No newline at end of file

From 574059c268ba2c261a67c13f5eec3ae35f28d4f4 Mon Sep 17 00:00:00 2001
From: Francesco Biscani <bluescarni@gmail.com>
Date: Thu, 2 Nov 2023 13:15:05 +0100
Subject: [PATCH 12/14] A first round of cleanups, missing headers, etc.

---
 include/heyoka/model/ffnn.hpp | 17 +++++++++++------
 include/heyoka/models.hpp     |  1 +
 src/model/ffnn.cpp            | 22 ++++++++++++++--------
 3 files changed, 26 insertions(+), 14 deletions(-)

diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
index a358e5f4a..62cdedbcc 100644
--- a/include/heyoka/model/ffnn.hpp
+++ b/include/heyoka/model/ffnn.hpp
@@ -9,13 +9,17 @@
 #ifndef HEYOKA_MODEL_FFNN_HPP
 #define HEYOKA_MODEL_FFNN_HPP
 
-#include <atomic>
+#include <cstdint>
+#include <functional>
 #include <tuple>
 #include <utility>
 #include <vector>
 
+#include <boost/numeric/conversion/cast.hpp>
+
 #include <heyoka/config.hpp>
 #include <heyoka/detail/igor.hpp>
+#include <heyoka/detail/type_traits.hpp>
 #include <heyoka/detail/visibility.hpp>
 #include <heyoka/expression.hpp>
 
@@ -25,8 +29,9 @@ namespace model
 {
 namespace detail
 {
+
 template <typename... KwArgs>
-auto ffnn_common_opts(KwArgs &&...kw_args)
+auto ffnn_common_opts(const KwArgs &...kw_args)
 {
     igor::parser p{kw_args...};
 
@@ -37,7 +42,7 @@ auto ffnn_common_opts(KwArgs &&...kw_args)
         if constexpr (p.has(kw::inputs)) {
             return std::vector<expression>{p(kw::inputs)};
         } else {
-            static_assert(::heyoka::detail::always_false_v<KwArgs...>,
+            static_assert(heyoka::detail::always_false_v<KwArgs...>,
                           "The 'inputs' keyword argument is necessary but it was not provided");
         }
     }();
@@ -47,7 +52,7 @@ auto ffnn_common_opts(KwArgs &&...kw_args)
         if constexpr (p.has(kw::nn_hidden)) {
             return std::vector<std::uint32_t>{p(kw::nn_hidden)};
         } else {
-            static_assert(::heyoka::detail::always_false_v<KwArgs...>,
+            static_assert(heyoka::detail::always_false_v<KwArgs...>,
                           "The 'nn_hidden' keyword argument is necessary but it was not provided");
         }
     }();
@@ -57,7 +62,7 @@ auto ffnn_common_opts(KwArgs &&...kw_args)
         if constexpr (p.has(kw::n_out)) {
             return std::uint32_t{p(kw::n_out)};
         } else {
-            static_assert(::heyoka::detail::always_false_v<KwArgs...>,
+            static_assert(heyoka::detail::always_false_v<KwArgs...>,
                           "The 'n_out' keyword argument is necessary but it was not provided");
         }
     }();
@@ -67,7 +72,7 @@ auto ffnn_common_opts(KwArgs &&...kw_args)
         if constexpr (p.has(kw::activations)) {
             return std::vector<std::function<expression(const expression &)>>{p(kw::activations)};
         } else {
-            static_assert(::heyoka::detail::always_false_v<KwArgs...>,
+            static_assert(heyoka::detail::always_false_v<KwArgs...>,
                           "The 'activations' keyword argument is necessary but it was not provided");
         }
     }();
diff --git a/include/heyoka/models.hpp b/include/heyoka/models.hpp
index 12b11a67f..545034687 100644
--- a/include/heyoka/models.hpp
+++ b/include/heyoka/models.hpp
@@ -10,6 +10,7 @@
 #define HEYOKA_MODELS_HPP
 
 #include <heyoka/model/cr3bp.hpp>
+#include <heyoka/model/ffnn.hpp>
 #include <heyoka/model/fixed_centres.hpp>
 #include <heyoka/model/mascon.hpp>
 #include <heyoka/model/nbody.hpp>
diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
index 1a5120765..93e530645 100644
--- a/src/model/ffnn.cpp
+++ b/src/model/ffnn.cpp
@@ -7,12 +7,15 @@
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #include <algorithm>
-#include <cmath>
-#include <iostream>
+#include <cassert>
+#include <cstdint>
+#include <functional>
+#include <stdexcept>
 #include <vector>
 
+#include <boost/numeric/conversion/cast.hpp>
+
 #include <fmt/core.h>
-#include <fmt/ranges.h>
 
 #include <heyoka/config.hpp>
 #include <heyoka/expression.hpp>
@@ -22,19 +25,20 @@ HEYOKA_BEGIN_NAMESPACE
 
 namespace model::detail
 {
+
 std::vector<expression> compute_layer(std::uint32_t layer_id, const std::vector<expression> &inputs,
                                       const std::vector<std::uint32_t> &n_neurons,
                                       const std::function<expression(const expression &)> &activation,
                                       const std::vector<expression> &nn_wb, std::uint32_t n_net_w,
                                       std::uint32_t &wcounter, std::uint32_t &bcounter)
 {
-    assert(layer_id > 0);
+    assert(layer_id > 0u);
     auto n_neurons_prev_layer = boost::numeric_cast<std::uint32_t>(inputs.size());
     auto n_neurons_curr_layer = n_neurons[layer_id];
 
     std::vector<expression> retval(n_neurons_curr_layer, 0_dbl);
-    for (std::uint32_t i = 0u; i < n_neurons_curr_layer; ++i) {
-        for (std::uint32_t j = 0u; j < n_neurons_prev_layer; ++j) {
+    for (std::uint32_t i = 0; i < n_neurons_curr_layer; ++i) {
+        for (std::uint32_t j = 0; j < n_neurons_prev_layer; ++j) {
 
             // Add the weight and update the weight counter
             retval[i] += nn_wb[wcounter] * inputs[j];
@@ -69,7 +73,7 @@ std::vector<expression> ffnn_impl(const std::vector<expression> &in, const std::
     if (n_out == 0) {
         throw std::invalid_argument("The number of network outputs cannot be zero.");
     }
-    if (!std::all_of(nn_hidden.begin(), nn_hidden.end(), [](std::uint32_t item) { return item > 0; })) {
+    if (!std::all_of(nn_hidden.begin(), nn_hidden.end(), [](std::uint32_t item) { return item > 0u; })) {
         throw std::invalid_argument("The number of neurons for each hidden layer must be greater than zero!");
     }
 
@@ -108,5 +112,7 @@ std::vector<expression> ffnn_impl(const std::vector<expression> &in, const std::
     }
     return retval;
 }
+
 } // namespace model::detail
-HEYOKA_END_NAMESPACE
\ No newline at end of file
+
+HEYOKA_END_NAMESPACE

From 7d474ec12e5b1494911e5872d0c04e707b702198 Mon Sep 17 00:00:00 2001
From: Francesco Biscani <bluescarni@gmail.com>
Date: Thu, 2 Nov 2023 15:36:31 +0100
Subject: [PATCH 13/14] Review/improvements to the ffn model.

---
 include/heyoka/model/ffnn.hpp | 116 ++++++++++++++++++++--------------
 src/model/ffnn.cpp            |  69 ++++++++++++--------
 test/model_ffnn.cpp           |  37 ++++++++---
 3 files changed, 138 insertions(+), 84 deletions(-)

diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
index 62cdedbcc..ea9560c26 100644
--- a/include/heyoka/model/ffnn.hpp
+++ b/include/heyoka/model/ffnn.hpp
@@ -16,12 +16,14 @@
 #include <vector>
 
 #include <boost/numeric/conversion/cast.hpp>
+#include <boost/safe_numerics/safe_integer.hpp>
 
 #include <heyoka/config.hpp>
 #include <heyoka/detail/igor.hpp>
 #include <heyoka/detail/type_traits.hpp>
 #include <heyoka/detail/visibility.hpp>
 #include <heyoka/expression.hpp>
+#include <heyoka/kw.hpp>
 
 HEYOKA_BEGIN_NAMESPACE
 
@@ -37,75 +39,91 @@ auto ffnn_common_opts(const KwArgs &...kw_args)
 
     static_assert(!p.has_unnamed_arguments(), "This function accepts only named arguments");
 
-    // Network inputs. Mandatory
-    auto inputs = [&p]() {
-        if constexpr (p.has(kw::inputs)) {
-            return std::vector<expression>{p(kw::inputs)};
-        } else {
-            static_assert(heyoka::detail::always_false_v<KwArgs...>,
-                          "The 'inputs' keyword argument is necessary but it was not provided");
+    // Network inputs. Mandatory.
+    // The kw::inputs must be a range of values from which
+    // an expression can be constructed.
+    std::vector<expression> inputs;
+    if constexpr (p.has(kw::inputs)) {
+        for (const auto &val : p(kw::inputs)) {
+            inputs.emplace_back(val);
         }
-    }();
-
-    // Number of hidden neurons per hidden layer. Mandatory
-    auto nn_hidden = [&p]() {
-        if constexpr (p.has(kw::nn_hidden)) {
-            return std::vector<std::uint32_t>{p(kw::nn_hidden)};
-        } else {
-            static_assert(heyoka::detail::always_false_v<KwArgs...>,
-                          "The 'nn_hidden' keyword argument is necessary but it was not provided");
+    } else {
+        static_assert(heyoka::detail::always_false_v<KwArgs...>,
+                      "The 'inputs' keyword argument is necessary but it was not provided");
+    }
+
+    // Number of hidden neurons per hidden layer. Mandatory.
+    // The kw::nn_hidden argument must be a range containing
+    // integral values.
+    std::vector<std::uint32_t> nn_hidden;
+    if constexpr (p.has(kw::nn_hidden)) {
+        for (const auto &nval : p(kw::nn_hidden)) {
+            nn_hidden.push_back(boost::numeric_cast<std::uint32_t>(nval));
         }
-    }();
+    } else {
+        static_assert(heyoka::detail::always_false_v<KwArgs...>,
+                      "The 'nn_hidden' keyword argument is necessary but it was not provided");
+    }
 
-    // Number of network outputs. Mandatory
+    // Number of network outputs. Mandatory.
+    // The kw::n_out argument must be of integral type.
     auto n_out = [&p]() {
         if constexpr (p.has(kw::n_out)) {
-            return std::uint32_t{p(kw::n_out)};
+            return boost::numeric_cast<std::uint32_t>(p(kw::n_out));
         } else {
             static_assert(heyoka::detail::always_false_v<KwArgs...>,
                           "The 'n_out' keyword argument is necessary but it was not provided");
         }
     }();
 
-    // Network activation functions. Mandatory
-    auto activations = [&p]() {
-        if constexpr (p.has(kw::activations)) {
-            return std::vector<std::function<expression(const expression &)>>{p(kw::activations)};
-        } else {
-            static_assert(heyoka::detail::always_false_v<KwArgs...>,
-                          "The 'activations' keyword argument is necessary but it was not provided");
+    // Network activation functions. Mandatory.
+    // The kw::activations argument must be a range containing values
+    // from which a std::function can be constructed.
+    std::vector<std::function<expression(const expression &)>> activations;
+    if constexpr (p.has(kw::activations)) {
+        for (const auto &f : p(kw::activations)) {
+            activations.emplace_back(f);
         }
-    }();
+    } else {
+        static_assert(heyoka::detail::always_false_v<KwArgs...>,
+                      "The 'activations' keyword argument is necessary but it was not provided");
+    }
+
+    // Network weights and biases. Optional, defaults to heyoka parameters.
+    // The kw::nn_wb argument, if present, must be a range of values from which
+    // expressions can be constructed.
+    std::vector<expression> nn_wb;
+    if constexpr (p.has(kw::nn_wb)) {
+        for (const auto &val : p(kw::nn_wb)) {
+            nn_wb.emplace_back(val);
+        }
+    } else {
+        // Safe counterpart to std::uint32_t in order to avoid
+        // overflows when manipulating indices and sizes.
+        using su32 = boost::safe_numerics::safe<std::uint32_t>;
 
-    // Network weights and biases. Defaults to heyoka parameters.
-    auto nn_wb = [&p, &nn_hidden, &inputs, n_out]() {
         // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons)
-        auto n_hidden_layers = boost::numeric_cast<std::uint32_t>(nn_hidden.size());
+        auto n_hidden_layers = su32(nn_hidden.size());
         // Number of neuronal layers (counting input and output)
         auto n_layers = n_hidden_layers + 2;
         // Number of inputs
-        auto n_in = boost::numeric_cast<std::uint32_t>(inputs.size());
+        auto n_in = su32(inputs.size());
         // Number of neurons per neuronal layer
-        std::vector<std::uint32_t> n_neurons = nn_hidden;
-        n_neurons.insert(n_neurons.begin(), n_in);
+        std::vector<su32> n_neurons{n_in};
+        n_neurons.insert(n_neurons.end(), nn_hidden.begin(), nn_hidden.end());
         n_neurons.insert(n_neurons.end(), n_out);
-        std::vector<expression> retval{};
-        if constexpr (p.has(kw::nn_wb)) {
-            retval = p(kw::nn_wb);
-        } else {
-            // Number of network parameters (wb: weights and biases, w: only weights)
-            std::uint32_t n_wb = 0u;
-            for (std::uint32_t i = 1u; i < n_layers; ++i) {
-                n_wb += n_neurons[i - 1] * n_neurons[i];
-                n_wb += n_neurons[i];
-            }
-            retval.resize(n_wb);
-            for (decltype(retval.size()) i = 0u; i < retval.size(); ++i) {
-                retval[i] = heyoka::par[boost::numeric_cast<std::uint32_t>(i)];
-            }
+
+        // Number of network parameters (wb: weights and biases, w: only weights)
+        su32 n_wb = 0;
+        for (su32 i = 1; i < n_layers; ++i) {
+            n_wb += n_neurons[i - 1] * n_neurons[i];
+            n_wb += n_neurons[i];
         }
-        return retval;
-    }();
+        nn_wb.resize(n_wb);
+        for (decltype(nn_wb.size()) i = 0; i < nn_wb.size(); ++i) {
+            nn_wb[i] = par[boost::numeric_cast<std::uint32_t>(i)];
+        }
+    }
 
     return std::tuple{std::move(inputs), std::move(nn_hidden), std::move(n_out), std::move(activations),
                       std::move(nn_wb)};
diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
index 93e530645..88ce78d2b 100644
--- a/src/model/ffnn.cpp
+++ b/src/model/ffnn.cpp
@@ -13,7 +13,7 @@
 #include <stdexcept>
 #include <vector>
 
-#include <boost/numeric/conversion/cast.hpp>
+#include <boost/safe_numerics/safe_integer.hpp>
 
 #include <fmt/core.h>
 
@@ -26,19 +26,26 @@ HEYOKA_BEGIN_NAMESPACE
 namespace model::detail
 {
 
-std::vector<expression> compute_layer(std::uint32_t layer_id, const std::vector<expression> &inputs,
-                                      const std::vector<std::uint32_t> &n_neurons,
+namespace
+{
+
+// Safe counterpart to std::uint32_t in order to avoid
+// overflows when manipulating indices and sizes.
+using su32 = boost::safe_numerics::safe<std::uint32_t>;
+
+std::vector<expression> compute_layer(su32 layer_id, const std::vector<expression> &inputs,
+                                      const std::vector<su32> &n_neurons,
                                       const std::function<expression(const expression &)> &activation,
-                                      const std::vector<expression> &nn_wb, std::uint32_t n_net_w,
-                                      std::uint32_t &wcounter, std::uint32_t &bcounter)
+                                      const std::vector<expression> &nn_wb, su32 n_net_w, su32 &wcounter,
+                                      su32 &bcounter)
 {
     assert(layer_id > 0u);
-    auto n_neurons_prev_layer = boost::numeric_cast<std::uint32_t>(inputs.size());
+    auto n_neurons_prev_layer = su32(inputs.size());
     auto n_neurons_curr_layer = n_neurons[layer_id];
 
-    std::vector<expression> retval(n_neurons_curr_layer, 0_dbl);
-    for (std::uint32_t i = 0; i < n_neurons_curr_layer; ++i) {
-        for (std::uint32_t j = 0; j < n_neurons_prev_layer; ++j) {
+    std::vector<expression> retval(static_cast<std::vector<expression>::size_type>(n_neurons_curr_layer), 0_dbl);
+    for (su32 i = 0; i < n_neurons_curr_layer; ++i) {
+        for (su32 j = 0; j < n_neurons_prev_layer; ++j) {
 
             // Add the weight and update the weight counter
             retval[i] += nn_wb[wcounter] * inputs[j];
@@ -54,13 +61,18 @@ std::vector<expression> compute_layer(std::uint32_t layer_id, const std::vector<
     return retval;
 }
 
+} // namespace
+
 std::vector<expression> ffnn_impl(const std::vector<expression> &in, const std::vector<std::uint32_t> &nn_hidden,
                                   std::uint32_t n_out,
                                   const std::vector<std::function<expression(const expression &)>> &activations,
                                   const std::vector<expression> &nn_wb)
 {
     // Sanity checks
-    if (nn_hidden.size() + 1 != activations.size()) {
+    if (activations.empty()) {
+        throw std::invalid_argument("Cannot create a FFNN with an empty list of activation functions");
+    }
+    if (nn_hidden.size() != activations.size() - 1u) {
         throw std::invalid_argument(fmt::format(
             "The number of hidden layers, as detected from the inputs, was {}, while "
             "the number of activation function supplied was {}. A FFNN needs exactly one more activation function "
@@ -70,29 +82,35 @@ std::vector<expression> ffnn_impl(const std::vector<expression> &in, const std::
     if (in.empty()) {
         throw std::invalid_argument("The inputs provided to the FFNN is an empty vector.");
     }
-    if (n_out == 0) {
+    if (n_out == 0u) {
         throw std::invalid_argument("The number of network outputs cannot be zero.");
     }
-    if (!std::all_of(nn_hidden.begin(), nn_hidden.end(), [](std::uint32_t item) { return item > 0u; })) {
+    if (!std::all_of(nn_hidden.begin(), nn_hidden.end(), [](auto item) { return item > 0u; })) {
         throw std::invalid_argument("The number of neurons for each hidden layer must be greater than zero!");
     }
+    if (std::any_of(activations.begin(), activations.end(), [](const auto &func) { return !func; })) {
+        throw std::invalid_argument("The list of activation functions cannot contain empty functions");
+    }
+
+    // From now on, always use safe arithmetics to compute/manipulate
+    // indices and sizes.
+    using detail::su32;
 
     // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons)
-    auto n_hidden_layers = boost::numeric_cast<std::uint32_t>(nn_hidden.size());
+    auto n_hidden_layers = su32(nn_hidden.size());
     // Number of neuronal layers (counting input and output)
     auto n_layers = n_hidden_layers + 2;
     // Number of inputs
-    auto n_in = boost::numeric_cast<std::uint32_t>(in.size());
+    auto n_in = su32(in.size());
     // Number of neurons per neuronal layer
-    std::vector<std::uint32_t> n_neurons = nn_hidden;
-    n_neurons.insert(n_neurons.begin(), n_in);
+    std::vector<su32> n_neurons{n_in};
+    n_neurons.insert(n_neurons.end(), nn_hidden.begin(), nn_hidden.end());
     n_neurons.insert(n_neurons.end(), n_out);
     // Number of network parameters (wb: weights and biases, w: only weights)
-    std::uint32_t n_net_wb = 0u;
-    std::uint32_t n_net_w = 0u;
-    for (std::uint32_t i = 1u; i < n_layers; ++i) {
-        n_net_wb += n_neurons[i - 1] * n_neurons[i];
-        n_net_w += n_neurons[i - 1] * n_neurons[i];
+    su32 n_net_wb = 0, n_net_w = 0;
+    for (su32 i = 1; i < n_layers; ++i) {
+        n_net_wb += n_neurons[i - 1u] * n_neurons[i];
+        n_net_w += n_neurons[i - 1u] * n_neurons[i];
         n_net_wb += n_neurons[i];
     }
     // Sanity check
@@ -100,15 +118,14 @@ std::vector<expression> ffnn_impl(const std::vector<expression> &in, const std::
         throw std::invalid_argument(fmt::format(
             "The number of network parameters, detected from its structure to be {}, does not match the size of "
             "the corresponding expressions: {}.",
-            n_net_wb, nn_wb.size()));
+            static_cast<std::uint32_t>(n_net_wb), nn_wb.size()));
     }
 
     // Now we build the expressions recursively transvering from layer to layer (L = f(Wx+b)))
     std::vector<expression> retval = in;
-    std::uint32_t wcounter = 0;
-    std::uint32_t bcounter = 0;
-    for (std::uint32_t i = 1u; i < n_layers; ++i) {
-        retval = detail::compute_layer(i, retval, n_neurons, activations[i - 1], nn_wb, n_net_w, wcounter, bcounter);
+    su32 wcounter = 0, bcounter = 0;
+    for (su32 i = 1; i < n_layers; ++i) {
+        retval = detail::compute_layer(i, retval, n_neurons, activations[i - 1u], nn_wb, n_net_w, wcounter, bcounter);
     }
     return retval;
 }
diff --git a/test/model_ffnn.cpp b/test/model_ffnn.cpp
index 9cd2fe141..4a5bd8c25 100644
--- a/test/model_ffnn.cpp
+++ b/test/model_ffnn.cpp
@@ -6,21 +6,29 @@
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
+#include <cstdint>
+#include <functional>
+#include <limits>
+#include <stdexcept>
+
 #include <fmt/core.h>
 #include <fmt/ranges.h>
 
 #include <heyoka/config.hpp>
 #include <heyoka/expression.hpp>
+#include <heyoka/kw.hpp>
 #include <heyoka/math.hpp>
 #include <heyoka/model/ffnn.hpp>
+#include <system_error>
 
 #include "catch.hpp"
-#include "heyoka/kw.hpp"
 
 using namespace heyoka;
 
 TEST_CASE("impl")
 {
+    using Catch::Matchers::Message;
+
     // A linear layer, just because
     auto linear = [](expression ret) -> expression { return ret; };
     // We also define a few symbols
@@ -45,6 +53,14 @@ TEST_CASE("impl")
     REQUIRE_THROWS_AS(
         model::detail::ffnn_impl({x}, {1}, 1, {heyoka::tanh, heyoka::tanh}, {1_dbl, 2_dbl, 3_dbl, 5_dbl, 6_dbl}),
         std::invalid_argument);
+    // 6 - empty activations list
+    REQUIRE_THROWS_MATCHES(model::detail::ffnn_impl({x}, {1}, 1, {}, {1_dbl, 2_dbl, 3_dbl, 5_dbl, 6_dbl}),
+                           std::invalid_argument,
+                           Message("Cannot create a FFNN with an empty list of activation functions"));
+    // 7 - empty activation functions
+    REQUIRE_THROWS_MATCHES(
+        model::detail::ffnn_impl({x}, {}, 1, {std::function<expression(const expression &)>{}}, {1_dbl, 2_dbl}),
+        std::invalid_argument, Message("The list of activation functions cannot contain empty functions"));
 
     // We now check some hand coded networks
     {
@@ -68,28 +84,31 @@ TEST_CASE("impl")
                                                {1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl});
         REQUIRE(my_net[0] == expression(heyoka::cos(1_dbl + 2_dbl * heyoka::sin(1_dbl + x + y))));
     }
+
+    // Check overflow detection.
+    REQUIRE_THROWS_AS(
+        model::detail::ffnn_impl({x}, {}, std::numeric_limits<std::uint32_t>::max(), {linear}, {1_dbl, 2_dbl}),
+        std::system_error);
 }
 
 TEST_CASE("igor_iface")
 {
     auto [x, y, z] = make_vars("x", "y", "z");
     {
-        auto igor_v = model::ffnn(
-            kw::inputs = {x, y}, kw::nn_hidden = std::vector<std::uint32_t>{2u}, kw::n_out = 1u,
-            kw::activations = std::vector<std::function<expression(const expression &)>>{heyoka::sin, heyoka::cos},
-            kw::nn_wb = {1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl});
+        auto igor_v = model::ffnn(kw::inputs = {x, y}, kw::nn_hidden = {2u}, kw::n_out = 1u,
+                                  kw::activations = {heyoka::sin, heyoka::cos},
+                                  kw::nn_wb = {1., 1., 1., 1., 1., 1., 1., 1., 1.});
         auto vanilla_v = model::detail::ffnn_impl({x, y}, {2u}, 1u, {heyoka::sin, heyoka::cos},
                                                   {1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl, 1_dbl});
         REQUIRE(igor_v == vanilla_v);
     }
     // We test the expected setting for the default weights+biases expressions to par[i].
     {
-        auto igor_v = model::ffnn(
-            kw::inputs = {x, y}, kw::nn_hidden = std::vector<std::uint32_t>{2u}, kw::n_out = 1u,
-            kw::activations = std::vector<std::function<expression(const expression &)>>{heyoka::sin, heyoka::cos});
+        auto igor_v = model::ffnn(kw::inputs = {x, y}, kw::nn_hidden = {2u}, kw::n_out = 1u,
+                                  kw::activations = {heyoka::sin, heyoka::cos});
         auto vanilla_v
             = model::detail::ffnn_impl({x, y}, {2u}, 1u, {heyoka::sin, heyoka::cos},
                                        {par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8]});
         REQUIRE(igor_v == vanilla_v);
     }
-}
\ No newline at end of file
+}

From b93e68cade790d44c2b128350da4c5671b570dc6 Mon Sep 17 00:00:00 2001
From: Francesco Biscani <bluescarni@gmail.com>
Date: Thu, 2 Nov 2023 15:41:50 +0100
Subject: [PATCH 14/14] Minutiae.

---
 include/heyoka/model/ffnn.hpp | 12 ++++++------
 src/model/ffnn.cpp            | 22 +++++++++++-----------
 2 files changed, 17 insertions(+), 17 deletions(-)

diff --git a/include/heyoka/model/ffnn.hpp b/include/heyoka/model/ffnn.hpp
index ea9560c26..e0141d763 100644
--- a/include/heyoka/model/ffnn.hpp
+++ b/include/heyoka/model/ffnn.hpp
@@ -40,7 +40,7 @@ auto ffnn_common_opts(const KwArgs &...kw_args)
     static_assert(!p.has_unnamed_arguments(), "This function accepts only named arguments");
 
     // Network inputs. Mandatory.
-    // The kw::inputs must be a range of values from which
+    // The kw::inputs argument must be a range of values from which
     // an expression can be constructed.
     std::vector<expression> inputs;
     if constexpr (p.has(kw::inputs)) {
@@ -102,18 +102,18 @@ auto ffnn_common_opts(const KwArgs &...kw_args)
         // overflows when manipulating indices and sizes.
         using su32 = boost::safe_numerics::safe<std::uint32_t>;
 
-        // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons)
+        // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons).
         auto n_hidden_layers = su32(nn_hidden.size());
-        // Number of neuronal layers (counting input and output)
+        // Number of neuronal layers (counting input and output).
         auto n_layers = n_hidden_layers + 2;
-        // Number of inputs
+        // Number of inputs.
         auto n_in = su32(inputs.size());
-        // Number of neurons per neuronal layer
+        // Number of neurons per neuronal layer.
         std::vector<su32> n_neurons{n_in};
         n_neurons.insert(n_neurons.end(), nn_hidden.begin(), nn_hidden.end());
         n_neurons.insert(n_neurons.end(), n_out);
 
-        // Number of network parameters (wb: weights and biases, w: only weights)
+        // Number of network parameters (wb: weights and biases, w: only weights).
         su32 n_wb = 0;
         for (su32 i = 1; i < n_layers; ++i) {
             n_wb += n_neurons[i - 1] * n_neurons[i];
diff --git a/src/model/ffnn.cpp b/src/model/ffnn.cpp
index 88ce78d2b..4e5358021 100644
--- a/src/model/ffnn.cpp
+++ b/src/model/ffnn.cpp
@@ -47,15 +47,15 @@ std::vector<expression> compute_layer(su32 layer_id, const std::vector<expressio
     for (su32 i = 0; i < n_neurons_curr_layer; ++i) {
         for (su32 j = 0; j < n_neurons_prev_layer; ++j) {
 
-            // Add the weight and update the weight counter
+            // Add the weight and update the weight counter.
             retval[i] += nn_wb[wcounter] * inputs[j];
             ++wcounter;
         }
 
-        // Add the bias and update the counter
+        // Add the bias and update the counter.
         retval[i] += nn_wb[bcounter + n_net_w];
         ++bcounter;
-        // Activation function
+        // Activation function.
         retval[i] = activation(retval[i]);
     }
     return retval;
@@ -68,7 +68,7 @@ std::vector<expression> ffnn_impl(const std::vector<expression> &in, const std::
                                   const std::vector<std::function<expression(const expression &)>> &activations,
                                   const std::vector<expression> &nn_wb)
 {
-    // Sanity checks
+    // Sanity checks.
     if (activations.empty()) {
         throw std::invalid_argument("Cannot create a FFNN with an empty list of activation functions");
     }
@@ -96,24 +96,24 @@ std::vector<expression> ffnn_impl(const std::vector<expression> &in, const std::
     // indices and sizes.
     using detail::su32;
 
-    // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons)
+    // Number of hidden layers (defined as all neuronal columns that are nor input nor output neurons).
     auto n_hidden_layers = su32(nn_hidden.size());
-    // Number of neuronal layers (counting input and output)
+    // Number of neuronal layers (counting input and output).
     auto n_layers = n_hidden_layers + 2;
-    // Number of inputs
+    // Number of inputs.
     auto n_in = su32(in.size());
-    // Number of neurons per neuronal layer
+    // Number of neurons per neuronal layer.
     std::vector<su32> n_neurons{n_in};
     n_neurons.insert(n_neurons.end(), nn_hidden.begin(), nn_hidden.end());
     n_neurons.insert(n_neurons.end(), n_out);
-    // Number of network parameters (wb: weights and biases, w: only weights)
+    // Number of network parameters (wb: weights and biases, w: only weights).
     su32 n_net_wb = 0, n_net_w = 0;
     for (su32 i = 1; i < n_layers; ++i) {
         n_net_wb += n_neurons[i - 1u] * n_neurons[i];
         n_net_w += n_neurons[i - 1u] * n_neurons[i];
         n_net_wb += n_neurons[i];
     }
-    // Sanity check
+    // Sanity check.
     if (nn_wb.size() != n_net_wb) {
         throw std::invalid_argument(fmt::format(
             "The number of network parameters, detected from its structure to be {}, does not match the size of "
@@ -121,7 +121,7 @@ std::vector<expression> ffnn_impl(const std::vector<expression> &in, const std::
             static_cast<std::uint32_t>(n_net_wb), nn_wb.size()));
     }
 
-    // Now we build the expressions recursively transvering from layer to layer (L = f(Wx+b)))
+    // Now we build the expressions recursively transvering from layer to layer (L = f(Wx+b))).
     std::vector<expression> retval = in;
     su32 wcounter = 0, bcounter = 0;
     for (su32 i = 1; i < n_layers; ++i) {