diff --git a/include/barry/barraydense-meat.hpp b/include/barry/barraydense-meat.hpp
index 29eb109f2..6df606996 100644
--- a/include/barry/barraydense-meat.hpp
+++ b/include/barry/barraydense-meat.hpp
@@ -449,10 +449,10 @@ template<typename Cell_Type, typename Data_Type> inline void BArrayDense<Cell_Ty
     if (check_bounds) 
         out_of_range(0u, i);
 
+    #pragma code_align 32
     #if defined(__OPENMP) || defined(_OPENMP)
     #pragma omp simd
     #endif
-    #pragma code_align 32
     for (size_t j = 0u; j < N; ++j) 
         x->operator[](j) = el[POS(j, i)];//this->get_cell(iter->first, i, false);
     
diff --git a/include/barry/model-bones.hpp b/include/barry/model-bones.hpp
index 4f32a7b23..69a43b49c 100644
--- a/include/barry/model-bones.hpp
+++ b/include/barry/model-bones.hpp
@@ -59,9 +59,11 @@ class Model {
      * - term k
      */
     ///@{
-    std::vector< std::vector< double > > stats_support;          ///< Sufficient statistics of the model (support)
-    std::vector< size_t >                stats_support_n_arrays; ///< Number of arrays included per support.
-    std::vector< std::vector< double > > stats_target;           ///< Target statistics of the model
+    std::vector< double >                stats_support;           ///< Sufficient statistics of the model (support)
+    std::vector< size_t >                stats_support_sizes;     ///< Number of vectors included in the support.
+    std::vector< size_t >                stats_support_sizes_acc; ///< Accumulated number of vectors included in the support.
+    std::vector< size_t >                stats_support_n_arrays;  ///< Number of arrays included per support.
+    std::vector< std::vector< double > > stats_target;            ///< Target statistics of the model
     std::vector< size_t >                arrays2support;
     ///@}
 
@@ -378,7 +380,7 @@ class Model {
      */
     ///@{
     std::vector< std::vector< double > > * get_stats_target();
-    std::vector< std::vector< double > > * get_stats_support();
+    std::vector< double > * get_stats_support();
     std::vector< size_t > * get_arrays2support();
     std::vector< std::vector< Array_Type > > * get_pset_arrays();
     std::vector< std::vector<double> > * get_pset_stats();  ///< Statistics of the support(s)
diff --git a/include/barry/model-meat.hpp b/include/barry/model-meat.hpp
index 325e10b2f..6d2ba9635 100644
--- a/include/barry/model-meat.hpp
+++ b/include/barry/model-meat.hpp
@@ -136,24 +136,29 @@ inline void Model<Array_Type, Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_T
     size_t ncores
 ) {
 
+    const size_t n = stats_support_sizes.size();
+
     // Barrier to make sure paralelization makes sense
-    if ((ncores > 1u) && (stats_support.size() < 1000u))
+    if ((ncores > 1u) && (n < 128u))
         ncores = 1u;
     
     #if defined(__OPENMP) || defined(_OPENMP)
     #pragma omp parallel for firstprivate(params) num_threads(ncores) \
-        shared(stats_support, normalizing_constants, first_calc_done) \
+        shared(n, normalizing_constants, first_calc_done, \
+            stats_support_sizes, stats_support_sizes_acc) \
         default(none)
     #endif
-    for (size_t i = 0u; i < stats_support.size(); ++i)
+    for (size_t i = 0u; i < n; ++i)
     {
 
         size_t k = params.size() + 1u;
-        size_t n = stats_support[i].size() / k;
+        size_t n = stats_support_sizes[i];
 
         first_calc_done[i] = true;
         normalizing_constants[i] = update_normalizing_constant(
-            params, &stats_support[i][0u], k, n
+            params, &stats_support[
+                stats_support_sizes_acc[i * k]
+                ], k, n
         );
 
     }
@@ -170,6 +175,8 @@ template <
     >
 inline Model<Array_Type,Data_Counter_Type,Data_Rule_Type, Data_Rule_Dyn_Type>::Model() :
     stats_support(0u),
+    stats_support_sizes(0u),
+    stats_support_sizes_acc(0u),
     stats_support_n_arrays(0u),
     stats_target(0u), arrays2support(0u),
     keys2support(0u),
@@ -206,6 +213,8 @@ inline Model<Array_Type,Data_Counter_Type,Data_Rule_Type, Data_Rule_Dyn_Type>::M
     size_t size_
     ) :
     stats_support(0u),
+    stats_support_sizes(0u),
+    stats_support_sizes_acc(0u),
     stats_support_n_arrays(0u),
     stats_target(0u), arrays2support(0u), keys2support(0u), 
     pset_arrays(0u), pset_stats(0u),
@@ -244,6 +253,8 @@ inline Model<Array_Type,Data_Counter_Type,Data_Rule_Type, Data_Rule_Dyn_Type>::M
     const Model<Array_Type,Data_Counter_Type,Data_Rule_Type, Data_Rule_Dyn_Type> & Model_
     ) : 
     stats_support(Model_.stats_support),
