Merge pull request #68 from virocon-organization/direct-sampling

Direct sampling contour

README.rst

@@ -28,11 +28,12 @@ The following methods are implemented in viroconcom:
 - Computing an environmental contour using either the
 
   - inverse first-order reliability method (IFORM),
-  - inverse second-order reliability method (ISORM) or the
+  - inverse second-order reliability method (ISORM),
+  - the direct sampling contour method or the
   - highest density contour method.
 
 
-ViroCon is written in Python 3.6.4. The software is seperated in two
+ViroCon is written in Python 3.7.4. The software is seperated in two
 main packages, viroconweb and viroconcom. This is the repository of
 viroconcom, which is the numerical core. It handles the statistical
 computations. viroconweb builds on viroconcom and is a web-based

examples/compare_contour_methods.py

@@ -1,7 +1,9 @@
 from viroconcom.params import ConstantParam, FunctionParam
 from viroconcom.distributions import WeibullDistribution, LognormalDistribution, \
     MultivariateDistribution
-from viroconcom.contours import IFormContour, ISormContour, HighestDensityContour
+from viroconcom.contours import IFormContour, ISormContour, \
+    HighestDensityContour, DirectSamplingContour, sort_points_to_form_continous_line
+from viroconcom.plot import plot_contour
 import matplotlib.pyplot as plt
 
 # Define the multivariate distribution given in the paper by Vanem and
@@ -19,43 +21,33 @@
 dependencies = [dep0, dep1]
 mul_dist = MultivariateDistribution(distributions, dependencies)
 
-# Compute an IFORM, an ISORM and a highest density contour.
+# Compute an IFORM, an ISORM, a direct sampling and a highest density contour.
 return_period = 50 # In years
 sea_state_duration = 6 # In hours
 iform_contour = IFormContour(mul_dist, return_period, sea_state_duration, 100)
 isorm_contour = ISormContour(mul_dist, return_period, sea_state_duration, 100)
+ds_contour = DirectSamplingContour(
+    mul_dist, return_period, sea_state_duration, 5000000)
 limits = [(0, 20), (0, 20)] # Limits of the computational domain
 deltas = [0.005, 0.005] # Dimensions of the grid cells
 hdens_contour = HighestDensityContour(
     mul_dist, return_period, sea_state_duration, limits, deltas)
+hdc_coordinates = sort_points_to_form_continous_line(
+    hdens_contour.coordinates[0][0], hdens_contour.coordinates[0][1])
 
-# Plot the three contours.
-plt.scatter(hdens_contour.coordinates[0][0], hdens_contour.coordinates[0][1],
-            label='highest density contour')
-plt.scatter(iform_contour.coordinates[0][0], iform_contour.coordinates[0][1],
-            label='IFORM contour')
-plt.scatter(isorm_contour.coordinates[0][0], isorm_contour.coordinates[0][1],
-            label='ISORM contour')
-plt.xlabel('significant wave height [m]')
-plt.ylabel('zero-upcrossing period [s]')
+# Plot the four contours (a similar plot was presented in the paper by
+# Haselsteiner et al. (2017; 10.1016/j.coastaleng.2017.03.002), Fig. 8 c).
+fig = plt.figure(figsize=(5, 5), dpi=150)
+ax = fig.add_subplot(111)
+plot_contour(iform_contour.coordinates[0][1], iform_contour.coordinates[0][0],
+             ax, line_style='b-', contour_label='IFORM contour')
+plot_contour(isorm_contour.coordinates[0][1], isorm_contour.coordinates[0][0],
+             ax, line_style='r-', contour_label='ISORM contour')
+plot_contour(ds_contour.coordinates[1], ds_contour.coordinates[0], ax,
+             line_style='g-', contour_label='Direct sampling contour')
+plot_contour(hdc_coordinates[1], hdc_coordinates[0], ax,
+             line_style='k-', contour_label='Highest density contour')
+plt.xlabel('Zero-up-crossing period (s)')
+plt.ylabel('Significant wave height (m)')
 plt.legend()
 plt.show()
-
-
-# To evalute viroconcom, we compare the maximum values of the contour with the
-# results from Haselsteiner et al. (2017; doi: 10.1016/j.coastaleng.2017.03.002).
-#
-# Results in Haselsteiner et al. (2017) with alpha = 1.37 * 10^-5 are:
-# Method           maximum Hs (m)     maximum Tz (s)
-# IFORM contour    15.23              13.96
-# ISORM contour    ca. 17.4           ca. 14.9
-# HDC              16.79              14.64
-print('Maximum values for the IFORM contour: ' +
-      '{:.2f}'.format(max(iform_contour.coordinates[0][0])) + ' m, '
-      + '{:.2f}'.format(max(iform_contour.coordinates[0][1])) + ' s')
-print('Maximum values for the ISORM contour: ' +
-      '{:.2f}'.format(max(isorm_contour.coordinates[0][0])) + ' m, '
-      + '{:.2f}'.format(max(isorm_contour.coordinates[0][1])) + ' s')
-print('Maximum values for the HDC: ' +
-      '{:.2f}'.format(max(hdens_contour.coordinates[0][0])) + ' m, '
-      + '{:.2f}'.format(max(hdens_contour.coordinates[0][1])) + ' s')

