Merge pull request #197 from GeoStat-Framework/refactor_transform

Refactor: transform submodule; Field storage control; CondSRF optimization

CHANGELOG.md

@@ -2,6 +2,35 @@
 
 All notable changes to **GSTools** will be documented in this file.
 
+
+## [1.3.3] - Pure Pink - ?
+
+### Enhancements
+See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197)
+- `gstools.transform`:
+  - add keywords `field`, `store` and `process` keyword to all transformations to control storage and respect `normalizer`
+  - added `apply_function` transformation
+  - added `apply` as wrapper for all transformations
+  - added `transform` method to all `Field` (sub)classes as interface to `transform.apply`
+  - added checks for normal fields to work smoothly with recently added `normalizer` submodule
+- `Field`:
+  - allow naming fields when generating and control storage with `store` keyword
+  - all subclasses now have the `post_process` keyword (apply mean, normalizer, trend)
+  - added subscription to access fields by name (`Field["field"]`)
+  - added `set_pos` method to set position tuple
+  - allow reusing present `pos` tuple
+  - added `pos`, `mesh_type`, `field_names`, `field_shape`, `all_fields` properties
+- `CondSRF`:
+  - memory optimization by forwarding `pos` from underlying `krige` instance
+  - only recalculate kriging field if `pos` tuple changed (optimized ensemble generation)
+- performance improvement by using `np.asarray` instead of `np.array` where possible
+- updated examples to use new features
+- added incomplete lower gamma function `inc_gamma_low` (for TPLGaussian spectral density)
+
+### Bugfixes
+- `inc_gamma` was defined wrong for integer `s < 0`
+
+
 ## [1.3.2] - Pure Pink - 2021-07
 
 ### Bugfixes
@@ -280,7 +309,7 @@ All notable changes to **GSTools** will be documented in this file.
 First release of GSTools.
 
 
-[Unreleased]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.2...HEAD
+[1.3.3]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.2...HEAD
 [1.3.2]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.1...v1.3.2
 [1.3.1]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.0...v1.3.1
 [1.3.0]: https://github.com/GeoStat-Framework/gstools/compare/v1.2.1...v1.3.0

docs/source/random.rst

@@ -2,8 +2,6 @@ gstools.random
 ==============
 
 .. automodule:: gstools.random
-   :members:
-   :undoc-members:
 
 .. raw:: latex
 

docs/source/tools.rst

@@ -2,8 +2,6 @@ gstools.tools
 =============
 
 .. automodule:: gstools.tools
-   :members:
-   :undoc-members:
 
 .. raw:: latex
 

docs/source/transform.rst

@@ -2,8 +2,6 @@ gstools.transform
 =================
 
 .. automodule:: gstools.transform
-   :members:
-   :undoc-members:
 
 .. raw:: latex
 

examples/01_random_field/01_srf_ensemble.py

@@ -4,6 +4,8 @@
 
 Creating an ensemble of random fields would also be
 a great idea. Let's reuse most of the previous code.
