more on PD and PI

Author: conradhuebler

src/capabilities/confscan.cpp

@@ -56,12 +56,12 @@ static const json ConfScanJson = {
     { "preventreorder", false },
     { "scaleLoose", 1.5 },
     { "scaleTight", 0.1 },
-    { "scaleLooseEnergy", 1.2 },
-    { "scaleTightEnergy", 0.1 },
-    { "scaleLooseRotational", 1.2 },
-    { "scaleTightRotational", 0.1 },
-    { "scaleLooseRipser", 1.2 },
-    { "scaleTightRipser", 0.1 },
+    { "sLE", 1.2 },
+    { "sTE", 0.1 },
+    { "sLI", 1.2 },
+    { "sTI", 0.1 },
+    { "sLH", 1.2 },
+    { "sTH", 0.1 },
     { "skip", 0 },
     { "allxyz", false },
     { "update", false },
@@ -315,13 +315,13 @@ class ConfScan : public CurcumaMethod {
     std::map<double, int> m_ordered_list;
     std::vector<std::pair<std::string, Molecule*>> m_molecules;
     double m_rmsd_threshold = 1.0, m_nearly_missed = 0.8, m_energy_cutoff = -1, m_reference_last_energy = 0, m_target_last_energy = 0, m_lowest_energy = 1, m_current_energy = 0;
-    double m_scaleTightEnergy = 0.1, m_scaleLooseEnergy = 1.5;
-    double m_scaleTightRotational = 0.1, m_scaleLooseRotational = 1.5;
-    double m_scaleTightRipser = 0.1, m_scaleLooseRipser = 1.5;
+    double m_sTE = 0.1, m_sLE = 1.5;
+    double m_sTI = 0.1, m_sLI = 1.5;
+    double m_sTH = 0.1, m_sLH = 1.5;
 
     double m_reference_restored_energy = -1e10, m_target_restored_energy = -1e10;
-    double m_diff_rot_threshold_loose = 0.0, m_diff_ripser_threshold_loose = 0.0, m_diff_energy_threshold_loose = 0.0;
-    double m_diff_rot_threshold_tight = 0.0, m_diff_ripser_threshold_tight = 0.0, m_diff_energy_threshold_tight = 0.0;
+    double m_dLI = 0.0, m_dLH = 0.0, m_dLE = 0.0;
+    double m_dTI = 0.0, m_dTH = 0.0, m_dTE = 0.0;
 
     std::vector<Molecule*> m_result, m_rejected_structures, m_stored_structures, m_previously_accepted;
     std::vector<const Molecule*> m_threshold;
@@ -336,7 +336,7 @@ class ConfScan : public CurcumaMethod {
     std::string m_molalign = "molalign";
 
     double m_domolalign = -1;
-    double m_last_diff = 0.0, m_last_ripser = 0.0, m_last_dE = -1, m_dE = -1, m_damping = 0.8;
+    double m_lastDI = 0.0, m_lastDH = 0.0, m_lastdE = -1, m_dE = -1, m_damping = 0.8;
     int m_maxmol = 0;
     int m_maxrank = 10000;
     int m_maxParam = -1;

src/capabilities/confscan.h

@@ -29,7 +29,7 @@
 
 using json = nlohmann::json;
 
-PersistentDiagram::PersistentDiagram(const json config)
+PersistentDiagram::PersistentDiagram(const json& config)
 {
     json j = MergeJson(RipserJson, config);
     // std::cout << j << std::endl;
@@ -43,13 +43,15 @@ PersistentDiagram::PersistentDiagram(const json config)
     m_std_x = j["ripser_stdx"];
     m_std_y = j["ripser_stdy"];
     m_dimension = j["ripser_dimension"];
+    m_epsilon = j["ripser_epsilon"];
+
     m_threshold = std::numeric_limits<float>::max();
 }
 
 dpairs PersistentDiagram::generatePairs()
 {
     coefficient_t modulus = 2;
-
+    std::vector<float> distance = m_compressed_lower_distance_matrix;
     compressed_lower_distance_matrix dist(std::move(m_compressed_lower_distance_matrix));
 
     value_t enclosing_radius = std::numeric_limits<value_t>::infinity();
@@ -75,6 +77,68 @@ dpairs PersistentDiagram::generatePairs()
     return final;
 }
 
+triples PersistentDiagram::generateTriples()
+{
+    coefficient_t modulus = 2;
+    std::vector<float> distance = m_compressed_lower_distance_matrix;
+    compressed_lower_distance_matrix dist(std::move(m_compressed_lower_distance_matrix));
+
+    value_t enclosing_radius = std::numeric_limits<value_t>::infinity();
+    if (m_threshold == std::numeric_limits<value_t>::max()) {
+        for (size_t i = 0; i < dist.size(); ++i) {
+            value_t r_i = -std::numeric_limits<value_t>::infinity();
+            for (size_t j = 0; j < dist.size(); ++j)
+                r_i = std::max(r_i, dist(i, j));
+            enclosing_radius = std::min(enclosing_radius, r_i);
+        }
+    }
+
+    auto result = ripser<compressed_lower_distance_matrix>(std::move(dist), m_dimension, enclosing_radius,
+        m_ratio, modulus)
+                      .compute_barcodes();
+
+    triples final;
+    for (const auto& a : result)
+        if (a.first != 0) {
+            for (const auto& b : a.second) {
+                triple t;
+                t.start = b.first;
+                t.end = b.second;
+                final.push_back(t);
+            }
+        } else {
+            int counter = 0;
+            std::vector<int> indices;
+            for (const auto& b : a.second)
+
+            {
+                bool found = false;
+                for (int i = 0; i < distance.size(); ++i) {
+                    if (std::abs(distance[i] - b.second) < 1e-6 && std::find(indices.begin(), indices.end(), i) == indices.end()) {
+                        //    std::cout << i << " " << distance[i] << " " << m_en_scaling[i] << " " << counter++ << std::endl;
+                        indices.push_back(i);
+                        triple t;
+                        t.start = b.first;
+                        t.end = b.second;
+                        t.scaling = m_en_scaling[i] + m_epsilon;
+                        final.push_back(t);
+
+                        found = true;
+                        break;
+                    }
+                }
+                if (found == false) {
+                    triple t;
+                    t.start = b.first;
+                    t.end = b.second;
+                    final.push_back(t);
+                }
+            }
+        }
+
+    return final;
+}
+
 Eigen::MatrixXd PersistentDiagram::generateImage(const dpairs& pairs)
 {
     Eigen::MatrixXd matrix = Eigen::MatrixXd::Zero(m_bins, m_bins);
@@ -86,7 +150,6 @@ Eigen::MatrixXd PersistentDiagram::generateImage(const dpairs& pairs)
         for (int j = 0; j < matrix.rows(); ++j) // iteration over x values
         {
             double cur_x = (m_xmax - m_xmin) * i_bins * (j);
-
             for (const auto& pair : pairs) {
                 double x = pair.first;
                 double y = pair.second;
@@ -96,3 +159,25 @@ Eigen::MatrixXd PersistentDiagram::generateImage(const dpairs& pairs)
     }
     return matrix;
 }
+
+Eigen::MatrixXd PersistentDiagram::generateImage(const triples& pairs)
+{
+    Eigen::MatrixXd matrix = Eigen::MatrixXd::Zero(m_bins, m_bins);
+
+    double i_bins = 1 / double(m_bins);
+    for (int i = 0; i < matrix.cols(); ++i) // iteration over y values
+    {
+        double cur_y = m_ymax - (m_ymax - m_ymin) * i_bins * (i);
+        for (int j = 0; j < matrix.rows(); ++j) // iteration over x values
+        {
+            double cur_x = (m_xmax - m_xmin) * i_bins * (j);
+            for (const auto& triple : pairs) {
+                double x = triple.start;
+                double y = triple.end;
+                double std = triple.scaling;
+                matrix(i, j) += m_scaling * exp(-((cur_x - x) * (cur_x - x) * m_std_x / std)) * exp((-(cur_y - y) * (cur_y - y) * m_std_y / std));
+            }
+        }
+    }
+    return matrix;
+}

src/capabilities/persistentdiagram.cpp

@@ -26,9 +26,17 @@
 
 #include <Eigen/Dense>
 
+struct triple {
+    double start = 0;
+    double end = 0;
+    double scaling = 1;
+};
+
 typedef std::pair<double, double> dpair;
 
 typedef std::vector<dpair> dpairs;
+typedef std::vector<triple> triples;
+
 using json = nlohmann::json;
 
 static const json RipserJson = {
@@ -41,18 +49,24 @@ static const json RipserJson = {
     { "ripser_stdx", 10 },
     { "ripser_stdy", 10 },
     { "ripser_ratio", 1 },
-    { "ripser_dimension", 2 }
+    { "ripser_dimension", 2 },
+    { "ripser_epsilon", 0.4 }
 };
 
 class PersistentDiagram {
 public:
-    PersistentDiagram(const json config = RipserJson);
+    PersistentDiagram(const json& config = RipserJson);
 
     inline void setDistanceMatrix(const std::vector<float>& vector)
     {
         m_compressed_lower_distance_matrix = vector;
     }
 
+    inline void setENScaling(const std::vector<double>& en_scaling)
+    {
+        m_en_scaling = en_scaling;
+    }
+
     inline void setDimension(int dim)
     {
         m_dimension = dim;
@@ -69,6 +83,7 @@ class PersistentDiagram {
     }
 
     dpairs generatePairs();
+    triples generateTriples();
 
     inline void setXRange(double xmin, double xmax)
     {
@@ -83,6 +98,7 @@ class PersistentDiagram {
     }
 
     Eigen::MatrixXd generateImage(const dpairs& pairs);
+    Eigen::MatrixXd generateImage(const triples& pairs);
 
     void setScaling(double scaling) { m_scaling = scaling; }
     void setBins(int bins) { m_bins = bins; }
@@ -98,5 +114,7 @@ class PersistentDiagram {
     double m_scaling = 0.1;
     double m_std_x = 10;
     double m_std_y = 10;
+    double m_epsilon = 0.4;
     std::vector<float> m_compressed_lower_distance_matrix;
+    std::vector<double> m_en_scaling;
 };

src/capabilities/persistentdiagram.h

@@ -486,6 +486,95 @@ static const std::vector<double> CovalentRadius = {
     1.46 // At - taken from Bi
 };
 
+static const std::vector<double> PaulingEN = {
+    -1, // leading index
+    2.2, // H
+    0, // He
+    0.98, // Li
+    1.57, // Be
+    2.04, // B
+    2.55, // C
+    3.04, // N
+    3.44, // O
+    3.98, // F
+    0, // Ne
+    0.93, // Na
+    1.31, // Mg
+    1.61, // Al
+    1.9, // S
+    2.19, // P
+    2.58, // S
+    3.16, // Cl
+    0, // Ar
+    0.82, // K
+    1.0, // Ca
+    1.36, // Sc
+    1.54, // Ti
+    1.63, // V
+    1.66, // Cr
+    1.55, // Mn
+    1.83, // Fe
+    1.88, // Co
+    1.91, // Ni
+    1.9, // Cu
+    1.65, // Zn
+    1.81, // Ga
+    2.01, // Ge
+    2.18, // As
+    2.55, // Se
+    2.96, // Br
+    0, // Kr
+    0.82, // Rb
+    0.95, // Sr
+    1.22, // Y
+    1.33, // Zr
+    1.6, // Nb
+    2.16, // Mo
+    1.9, // Tc
+    2.2, // Ru
+    2.28, // Rh
+    2.20, // Pd
+    1.93, // Ag
+    1.69, // Cd
+    1.78, // In
+    1.96, // Sn
+    2.05, // Sb
+    2.1, // Te
+    2.66, // I
+    0, // Xe
+    0.79, // Cs
+    0.89, // Ba
+    1.1, // La
+    1.12, // Ce - taken from La
+    1.13, // Pr - taken from La
+    1.14, // Nd - taken from La
+    1.13, // Pm - taken from La
+    1.17, // Sm - taken from La
+    1.2, // Eu - taken from La
+    1.2, // Gd - taken from La
+    1.2, // Tb - taken from La
+    1.22, // Dy - taken from La
+    1.23, // Ho - taken from La
+    1.24, // Er - taken from La
+    1.25, // Tm - taken from La
+    1.1, // Yb - taken from La
+    1.27, // Lu
+    1.3, // Hf
+    1.5, // Ta
+    2.36, // W
+    1.9, // Re
+    2.2, // Os
+    2.2, // Ir
+    2.28, // Pt
+    2.54, // Au
+    2, // Hg
+    2.04, // Tl
+    2.33, // Pb
+    2.02, // Bi
+    2.0, // Po - taken from Bi
+    2.2 // At - taken from Bi
+};
+
 // The old ones
 /* for now, taken from here
  * https://de.wikipedia.org/wiki/Van-der-Waals-Radius */

src/core/elements.h

@@ -769,6 +769,17 @@ std::vector<float> Molecule::LowerDistanceVector() const
     return vector;
 }
 
+std::vector<double> Molecule::DeltaEN() const
+{
+    std::vector<double> vector;
+    for (int i = 0; i < AtomCount(); ++i) {
+        for (int j = 0; j < i; ++j) {
+            vector.push_back(std::abs(Elements::PaulingEN[m_atoms[i]] - Elements::PaulingEN[m_atoms[j]]));
+        }
+    }
+    return vector;
+}
+
 Position Molecule::Centroid(bool protons, int fragment) const
 {
     return GeometryTools::Centroid(getGeometryByFragment(fragment, protons));

src/core/molecule.cpp

@@ -101,6 +101,7 @@ class Molecule
 
     std::string LowerDistanceMatrix() const;
     std::vector<float> LowerDistanceVector() const;
+    std::vector<double> DeltaEN() const;
 
     bool setGeometry(const Geometry &geometry);
     bool setGeometryByFragment(const Geometry& geometry, int fragment, bool protons = true);