examples/direct_sampling_contour.py

@@ -0,0 +1,37 @@
+from viroconcom.params import ConstantParam, FunctionParam
+from viroconcom.distributions import WeibullDistribution, LognormalDistribution, \
+    MultivariateDistribution
+from viroconcom.contours import DirectSamplingContour
+import matplotlib.pyplot as plt
+
+# Define a Weibull distribution representing significant wave height.
+shape = ConstantParam(1.471)
+loc = ConstantParam(0.8888)
+scale = ConstantParam(2.776)
+dist0 = WeibullDistribution(shape, loc, scale)
+dep0 = (None, None, None)  # All three parameters are independent.
+
+# Define a lognormal distribution representing spectral peak period.
+my_sigma = FunctionParam("exp3", 0.0400, 0.1748, -0.2243)
+my_mu = FunctionParam("power3", 0.1, 1.489, 0.1901)
+dist1 = LognormalDistribution(sigma=my_sigma, mu=my_mu)
+dep1 = (0, None, 0)  # Parameter one and three depend on dist0.
+
+# Create a multivariate distribution by bundling the two distributions.
+distributions = [dist0, dist1]
+dependencies = [dep0, dep1]
+mul_dist = MultivariateDistribution(distributions, dependencies)
+
+# Compute a 1-year direct sampling contour based on drawing 10^6 observations.
+contour = DirectSamplingContour(mul_var_dist=mul_dist, return_period=1,
+                                state_duration=6, n=1000000, deg_step=6)
+
+# Plot the contour and the sample.
+plt.scatter(contour.sample[1], contour.sample[0], marker='.')
+plt.plot(contour.coordinates[1], contour.coordinates[0], color='red')
+plt.plot([contour.coordinates[1][-1], contour.coordinates[1][0]],
+         [contour.coordinates[0][-1], contour.coordinates[0][0]], color='red')
+plt.title('1-year direct sampling contour')
+plt.ylabel('Significant wave height (m)')
+plt.xlabel('Zero-up-crossing period (s)')
+plt.show()

tests/test_contours.py

@@ -17,7 +17,8 @@
 from viroconcom.distributions import (
     WeibullDistribution, ExponentiatedWeibullDistribution, LognormalDistribution,
     NormalDistribution, MultivariateDistribution)
-from viroconcom.contours import IFormContour, ISormContour, HighestDensityContour
+from viroconcom.contours import IFormContour, ISormContour, \
+    HighestDensityContour, DirectSamplingContour
 
 
 _here = os.path.dirname(__file__)
@@ -711,5 +712,76 @@ def test_setup_HDC_limits_Tuple_length(self):
     #     ax2.title.set_text('Sorted')
     #     #plt.show()
 
+
+class DirectSamplingTest(unittest.TestCase):
+    def _setup(self,
+               dep1=(None, None, None),
+               dep2=(0, None, 0),
+               par1=(ConstantParam(1.471), ConstantParam(0.8888),
+                     ConstantParam(2.776)),
+               par2=(FunctionParam('exp3', 0.0400, 0.1748, -0.2243), None,
+                     FunctionParam('power3', 0.1, 1.489, 0.1901))
+               ):
+        """
+        Create a 1-year contour (same as in DOI: 10.1016/j.oceaneng.2012.12.034).
+        """
+        return_period = 1 # years
+        state_duration = 6 # hours
+
+        # Define dependency tuple.
+        self.dep1 = dep1
+        self.dep2 = dep2
+
+        # Define parameters.
+        self.par1 = par1
+        self.par2 = par2
+
+        # Create distributions.
+        dist1 = WeibullDistribution(*par1)
+        dist2 = LognormalDistribution(sigma=par2[0], mu=par2[2])
+
+        distributions = [dist1, dist2]
+        dependencies = [dep1, dep2]
+
+        mul_dist = MultivariateDistribution(distributions, dependencies)
+
+        # Calculate contour.
+        ds_contour = DirectSamplingContour(
+            mul_dist, return_period, state_duration, 500000, 6)
+        return ds_contour
+
+    def test_direct_sampling_contour(self):
+        """
+        Computes a direct sampling contour and compares it with results
+        presented in the literature (DOI: 10.1016/j.oceaneng.2012.12.034).
+        """
+        contour = self._setup()
+        ref_contour_hs_1 = \
+            [9.99, 10.65, 10.99, 11.25, 11.25, 11.41, 11.42, 11.46, 11.48,
+             11.54, 11.57, 11.56, 11.58, 11.59, 11.59, 11.60, 11.60, 11.59,
+             11.59, 11.56, 11.53, 11.46, 11.26, 10.88, 7.44, 2.05]
+        ref_contour_tz_1 = \
+            [12.34, 12.34, 12.31, 12.25, 12.25, 12.18, 12.17, 12.15, 12.13,
+             12.06, 12.02, 12.03, 12.00, 11.96, 11.95, 11.86, 11.84, 11.77,
+             11.76, 11.67, 11.60, 11.47, 11.20, 10.77, 7.68, 3.76]
+        np.testing.assert_allclose(contour.coordinates[0][0:26], ref_contour_hs_1, atol=0.5)
+        np.testing.assert_allclose(contour.coordinates[1][0:26], ref_contour_tz_1, atol=0.5)
+
+    def test_3d_ds_contour(self):
+        """
+        Tests whether calculating a 3D DS contour raises the right error.
+        """
+        par = (ConstantParam(1.471), ConstantParam(0.8888),
+                ConstantParam(2.776))
+        dist = WeibullDistribution(*par)
+        dep = (None, None, None)
+        distributions = (dist, dist, dist)
+        dependencies = (dep, dep, dep)
+        mul_dist = MultivariateDistribution(distributions, dependencies)
+
+        # Calculate a 3D contour, which should raise an error.
+        self.assertRaises(NotImplementedError,
+                          DirectSamplingContour, mul_dist, 1, 3, 500000, 6)
+
 if __name__ == '__main__':