+
+We will set the position tuple `pos` before generation to reuse it afterwards.
 """
 
 import numpy as np
@@ -14,26 +16,26 @@
 
 model = gs.Gaussian(dim=2, var=1, len_scale=10)
 srf = gs.SRF(model)
+srf.set_pos([x, y], "structured")
 
 ###############################################################################
 # This time, we did not provide a seed to :any:`SRF`, as the seeds will used
 # during the actual computation of the fields. We will create four ensemble
-# members, for better visualisation and save them in a list and in a first
+# members, for better visualisation, save them in to srf class and in a first
 # step, we will be using the loop counter as the seeds.
 
-
 ens_no = 4
-field = []
 for i in range(ens_no):
-    field.append(srf.structured([x, y], seed=i))
+    srf(seed=i, store=f"field{i}")
 
 ###############################################################################
-# Now let's have a look at the results:
+# Now let's have a look at the results. We can access the fields by name or
+# index:
 
 fig, ax = pt.subplots(2, 2, sharex=True, sharey=True)
 ax = ax.flatten()
 for i in range(ens_no):
-    ax[i].imshow(field[i].T, origin="lower")
+    ax[i].imshow(srf[i].T, origin="lower")
 pt.show()
 
 ###############################################################################
@@ -44,9 +46,8 @@
 # provides a seed generator :any:`MasterRNG`. The loop, in which the fields are
 # generated would then look like
 
-
 from gstools.random import MasterRNG
 
 seed = MasterRNG(20170519)
 for i in range(ens_no):
-    field.append(srf.structured([x, y], seed=seed()))
+    srf(seed=seed(), store=f"better_field{i}")

examples/06_conditioned_fields/00_condition_ensemble.py

@@ -26,15 +26,20 @@
 model = gs.Gaussian(dim=1, var=0.5, len_scale=1.5)
 krige = gs.krige.Ordinary(model, cond_pos, cond_val)
 cond_srf = gs.CondSRF(krige)
+cond_srf.set_pos(gridx)
 
 ###############################################################################
+# To generate the ensemble we will use a seed-generator.
+# We can specify individual names for each field by the keyword `store`:
 
-fields = []
+seed = gs.random.MasterRNG(20170519)
 for i in range(100):
-    fields.append(cond_srf(gridx, seed=i))
+    cond_srf(seed=seed(), store=f"f{i}")
     label = "Conditioned ensemble" if i == 0 else None
-    plt.plot(gridx, fields[i], color="k", alpha=0.1, label=label)
-plt.plot(gridx, cond_srf.krige(gridx, only_mean=True), label="estimated mean")
+    plt.plot(gridx, cond_srf[f"f{i}"], color="k", alpha=0.1, label=label)
+
+fields = [cond_srf[f"f{i}"] for i in range(100)]
+plt.plot(gridx, cond_srf.krige(only_mean=True), label="estimated mean")
 plt.plot(gridx, np.mean(fields, axis=0), linestyle=":", label="Ensemble mean")
 plt.plot(gridx, cond_srf.krige.field, linestyle="dashed", label="kriged field")
 plt.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")

examples/06_conditioned_fields/01_2D_condition_ensemble.py

@@ -20,34 +20,40 @@
 model = gs.Gaussian(dim=2, var=0.5, len_scale=5, anis=0.5, angles=-0.5)
 krige = gs.Krige(model, cond_pos=cond_pos, cond_val=cond_val)
 cond_srf = gs.CondSRF(krige)
+cond_srf.set_pos([x, y], "structured")
 
 ###############################################################################
-# We create a list containing the generated conditioned fields.
+# To generate the ensemble we will use a seed-generator.
+# By specifying ``store=[f"fld{i}", False, False]``, only the conditioned field
+# is stored with the specified name. The raw random field and the raw kriging
+# field is not stored. This way, we can access each conditioned field by index
+# ``cond_srf[i]``:
 
+seed = gs.random.MasterRNG(20170519)
 ens_no = 4
-field = []
 for i in range(ens_no):
-    field.append(cond_srf.structured([x, y], seed=i))
+    cond_srf(seed=seed(), store=[f"fld{i}", False, False])
 
 ###############################################################################
 # Now let's have a look at the pairwise differences between the generated
 # fields. We will see, that they coincide at the given conditions.
 
 fig, ax = plt.subplots(ens_no + 1, ens_no + 1, figsize=(8, 8))
 # plotting kwargs for scatter and image
-sc_kwargs = dict(c=cond_val, edgecolors="k", vmin=0, vmax=np.max(field))
-im_kwargs = dict(extent=2 * [0, 5], origin="lower", vmin=0, vmax=np.max(field))
+vmax = np.max(cond_srf.all_fields)
+sc_kw = dict(c=cond_val, edgecolors="k", vmin=0, vmax=vmax)
+im_kw = dict(extent=2 * [0, 5], origin="lower", vmin=0, vmax=vmax)
 for i in range(ens_no):
     # conditioned fields and conditions
-    ax[i + 1, 0].imshow(field[i].T, **im_kwargs)
-    ax[i + 1, 0].scatter(*cond_pos, **sc_kwargs)
-    ax[i + 1, 0].set_ylabel(f"Field {i+1}", fontsize=10)
-    ax[0, i + 1].imshow(field[i].T, **im_kwargs)
-    ax[0, i + 1].scatter(*cond_pos, **sc_kwargs)
-    ax[0, i + 1].set_title(f"Field {i+1}", fontsize=10)
+    ax[i + 1, 0].imshow(cond_srf[i].T, **im_kw)
+    ax[i + 1, 0].scatter(*cond_pos, **sc_kw)
+    ax[i + 1, 0].set_ylabel(f"Field {i}", fontsize=10)
+    ax[0, i + 1].imshow(cond_srf[i].T, **im_kw)
+    ax[0, i + 1].scatter(*cond_pos, **sc_kw)
+    ax[0, i + 1].set_title(f"Field {i}", fontsize=10)
     # absolute differences
     for j in range(ens_no):
-        ax[i + 1, j + 1].imshow(np.abs(field[i] - field[j]).T, **im_kwargs)
+        ax[i + 1, j + 1].imshow(np.abs(cond_srf[i] - cond_srf[j]).T, **im_kw)
 
 # beautify plots
 ax[0, 0].axis("off")
@@ -56,3 +62,9 @@
     a.set_xticks([]), a.set_yticks([])
 fig.subplots_adjust(wspace=0, hspace=0)
 fig.show()
+
+###############################################################################
+# To check if the generated fields are correct, we can have a look at their
+# names:
+
+print(cond_srf.field_names)

examples/07_transformations/00_log_normal.py

@@ -3,6 +3,8 @@
 -----------------
 
 Here we transform a field to a log-normal distribution:
+
+See :any:`transform.normal_to_lognormal`
 """
 import gstools as gs
 
