creted methods for carrying out the comparison between the polynomial fits for...

creted methods for carrying out the comparison between the polynomial fits for different solvation energies from different solvents. this is using both the RMSD between the curves, and the Hausdorff distance.

solventmapcreator/polynomialanalysis/polynomialcomparison.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script for the comparison of the polynomial fits of different solvents.
+
+@author: mark
+"""
+
+import logging
+import numpy as np
+import numpy.polynomial.polynomial as poly
+import scipy.spatial.distance as sciptdist
+
+logging.basicConfig()
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.WARN)
+
+def calculate_hausdorff_distance_matrix(x_values, polynomial_dict_by_solvent_id):
+    """This calculates the hausdorff distance between all polynomial curves.
+    """
+    sorted_id_list = get_sorted_solvent_id_list(polynomial_dict_by_solvent_id)
+    hausdorff_matrix = np.zeros((len(sorted_id_list), len(sorted_id_list)))
+    for i in range(len(sorted_id_list)):
+        poly_values_i = polynomial_dict_by_solvent_id[sorted_id_list[i]]
+        for j in range(i):
+            poly_values_j = polynomial_dict_by_solvent_id[sorted_id_list[j]]
+            hausdorff_value = calculate_hausdorff_distance(x_values, poly_values_i, poly_values_j)
+            if i == j:
+                hausdorff_matrix[i][i] = hausdorff_value
+            else:
+                hausdorff_matrix[i][j] = hausdorff_value
+                hausdorff_matrix[j][i] = hausdorff_value
+    return hausdorff_matrix
+
+def calculate_rmsd_matrix(polynomial_dict_by_solvent_id):
+    """This calculates the rmsd between all curves, and outputs a matrix.
+    """
+    sorted_id_list = get_sorted_solvent_id_list(polynomial_dict_by_solvent_id)
+    rmsd_matrix = np.zeros((len(sorted_id_list), len(sorted_id_list)))
+    for i in range(len(sorted_id_list)):
+        poly_values_i = polynomial_dict_by_solvent_id[sorted_id_list[i]]
+        for j in range(i):
+            poly_values_j = polynomial_dict_by_solvent_id[sorted_id_list[j]]
+            rmsd_value = calculate_rmsd_between_curves(poly_values_i, poly_values_j)
+            #Noting that the matrix is symmetric
+            if i == j:
+                rmsd_matrix[i][i] = rmsd_value
+            else:
+                rmsd_matrix[i][j] = rmsd_value
+                rmsd_matrix[j][i] = rmsd_value
+    return rmsd_matrix
+
+def get_sorted_solvent_id_list(polynomial_dict_by_solvent_id):
+    """This returns a sorted list ofthe solvent IDs.
+    """
+    return sorted(polynomial_dict_by_solvent_id.keys())
+
+def calculate_hausdorff_distance(x_values, poly_values_1, poly_values_2):
+    """This calculates the Hausdorff distance between the 2 polynomial curves.
+    """
+    poly_values_1_with_x = create_array_for_hausdorff_analysis(x_values, poly_values_1)
+    poly_values_2_with_x = create_array_for_hausdorff_analysis(x_values, poly_values_2)
+    return max(sciptdist.directed_hausdorff(poly_values_1_with_x, poly_values_2_with_x)[0],
+               sciptdist.directed_hausdorff(poly_values_2_with_x, poly_values_1_with_x)[0])
+
+def calculate_rmsd_between_curves(poly_values_1, poly_values_2):
+    """This calculates the RMSD between the Assoicaiton energies along
+    the both curves. Note that both arrays must have the same dimensionality
+    for this analysis.
+    """
+    return np.sqrt(((poly_values_1 - poly_values_2)**2).mean())
+
+def create_array_for_hausdorff_analysis(x_values, polynomial_values):
+    """This creates a array of the correct dimensionality for hausdorff analysis.
+    """
+    array_size_2n = np.array([x_values, polynomial_values])
+    array_size_n2 = array_size_2n.transpose()
+    return array_size_n2
+
+def calculate_polynomial_values_all_solvents(x_values, polynomial_order,
+                                             polynomial_data_by_solvent_id):
+    """This calculates the polynomial vlaues for each solvent ID, and creates
+    a new dictionary, where the key is the solvent ID, and te value is the
+    polynomial values at the given x points.
+    """
+    poly_values_by_solvent_id = {}
+    for solvent_id, solvent_info_dict in polynomial_data_by_solvent_id.items():
+        poly_values = calculate_polynomial_values_from_dict(x_values, solvent_info_dict, polynomial_order)
+        poly_values_by_solvent_id[solvent_id] = poly_values
+    return poly_values_by_solvent_id
+
+def calculate_polynomial_values_from_dict(x_values, solvent_info_dict, polynomial_order):
+    """This calculates the polynomial values for the given info dict.
+    """
+    polynomial_coefficients = get_polynomial_coefficients(solvent_info_dict, polynomial_order)
+    return calculate_polynomial_values(x_values, polynomial_coefficients)
+
+def get_polynomial_coefficients(solvent_info_dict, polynomial_order):
+    """This returns the polynomial coefficients for the selected order.
+    """
+    return solvent_info_dict[polynomial_order]["coefficients"]
+
+
+def calculate_polynomial_values(x_values, polynomial_coefficients):
+    """This returns the values for the given polynomial function coefficients.
+    """
+    return poly.polyval(x_values, polynomial_coefficients)

solventmapcreator/test/polynomialanalysistest/polynomialcomparisontest.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script for tests of the polynomial comparison.
+
+@author: mark
+"""
+
+import logging
+import unittest
+import numpy as np
+import solventmapcreator.polynomialanalysis.polynomialcomparison as polynomialcomparison
+
+logging.basicConfig()
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.WARN)
+
+class PolynomialComparisonTestCase(unittest.TestCase):
+    """Test case for the polynomialcomparison script.
+    """
+    def setUp(self):
+        """Set up before tests
+        """
+        self.poly_data_by_solvent_id = {"poly_fit":{2:{"coefficients":np.array([0.0, 
+                                                                                             0.0, 1.0]),
+                                                                "RMSE":0.0254303,
+                                                                "order":1,
+                                                                "covar":0.0008143}},
+                                        "poly_fit2":{2:{"coefficients":np.array([1.0, 
+                                                                                          0.0, 1.0]),
+                                                                 "RMSE":0.0354303,
+                                                                 "order":1,
+                                                                 "covar":0.0008143}},
+                                        "poly_fit3":{2:{"coefficients":np.array([0.0, 
+                                                                                          0.0, 2.0]),
+                                                                 "RMSE":0.0354303,
+                                                                 "order":1,
+                                                                 "covar":0.0008143}}}
+        self.x_values = np.array([-2.0, -1.0, 0.0, 1.0, 2.0])
+        self.poly_values_dict = {"poly_fit":np.array([4.0, 1.0, 0.0, 1.0, 4.0]),
+                                 "poly_fit2":np.array([5.0, 2.0, 1.0, 2.0, 5.0]), 
+                                 "poly_fit3":np.array([8.0, 2.0, 0.0, 2.0, 8.0])}
+    def tearDown(self):
+        """Tear down after tests.
+        """
+        del self.poly_data_by_solvent_id
+        del self.x_values
+        del self.poly_values_dict
+    def test_calculate_hausdorff_distance_matrix(self):
+        """
+        """
+        expected_matrix = np.array([[0.0, 1.0, 4.0],
+                                    [1.0, 0.0, 3.0],
+                                    [4.0, 3.0, 0.0]])
+        actual_matrix = polynomialcomparison.calculate_hausdorff_distance_matrix(self.x_values,
+                                                                                 self.poly_values_dict)
+        np.testing.assert_array_almost_equal(expected_matrix, actual_matrix)
+    def test_calculate_rmsd_matrix(self):
+        """Test to see if expected RMSD matrix is produced.
+        """
+        expected_matrix = np.array([[0.0, 1.0, 2.607681],
+                                    [1.0, 0.0, 1.949359],
+                                    [2.607681, 1.949359, 0.0]])
+        actual_matrix = polynomialcomparison.calculate_rmsd_matrix(self.poly_values_dict)
+        np.testing.assert_array_almost_equal(expected_matrix, actual_matrix)
+    def test_get_sorted_solvent_id_list(self):
+        """Test to see if expected list is returned.
+        """
+        expected_list = ["poly_fit", "poly_fit2", "poly_fit3"]
+        actual_list = polynomialcomparison.get_sorted_solvent_id_list(self.poly_values_dict)
+        self.assertListEqual(expected_list, actual_list)
+    def test_calculate_hausdorff_distance(self):
+        """Test to see if expected Hausdorff distance is calculated.
+        """
+        expected_hausdorff_distance = 1.0
+        actual_hausdorff_distance = polynomialcomparison.calculate_hausdorff_distance(self.x_values, np.array([4.0, 1.0, 0.0, 1.0, 4.0]),np.array([5.0, 2.0, 1.0, 2.0, 5.0]))
+        self.assertAlmostEqual(expected_hausdorff_distance, actual_hausdorff_distance)
+    def test_calculate_rmsd_between_curves(self):
+        """test to see if expected RMSD is calculated.
+        """
+        expected_rmsd = 1.0
+        actual_rmsd = polynomialcomparison.calculate_rmsd_between_curves(np.array([4.0, 1.0, 0.0, 1.0, 4.0]),np.array([5.0, 2.0, 1.0, 2.0, 5.0]))
+        self.assertAlmostEqual(expected_rmsd, actual_rmsd)
+    def test_create_array_for_hausdorff_analysis(self):
+        """Test to see if expected array is produced.
+        """
+        expected_array = np.array([[-2.0, 4.0], [-1.0, 1.0], [0.0, 0.0], [1.0, 1.0], [2.0, 4.0]])
+        actual_array = polynomialcomparison.create_array_for_hausdorff_analysis(self.x_values, np.array([4.0, 1.0, 0.0, 1.0, 4.0]))
+        np.testing.assert_array_almost_equal(expected_array, actual_array)
+    def test_calculate_polynomial_values_all_solvents(self):
+        """Test to see if expected dictionary is produced.
+        """
+        expected_dict = {"poly_fit":np.array([4.0, 1.0, 0.0, 1.0, 4.0]),
+                         "poly_fit2":np.array([5.0, 2.0, 1.0, 2.0, 5.0]), 
+                         "poly_fit3":np.array([8.0, 2.0, 0.0, 2.0, 8.0])}
+        actual_dict = polynomialcomparison.calculate_polynomial_values_all_solvents(self.x_values, 2, self.poly_data_by_solvent_id)
+        self.assertListEqual(sorted(expected_dict.keys()), sorted(actual_dict.keys()))
+        for key in expected_dict.keys():
+            np.testing.assert_array_almost_equal(expected_dict[key], actual_dict[key])
+    def test_calculate_polynomial_values_from_dict(self):
+        """Test to see if expected array is produced.
+        """
+        expected_array = np.array([4.0, 1.0, 0.0, 1.0, 4.0])
+        actual_array = polynomialcomparison.calculate_polynomial_values_from_dict(self.x_values, self.poly_data_by_solvent_id["poly_fit"], 2)
+        np.testing.assert_array_almost_equal(expected_array, actual_array)
+    def test_get_polynomial_coefficients(self):
+        """Test to see if expected coefficients are returned.
+        """
+        expected_coefficients = np.array([0.0, 0.0, 1.0])
+        actual_coefficients = polynomialcomparison.get_polynomial_coefficients(self.poly_data_by_solvent_id["poly_fit"], 2)
+        np.testing.assert_array_almost_equal(expected_coefficients, actual_coefficients)
+    def test_calculate_polynomial_values(self):
+        """Test to see if expected values are returned.
+        """
+        expected_values = np.array([4.0, 1.0, 0.0, 1.0, 4.0])
+        poly_coefficients = self.poly_data_by_solvent_id["poly_fit"][2]["coefficients"]
+        actual_values = polynomialcomparison.calculate_polynomial_values(self.x_values, poly_coefficients)
+        np.testing.assert_array_almost_equal(expected_values, actual_values)
\ No newline at end of file

