minimal changes needed by FitSNAP PR #278

.gitignore

@@ -17,3 +17,8 @@ a.out
 __pycache__
 
 /ML-PACE/cnpy
+.DS_Store
+*.egg-info
+*.so
+*.pyc
+/python/build

ML-PACE/ace-evaluator/ace_array2dlm.h

@@ -105,7 +105,7 @@ class Array2DLM : public ContiguousArrayND<T> {
     void init(LS_TYPE lmax, string array_name = "Array2DLM") {
         if (is_proxy) {
             char s[1024];
-            sprintf(s, "Could not re-initialize proxy-array %s\n", this->array_name.c_str());
+            snprintf(s, sizeof(s), "Could not re-initialize proxy-array %s\n", this->array_name.c_str());
             throw logic_error(s);
         }
         this->lmax = lmax;

ML-PACE/ace-evaluator/ace_arraynd.h

@@ -161,7 +161,7 @@ template <class T> class Array1D: public ContiguousArrayND<T>  {
 
         if((i0>=dim[0])){ //(i0<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
+            snprintf(buf, sizeof(buf), "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
             throw std::out_of_range(buf);
         }
 
@@ -416,13 +416,13 @@ template <class T> class Array2D: public ContiguousArrayND<T>  {
 
         if((i0>=dim[0])){ //(i0<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
+            snprintf(buf, sizeof(buf), "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
             throw std::out_of_range(buf);
         }
 
         if((i1>=dim[1])){ //(i1<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1]-1);
+            snprintf(buf, sizeof(buf), "%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1]-1);
             throw std::out_of_range(buf);
         }
 
@@ -700,19 +700,19 @@ template <class T> class Array3D: public ContiguousArrayND<T>  {
 
         if((i0>=dim[0])){ //(i0<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
+            snprintf(buf, sizeof(buf), "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
             throw std::out_of_range(buf);
         }
 
         if((i1>=dim[1])){ //(i1<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1]-1);
+            snprintf(buf, sizeof(buf), "%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1]-1);
             throw std::out_of_range(buf);
         }
 
         if((i2>=dim[2])){ //(i2<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2]-1);
+            snprintf(buf, sizeof(buf), "%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2]-1);
             throw std::out_of_range(buf);
         }
 
@@ -1004,25 +1004,25 @@ template <class T> class Array4D: public ContiguousArrayND<T>  {
 
         if((i0>=dim[0])){ //(i0<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
+            snprintf(buf, sizeof(buf), "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
             throw std::out_of_range(buf);
         }
 
         if((i1>=dim[1])){ //(i1<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1]-1);
+            snprintf(buf, sizeof(buf), "%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1]-1);
             throw std::out_of_range(buf);
         }
 
         if((i2>=dim[2])){ //(i2<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2]-1);
+            snprintf(buf, sizeof(buf), "%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2]-1);
             throw std::out_of_range(buf);
         }
 
         if((i3>=dim[3])){ //(i3<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i3=%ld out of range (0, %ld)\n", array_name.c_str(), i3, dim[3]-1);
+            snprintf(buf, sizeof(buf), "%s: index i3=%ld out of range (0, %ld)\n", array_name.c_str(), i3, dim[3]-1);
             throw std::out_of_range(buf);
         }
 
@@ -1328,31 +1328,31 @@ template <class T> class Array5D: public ContiguousArrayND<T>  {
 
         if((i0>=dim[0])){ //(i0<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
+            snprintf(buf, sizeof(buf), "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
             throw std::out_of_range(buf);
         }
 
         if((i1>=dim[1])){ //(i1<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1]-1);
+            snprintf(buf, sizeof(buf), "%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1]-1);
             throw std::out_of_range(buf);
         }
 
         if((i2>=dim[2])){ //(i2<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2]-1);
+            snprintf(buf, sizeof(buf), "%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2]-1);
             throw std::out_of_range(buf);
         }
 
         if((i3>=dim[3])){ //(i3<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i3=%ld out of range (0, %ld)\n", array_name.c_str(), i3, dim[3]-1);
+            snprintf(buf, sizeof(buf), "%s: index i3=%ld out of range (0, %ld)\n", array_name.c_str(), i3, dim[3]-1);
             throw std::out_of_range(buf);
         }
 
         if((i4>=dim[4])){ //(i4<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i4=%ld out of range (0, %ld)\n", array_name.c_str(), i4, dim[4]-1);
+            snprintf(buf, sizeof(buf), "%s: index i4=%ld out of range (0, %ld)\n", array_name.c_str(), i4, dim[4]-1);
             throw std::out_of_range(buf);
         }
 
@@ -1672,37 +1672,37 @@ template <class T> class Array6D: public ContiguousArrayND<T>  {
 
         if((i0>=dim[0])){ //(i0<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
+            snprintf(buf, sizeof(buf), "%s: index i0=%ld out of range (0, %ld)\n", array_name.c_str(), i0, dim[0]-1);
             throw std::out_of_range(buf);
         }
 
         if((i1>=dim[1])){ //(i1<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1]-1);
+            snprintf(buf, sizeof(buf), "%s: index i1=%ld out of range (0, %ld)\n", array_name.c_str(), i1, dim[1]-1);
             throw std::out_of_range(buf);
         }
 
         if((i2>=dim[2])){ //(i2<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2]-1);
+            snprintf(buf, sizeof(buf), "%s: index i2=%ld out of range (0, %ld)\n", array_name.c_str(), i2, dim[2]-1);
             throw std::out_of_range(buf);
         }
 
         if((i3>=dim[3])){ //(i3<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i3=%ld out of range (0, %ld)\n", array_name.c_str(), i3, dim[3]-1);
+            snprintf(buf, sizeof(buf), "%s: index i3=%ld out of range (0, %ld)\n", array_name.c_str(), i3, dim[3]-1);
             throw std::out_of_range(buf);
         }
 
         if((i4>=dim[4])){ //(i4<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i4=%ld out of range (0, %ld)\n", array_name.c_str(), i4, dim[4]-1);
+            snprintf(buf, sizeof(buf), "%s: index i4=%ld out of range (0, %ld)\n", array_name.c_str(), i4, dim[4]-1);
             throw std::out_of_range(buf);
         }
 
         if((i5>=dim[5])){ //(i5<0)||
             char buf[1024];            
-            sprintf(buf,"%s: index i5=%ld out of range (0, %ld)\n", array_name.c_str(), i5, dim[5]-1);
+            snprintf(buf, sizeof(buf), "%s: index i5=%ld out of range (0, %ld)\n", array_name.c_str(), i5, dim[5]-1);
             throw std::out_of_range(buf);
         }
 

ML-PACE/ace-evaluator/ace_c_basis.cpp

@@ -1369,89 +1369,108 @@ void ACECTildeBasisSet::load_yaml(const string &yaml_file_name) {
         if (mu > nelements - 1)
             throw invalid_argument("yace::functions has species type key larger than nelements");
 
-        total_basis_size_rank1[mu] = 0;
-        total_basis_size[mu] = 0;
-
         auto ctildefunc_vec_yaml = p.second;
 
-        vector<ACECTildeBasisFunction> ctildefunc_vec;//TODO: read write is_half_ms_basis
+        // -------------------------------------------------------------------------
+        // REFACTOR START: Pre-calculate sizes to avoid vectors and stack objects
+        // -------------------------------------------------------------------------
+
+        int count_rank1 = 0;
+        int count_other = 0;
+
+        // Pass 1: Count how many basis functions of each type exist
+        for (const auto &ctildefunc_yaml : ctildefunc_vec_yaml) {
+            SHORT_INT_TYPE rank = ctildefunc_yaml["rank"].as<SHORT_INT_TYPE>();
+            if (rank == 1) count_rank1++;
+            else count_other++;
+        }
+
+        total_basis_size_rank1[mu] = count_rank1;
+        total_basis_size[mu] = count_other;
+
+        // Allocate the final arrays immediately
+        basis_rank1[mu] = new ACECTildeBasisFunction[count_rank1];
+        basis[mu] = new ACECTildeBasisFunction[count_other];
+
+        int func_ind_rank1 = 0;
+        int func_ind_other = 0;
+
+        // Pass 2: Fill data directly into the heap memory
         for (const auto &ctildefunc_yaml: ctildefunc_vec_yaml) {
+
+            SHORT_INT_TYPE rank = ctildefunc_yaml["rank"].as<SHORT_INT_TYPE>();
+
+            // Pointer to the actual slot in the array we want to fill
+            ACECTildeBasisFunction* target;
+
+            if (rank == 1) {
+                target = &basis_rank1[mu][func_ind_rank1];
+                func_ind_rank1++;
+            } else {
+                target = &basis[mu][func_ind_other];
+                func_ind_other++;
+            }
 
-            ACECTildeBasisFunction ctildefunc;
+            // --- Fill the target object directly (No Copying!) ---
 
-            ctildefunc.mu0 = ctildefunc_yaml["mu0"].as<SHORT_INT_TYPE>();
-            ctildefunc.rank = ctildefunc_yaml["rank"].as<SHORT_INT_TYPE>();
-            ctildefunc.ndensity = ctildefunc_yaml["ndensity"].as<SHORT_INT_TYPE>();
-            ctildefunc.num_ms_combs = ctildefunc_yaml["num_ms_combs"].as<SHORT_INT_TYPE>();
+            target->mu0 = ctildefunc_yaml["mu0"].as<SHORT_INT_TYPE>();
+            target->rank = rank;
+            target->ndensity = ctildefunc_yaml["ndensity"].as<SHORT_INT_TYPE>();
+            target->num_ms_combs = ctildefunc_yaml["num_ms_combs"].as<SHORT_INT_TYPE>();
 
+            // Aggregate rankmax
+            if (this->rankmax < target->rank)
+                this->rankmax = target->rank;
+
+            // Direct Allocation: MUS
             int_vec = ctildefunc_yaml["mus"].as<vector<int>>();
-            if (int_vec.size() != ctildefunc.rank)
+            if (int_vec.size() != target->rank)
                 throw invalid_argument("mus:: not sufficient number of values");
-            ctildefunc.mus = new SPECIES_TYPE[ctildefunc.rank];
-            for (int r = 0; r < ctildefunc.rank; r++)
-                ctildefunc.mus[r] = int_vec.at(r);
+            target->mus = new SPECIES_TYPE[target->rank];
+            for (int r = 0; r < target->rank; r++)
+                target->mus[r] = int_vec.at(r);
 
+            // Direct Allocation: NS
             int_vec = ctildefunc_yaml["ns"].as<vector<int>>();
-            if (int_vec.size() != ctildefunc.rank)
+            if (int_vec.size() != target->rank)
                 throw invalid_argument("ns:: not sufficient number of values");
-            ctildefunc.ns = new NS_TYPE[ctildefunc.rank];
-            for (int r = 0; r < ctildefunc.rank; r++)
-                ctildefunc.ns[r] = int_vec.at(r);
-
+            target->ns = new NS_TYPE[target->rank];
+            for (int r = 0; r < target->rank; r++)
+                target->ns[r] = int_vec.at(r);
 
+            // Direct Allocation: LS
             int_vec = ctildefunc_yaml["ls"].as<vector<int>>();
-            if (int_vec.size() != ctildefunc.rank)
+            if (int_vec.size() != target->rank)
                 throw invalid_argument("ls:: not sufficient number of values");
-            ctildefunc.ls = new LS_TYPE[ctildefunc.rank];
-            for (int r = 0; r < ctildefunc.rank; r++)
-                ctildefunc.ls[r] = int_vec.at(r);
+            target->ls = new LS_TYPE[target->rank];
+            for (int r = 0; r < target->rank; r++)
+                target->ls[r] = int_vec.at(r);
 
-            //this->ms_combs; //[num_ms_combs * rank]
+            // Direct Allocation: MS_COMBS
             int_vec = ctildefunc_yaml["ms_combs"].as<vector<int>>();
-            if (int_vec.size() != ctildefunc.rank * ctildefunc.num_ms_combs)
+            if (int_vec.size() != target->rank * target->num_ms_combs)
                 throw invalid_argument("ms_combs:: not sufficient number of values");
-            ctildefunc.ms_combs = new MS_TYPE[ctildefunc.rank * ctildefunc.num_ms_combs];
-            for (int r = 0; r < ctildefunc.rank * ctildefunc.num_ms_combs; r++)
-                ctildefunc.ms_combs[r] = int_vec.at(r);
+            target->ms_combs = new MS_TYPE[target->rank * target->num_ms_combs];
+            for (int r = 0; r < target->rank * target->num_ms_combs; r++)
+                target->ms_combs[r] = int_vec.at(r);
 
-
-            // this->ctildes; //[num_of_ms_combs * ndensity]
+            // Direct Allocation: CTILDES
             double_vec = ctildefunc_yaml["ctildes"].as<vector<DOUBLE_TYPE >>();
-            if (double_vec.size() != ctildefunc.ndensity * ctildefunc.num_ms_combs)
+            if (double_vec.size() != target->ndensity * target->num_ms_combs)
                 throw invalid_argument("ctildes:: not sufficient number of values");
-            ctildefunc.ctildes = new DOUBLE_TYPE[ctildefunc.ndensity * ctildefunc.num_ms_combs];
-            for (int r = 0; r < ctildefunc.ndensity * ctildefunc.num_ms_combs; r++)
-                ctildefunc.ctildes[r] = double_vec.at(r);
-
-            ctildefunc_vec.emplace_back(ctildefunc);
-
-            if (ctildefunc.rank == 1)
-                total_basis_size_rank1[mu]++;
-            else
-                total_basis_size[mu]++;
+            target->ctildes = new DOUBLE_TYPE[target->ndensity * target->num_ms_combs];
+            for (int r = 0; r < target->ndensity * target->num_ms_combs; r++)
+                target->ctildes[r] = double_vec.at(r);
 
-        } // end for over ctildefunc_vec_yaml
+            // Set ownership flag (just in case)
+            target->is_proxy = false;
 
-//        cout << "total_basis_size_rank1[mu]=" << total_basis_size_rank1[mu] << endl;
-//        cout << "total_basis_size[mu]=" << total_basis_size[mu] << endl;
-
-        basis_rank1[mu] = new ACECTildeBasisFunction[total_basis_size_rank1[mu]];
-        basis[mu] = new ACECTildeBasisFunction[total_basis_size[mu]];
-        int func_ind_rank1 = 0, func_ind = 0;
-        for (const ACECTildeBasisFunction ctildefunc: ctildefunc_vec) {
-            if (ctildefunc.rank == 1) {
-                basis_rank1[mu][func_ind_rank1] = ctildefunc;
-                func_ind_rank1++;
-            } else {
-                basis[mu][func_ind] = ctildefunc;
-                func_ind++;
-            }
-            //aggregate rankmax
-            if (this->rankmax < ctildefunc.rank)
-                this->rankmax = ctildefunc.rank;
-        }
-    }
+        } // end loop over yaml functions
+        // -------------------------------------------------------------------------
+        // REFACTOR END
+        // -------------------------------------------------------------------------
 
+    } // end loop over species (p)
 
     pack_flatten_basis();
 }

ML-PACE/ace-evaluator/ace_c_basisfunction.h

@@ -56,7 +56,7 @@ struct ACEAbstractBasisFunction {
     /**
     * flattened array of  computed combinations of (m1, m2, ..., m_rank)
     * which have non-zero general Clebsch-Gordan coefficient:
-    *  \f$ \mathbf{m}_1, \dots, \mathbf{m}_\mathrm{num ms combs}\f$ =
+    * \f$ \mathbf{m}_1, \dots, \mathbf{m}_\mathrm{num ms combs}\f$ =
     * \f$ (m_1, m_2, \dots, m_{rank})_1,  \dots, (m_1, m_2, \dots, m_{rank})_{\mathrm{num ms combs}} \f$,
     * size =  num_ms_combs * rank,
     * effective shape: [num_ms_combs][rank]
@@ -211,8 +211,10 @@ struct ACECTildeBasisFunction : public ACEAbstractBasisFunction {
     /**
      * Copy constructor, to fulfill the Rule of Three.
      * Always copy the dynamic memory, even if the source is a proxy object.
+     * FIX: Explicitly initialize ctildes to nullptr to avoid deleting garbage in _copy_from
      */
-    ACECTildeBasisFunction(const ACECTildeBasisFunction &other) {
+    ACECTildeBasisFunction(const ACECTildeBasisFunction &other) 
+        : ACEAbstractBasisFunction(), ctildes(nullptr) {
         _copy_from(other);
     }
 

ML-PACE/ace-evaluator/ace_radial.cpp

@@ -183,7 +183,7 @@ Function that computes Chebyshev polynomials of first and second kind
 void ACERadialFunctions::calcCheb(NS_TYPE n, DOUBLE_TYPE x) {
     if (n < 0) {
         char s[1024];
-        sprintf(s, "The order n of the polynomials should be positive %d\n", n);
+        snprintf(s, sizeof(s), "The order n of the polynomials should be positive %d\n", n);
         throw std::invalid_argument(s);
     }
     DOUBLE_TYPE twox = 2.0 * x;

ML-PACE/ace/ace_b_basis.cpp

@@ -512,7 +512,7 @@ void order_and_compress_b_basis_function(ACEBBasisFunction &func) {
                 s << "->>sorted XS-ns-ls-ms combinations: {\n";
                 char buf[1024];
                 for (const auto &tup: v) {
-                    sprintf(buf, "(%d, %d, %d, %d)\n", get<0>(tup), get<1>(tup), get<2>(tup), get<3>(tup));
+                    snprintf(buf, sizeof(buf), "(%d, %d, %d, %d)\n", get<0>(tup), get<1>(tup), get<2>(tup), get<3>(tup));
                     s << buf;
                 }
                 s << "}";

ML-PACE/ace/ace_b_basis.h

@@ -242,7 +242,7 @@ class ACEBBasisSet : public ACEFlattenBasisSet {
 
     void initialize_basis(BBasisConfiguration &basisSetup);
 
-    void _clean_contiguous_arrays();
+    void _clean_contiguous_arrays() override;
 
 
     vector<DOUBLE_TYPE> get_all_coeffs() const override;