@@ -11,5 +13,5 @@
 model = gs.Gaussian(dim=2, var=1, len_scale=10)
 srf = gs.SRF(model, seed=20170519)
 srf.structured([x, y])
-gs.transform.normal_to_lognormal(srf)
+srf.transform("normal_to_lognormal")  # also "lognormal" works
 srf.plot()

examples/07_transformations/01_binary.py

@@ -5,6 +5,8 @@
 Here we transform a field to a binary field with only two values.
 The dividing value is the mean by default and the upper and lower values
 are derived to preserve the variance.
+
+See :any:`transform.binary`
 """
 import gstools as gs
 
@@ -13,5 +15,5 @@
 model = gs.Gaussian(dim=2, var=1, len_scale=10)
 srf = gs.SRF(model, seed=20170519)
 srf.structured([x, y])
-gs.transform.binary(srf)
+srf.transform("binary")
 srf.plot()

examples/07_transformations/02_discrete.py

@@ -6,32 +6,38 @@
 If we do not give thresholds, the pairwise means of the given
 values are taken as thresholds.
 If thresholds are given, arbitrary values can be applied to the field.
+
+See :any:`transform.discrete`
 """
 import numpy as np
 import gstools as gs
 
-# structured field with a size of 100x100 and a grid-size of 0.5x0.5
+# Structured field with a size of 100x100 and a grid-size of 0.5x0.5
 x = y = np.arange(200) * 0.5
 model = gs.Gaussian(dim=2, var=1, len_scale=5)
 srf = gs.SRF(model, seed=20170519)
-
-# create 5 equidistanly spaced values, thresholds are the arithmetic means
 srf.structured([x, y])
-discrete_values = np.linspace(np.min(srf.field), np.max(srf.field), 5)
-gs.transform.discrete(srf, discrete_values)
-srf.plot()
 
-# calculate thresholds for equal shares
+###############################################################################
+# Create 5 equidistanly spaced values, thresholds are the arithmetic means
+
+values1 = np.linspace(np.min(srf.field), np.max(srf.field), 5)
+srf.transform("discrete", store="f1", values=values1)
+srf.plot("f1")
+
+###############################################################################
+# Calculate thresholds for equal shares
 # but apply different values to the separated classes
-discrete_values2 = [0, -1, 2, -3, 4]
-srf.structured([x, y])
-gs.transform.discrete(srf, discrete_values2, thresholds="equal")
-srf.plot()
 
