Repository: darioizzo/dcgp (language: python)

Issue:
[BUG] n_eph >= 5 causes crash for symbolic regression
**Describe the bug**
Large numbers of ephemeral constants causes Python to crash when doing symbolic regression.

**To Reproduce**
Steps to reproduce the behavior:
1. Copy the code from here: http://darioizzo.github.io/dcgp/notebooks/symbolic_regression_3.html
2. Choose `n_eph` to be 5 or larger in the argument to `dcgpy.symbolic_regression`.
3. Run the code.
This causes a Python crash.

**Screenshots**
Here is the log, with verbosity=1000 (I'm not sure what the highest verbosity level is)

```
2aab23143000-2aab23144000 ---p 00000000 00:00 0 
2aab23144000-2aab23384000 rw-p 00000000 00:00 0 
2aab24000000-2aab24021000 rw-p 00000000 00:00 0 
2aab24021000-2aab28000000 ---p 00000000 00:00 0 
2aab28000000-2aab28021000 rw-p 00000000 00:00 0 
2aab28021000-2aab2c000000 ---p 00000000 00:00 0 
555555554000-5555555af000 r--p 00000000 00:2c 4221074                    /mnt/home/mcranmer/miniconda3/envs/dcgpy/bin/python3.7
5555555af000-555555788000 r-xp 0005b000 00:2c 4221074                    /mnt/home/mcranmer/miniconda3/envs/dcgpy/bin/python3.7
555555788000-55555582f000 r--p 00234000 00:2c 4221074                    /mnt/home/mcranmer/miniconda3/envs/dcgpy/bin/python3.7
55555582f000-555555832000 r--p 002da000 00:2c 4221074                    /mnt/home/mcranmer/miniconda3/envs/dcgpy/bin/python3.7
555555832000-55555589b000 rw-p 002dd000 00:2c 4221074                    /mnt/home/mcranmer/miniconda3/envs/dcgpy/bin/python3.7
55555589b000-555559556000 rw-p 00000000 00:00 0                          [heap]
7ffffffdc000-7ffffffff000 rw-p 00000000 00:00 0                          [stack]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0                  [vsyscall]

   Gen:        Fevals:     Best loss: Ndf size:  Compl.:
      0              0        1.41375         4         7
[I 22:01:59.528 NotebookApp] KernelRestarter: restarting kernel (1/5), keep random ports
```

**Environment (please complete the following information):**
 - OS: CentOS Linux release 7.8.2003 (Core)
 - Installation method: `conda create --name dcgpy --channel conda-forge python=3.7 numpy scipy matplotlib cython h5py pyfftw notebook tqdm jupyter jupyterlab mpi4py openmpi astropy ipywidgets dcgp-python -y`

Versions:

``` 
python: 3.7.5
dcgpy: 1.5
pygmo: 2.15.0
```



Discussion / hints:
I can confirm this bug. Here is another piece of code that can be used to investigate:
```python
import dcgpy
import pygmo as pg
import numpy as np

X = np.random.random( (5795, 6) )
Y = np.random.random( (5795, 1) )
ss = dcgpy.kernel_set_double(["sum", "diff", "mul", "sin", "cos"])

udp = dcgpy.symbolic_regression(
    points = X, labels = Y, kernels=ss(),
    rows = 8,
    cols = 50,
    n_eph = 6,
    levels_back = 40,
    multi_objective=True)

prob = pg.problem(udp)

uda  = dcgpy.momes4cgp(gen = 250, max_mut = 4)
alg = pg.algorithm(uda)
alg.set_verbosity(10)

pop = pg.population(prob, 20)
pop = alg.evolve(pop)
```
setting `n_eph` to any value greater than 4 results in an interpreter crash.

Fix patch:
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 11705e4e..5c9569e7 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -18,6 +18,8 @@ ADD_EXAMPLE(symbolic_regression_4)
 # UDAs
 ADD_EXAMPLE(es4cgp)
 ADD_EXAMPLE(moes4cgp)
+ADD_EXAMPLE(mes4cgp)
+ADD_EXAMPLE(momes4cgp)
 
 
 # Solving differential equations examples
diff --git a/examples/es4cgp.cpp b/examples/es4cgp.cpp
index 5de7837f..a7544fa3 100644
--- a/examples/es4cgp.cpp
+++ b/examples/es4cgp.cpp
@@ -40,7 +40,7 @@ int main(int ac, char *av[])
     bool bfe;
     po::options_description desc("Allowed options");
     desc.add_options()("help, h", "produce help message")(
-        "generations, g", po::value<unsigned>(&gen)->default_value(10000u), "number of generations")(
+        "gen, g", po::value<unsigned>(&gen)->default_value(10000u), "number of generations")(
         "max_mut, m", po::value<unsigned>(&max_mut)->default_value(15u), "maximum number of mutations allowed")(
         "n_eph, c", po::value<unsigned>(&n_eph)->default_value(3u), "number of constants in the expression")(
         "verbosity, v", po::value<unsigned>(&verbosity)->default_value(500u), "the screen frequency")(
diff --git a/examples/mes4cgp.cpp b/examples/mes4cgp.cpp
new file mode 100644
index 00000000..38472d08
--- /dev/null
+++ b/examples/mes4cgp.cpp
@@ -0,0 +1,89 @@
+#include <boost/algorithm/string.hpp>
+#include <boost/program_options.hpp>
+
+#include <iostream>
+#include <pagmo/algorithm.hpp>
+#include <pagmo/io.hpp>
+#include <pagmo/population.hpp>
+#include <pagmo/problem.hpp>
+#include <vector>
+
+// patch to make this compile in clang-cl
+#if defined(_MSC_VER) && defined(__clang__)
+#define and &&
+#define or ||
+#define not !
+#endif
+
+#include <symengine/expression.h>
+
+#if defined(_MSC_VER) && defined(__clang__)
+#undef and
+#undef or
+#undef not
+#endif
+
+#include <dcgp/algorithms/mes4cgp.hpp>
+#include <dcgp/gym.hpp>
+#include <dcgp/kernel_set.hpp>
+#include <dcgp/problems/symbolic_regression.hpp>
+
+namespace po = boost::program_options;
+using namespace dcgp;
+using namespace boost::algorithm;
+
+int main(int ac, char *av[])
+{
+    // PREAMBLE ----------------------------------------------------------------
+    // We deal with the command line options as to create a nicely commented executable
+    unsigned gen, max_mut, n_eph, pop_size, verbosity;
+    bool bfe;
+    po::options_description desc("Allowed options");
+    desc.add_options()("help, h", "produce help message")(
+        "gen, g", po::value<unsigned>(&gen)->default_value(10000u), "number of generations")(
+        "max_mut, m", po::value<unsigned>(&max_mut)->default_value(15u), "maximum number of mutations allowed")(
+        "n_eph, c", po::value<unsigned>(&n_eph)->default_value(3u), "number of constants in the expression")(
+        "verbosity, v", po::value<unsigned>(&verbosity)->default_value(500u), "the screen frequency")(
+        "pop_size, p", po::value<unsigned>(&pop_size)->default_value(4u), "population size");
+
+    po::variables_map vm;
+    po::store(po::parse_command_line(ac, av, desc), vm);
+    po::notify(vm);
+    if (vm.count("help")) {
+        std::cout << desc << "\n";
+        return 0;
+    }
+    // ---------------------------------------------------------------------------
+
+    // We load the data
+    std::vector<std::vector<double>> X, Y;
+    gym::generate_P1(X, Y);
+    // We instantiate a symbolic regression problem with no ephemeral constants.
+    symbolic_regression udp(X, Y, 1u, 16u, 17u, 2u, kernel_set<double>({"sum", "diff", "mul", "pdiv"})(), n_eph, false,
+                            0u, "MSE");
+
+    // We init the population
+    pagmo::population pop{udp, pop_size};
+
+    // And we define an evolutionary startegy
+    dcgp::mes4cgp uda(gen, max_mut, 1e-8);
+    pagmo::algorithm algo{uda};
+    algo.set_verbosity(verbosity);
+
+    // We evolve the population
+    pop = algo.evolve(pop);
+
+    // We print on screen the best found
+    auto idx = pop.best_idx();
+    auto prettier = udp.prettier(pop.get_x()[idx]);
+    trim_left_if(prettier, is_any_of("["));
+    trim_right_if(prettier, is_any_of("]"));
+    pagmo::print("\nBest fitness: ", pop.get_f()[idx], "\n");
+    pagmo::print("Chromosome: ", pop.get_x()[idx], "\n");
+    pagmo::print("Pretty Formula: ", udp.pretty(pop.get_x()[idx]), "\n");
+    pagmo::print("Prettier Formula: ", prettier, "\n");
+    pagmo::print("Expanded Formula: ", SymEngine::expand(SymEngine::Expression(prettier)), "\n");
+
+    return 0;
+}
+
diff --git a/examples/moes4cgp.cpp b/examples/moes4cgp.cpp
index 7b531cd9..3ee40a54 100644
--- a/examples/moes4cgp.cpp
+++ b/examples/moes4cgp.cpp
@@ -40,7 +40,7 @@ int main(int ac, char *av[])
     bool bfe;
     po::options_description desc("Allowed options");
     desc.add_options()("help, h", "produce help message")(
-        "generations, g", po::value<unsigned>(&gen)->default_value(10000u), "number of generations")(
+        "gen, g", po::value<unsigned>(&gen)->default_value(10000u), "number of generations")(
         "max_mut, m", po::value<unsigned>(&max_mut)->default_value(15u), "maximum number of mutations allowed")(
         "n_eph, c", po::value<unsigned>(&n_eph)->default_value(3u), "number of constants in the expression")(
         "verbosity, v", po::value<unsigned>(&verbosity)->default_value(500u), "the screen frequency")(
diff --git a/examples/momes4cgp.cpp b/examples/momes4cgp.cpp
new file mode 100644
index 00000000..c0e1712a
--- /dev/null
+++ b/examples/momes4cgp.cpp
@@ -0,0 +1,99 @@
+#include <boost/algorithm/string.hpp>
+#include <boost/program_options.hpp>
+
+#include <iostream>
+#include <pagmo/algorithm.hpp>
+#include <pagmo/io.hpp>
+#include <pagmo/population.hpp>
+#include <pagmo/problem.hpp>
+#include <vector>
+
+// patch to make this compile in clang-cl
+#if defined(_MSC_VER) && defined(__clang__)
+#define and &&
+#define or ||
+#define not !
+#endif
+
+#include <symengine/expression.h>
+
+#if defined(_MSC_VER) && defined(__clang__)
+#undef and
+#undef or
+#undef not
+#endif
+
+#include <dcgp/algorithms/momes4cgp.hpp>
+#include <dcgp/gym.hpp>
+#include <dcgp/kernel_set.hpp>
+#include <dcgp/problems/symbolic_regression.hpp>
+
+namespace po = boost::program_options;
+using namespace dcgp;
+using namespace boost::algorithm;
+
+int main(int ac, char *av[])
+{
+    // PREAMBLE ----------------------------------------------------------------
+    // We deal with the command line options as to create a nicely commented executable
+    unsigned gen, max_mut, n_eph, pop_size, verbosity;
+    bool bfe;
+    po::options_description desc("Allowed options");
+    desc.add_options()("help, h", "produce help message")(
+        "gen, g", po::value<unsigned>(&gen)->default_value(10000u), "number of generations")(
+        "max_mut, m", po::value<unsigned>(&max_mut)->default_value(15u), "maximum number of mutations allowed")(
+        "n_eph, c", po::value<unsigned>(&n_eph)->default_value(3u), "number of constants in the expression")(
+        "verbosity, v", po::value<unsigned>(&verbosity)->default_value(500u), "the screen frequency")(
+        "pop_size, p", po::value<unsigned>(&pop_size)->default_value(4u), "population size");
+
+    po::variables_map vm;
+    po::store(po::parse_command_line(ac, av, desc), vm);
+    po::notify(vm);
+    if (vm.count("help")) {
+        std::cout << desc << "\n";
+        return 0;
+    }
+    // ---------------------------------------------------------------------------
+
+    // We load the data
+    std::vector<std::vector<double>> X, Y;
+    gym::generate_P1(X, Y);
+    // We instantiate a symbolic regression problem with no ephemeral constants.
+    symbolic_regression udp(X, Y, 1u, 16u, 17u, 2u, kernel_set<double>({"sum", "diff", "mul", "pdiv"})(), n_eph, true,
+                            0u, "MSE");
+
+    // We init the population
+    pagmo::population pop{udp, pop_size};
+
+    for (decltype(pop_size) i = 0u; i < pop_size; ++i) {
+        auto prettier = udp.prettier(pop.get_x()[i]);
+        trim_left_if(prettier, is_any_of("["));
+        trim_right_if(prettier, is_any_of("]"));
+        pagmo::print(std::setw(2), i + 1, " - Loss: ", std::setw(13), std::left, pop.get_f()[i][0], std::setw(15),
+                     "Complexity: ", std::left, std::setw(5), pop.get_f()[i][1], std::setw(10), "Formula: ", prettier,
+                     "\n");
+    }
+
+    // And we define an evolutionary startegy
+    dcgp::momes4cgp uda(gen, max_mut, 1e-8);
+    pagmo::algorithm algo{uda};
+    algo.set_verbosity(verbosity);
+
+    // We evolve the population
+    pop = algo.evolve(pop);
+
+    // Finally we print on screen the non dominated front.
+    pagmo::print("\nNon dominated Front at the end:\n");
+    auto ndf = pagmo::non_dominated_front_2d(pop.get_f());
+    for (decltype(ndf.size()) i = 0u; i < ndf.size(); ++i) {
+        auto idx = ndf[i];
+        auto prettier = udp.prettier(pop.get_x()[idx]);
+        trim_left_if(prettier, is_any_of("["));
+        trim_right_if(prettier, is_any_of("]"));
+        pagmo::print(std::setw(2), i + 1, " - Loss: ", std::setw(13), std::left, pop.get_f()[idx][0], std::setw(15),
+                     "Complexity: ", std::left, std::setw(5), pop.get_f()[idx][1], std::setw(10), "Formula: ", prettier,
+                     "\n");
+    }
+
+    return 0;
+}
diff --git a/include/dcgp/algorithms/mes4cgp.hpp b/include/dcgp/algorithms/mes4cgp.hpp
index 0da05e86..3a218500 100644
--- a/include/dcgp/algorithms/mes4cgp.hpp
+++ b/include/dcgp/algorithms/mes4cgp.hpp
@@ -140,6 +140,8 @@ class mes4cgp
                        [](double a) { return boost::numeric_cast<unsigned>(a); });
         // The Hessian sparsity is defined here and it will not change.
         auto hs = prob.hessians_sparsity();
+        // Hessian 
+        Eigen::MatrixXd fullH;
         // (active) Hessian
         Eigen::MatrixXd H;
         // (active) Gradient
@@ -173,7 +175,7 @@ class mes4cgp
             std::vector<pagmo::vector_double> mutated_f(NP, best_f);
             for (decltype(NP) i = 0u; i < NP; ++i) {
                 cgp.set(best_xu);
-                cgp.mutate_active(dis(m_e));
+                cgp.mutate_random(dis(m_e));
                 std::vector<unsigned> mutated_xu = cgp.get();
                 std::transform(mutated_xu.begin(), mutated_xu.end(), mutated_x[i].data() + n_eph,
                                [](unsigned a) { return boost::numeric_cast<double>(a); });
@@ -201,9 +203,10 @@ class mes4cgp
                     std::vector<pagmo::vector_double::size_type> non_zero, zero;
                     for (decltype(grad.size()) j = 0u; j < grad.size(); ++j) {
                         if (grad[j] != 0.) {
-                            non_zero.emplace_back(j);
+                            non_zero.push_back(j);
                         } else {
-                            zero.emplace_back(j);
+                            // if the constant is inactive we reset it randomly in the best chromosome
+                            best_x[j] = std::uniform_real_distribution<double>(-10., 10.)(m_e);
                         }
                     }
                     auto n_non_zero = non_zero.size();
@@ -216,24 +219,32 @@ class mes4cgp
                                 = mutated_x[i][non_zero[0]] - grad[non_zero[0]] / hess[0][non_zero[0]];
                         }
                     } else if (n_non_zero > 1u) {
-                        // The (active) Hessian stored in a square Eigen matrix. Rows and Cols
-                        // will be removed later so that size will be n_non_zero x n_non_zero
-                        H = Eigen::MatrixXd::Zero(_(n_eph), _(n_eph));
+                        // The full Hessian
+                        fullH = Eigen::MatrixXd::Zero(_(n_eph), _(n_eph));
+                        // (active) Hessian
+                        H = Eigen::MatrixXd::Zero(_(n_non_zero), _(n_non_zero));
                         // (active) Gradient stored in a column vector
                         G = Eigen::MatrixXd::Zero(_(n_non_zero), 1u);
                         // (active) Constants stored in a column vector
                         C = Eigen::MatrixXd::Zero(_(n_non_zero), 1u);
                         // We construct the hessian including all the zero elements
                         for (decltype(hess[0].size()) j = 0u; j < hess[0].size(); ++j) {
-                            H(_(hs[0][j].first), _(hs[0][j].second)) = hess[0][j];
-                            H(_(hs[0][j].second), _(hs[0][j].first)) = hess[0][j];
+                            fullH(_(hs[0][j].first), _(hs[0][j].second)) = hess[0][j];
+                            fullH(_(hs[0][j].second), _(hs[0][j].first)) = hess[0][j];
                         }
-                        // We remove the rows and columns corresponding to ephemeral constants not
-                        // in the expression
-                        for (decltype(zero.size()) j = 0u; j < zero.size(); ++j) {
-                            remove_row(H, non_zero[j] - j);
-                            remove_column(H, non_zero[j] - j);
+                        // The following nested fors loops copy into the active hessian all cols and rows tht are non zero, i.e.
+                        // corresponding to eph constants that actually are in the expression
+                        // Probably this can be coded more efficiently using operations on rows and columns instead.
+                        decltype(non_zero.size()) new_row_idx = 0u;
+                        for (decltype(non_zero.size()) row_idx = 0u; row_idx < non_zero.size(); ++row_idx) {
+                            decltype(non_zero.size()) new_col_idx = 0u;
+                            for (decltype(non_zero.size()) col_idx = 0u; col_idx < non_zero.size(); ++col_idx) {
+                                H(_(new_row_idx), _(new_col_idx)) = fullH(_(non_zero[row_idx]), _(non_zero[col_idx]));
+                                new_col_idx++;
+                            }
+                            new_row_idx++;
                         }
+
                         // We now construct the (active) gradient and the (active) constant column vectors.
                         for (decltype(n_non_zero) j = 0u; j < n_non_zero; ++j) {
                             G(_(j), 0) = grad[non_zero[j]];
@@ -399,29 +410,6 @@ class mes4cgp
     }
 
 private:
-    static inline void remove_row(Eigen::MatrixXd &matrix, pagmo::vector_double::size_type row_to_remove)
-    {
-        Eigen::Index n_rows = matrix.rows() - 1u;
-        Eigen::Index n_cols = matrix.cols();
-
-        if (_(row_to_remove) < n_rows)
-            matrix.block(_(row_to_remove), 0, n_rows - _(row_to_remove), n_cols)
-                = matrix.block(_(row_to_remove) + 1, 0, n_rows - _(row_to_remove), n_cols);
-
-        matrix.conservativeResize(n_rows, n_cols);
-    }
-
-    static inline void remove_column(Eigen::MatrixXd &matrix, pagmo::vector_double::size_type col_to_remove)
-    {
-        Eigen::Index n_rows = matrix.rows();
-        Eigen::Index n_cols = matrix.cols() - 1u;
-
-        if (_(col_to_remove) < n_cols)
-            matrix.block(0, _(col_to_remove), n_rows, n_cols - _(col_to_remove))
-                = matrix.block(0, _(col_to_remove) + 1, n_rows, n_cols - _(col_to_remove));
-
-        matrix.conservativeResize(n_rows, n_cols);
-    }
     // Eigen stores indexes and sizes as signed types, while dCGP (and pagmo)
     // uses STL containers thus sizes and indexes are unsigned. To
     // make the conversion as painless as possible this template is provided
diff --git a/include/dcgp/algorithms/momes4cgp.hpp b/include/dcgp/algorithms/momes4cgp.hpp
index e3530a2f..34e37e58 100644
--- a/include/dcgp/algorithms/momes4cgp.hpp
+++ b/include/dcgp/algorithms/momes4cgp.hpp
@@ -125,12 +125,19 @@ class momes4cgp
         auto cgp = udp_ptr->get_cgp();
         // How many ephemeral constants?
         auto n_eph = prob.get_ncx();
-        // The hessian will be stored in a square Eigen for inversion
-        Eigen::MatrixXd H = Eigen::MatrixXd::Zero(_(n_eph), _(n_eph));
-        // Gradient and Constants will be stored in these column vectors
-        Eigen::MatrixXd G = Eigen::MatrixXd::Zero(_(n_eph), 1);
-        Eigen::MatrixXd C = Eigen::MatrixXd::Zero(_(n_eph), 1);
+        // The Hessian sparsity is defined here and it will not change.
         auto hs = prob.hessians_sparsity();
+        // Hessian 
+        Eigen::MatrixXd fullH;
+        // (active) Hessian
+        Eigen::MatrixXd H;
+        // (active) Gradient
+        Eigen::MatrixXd G;
+        // (active) Constants
+        Eigen::MatrixXd C;
+        // Full pivoting LU decomposition to check that H is invertible, find the inverse
+        // and see if H is positive definite
+        Eigen::FullPivLU<Eigen::MatrixXd> fullpivlu;
 
         // Main loop
         for (decltype(m_gen) gen = 1u; gen <= m_gen; ++gen) {
@@ -178,7 +185,7 @@ class momes4cgp
                 // Use it to set the CGP
                 cgp.set(mutated_xu);
                 // Mutate the expression
-                cgp.mutate_active(n_active_mutations[i]);
+                cgp.mutate_random(n_active_mutations[i]);
                 mutated_xu = cgp.get();
                 // Put it back
                 std::transform(mutated_xu.begin(), mutated_xu.end(), mutated_x[i].data() + n_eph,
@@ -203,9 +210,10 @@ class momes4cgp
                     std::vector<pagmo::vector_double::size_type> non_zero, zero;
                     for (decltype(grad.size()) j = 0u; j < grad.size(); ++j) {
                         if (grad[j] != 0.) {
-                            non_zero.emplace_back(j);
+                            non_zero.push_back(j);
                         } else {
-                            zero.emplace_back(j);
+                            // if the constant is inactive we reset it randomly in the mutant chromosome
+                            mutated_x[i][j] = std::uniform_real_distribution<double>(-10., 10.)(m_e);
                         }
                     }
                     auto n_non_zero = non_zero.size();
@@ -218,27 +226,32 @@ class momes4cgp
                                 = mutated_x[i][non_zero[0]] - grad[non_zero[0]] / hess[0][non_zero[0]];
                         }
                     } else if (n_non_zero > 1u) {
-                        // The (active) Hessian stored in a square Eigen matrix. Rows and Cols
-                        // will be removed later so that size will be n_non_zero x n_non_zero
-                        H = Eigen::MatrixXd::Zero(_(n_eph), _(n_eph));
+                        // The full Hessian
+                        fullH = Eigen::MatrixXd::Zero(_(n_eph), _(n_eph));
+                        // (active) Hessian
+                        H = Eigen::MatrixXd::Zero(_(n_non_zero), _(n_non_zero));
                         // (active) Gradient stored in a column vector
                         G = Eigen::MatrixXd::Zero(_(n_non_zero), 1u);
                         // (active) Constants stored in a column vector
                         C = Eigen::MatrixXd::Zero(_(n_non_zero), 1u);
-                        // Full pivoting LU decomposition to check that H is invertible, find the inverse
-                        // and see if H is positive definite
-                        Eigen::FullPivLU<Eigen::MatrixXd> fullpivlu;
                         // We construct the hessian including all the zero elements
                         for (decltype(hess[0].size()) j = 0u; j < hess[0].size(); ++j) {
-                            H(_(hs[0][j].first), _(hs[0][j].second)) = hess[0][j];
-                            H(_(hs[0][j].second), _(hs[0][j].first)) = hess[0][j];
+                            fullH(_(hs[0][j].first), _(hs[0][j].second)) = hess[0][j];
+                            fullH(_(hs[0][j].second), _(hs[0][j].first)) = hess[0][j];
                         }
-                        // We remove the rows and columns corresponding to ephemeral constants not
-                        // in the expression
-                        for (decltype(zero.size()) j = 0u; j < zero.size(); ++j) {
-                            remove_row(H, non_zero[j] - j);
-                            remove_column(H, non_zero[j] - j);
+                        // The following nested fors loops copy into the active hessian all cols and rows tht are non zero, i.e.
+                        // corresponding to eph constants that actually are in the expression
+                        // Probably this can be coded more efficiently using operations on rows and columns instead.
+                        decltype(non_zero.size()) new_row_idx = 0u;
+                        for (decltype(non_zero.size()) row_idx = 0u; row_idx < non_zero.size(); ++row_idx) {
+                            decltype(non_zero.size()) new_col_idx = 0u;
+                            for (decltype(non_zero.size()) col_idx = 0u; col_idx < non_zero.size(); ++col_idx) {
+                                H(_(new_row_idx), _(new_col_idx)) = fullH(_(non_zero[row_idx]), _(non_zero[col_idx]));
+                                new_col_idx++;
+                            }
+                            new_row_idx++;
                         }
+
                         // We now construct the (active) gradient and the (active) constant column vectors.
                         for (decltype(n_non_zero) j = 0u; j < n_non_zero; ++j) {
                             G(_(j), 0) = grad[non_zero[j]];
diff --git a/include/dcgp/expression.hpp b/include/dcgp/expression.hpp
index aaf5647e..930773bc 100644
--- a/include/dcgp/expression.hpp
+++ b/include/dcgp/expression.hpp
@@ -138,11 +138,11 @@ class expression
         for (auto i = 0u; i < m_x.size(); ++i) {
             m_x[i] = std::uniform_int_distribution<unsigned>(m_lb[i], m_ub[i])(m_e);
         }
-        // We init the ephemeral constants with 1, 2, 3, 4, 5 ...
-        for (auto i = 1u; i <= n_eph; ++i) {
-            m_eph_val.push_back(static_cast<T>(i));
+        // We init the ephemeral constants in [-10, 10]
+        for (auto i = 0u; i < n_eph; ++i) {
+            m_eph_val.push_back(static_cast<T>(std::uniform_real_distribution<double>(-10., 10.)(m_e)));
         }
-        // We init the ephemeral constants with c1, c2, c3, c4, c5 ...
+        // We init the ephemeral constants with c1,c2,c3,c4,c5 ...
         for (auto i = 1u; i <= n_eph; ++i) {
             m_eph_symb.emplace_back("c" + std::to_string(i));
         }


Test patch:
diff --git a/dcgpy/test.py b/dcgpy/test.py
index f2bbb2cd..1785d71b 100644
--- a/dcgpy/test.py
+++ b/dcgpy/test.py
@@ -226,7 +226,6 @@ def test_double(self):
                         n_eph=2,
                         seed=33)
         # Ephemeral value attributes tests
-        self.assertEqual(ex.eph_val, [1, 2])
         self.assertEqual(ex.eph_symb, ["c1", "c2"])
         ex.eph_val = [-0.2, 0.3]
         self.assertEqual(ex.eph_val, [-0.2, 0.3])
diff --git a/test/expression.cpp b/test/expression.cpp
index a71d40af..a1ffc887 100644
--- a/test/expression.cpp
+++ b/test/expression.cpp
@@ -372,7 +372,6 @@ BOOST_AUTO_TEST_CASE(ephemeral_constants_test)
     {
         kernel_set<double> basic_set({"sum", "diff", "mul", "div"});
         expression<double> ex(2, 2, 1, 10, 11, 2, basic_set(), 2u, 123u);
-        BOOST_CHECK((ex.get_eph_val() == std::vector<double>{1., 2.}));
         ex.set_eph_val({3.1, -1.2});
         BOOST_CHECK((ex.get_eph_val() == std::vector<double>{3.1, -1.2}));
         BOOST_CHECK((ex.get_eph_symb() == std::vector<std::string>{"c1", "c2"}));