-    unittest.main()
+    unittest.main()

tests/test_distributions.py

@@ -4,13 +4,12 @@
 import numpy as np
 import scipy.stats as sts
 
-from .context import viroconcom
-
 from viroconcom.params import ConstantParam, FunctionParam
 from viroconcom.distributions import (
     WeibullDistribution, ExponentiatedWeibullDistribution,
     LognormalDistribution, NormalDistribution,
     MultivariateDistribution)
+from viroconcom.fitting import Fit
 
 class MultivariateDistributionTest(unittest.TestCase):
     """
@@ -97,8 +96,6 @@ def test_add_distribution_dependencies_value(self):
         with self.assertRaises(ValueError):
             MultivariateDistribution(self.distributions, dependencies)
 
-
-
     def test_add_distribution_not_iterable(self):
         """
         Tests the function when both distributions and dependencies
@@ -110,6 +107,61 @@ def test_add_distribution_not_iterable(self):
         with self.assertRaises(ValueError):
             MultivariateDistribution(distributions, dependencies)
 
+    def test_multivariate_draw_sample(self):
+        """
+        Tests the draw_sample() function of MulvariateDistribution.
+        """
+
+        # Create a MultivariateDistribution, the joint distribution for Hs-Tz
+        # presented in Vanem et al. (2012).
+        # Start with a Weibull distribution for wave height, Hs.
+        shape = ConstantParam(1.471)
+        loc = ConstantParam(0.8888)
+        scale = ConstantParam(2.776)
+        par0 = (shape, loc, scale)
+        dist0 = WeibullDistribution(*par0)
+
+        # Conditonal lognormal distribution for period, Tz.
+        dist1_shape = FunctionParam('exp3', 0.0400, 0.1748, -0.2243)
+        dist1_loc = None
+        dist1_scale = FunctionParam('power3', 0.1, 1.489, 0.1901)
+        par1 = (dist1_shape, dist1_loc, dist1_scale)
+        dist1 = LognormalDistribution(*par1)
+
+        distributions = [dist0, dist1]
+
+        dep0 = (None, None, None)
+        dep1 = (0, None, 0)
+        dependencies = [dep0, dep1]
+
+        mul_var_dist = MultivariateDistribution(distributions, dependencies)
+
+        n=10000
+        sample = mul_var_dist.draw_sample(n)
+        self.assertEqual(n, sample[0].size)
+        self.assertEqual(n, sample[1].size)
+
+        # Fit the sample to the correct model structure and compare the
+        # estimated parameters with the true parameters.
+        dist_description_0 = {'name': 'Weibull',
+                              'dependency': (None, None, None),
+                              'width_of_intervals': 0.5}
+        dist_description_1 = {'name': 'Lognormal',
+                              'dependency': (0, None, 0),
+                              'functions': ('exp3', None, 'power3')}
+        fit = Fit(sample, [dist_description_0, dist_description_1])
+        fitted_dist0 = fit.mul_var_dist.distributions[0]
+        fitted_dist1 = fit.mul_var_dist.distributions[1]
+        self.assertAlmostEqual(fitted_dist0.shape(0), shape(0), delta=0.1)
+        self.assertAlmostEqual(fitted_dist0.loc(0), loc(0), delta=0.1)
+        self.assertAlmostEqual(fitted_dist0.scale(0), scale(0), delta=0.1)
+        self.assertAlmostEqual(fitted_dist1.shape.a, 0.04, delta=0.1)
+        self.assertAlmostEqual(fitted_dist1.shape.b, 0.1748, delta=0.1)
+        self.assertAlmostEqual(fitted_dist1.shape.c, -0.2243, delta=0.15)
+        self.assertAlmostEqual(fitted_dist1.scale(1), dist1_scale(1), delta=0.1)
+        self.assertAlmostEqual(fitted_dist1.scale(2), dist1_scale(2), delta=0.1)
+
+
     def test_latex_representation(self):
         """
         Tests if the latex representation is correct.