Repository: earth-chris/elapid

Issue:
Create a niche envelope model class
This is a presence-only empirical model (no absence/pseudoabsence required). It would find the min/max range of all covariates in the training data, and use boolean logic to flag all predictions to be inside (`True`) or outside (`False`) the range of suitability for ech covariate. Rows that match all conditions (the `union` of all logic checks, or `np.all()`) get flagged as `1` (inside the envelope of suitability).

Some features to include:

- allow setting a percentile instead of the min/max to filter outliers.
- allow the user to pass a set of `mins` or `maxs`, which would kinda defeat the purpose of the 'model' but allow the rapid application of model logic.
- compute the `intersection` of suitability instead of the `union`, returning a continuous 0.0-1.0 range of suitability based on the product of all covariates that are in range.

Discussion:

Patch:
diff --git a/elapid/__version__.py b/elapid/__version__.py
index 7a13148..02ed9ef 100644
--- a/elapid/__version__.py
+++ b/elapid/__version__.py
@@ -1,1 +1,1 @@
-"0.3.2"
+"0.3.3"
diff --git a/elapid/config.py b/elapid/config.py
index 16f2178..04e36f9 100644
--- a/elapid/config.py
+++ b/elapid/config.py
@@ -1,5 +1,5 @@
 """SDM model configuration parameters."""
-from typing import Union
+from typing import Tuple, Union
 
 
 class MaxentConfig:
@@ -49,3 +49,7 @@ class RegularizationConfig:
     hinge: list = [[0, 1], [0.5, 0.5]]
     threshold: list = [[0, 100], [2, 1]]
     categorical: list = [[0, 10, 17], [0.65, 0.5, 0.25]]
+
+
+class NicheEnvelopeConfig:
+    percentile_range: Tuple[float, float] = [2.5, 97.5]
diff --git a/elapid/features.py b/elapid/features.py
index 1a13281..cb32e79 100644
--- a/elapid/features.py
+++ b/elapid/features.py
@@ -346,6 +346,8 @@ class MaxentFeatureTransformer(BaseEstimator):
     n_threshold_features_: int = None
     categorical_: list = None
     continuous_: list = None