+    stats_support_sizes(Model_.stats_support_sizes),
+    stats_support_sizes_acc(Model_.stats_support_sizes_acc),
     stats_support_n_arrays(Model_.stats_support_n_arrays),
     stats_target(Model_.stats_target),
     arrays2support(Model_.arrays2support),
@@ -301,6 +312,8 @@ inline Model<Array_Type,Data_Counter_Type,Data_Rule_Type, Data_Rule_Dyn_Type> &
             delete rules_dyn;
         
         stats_support              = Model_.stats_support;
+        stats_support_sizes        = Model_.stats_support_sizes;
+        stats_support_sizes_acc    = Model_.stats_support_sizes_acc;
         stats_support_n_arrays     = Model_.stats_support_n_arrays;
         stats_target               = Model_.stats_target;
         arrays2support             = Model_.arrays2support;
@@ -498,11 +511,14 @@ inline size_t Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
     MapVec_type< double, size_t >::const_iterator locator = keys2support.find(key);
     if (force_new | (locator == keys2support.end()))
     {
+
+        // Current size of the support stats
+        size_t stats_support_size = stats_support.size();
         
         // Adding to the map
-        keys2support[key] = stats_support.size();
+        keys2support[key] = stats_support_sizes.size();
         stats_support_n_arrays.push_back(1u);       // How many elements now
-        arrays2support.push_back(stats_support.size()); // Map of the array id to the support
+        arrays2support.push_back(stats_support_sizes.size()); // Map of the array id to the support
         
         // Computing support using the counters included in the model
         support_fun.reset_array(Array_);
@@ -583,16 +599,28 @@ inline size_t Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
 
             }
 
-            stats_support.push_back(s_new);
+            for (auto & s : s_new)
+                stats_support.push_back(s);
+            
 
-        } else 
-            stats_support.push_back(support_fun.get_counts());
+        } else {
+            for (const auto & s: support_fun.get_counts())
+                stats_support.push_back(s);
+        }
         
         // Making room for the previous parameters. This will be used to check if
         // the normalizing constant has been updated or not.
         params_last.push_back(stats_target[0u]);
         normalizing_constants.push_back(0.0);
         first_calc_done.push_back(false);
+
+        // Incrementing the size of the support set
+        stats_support_sizes.push_back(
+            
+            (stats_support.size() - stats_support_size)/
+                (counter_fun.size() + 1u)
+
+            );
         
         return arrays2support.size() - 1u;
         
@@ -622,8 +650,18 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
 
     size_t idx = arrays2support[i];
 
+    // Cummulative sum of the stats_support_sizes
+    std::vector< size_t > stats_support_sizes_acc;
+    stats_support_sizes_acc.reserve(stats_support_sizes.size());
+    size_t acc = 0u;
+    for (size_t i = 0u; i < stats_support_sizes.size(); ++i)
+    {
+        acc += stats_support_sizes[i];
+        stats_support_sizes_acc.push_back(acc);
+    }
+
     // Checking if this actually has a change of happening
-    if (this->stats_support[idx].size() == 0u)
+    if (this->stats_support_sizes[idx] == 0u)
         return as_log ? -std::numeric_limits<double>::infinity() : 0.0;
     
     // Checking if we have updated the normalizing constant or not
@@ -633,10 +671,12 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
         first_calc_done[idx] = true;
         
         size_t k = params.size() + 1u;
