Custoim weights included in point cloud

Author: SimonMolinsky

src/pyinterpolate/semivariogram/experimental/experimental_semivariogram.py

@@ -14,6 +14,10 @@
 from pyinterpolate.semivariogram.experimental.functions.semivariance import \
     semivariance_fn
 from pyinterpolate.semivariogram.lags.lags import get_lags
+from semivariogram.experimental.functions.directional import \
+    from_ellipse_weighted_cloud
+from semivariogram.experimental.functions.general import \
+    omnidirectional_weighted_semivariogram_cloud
 
 
 def calculate_semivariance(ds: Union[np.ndarray, VariogramPoints],
@@ -296,6 +300,77 @@ def directional_semivariance_cloud(points: np.ndarray,
     return output_semivariances
 
 
+def directional_weighted_semivariance_cloud(points: np.ndarray,
+                                            lags: Union[List, np.ndarray],
+                                            custom_weights: np.ndarray,
+                                            direction: float,
+                                            tolerance: float,
+                                            method: str) -> Dict:
+    """
+    Function calculates directional semivariances.
+
+    Parameters
+    ----------
+    points : numpy array
+        ``[x, y, value]``
+
+    lags : numpy array
+
+    custom_weights : numpy array
+
+    direction : float, optional
+        Direction of semivariogram, values from 0 to 360 degrees:
+
+        - 0 or 180: is E-W,
+        - 90 or 270 is N-S,
+        - 45 or 225 is NE-SW,
+        - 135 or 315 is NW-SE.
+
+    tolerance : float, optional
+        If ``tolerance`` is 0 then points must be placed at a single line with
+        the beginning in the origin of the coordinate system and the
+        direction given by y axis and direction parameter.
+        If ``tolerance`` is ``> 0`` then the bin is selected as an elliptical
+        area with major axis pointed in the same direction as the line for
+        ``0`` tolerance.
+
+        * The major axis size == ``step_size``.
+        * The minor axis size is ``tolerance * step_size``
+        * The baseline point is at a center of the ellipse.
+        * The ``tolerance == 1`` creates an omnidirectional semivariogram.
+
+    method : str
+        Neighbors selection in a given direction. Available methods:
+
+        * "triangle" or "t", default method where point neighbors are
+          selected from a triangular area,
+        * "ellipse" or "e", more accurate method but also slower.
+
+    Returns
+    -------
+    : Dict
+        ``{lag: [semivariances]}``
+    """
+    output_semivariances = dict()
+
+    if method == "e" or method == "ellipse":
+        output_semivariances = from_ellipse_weighted_cloud(
+            points,
+            lags,
+            custom_weights,
+            direction,
+            tolerance)
+    elif method == "t" or method == "triangle":
+        output_semivariances = from_triangle_cloud(
+            points=points,
+            lags=lags,
+            direction=direction,
+            tolerance=tolerance
+        )
+
+    return output_semivariances
+
+
 def directional_semivariance(points: np.ndarray,
                              lags: Union[List, np.ndarray],
                              direction: float,
@@ -352,7 +427,7 @@ def directional_semivariance(points: np.ndarray,
 
     if custom_weights is None:
         if (dir_neighbor_selection_method == "e" or
-            dir_neighbor_selection_method == "ellipse"):
+                dir_neighbor_selection_method == "ellipse"):
             output_semivariances = from_ellipse(
                 semivariance_fn,
                 points,
@@ -408,13 +483,16 @@ def omnidirectional_semivariance(points: np.ndarray,
     else:
 
         if custom_weights is not None:
-            # todo warning message - custom_weights not included
-            pass
-
-        sorted_semivariances = omnidirectional_semivariogram_cloud(
-            points=points,
-            lags=lags
-        )
+            sorted_semivariances = omnidirectional_weighted_semivariogram_cloud(
+                points=points,
+                lags=lags,
+                custom_weights=custom_weights
+            )
+        else:
+            sorted_semivariances = omnidirectional_semivariogram_cloud(
+                points=points,
+                lags=lags
+            )
 
     return sorted_semivariances
 
@@ -592,16 +670,22 @@ def point_cloud_semivariance(ds: Union[np.ndarray, VariogramPoints],
     if direction is not None and tolerance is not None:
 
         if custom_weights is not None:
-            # todo warning message - custom_weights not included
-            pass
-
-        experimental_semivariances = directional_semivariance_cloud(
-            ds.points,
-            lags,
-            direction,
-            tolerance,
-            dir_neighbors_selection_method
-        )
+            experimental_semivariances = directional_weighted_semivariance_cloud(
+                ds.points,
+                lags,
+                custom_weights,
+                direction,
+                tolerance,
+                dir_neighbors_selection_method
+            )
+        else:
+            experimental_semivariances = directional_semivariance_cloud(
+                ds.points,
+                lags,
+                direction,
+                tolerance,
+                dir_neighbors_selection_method
+            )
     else:
         experimental_semivariances = omnidirectional_semivariance(
             ds.points, lags, custom_weights, as_point_cloud=True

src/pyinterpolate/semivariogram/experimental/functions/directional.py

@@ -134,6 +134,65 @@ def from_ellipse_cloud(
     return point_cloud
 
 
+def from_ellipse_weighted_cloud(
+        points: np.ndarray,
+        lags: Union[List, np.ndarray],
+        custom_weights: np.ndarray,
+        direction: float,
+        tolerance: float):
+    """
+    Function calculates weighted point cloud semivariances from the elliptical
+    neighborhood.
+
+    Parameters
+    ----------
+    points : numpy array
+        Coordinates and their values: ``[x, y, value]``.
+
+    lags : numpy array
+        Bins array.
+
+    custom_weights : numpy array
+        Custom weights applied to points.
+
+    direction : float
+        Direction of semivariogram, values from 0 to 360 degrees:
+
+        - 0 or 180: is E-W,
+        - 90 or 270 is N-S,
+        - 45 or 225 is NE-SW,
+        - 135 or 315 is NW-SE.
+
+    tolerance : float, default=1
+        Value in range (0-1] to calculate semi-minor axis length of the
+        search area. If tolerance is close to 0 then points must be placed
+        in a single line with beginning in the origin of coordinate system
+        and direction given by y axis and direction parameter.
+            * The major axis length == step_size,
+            * The minor axis size == tolerance * step_size.
+            * Tolerance == 1 creates the omnidirectional semivariogram.
+
+    Returns
+    -------
+    output : Dict
+        ``{lag: [semivariances]}``
+    """
+    point_cloud = dict()
+    w_matrix = define_whitening_matrix(theta=direction,
+                                       minor_axis_size=tolerance)
+    for idx in range(len(lags)):
+        values = _get_weighted_values_ellipse_cloud(
+            points=points,
+            lags=lags,
+            custom_weights=custom_weights,
+            lag_idx=idx,
+            w_matrix=w_matrix
+        )
+        point_cloud[lags[idx]] = values
+
+    return point_cloud
+
+
 def from_triangle(
         fn: Callable,
         points: np.ndarray,
@@ -533,6 +592,79 @@ def _get_values_ellipse_cloud(
     return np.array([])
 
 
+def _get_weighted_values_ellipse_cloud(
+        points: np.ndarray,
+        lags: Union[List, np.ndarray],
+        custom_weights: np.ndarray,
+        lag_idx: int,
+        w_matrix: np.ndarray) -> np.ndarray:
+    """
+    Function selects semivariances or covariances per lag in the ellipse.
+
+    Parameters
+    ----------
+    points : numpy array
+        ``[x, y, value]``
+
+    lags : list or numpy array
+        The list of lags.
+
+    custom_weights : numpy array
+        Custom weights applied to points.
+
+    lag_idx : int
+        Current lag index.
+
+    w_matrix : numpy array
+        Matrix used for masking values in ellipse.
+
+    Returns
+    -------
+    values : numpy array
+        Semivariances for a given lag.
+    """
+    coords = points[:, :-1]
+
+    current_lag, previous_lag = get_current_and_previous_lag(
+        lag_idx=lag_idx, lags=lags
+    )
+
+    step_size = current_lag - previous_lag
+    semivariances = []
+
+    for idx, point in enumerate(points):
+
+        coordinates = point[:-1]
+
+        mask = select_points_within_ellipse(
+            ellipse_center=coordinates,
+            other_points=coords,
+            lag=current_lag,
+            step_size=step_size,
+            w_matrix=w_matrix
+        )
+
+        points_in_range = points[mask, -1]
+
+        # Extend list for calculations
+        if len(points_in_range) > 0:
+            # Calculate semivariances
+            w_vals = custom_weights[mask]
+            w_of_single_point = custom_weights[idx]
+            w_mult = w_of_single_point * w_vals
+            w_sum = w_of_single_point + w_vals
+            w_h = w_mult / w_sum
+            z_h = (points_in_range - point[-1]) ** 2
+            nominator_z = z_h - np.average(points_in_range,
+                                           weights=w_vals)
+            nominator_z = np.sum(nominator_z * w_h)
+            gamma = nominator_z / np.sum(w_h)
+            average_semivariance = 0.5 * gamma
+            semivariances.append(average_semivariance)
+
+    return np.array(semivariances)
+
+
 def _get_values_triangle(
         fn: Callable,
         points: np.ndarray,