+    categorical_pd_: list = None
+    continuous_pd_: list = None
     labels_: list = None
     estimators_: dict = {
         "linear": None,
diff --git a/elapid/models.py b/elapid/models.py
index 9f4e213..fa97606 100644
--- a/elapid/models.py
+++ b/elapid/models.py
@@ -1,5 +1,5 @@
 """Classes for training species distribution models."""
-from typing import List, Union
+from typing import List, Tuple, Union
 
 import numpy as np
 import pandas as pd
@@ -7,13 +7,13 @@ from glmnet.logistic import LogitNet
 from sklearn.base import BaseEstimator
 
 from elapid import features as _features
-from elapid.config import MaxentConfig
+from elapid.config import MaxentConfig, NicheEnvelopeConfig
 from elapid.types import ArrayLike, Number, validate_feature_types
-from elapid.utils import n_cpus
+from elapid.utils import make_band_labels, n_cpus
 
 
 class MaxentModel(BaseEstimator):
-    """Creates a model estimator for Maxent-style species distribution models."""
+    """Model estimator for Maxent-style species distribution models."""
 
     feature_types_: list = MaxentConfig.feature_types
     tau_: float = MaxentConfig.tau
@@ -119,9 +119,6 @@ class MaxentModel(BaseEstimator):
             feature_labels: list of length n_features, with labels identifying each column's feature type
                 with options ["linear", "quadratic", "product", "threshold", "hinge", "categorical"]
                 must be set if `is_features=True`
-
-        Returns:
-            None: Updates the model object
         """
         # fir the feature transformer
         if is_features:
@@ -227,7 +224,7 @@ class MaxentModel(BaseEstimator):
         y: ArrayLike,
         categorical: list = None,
         labels: list = None,
-        transform: str = "logistic",
+        transform: str = "cloglog",
         is_features: bool = False,
         feature_labels: list = None,
     ) -> ArrayLike:
@@ -239,7 +236,7 @@ class MaxentModel(BaseEstimator):
             categorical: column indices indicating which columns are categorical
             labels: Covariate labels. Ignored if x is a pandas DataFrame
             transform: maxent model transformation type. select from
-                ["raw", "exponential", "logistic", "cloglog"].
+                ["raw", "logistic", "cloglog"].
             is_features: specify that x data has already been transformed from covariates to features
             feature_labels: list of length n_features, with labels identifying each column's feature type
                 with options ["linear", "quadratic", "product", "threshold", "hinge", "categorical"]
@@ -284,6 +281,178 @@ class MaxentModel(BaseEstimator):
         self.initialized_ = True
 
 
+class NicheEnvelopeModel(BaseEstimator):
+    """Model estimator for niche envelope-style models."""
+
+    percentile_range_: Tuple[float, float] = None
+    feature_mins_: np.ndarray = None
+    feature_maxs_: np.ndarray = None
+    categorical_estimator: BaseEstimator = None
+    categorical_: list = None
+    continuous_: list = None
+    categorical_pd_: list = None
+    continuous_pd_: list = None
+    in_categorical_: np.ndarray = None
+
+    def __init__(self, percentile_range: Tuple[float, float] = NicheEnvelopeConfig.percentile_range):
+        """Create a niche envelope model estimator.
+
+        Args:
+            percentile_range: covariate values within this range are flagged as suitable habitat
+                using a narrow range like [10, 90] drops more areas from suitability maps
+                while [0, 100] creates an envelope around the full range of observed
+                covariates at all y==1 locations.
+        """
+        self.percentile_range_ = percentile_range
+        self.categorical_estimator = _features.CategoricalTransformer()
+
+    def _format_covariate_data(self, x: ArrayLike) -> Tuple[np.array, np.array]:
+        """Reads input x data and formats it to consistent array dtypes.
+
+        Args:
+            x: array-like of shape (n_samples, n_features)
+
+        Returns:
+            (continuous, categorical) tuple of ndarrays with continuous and
+                categorical covariate data.
+        """
+        if isinstance(x, np.ndarray):
+            if self.categorical_ is None:
+                con = x
+                cat = None
+            else:
+                con = x[:, self.continuous_]
+                cat = x[:, self.categorical_]
+
+        elif isinstance(x, pd.DataFrame):
+            con = x[self.continuous_pd_].to_numpy()
+            if len(self.categorical_pd_) > 0:
+                cat = x[self.categorical_pd_].to_numpy()
+            else:
+                cat = None
+
+        else:
+            raise TypeError(f"Unsupported x dtype: {type(x)}. Must be pd.DataFrame or np.array")
+
+        return con, cat
+
+    def _format_labels_and_dtypes(self, x: ArrayLike, categorical: list = None, labels: list = None) -> None:
+        """Read input x data and lists of categorical data indices and band
+            labels to format and store this info for later indexing.
+
+        Args:
+            s: array-like of shape (n_samples, n_features)
+            categorical: indices indicating which x columns are categorical
+            labels: covariate column labels. ignored if x is a pandas DataFrame
+        """
+        if isinstance(x, np.ndarray):
+            nrows, ncols = x.shape
+            if categorical is None:
+                continuous = list(range(ncols))
+            else:
+                continuous = list(set(range(ncols)).difference(set(categorical)))
+            self.labels_ = labels or make_band_labels(ncols)
+            self.categorical_ = categorical
+            self.continuous_ = continuous
+
+        elif isinstance(x, pd.DataFrame):
+            x.drop(["geometry"], axis=1, errors="ignore", inplace=True)
+            self.labels_ = labels or list(x.columns)
+
+            # store both pandas and numpy indexing of these values
+            self.continuous_pd_ = list(x.select_dtypes(exclude="category").columns)
+            self.categorical_pd_ = list(x.select_dtypes(include="category").columns)
+
+            all_columns = list(x.columns)
+            self.continuous_ = [all_columns.index(item) for item in self.continuous_pd_ if item in all_columns]
+            if len(self.categorical_pd_) != 0:
+                self.categorical_ = [all_columns.index(item) for item in self.categorical_pd_ if item in all_columns]
+            else:
+                self.categorical_ = None
+
+    def fit(self, x: ArrayLike, y: ArrayLike, categorical: list = None, labels: list = None) -> None:
+        """Trains a maxent model using a set of covariates and presence/background points.
+
+        Args:
+            x: array-like of shape (n_samples, n_features) with covariate data
+            y: array-like of shape (n_samples,) with binary presence/background (1/0) values
+            categorical: indices for which columns are categorical
+            labels: covariate labels. ignored if x is a pandas DataFrame
+        """
+        # format the input x data
+        self._format_labels_and_dtypes(x, categorical=categorical, labels=labels)
+        con, cat = self._format_covariate_data(x)
+
+        # estimate the feature range of the continuous data
+        self.feature_mins_ = np.percentile(con[y == 1], self.percentile_range_[0], axis=0)
+        self.feature_maxs_ = np.percentile(con[y == 1], self.percentile_range_[1], axis=0)
+
+        # one-hot encode the categorical data and label the classes with
+        if cat is not None:
+            ohe = self.categorical_estimator.fit_transform(cat)
+            self.in_categorical_ = np.any(ohe[y == 1], axis=0)
+
+    def predict(self, x: ArrayLike, overlay: str = "average") -> np.ndarray:
+        """Applies a model to a set of covariates or features. Requires that a model has been fit.
+
+        Args:
+            x: array-like of shape (n_samples, n_features) with covariate data
+            overlay: niche envelope overlap type.
+                select from ["average", "intersection", "union"]
+
+        Returns:
+            array-like of shape (n_samples,) with model predictions
+        """
+        overlay = overlay.lower()
+        overlay_options = ["average", "intersection", "union"]
+        assert overlay in overlay_options, f"overlay must be one of {', '.join(overlay_options)}"
+
+        # format the input data
+        con, cat = self._format_covariate_data(x)
+        nrows, ncols = x.shape
+
+        # any value within the transformed range is considered within the envelope
+        in_range = (con >= self.feature_mins_) * (con <= self.feature_maxs_)
+
+        # map the class locations where the species has been observed
+        if cat is not None:
+            ohe = self.categorical_estimator.transform(cat)
+            should_be_here = ohe[:, self.in_categorical_].any(axis=1).reshape((nrows, 1))
+            shouldnt_be_here = (~ohe[:, ~self.in_categorical_].any(axis=1)).reshape((nrows, 1))
+            in_range = np.concatenate((in_range, should_be_here, shouldnt_be_here), axis=1)
+
+        # comput envelope based on the overlay method
+        if overlay == "average":
+            ypred = np.mean(in_range, axis=1, dtype="float32")
+
+        elif overlay == "intersection":
+            ypred = np.all(in_range, axis=1).astype("uint8")
+
+        elif overlay == "union":
+            ypred = np.any(in_range, axis=1).astype("uint8")
+
+        return ypred
+
+    def fit_predict(
+        self, x: ArrayLike, y: ArrayLike, categorical: list = None, labels: list = None, overlay: str = "average"
+    ) -> np.ndarray:
+        """Trains and applies a model to x/y data.
+
+        Args:
+            x: array-like of shape (n_samples, n_features) with covariate data
+            y: array-like of shape (n_samples,) with binary presence/background (1/0) values
+            categorical: column indices indicating which columns are categorical
+            labels: Covariate labels. Ignored if x is a pandas DataFrame
+            overlay: maxent model transformation type.
+                select from ["average", "intersection", "union"]
+
+        Returns:
+            array-like of shape (n_samples,) with model predictions
+        """
+        self.fit(x, y, categorical=categorical, labels=labels)
+        return self.predict(x, overlay=overlay)
+
+
 def maxent_alpha(raw: np.ndarray) -> float:
     """Compute the sum-to-one alpha maxent model parameter.
 


Test patch:
diff --git a/tests/test_models.py b/tests/test_models.py
index d8e418a..5beeb8d 100644
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -1,3 +1,4 @@
+import numpy as np
 from sklearn import metrics
 
 from elapid import models
@@ -11,8 +12,6 @@ def test_MaxentModel_flow():
     model.fit(x, y)
     ypred = model.predict(x)
     assert len(ypred) == len(y)
-    print(ypred)
-    print(ypred.max())
     assert ypred.max() <= 1.0
     assert ypred.min() >= 0.0
 
@@ -76,3 +75,24 @@ def test_MaxentModel_feature_types():
         ypred = model.predict(x)
         auc_score = metrics.roc_auc_score(y, ypred)
         assert 0.5 <= auc_score <= 1.0
+
+
+def test_NicheEnvelopeModel():
+    # test the full range of values, ensuring all y=1 points are included
+    ne = models.NicheEnvelopeModel(percentile_range=[0, 100])
+    ne.fit(x, y)
+    union = ne.predict(x, overlay="union")
+    intersection = ne.predict(x, overlay="intersection")
+    average = ne.predict(x, overlay="average")
+
+    assert np.min([union, intersection, average]) <= 1
+    assert np.max([union, intersection, average]) >= 0
+
+    assert y.sum() == union[y == 1].sum()
+    assert intersection[y == 1].mean() <= average[y == 1].mean() <= union[y == 1].mean()
+
+    narrow = models.NicheEnvelopeModel(percentile_range=[20, 80])
+    narrow.fit(x, y)
+
+    average_narrow = narrow.predict(x, overlay="average")
+    assert average_narrow[y == 1].sum() < average[y == 1].sum()