-        size_t n = stats_support[idx].size() / k;
+        size_t n = stats_support_sizes[idx];
 
         normalizing_constants[idx] = update_normalizing_constant(
-            params, &stats_support[idx][0u], k, n
+            params, &stats_support[
+                stats_support_sizes_acc[idx] * k
+                ], k, n
         );
         
         params_last[idx] = params;
@@ -691,7 +731,7 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
     }
 
     // Checking if this actually has a change of happening
-    if (this->stats_support[loc].size() == 0u)
+    if (this->stats_support_sizes[loc] == 0u)
         return as_log ? -std::numeric_limits<double>::infinity() : 0.0;
     
     // Counting stats_target
@@ -711,10 +751,10 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
         first_calc_done[loc] = true;
 
         size_t k = params.size() + 1u;
-        size_t n = stats_support[loc].size() / k;
+        size_t n = stats_support_sizes[loc];
         
         normalizing_constants[loc] = update_normalizing_constant(
-            params, &stats_support[loc][0u], k, n
+            params, &stats_support[loc], k, n
         );
         
         params_last[loc] = params;
@@ -768,7 +808,7 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
         
 
     // Checking if this actually has a change of happening
-    if (this->stats_support[loc].size() == 0u)
+    if (this->stats_support_sizes[loc] == 0u)
     {
         throw std::logic_error("The support set for this array is empty.");
     }
@@ -779,10 +819,12 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
         first_calc_done[loc] = true;
         
         size_t k = params.size() + 1u;
-        size_t n = stats_support[loc].size() / k;
+        size_t n = stats_support_sizes[loc];
 
         normalizing_constants[loc] = update_normalizing_constant(
-            params, &stats_support[loc][0u], k, n
+            params, &stats_support[
+                stats_support_sizes_acc[loc] * k
+                ], k, n
         );
         
         params_last[loc] = params;
@@ -838,7 +880,7 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
     }
 
     // Checking if this actually has a change of happening
-    if (this->stats_support[loc].size() == 0u)
+    if (this->stats_support_sizes[loc] == 0u)
     {
         throw std::logic_error("The support set for this array is empty.");
     }
@@ -849,10 +891,12 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
         first_calc_done[loc] = true;
         
         size_t k = params.size() + 1u;
-        size_t n = stats_support[loc].size() / k;
+        size_t n = stats_support_sizes[loc];
 
         normalizing_constants[loc] = update_normalizing_constant(
-            params, &stats_support[loc][0u], k, n
+            params, &stats_support[
+                stats_support_sizes_acc[loc] * k
+            ], k, n
         );
         
         params_last[loc] = params;
