Merge branch 'master' into doctest

.gitignore

@@ -23,6 +23,7 @@ var/
 *.egg-info/
 .installed.cfg
 *.egg
+benchmarks
 
 # PyInstaller
 #  Usually these files are written by a python script from a template

CHANGELOG.md

@@ -0,0 +1,11 @@
+
+# Changelog
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
+
+## [Unreleased]
+### Added
+	- volumetric_liquid_water is a new (and recommended) way to set the amount of water in a layer using make_snowpack.
+
+	- the 'emmodel' argument in make_model can now be a dict mapping a different emmodel for each sort of layer medium. Useful for snow + sea-ice for instance when the emmodels must be different for snow and ice.

README.md

@@ -5,7 +5,7 @@ Snow Microwave Radiative Transfer model
 [SMRT](https://www.smrt-model.science/) is a radiative transfer model to compute emission and backscatter from snowpack.
 
 Getting started is easy, follow the [instructions](https://www.smrt-model.science/getstarted.html) and explore the other repositories
-with examples in the ['smrt-model' github organization](https://github.com/smrt-model).
+with examples in the ['smrt-model' github organization](https://github.com/smrt-model) or read the detailed ['documentation'](https://smrt.readthedocs.io/en/latest/).
 
 If you want to try without installing anything on your computer, use free mybinder.org notenooks: [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/smrt-model/smrt/master?filepath=examples/iba_onelayer.ipynb)
 
@@ -18,7 +18,7 @@ WARRANTIES.
 
 DISCLAIMER: This version of SMRT is under peer review. Please use this software with caution, ask for assistance if needed, and let us know any feedback you may have.
 
-Copyright (c) 2016-2018 Ghislain Picard, Melody Sandells, Henning Löwe.
+Copyright (c) 2016-2022 Ghislain Picard, Melody Sandells, Henning Löwe.
 
 
 Other contributions
@@ -28,3 +28,5 @@ Other contributions
  - Ludovic Brucker
  - Marion Leduc-Leballeur
  - Mai Winstrup
+ - Carlo Marin
+

binder/requirements.txt

@@ -1,5 +1,4 @@
 numpy
-six
 pandas
 scipy
 xarray

smrt/core/error.py

@@ -2,13 +2,18 @@
 
 """
 
+import warnings
 
 
 class SMRTError(Exception):
     """Error raised by the model"""
     pass
 
 
-class SMRTWarning(Exception):
+class SMRTWarning(Warning):
     """Warning raised by the model"""
     pass
+
+
+def smrt_warn(message):
+    warnings.warn(message, category=SMRTWarning, stacklevel=2)

smrt/core/layer.py

@@ -38,7 +38,6 @@ class Layer(object):
     def __init__(self, thickness, microstructure_model=None,
                  temperature=FREEZING_POINT, permittivity_model=None, inclusion_shape=None,
                  **kwargs):
-
         """ Build a snow layer.
 
         :param thickness: thickness of snow layer in m
@@ -51,7 +50,9 @@ def __init__(self, thickness, microstructure_model=None,
 """
 
         self.thickness = thickness
-        self.microstructure_model = microstructure_model
+        # TODO Ghi: send a warning for non valid_args
+        if thickness == 0:
+            raise SMRTError("Layer with thickness = 0 (or even <~wavelength) is not recommended, part of the code does not support it.")
 
         self.temperature = temperature
         if temperature < 0:
@@ -60,25 +61,21 @@ def __init__(self, thickness, microstructure_model=None,
         self.permittivity_model = permittivity_model
         self.inclusion_shape = inclusion_shape
 
+        self.microstructure_model = microstructure_model
+
         # manage the microstructure parameters
         if microstructure_model is not None:
             valid_args = microstructure_model.valid_arguments()
-            params = {k: kwargs[k] for k in kwargs if k in valid_args}  # filter valid args
+            microstructure_params = {k: kwargs[k] for k in kwargs if k in valid_args}  # filter valid args
 
+            # make an instance of the micro-structure model
+            self.microstructure = microstructure_model(microstructure_params)
         else:
-            params = {}
-
-        # TODO Ghi: send a warning for non valid_args
-        if thickness == 0:
-            raise SMRTError("Layer with thickness = 0 (or even <~wavelength) is not recommended, part of the code does not support it.")
-
-        # make an instance of the micro-structure model
-        if microstructure_model is not None:
-            self.microstructure = microstructure_model(params)
+            microstructure_params = {}
 
         # other params
         for k in kwargs:
-            if k in params:
+            if k in microstructure_params:
                 continue
             setattr(self, k, kwargs[k])
 
@@ -126,7 +123,6 @@ def permittivity(self, i, frequency):
             # be called with the layer argument, but the initial permittivity_model function never heard about layers
             return self.permittivity_model[i](frequency, layer_to_inject=self)
 
-
         else:  # assume it is independent of the frequency.
             return self.permittivity_model[i]
 
@@ -172,18 +168,34 @@ def inverted_medium(self):
 
     def __setattr__(self, name, value):
 
-        if not hasattr(self, "read_only_attributes"):
-            super().__setattr__("read_only_attributes", set())
+        if hasattr(self, "read_only_attributes") and (name in self.read_only_attributes):
+            raise SMRTError(f"The attribute '{name}' is read-only, because setting its value requires recalculation."
+                            " In general, this is solved by using the update method.")
 
-        if name in self.read_only_attributes:
-            raise AttributeError("The attribute '%s' is read-only, setting its value would be subject to side effect. You need to create a new layer." % name)
         super().__setattr__(name, value)
 
+        # the callback system has been deactivated and replaced by the update method.
+        # See here for why the option __setattr__ has been chosen:
+        # https://stackoverflow.com/questions/17576009/python-class-property-use-setter-but-evade-getter
+        # if hasattr(self, "attributes_with_callback") and (name in self.attributes_with_callback):
+        #     # get the callback function
+        #     callback = self.attributes_with_callback[name]
+        #     # call the callback function
+        #     callback(self, name)
+
+    def update(self, **kwargs):
+        """update the attributes. This method is to be used when recalculation of the state of the object
+        is necessary. See for instance :py:class:`~smrt.inputs.make_medium.SnowLayer`.
+        """
+        for k, v in kwargs.items():
+            setattr(self, k, v)
+
 
 def get_microstructure_model(modulename, classname=None):
     """return the class corresponding to the microstructure_model defined in modulename.
 
-    This function import the correct module if possible and return the class. It is used internally and should not be needed for normal usage.
+    This function import the correct module if possible and return the class.
+    It is used internally and should not be needed for normal usage.
 
     :param modulename: name of the python module in smrt/microstructure_model
     """

smrt/core/lib.py

@@ -5,6 +5,7 @@
 import numpy as np
 import pandas as pd
 import scipy.sparse
+from smrt.core.optional_numba import numba
 
 from smrt.core.error import SMRTError
 
@@ -48,11 +49,11 @@ def check_argument_size(x, n, name=None):
 def is_sequence(x):
     # maybe not the smartest way...
     return (
-            isinstance(x, Sequence) or \
-            isinstance(x, np.ndarray) or \
-            isinstance(x, pd.DataFrame) or \
-            isinstance(x, pd.Series)
-            ) and not isinstance(x, str)
+        isinstance(x, Sequence) or
+        isinstance(x, np.ndarray) or
+        isinstance(x, pd.DataFrame) or
+        isinstance(x, pd.Series)
+    ) and not isinstance(x, str)
 
 
 def len_atleast_1d(x):
@@ -143,7 +144,7 @@ def __isub__(self, other):
             self.diag -= other.diag
             return self
         else:
-        
+
             raise NotImplementedError
 
     def __getitem__(self, key):
@@ -262,7 +263,7 @@ def compress(self, mode=None, auto_reduce_npol=False):
 
         elif self.mtype == "diagonal4":
             if self.values.shape[0] == 3 and auto_reduce_npol and mode == 0:
-                ## 3pol->2pol
+                # 3pol->2pol
                 mat = self.values[0:2, :]
             else:
                 mat = self.values
@@ -354,12 +355,17 @@ def isnull(m):
         m = m.diagonal()
 
     return (m is 0) or \
-            (getattr(m, "mtype", None) == "0") or \
-             (not np.any(m))
+        (getattr(m, "mtype", None) == "0") or \
+        (not np.any(m))
 
 
-def abs2(c):
-    return c.real**2 + c.imag**2
+if numba:
+    @numba.vectorize([numba.float64(numba.complex128), numba.float32(numba.complex64)], cache=True)
+    def abs2(x):
+        return x.real**2 + x.imag**2
+else:
+    def abs2(x):
+        return x.real**2 + x.imag**2
 
 
 def generic_ft_even_matrix(phase_function, m_max, nsamples=None):
@@ -450,4 +456,3 @@ def set_max_numerical_threads(nthreads):
     os.environ['OMP_NUM_THREADS'] = nthreads
     os.environ['VECLIB_MAXIMUM_THREADS'] = nthreads
     os.environ['NUMEXPR_NUM_THREADS'] = nthreads
-