-# user defined thresholds
+values2 = [0, -1, 2, -3, 4]
+srf.transform("discrete", store="f2", values=values2, thresholds="equal")
+srf.plot("f2")
+
+###############################################################################
+# Create user defined thresholds
+# and apply different values to the separated classes
+
+values3 = [0, 1, 10]
 thresholds = [-1, 1]
-# apply different values to the separated classes
-discrete_values3 = [0, 1, 10]
-srf.structured([x, y])
-gs.transform.discrete(srf, discrete_values3, thresholds=thresholds)
-srf.plot()
+srf.transform("discrete", store="f3", values=values3, thresholds=thresholds)
+srf.plot("f3")

examples/07_transformations/03_zinn_harvey.py

@@ -6,6 +6,8 @@
 `Zinn & Harvey (2003) <https://www.researchgate.net/publication/282442995_zinnharvey2003>`__.
 With this transformation, one could overcome the restriction that in ordinary
 Gaussian random fields the mean values are the ones being the most connected.
+
+See :any:`transform.zinnharvey`
 """
 import gstools as gs
 
@@ -14,5 +16,5 @@
 model = gs.Gaussian(dim=2, var=1, len_scale=10)
 srf = gs.SRF(model, seed=20170519)
 srf.structured([x, y])
-gs.transform.zinnharvey(srf, conn="high")
+srf.transform("zinnharvey", conn="high")
 srf.plot()

examples/07_transformations/04_bimodal.py

@@ -1,13 +1,15 @@
 """
-bimodal fields
----------------
+Bimodal fields
+--------------
 
 We provide two transformations to obtain bimodal distributions:
 
 * `arcsin <https://en.wikipedia.org/wiki/Arcsine_distribution>`__.
 * `uquad <https://en.wikipedia.org/wiki/U-quadratic_distribution>`__.
 
 Both transformations will preserve the mean and variance of the given field by default.
+
+See: :any:`transform.normal_to_arcsin` and :any:`transform.normal_to_uquad`
 """
 import gstools as gs
 
@@ -16,5 +18,5 @@
 model = gs.Gaussian(dim=2, var=1, len_scale=10)
 srf = gs.SRF(model, seed=20170519)
 field = srf.structured([x, y])
-gs.transform.normal_to_arcsin(srf)
+srf.transform("normal_to_arcsin")  # also "arcsin" works
 srf.plot()

examples/07_transformations/05_combinations.py

@@ -9,21 +9,32 @@
 Then we apply the Zinn & Harvey transformation to connect the low values.
 Afterwards the field is transformed to a binary field and last but not least,
 we transform it to log-values.
+
+We can select the desired field by its name and we can define an output name
+to store the field.
+
+If you don't specify `field` and `store` everything happens inplace.
 """
+# sphinx_gallery_thumbnail_number = 1
 import gstools as gs
 
 # structured field with a size of 100x100 and a grid-size of 1x1
 x = y = range(100)
 model = gs.Gaussian(dim=2, var=1, len_scale=10)
 srf = gs.SRF(model, mean=-9, seed=20170519)
 srf.structured([x, y])
-gs.transform.normal_force_moments(srf)
-gs.transform.zinnharvey(srf, conn="low")
-gs.transform.binary(srf)
-gs.transform.normal_to_lognormal(srf)
-srf.plot()
+srf.transform("force_moments", field="field", store="f_forced")
+srf.transform("zinnharvey", field="f_forced", store="f_zinnharvey", conn="low")
+srf.transform("binary", field="f_zinnharvey", store="f_binary")
+srf.transform("lognormal", field="f_binary", store="f_result")
+srf.plot(field="f_result")
 
 ###############################################################################
 # The resulting field could be interpreted as a transmissivity field, where
 # the values of low permeability are the ones being the most connected
 # and only two kinds of soil exist.
+#
+# All stored fields can be accessed and plotted by name:
+
+print("Max binary value:", srf.f_binary.max())
+srf.plot(field="f_zinnharvey")