@@ -883,7 +927,8 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
     {
         #if defined(__OPENMP) || defined(_OPENMP)
         #pragma omp parallel for num_threads(ncores) \
-            shared(normalizing_constants, params_last, first_calc_done, stats_support) \
+            shared(normalizing_constants, params_last, first_calc_done, \
+                stats_support, stats_support_sizes, stats_support_sizes_acc) \
             firstprivate(params)
         #endif
         for (size_t i = 0u; i < params_last_size; ++i)
@@ -893,11 +938,13 @@ inline double Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
             {
 
                 size_t k = params.size() + 1u;
-                size_t n = stats_support[i].size() / k;
+                size_t n = stats_support_sizes[i];
                 
                 first_calc_done[i] = true;
                 normalizing_constants[i] = update_normalizing_constant(
-                    params, &stats_support[i][0u], k, n
+                    params, &stats_support[
+                        stats_support_sizes_acc[i] * k
+                    ], k, n
                 );
                 
                 params_last[i] = params;
@@ -988,20 +1035,30 @@ inline void Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_Ty
     if (i >= arrays2support.size())
         throw std::range_error("The requested support is out of range");
 
-    const auto & S = stats_support[arrays2support[i]];
+    // const auto & S = stats_support[arrays2support[i]];
+    size_t array_id = arrays2support[i];
 
     size_t k       = nterms();
-    size_t nunique = S.size() / (k + 1u);
+    size_t nunique = stats_support_sizes.size();
 
     for (size_t l = 0u; l < nunique; ++l)
     {
 
         printf_barry("% 5li ", l);
 
-        printf_barry("counts: %.0f motif: ", S[l * (k + 1u)]);
+        printf_barry("counts: %.0f motif: ", stats_support[
+            stats_support_sizes_acc[l] * (k + 1u) 
+            // l * (k + 1u)
+            ]);
         
         for (size_t j = 0u; j < k; ++j)
-            printf_barry("%.2f, ", S[l * (k + 1) + j + 1]);
+        {
+            printf_barry(
+                "%.2f, ",
+                stats_support_sizes_acc[
+                    stats_support_sizes_acc[l] * (k + 1) + j + 1
+                    ]);
+        }
 
         printf_barry("\n");
 
@@ -1023,14 +1080,14 @@ inline void Model<Array_Type,Data_Counter_Type,Data_Rule_Type, Data_Rule_Dyn_Typ
     int min_v = std::numeric_limits<int>::max();
     int max_v = 0;
 
-    for (const auto & stat : this->stats_support)
+    for (const auto & stat : this->stats_support_sizes)
     {
 
-        if (static_cast<int>(stat.size()) > max_v)
-            max_v = static_cast<int>(stat.size());
+        if (static_cast<int>(stat) > max_v)
+            max_v = static_cast<int>(stat);
         
-        if (static_cast<int>(stat.size()) < min_v)
-            min_v = static_cast<int>(stat.size());
+        if (static_cast<int>(stat) < min_v)
+            min_v = static_cast<int>(stat);
 
     }  
 
@@ -1085,7 +1142,7 @@ inline size_t Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
 {
 
     // INITIALIZED()
-    return this->stats_support.size();
+    return this->stats_support_sizes.size();
 
 } 
 
@@ -1121,11 +1178,12 @@ inline size_t Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_
 {
 
     // INITIALIZED()
-    size_t tot = 0u;
-    for (auto& a : stats_support)
-        tot += a.size();
+    return stats_support_sizes_acc.back();
+    // size_t tot = 0u;
+    // for (auto& a : stats_support)
+    //     tot += a.size();
 
-    return tot;
+    // return tot;
 
 }
 
@@ -1247,9 +1305,9 @@ inline Array_Type Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_
         // throw std::out_of_range("Sampling from an array that has no support in the model.");
 
         // Adding to the map
-        keys2support[key] = stats_support.size();
+        keys2support[key] = stats_support_sizes.size();
         stats_support_n_arrays.push_back(1u);       // How many elements now
-        arrays2support.push_back(stats_support.size()); // Map of the array id to the support
+        arrays2support.push_back(stats_support_sizes.size()); // Map of the array id to the support
         
         // Computing support using the counters included in the model
         support_fun.reset_array(Array_);
@@ -1311,10 +1369,16 @@ inline Array_Type Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_
 
             }
 
-            stats_support.push_back(s_new);
+            for (auto & s : s_new)
+                stats_support.push_back(s);
+            // stats_support.push_back(s_new);
+
+        } else {
+            for (auto & s : support_fun.get_counts())
+                stats_support.push_back(s);
 
-        } else 
-            stats_support.push_back(support_fun.get_counts());
+            // stats_support.push_back(support_fun.get_counts());
+        }
         
         // Making room for the previous parameters. This will be used to check if
         // the normalizing constant has been updated or not.
@@ -1461,7 +1525,7 @@ inline std::vector< std::vector< double > > * Model<Array_Type,Data_Counter_Type
 }
 
 template <typename Array_Type, typename Data_Counter_Type, typename Data_Rule_Type, typename Data_Rule_Dyn_Type>
-inline std::vector< std::vector< double > > *
+inline std::vector< double > *
 Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_Type>::get_stats_support()
 {
     return &stats_support;
@@ -1508,35 +1572,42 @@ Model<Array_Type,Data_Counter_Type, Data_Rule_Type, Data_Rule_Dyn_Type>::set_tra
 
     size_t k = counters->size(); 
 
+    auto stats_support_old = get_stats_support();
+
     // Applying over the support
-    for (auto & s : stats_support)
+    for (auto & n : stats_support_sizes)
     {
 
-        // Making room for the new support
-        std::vector< double > s_new(0u);
-        s_new.reserve(s.size());
-
-        size_t n = s.size() / (k + 1u);
-
         // Iterating through the unique sets
+        // size_t n = s;
         for (size_t i = 0; i < n; ++i)
         {
 
-            // Appending size
-            s_new.push_back(s[i * (k + 1u)]);
-
             // Applying transformation and adding to the new set
-            auto res = transform_model_fun(&s[i * (k + 1u) + 1u], k);
+            auto res = transform_model_fun(
+                &stats_support_old[
+                    stats_support_sizes_acc[i] * (k + 1u) + 1u
+                ],
+                k
+                );
 
             if (res.size() != transform_model_term_names.size())
                 throw std::length_error("The transform vector from -transform_model_fun- does not match the size of -transform_model_term_names-.");
 
-            std::copy(res.begin(), res.end(), std::back_inserter(s_new));
+            // Weigth
+            stats_support[
+                stats_support_sizes_acc[i] * (res.size() + 1u)
+                ] = stats_support_old[
+                    stats_support_sizes_acc[i] * (k + 1u) 
+                ];
 
-        }
+            // Copying the rest of the elements
+            for (size_t j = 0u; j < res.size(); ++j)
+                stats_support[
+                    stats_support_sizes_acc[i] * (res.size() + 1u) + j + 1u
+                    ] = res[j];
 
-        // Exchanging with the original
-        std::swap(s, s_new);
+        }
 
     }
 