solventmapcreator/test/solvationmapcreatortests.py

@@ -13,6 +13,7 @@ from solventmapcreator.test.iotest.polynomialdatawritertest import PolynomialDat
 from solventmapcreator.test.iotest.polynomialdatareadertest import PolynomialDataReaderTestCase
 from solventmapcreator.test.solvationcalculationtest.solvationcalculatortest import SolvationCalculatorTestCase
 from solventmapcreator.test.polynomialanalysistest.polynomialdataanalysistest import PolynomialDataAnalysisTestCase
+from solventmapcreator.test.polynomialanalysistest.polynomialcomparisontest import PolynomialComparisonTestCase
 
 logging.basicConfig()
 LOGGER = logging.getLogger(__name__)
@@ -23,7 +24,8 @@ IO_TEST_CASES = [SolvationEnergyReaderTestCase, PolynomialDataWriterTestCase,
 
 SOLVATION_CALCULATION_TEST_CASES = [SolvationCalculatorTestCase]
 
-POLYNOMIAL_ANALYSIS_TEST_CASES = [PolynomialDataAnalysisTestCase]
+POLYNOMIAL_ANALYSIS_TEST_CASES = [PolynomialDataAnalysisTestCase,
+                                  PolynomialComparisonTestCase]
 
 def test_suite():
     """Function creates a test suite and then loads all the tests from the