diff --git a/include/barry/models/geese/flock-bones.hpp b/include/barry/models/geese/flock-bones.hpp
index 425db5256..fb06d92a7 100644
--- a/include/barry/models/geese/flock-bones.hpp
+++ b/include/barry/models/geese/flock-bones.hpp
@@ -53,7 +53,7 @@ class Flock {
     // void add_geese(Geese x);
     PhyloCounters * get_counters();
     PhyloSupport *  get_support_fun();
-    std::vector< std::vector< double > > * get_stats_support();
+    std::vector< double > * get_stats_support();
     std::vector< std::vector< double > > * get_stats_target();
     PhyloModel *  get_model();
 
diff --git a/include/barry/models/geese/flock-meat.hpp b/include/barry/models/geese/flock-meat.hpp
index 3d5c1c173..c663484a5 100644
--- a/include/barry/models/geese/flock-meat.hpp
+++ b/include/barry/models/geese/flock-meat.hpp
@@ -114,7 +114,7 @@ inline PhyloSupport *  Flock::get_support_fun()
 
 }
 
-inline std::vector< std::vector< double > > *  Flock::get_stats_support()
+inline std::vector< double > *  Flock::get_stats_support()
 {
 
     return this->model.get_stats_support();
diff --git a/include/barry/support-bones.hpp b/include/barry/support-bones.hpp
index 726b665e1..a87180c83 100644
--- a/include/barry/support-bones.hpp
+++ b/include/barry/support-bones.hpp
@@ -74,9 +74,9 @@ class Support {
     size_t max_num_elements = BARRY_MAX_NUM_ELEMENTS;
     
     // Temp variables to reduce memory allocation
-    std::vector< double >                current_stats;
-    std::vector< size_t >                coordinates_free;
-    std::vector< size_t >                coordinates_locked;
+    std::vector< double > current_stats;
+    std::vector< size_t > coordinates_free;
+    std::vector< size_t > coordinates_locked;
     size_t coordiantes_n_free;
     size_t coordiantes_n_locked;
     std::vector< double > change_stats;
@@ -184,6 +184,7 @@ class Support {
     );
     
     std::vector< double > get_counts() const;
+    const std::vector< double > & get_counts() const;
     std::vector< double > * get_current_stats(); ///< List current statistics.
     void print() const;
     
diff --git a/include/barry/support-meat.hpp b/include/barry/support-meat.hpp
index 60ed03378..04eb6426c 100644
--- a/include/barry/support-meat.hpp
+++ b/include/barry/support-meat.hpp
@@ -561,6 +561,13 @@ inline std::vector< double > Support<Array_Type,Data_Counter_Type,Data_Rule_Type
     
 }
 
+template <typename Array_Type, typename Data_Counter_Type, typename Data_Rule_Type, typename Data_Rule_Dyn_Type>
+inline const std::vector< double > & Support<Array_Type,Data_Counter_Type,Data_Rule_Type, Data_Rule_Dyn_Type>::get_counts() const {
+    
+    return data.get_data(); 
+    
+}
+
 // template <typename Array_Type, typename Data_Counter_Type, typename Data_Rule_Type, typename Data_Rule_Dyn_Type>
 // inline const MapVec_type<> * Support<Array_Type,Data_Counter_Type,Data_Rule_Type, Data_Rule_Dyn_Type>::get_counts_ptr() const {