added the possibility to specify a custom loss function

Author: mathias-von-ottenbreit

API_REFERENCE_FOR_REGRESSION.md

@@ -1,6 +1,6 @@
 # APLRRegressor
 
-## class aplr.APLRRegressor(m:int=1000, v:float=0.1, random_state:int=0, loss_function:str="mse", link_function:str="identity", n_jobs:int=0, validation_ratio:float=0.2, bins:int=300, max_interaction_level:int=1, max_interactions:int=100000, min_observations_in_split:int=20, ineligible_boosting_steps_added:int=10, max_eligible_terms:int=5, verbosity:int=0, dispersion_parameter:float=1.5, validation_tuning_metric:str="default", quantile:float=0.5, calculate_custom_validation_error_function:Optional[Callable[[npt.ArrayLike, npt.ArrayLike, npt.ArrayLike, npt.ArrayLike], float]]=None)
+## class aplr.APLRRegressor(m:int=1000, v:float=0.1, random_state:int=0, loss_function:str="mse", link_function:str="identity", n_jobs:int=0, validation_ratio:float=0.2, bins:int=300, max_interaction_level:int=1, max_interactions:int=100000, min_observations_in_split:int=20, ineligible_boosting_steps_added:int=10, max_eligible_terms:int=5, verbosity:int=0, dispersion_parameter:float=1.5, validation_tuning_metric:str="default", quantile:float=0.5, calculate_custom_validation_error_function:Optional[Callable[[npt.ArrayLike, npt.ArrayLike, npt.ArrayLike, npt.ArrayLike], float]]=None, calculate_custom_loss_function:Optional[Callable[[npt.ArrayLike, npt.ArrayLike, npt.ArrayLike, npt.ArrayLike], float]]=None, calculate_custom_negative_gradient_function:Optional[Callable[[npt.ArrayLike, npt.ArrayLike, npt.ArrayLike], npt.ArrayLike]]=None)
 
 ### Constructor parameters
 
@@ -14,7 +14,7 @@ The learning rate. Must be greater than zero and not more than one. The higher t
 Used to randomly split training observations into training and validation if ***validation_set_indexes*** is not specified when fitting.
 
 #### loss_function (default = "mse")
-Determines the loss function used. Allowed values are "mse", "binomial", "poisson", "gamma", "tweedie", "group_mse", "mae", "quantile", "negative_binomial", "cauchy" and "weibull". This is used together with ***link_function***. When ***loss_function*** is "group_mse" then the "group" argument in the ***fit*** method must be provided. In the latter case APLR will try to minimize group MSE when training the model. The ***loss_function*** "quantile" is used together with the ***quantile*** constructor parameter.
+Determines the loss function used. Allowed values are "mse", "binomial", "poisson", "gamma", "tweedie", "group_mse", "mae", "quantile", "negative_binomial", "cauchy", "weibull" and "custom_function". This is used together with ***link_function***. When ***loss_function*** is "group_mse" then the "group" argument in the ***fit*** method must be provided. In the latter case APLR will try to minimize group MSE when training the model. The ***loss_function*** "quantile" is used together with the ***quantile*** constructor parameter. When ***loss_function*** is "custom_function" then the constructor parameters ***calculate_custom_loss_function*** and ***calculate_custom_negative_gradient_function***, both described below, must be provided.
 
 #### link_function (default = "identity")
 Determines how the linear predictor is transformed to predictions. Allowed values are "identity", "logit" and "log". For an ordinary regression model use ***loss_function*** "mse" and ***link_function*** "identity". For logistic regression use ***loss_function*** "binomial" and ***link_function*** "logit". For a multiplicative model use the "log" ***link_function***. The "log" ***link_function*** often works best with a "poisson", "gamma", "tweedie", "negative_binomial" or "weibull" ***loss_function***, depending on the data. The ***loss_function*** "poisson", "gamma", "tweedie", "negative_binomial" or "weibull" should only be used with the "log" ***link_function***. Inappropriate combinations of ***loss_function*** and ***link_function*** may result in a warning message when fitting the model and/or a poor model fit. Please note that values other than "identity" typically require a significantly higher ***m*** (or ***v***) in order to converge.
@@ -56,14 +56,32 @@ Specifies which metric to use for validating the model and tuning ***m***. Avail
 Specifies the quantile to use when ***loss_function*** is "quantile".
 
 #### calculate_custom_validation_error_function (default = None)
-An optional Python function that calculates validation error if ***validation_tuning_metric*** is "custom_function". Example:
+A Python function that calculates validation error if ***validation_tuning_metric*** is "custom_function". Example:
 
 ```
 def custom_validation_error_function(y, predictions, sample_weight, group):
     squared_errors = (y-predictions)**2
     return squared_errors.mean()
 ```
 
+#### calculate_custom_loss_function (default = None)
+A Python function that calculates loss if ***loss_function*** is "custom_function". Example:
+
+```
+def custom_loss_function(y, predictions, sample_weight, group):
+    squared_errors = (y-predictions)**2
+    return squared_errors.mean()
+```
+
+#### calculate_custom_negative_gradient_function (default = None)
+A Python function that calculates the negative gradient if ***loss_function*** is "custom_function". The negative gradient should be proportional to the negative of the first order differentiation of the custom loss function (***calculate_custom_loss_function***) with respect to the predictions. Example:
+
+```
+def custom_negative_gradient_function(y, predictions, group):
+    residuals = y-predictions
+    return residuals
+```
+
 ## Method: fit(X:npt.ArrayLike, y:npt.ArrayLike, sample_weight:npt.ArrayLike = np.empty(0), X_names:List[str]=[], validation_set_indexes:List[int]=[], prioritized_predictors_indexes:List[int]=[], monotonic_constraints:List[int]=[], group:npt.ArrayLike = np.empty(0), interaction_constraints:List[int]=[])
 
 ***This method fits the model to data.***

aplr/aplr.py

@@ -9,7 +9,9 @@ def __init__(self, m:int=1000, v:float=0.1, random_state:int=0, loss_function:st
                  validation_ratio:float=0.2, bins:int=300, max_interaction_level:int=1, max_interactions:int=100000, 
                  min_observations_in_split:int=20, ineligible_boosting_steps_added:int=10, max_eligible_terms:int=5, verbosity:int=0, 
                  dispersion_parameter:float=1.5, validation_tuning_metric:str="default", quantile:float=0.5,
-                 calculate_custom_validation_error_function:Optional[Callable[[npt.ArrayLike, npt.ArrayLike, npt.ArrayLike, npt.ArrayLike], float]]=None):
+                 calculate_custom_validation_error_function:Optional[Callable[[npt.ArrayLike, npt.ArrayLike, npt.ArrayLike, npt.ArrayLike], float]]=None,
+                 calculate_custom_loss_function:Optional[Callable[[npt.ArrayLike, npt.ArrayLike, npt.ArrayLike, npt.ArrayLike], float]]=None,
+                 calculate_custom_negative_gradient_function:Optional[Callable[[npt.ArrayLike, npt.ArrayLike, npt.ArrayLike], npt.ArrayLike]]=None):
         self.m=m
         self.v=v
         self.random_state=random_state
@@ -28,6 +30,8 @@ def __init__(self, m:int=1000, v:float=0.1, random_state:int=0, loss_function:st
         self.validation_tuning_metric=validation_tuning_metric
         self.quantile=quantile
         self.calculate_custom_validation_error_function=calculate_custom_validation_error_function
+        self.calculate_custom_loss_function=calculate_custom_loss_function
+        self.calculate_custom_negative_gradient_function=calculate_custom_negative_gradient_function
 
         #Creating aplr_cpp and setting parameters
         self.APLRRegressor=aplr_cpp.APLRRegressor()
@@ -53,6 +57,8 @@ def __set_params_cpp(self):
         self.APLRRegressor.validation_tuning_metric=self.validation_tuning_metric
         self.APLRRegressor.quantile=self.quantile
         self.APLRRegressor.calculate_custom_validation_error_function=self.calculate_custom_validation_error_function
+        self.APLRRegressor.calculate_custom_loss_function=self.calculate_custom_loss_function
+        self.APLRRegressor.calculate_custom_negative_gradient_function=self.calculate_custom_negative_gradient_function
 
     def fit(self, X:npt.ArrayLike, y:npt.ArrayLike, sample_weight:npt.ArrayLike = np.empty(0), X_names:List[str]=[], validation_set_indexes:List[int]=[], prioritized_predictors_indexes:List[int]=[], monotonic_constraints:List[int]=[], group:npt.ArrayLike = np.empty(0), interaction_constraints:List[int]=[]):
         self.__set_params_cpp()
@@ -123,7 +129,9 @@ def get_params(self, deep=True):
             "dispersion_parameter":self.dispersion_parameter,
             "validation_tuning_metric":self.validation_tuning_metric,
             "quantile":self.quantile,
-            "calculate_custom_validation_error_function":self.calculate_custom_validation_error_function
+            "calculate_custom_validation_error_function":self.calculate_custom_validation_error_function,
+            "calculate_custom_loss_function":self.calculate_custom_loss_function,
+            "calculate_custom_negative_gradient_function":self.calculate_custom_negative_gradient_function
         }
 
     #For sklearn

cpp/APLRRegressor.h

@@ -104,7 +104,7 @@ class APLRRegressor
     void throw_error_if_vector_contains_negative_values(const VectorXd &y, const std::string &error_message);
     void throw_error_if_response_is_not_greater_than_zero(const VectorXd &y, const std::string &error_message);
     void throw_error_if_dispersion_parameter_is_invalid();
-    VectorXd differentiate_predictions();
+    VectorXd differentiate_predictions_wrt_linear_predictor();
     void scale_training_observations_if_using_log_link_function();
     void revert_scaling_if_using_log_link_function();
     void cap_predictions_to_minmax_in_training(VectorXd &predictions);
@@ -144,13 +144,17 @@ class APLRRegressor
     std::string validation_tuning_metric;
     double quantile;
     std::function<double(const VectorXd &y, const VectorXd &predictions, const VectorXd &sample_weight, const VectorXi &group)> calculate_custom_validation_error_function;
+    std::function<double(const VectorXd &y, const VectorXd &predictions, const VectorXd &sample_weight, const VectorXi &group)> calculate_custom_loss_function;
+    std::function<VectorXd(const VectorXd &y, const VectorXd &predictions, const VectorXi &group)> calculate_custom_negative_gradient_function;
 
     APLRRegressor(size_t m=1000,double v=0.1,uint_fast32_t random_state=std::numeric_limits<uint_fast32_t>::lowest(),std::string loss_function="mse",
         std::string link_function="identity", size_t n_jobs=0, double validation_ratio=0.2,
         size_t reserved_terms_times_num_x=100, size_t bins=300,size_t verbosity=0,size_t max_interaction_level=1,size_t max_interactions=100000,
         size_t min_observations_in_split=20, size_t ineligible_boosting_steps_added=10, size_t max_eligible_terms=5,double dispersion_parameter=1.5,
         std::string validation_tuning_metric="default", double quantile=0.5,
-        const std::function<double(VectorXd,VectorXd,VectorXd,VectorXi)> &calculate_custom_validation_error_function={});
+        const std::function<double(VectorXd,VectorXd,VectorXd,VectorXi)> &calculate_custom_validation_error_function={},
+        const std::function<double(VectorXd,VectorXd,VectorXd,VectorXi)> &calculate_custom_loss_function={},
+        const std::function<VectorXd(VectorXd,VectorXd,VectorXi)> &calculate_custom_negative_gradient_function={});
     APLRRegressor(const APLRRegressor &other);
     ~APLRRegressor();
     void fit(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight=VectorXd(0),const std::vector<std::string> &X_names={},
@@ -177,15 +181,18 @@ APLRRegressor::APLRRegressor(size_t m,double v,uint_fast32_t random_state,std::s
     double validation_ratio,size_t reserved_terms_times_num_x,size_t bins,size_t verbosity,size_t max_interaction_level,
     size_t max_interactions,size_t min_observations_in_split,size_t ineligible_boosting_steps_added,size_t max_eligible_terms,double dispersion_parameter,
     std::string validation_tuning_metric, double quantile, 
-    const std::function<double(VectorXd,VectorXd,VectorXd,VectorXi)> &calculate_custom_validation_error_function):
+    const std::function<double(VectorXd,VectorXd,VectorXd,VectorXi)> &calculate_custom_validation_error_function,
+    const std::function<double(VectorXd,VectorXd,VectorXd,VectorXi)> &calculate_custom_loss_function,
+    const std::function<VectorXd(VectorXd,VectorXd,VectorXi)> &calculate_custom_negative_gradient_function):
         reserved_terms_times_num_x{reserved_terms_times_num_x},intercept{NAN_DOUBLE},m{m},v{v},
         loss_function{loss_function},link_function{link_function},validation_ratio{validation_ratio},n_jobs{n_jobs},random_state{random_state},
         bins{bins},verbosity{verbosity},max_interaction_level{max_interaction_level},intercept_steps{VectorXd(0)},
         max_interactions{max_interactions},interactions_eligible{0},validation_error_steps{VectorXd(0)},
         min_observations_in_split{min_observations_in_split},ineligible_boosting_steps_added{ineligible_boosting_steps_added},
         max_eligible_terms{max_eligible_terms},number_of_base_terms{0},dispersion_parameter{dispersion_parameter},min_training_prediction_or_response{NAN_DOUBLE},
         max_training_prediction_or_response{NAN_DOUBLE}, validation_tuning_metric{validation_tuning_metric},
-        validation_indexes{std::vector<size_t>(0)}, quantile{quantile}, calculate_custom_validation_error_function{calculate_custom_validation_error_function}
+        validation_indexes{std::vector<size_t>(0)}, quantile{quantile}, calculate_custom_validation_error_function{calculate_custom_validation_error_function},
+        calculate_custom_loss_function{calculate_custom_loss_function},calculate_custom_negative_gradient_function{calculate_custom_negative_gradient_function}
 {
 }
 
@@ -201,7 +208,8 @@ APLRRegressor::APLRRegressor(const APLRRegressor &other):
     feature_importance{other.feature_importance},dispersion_parameter{other.dispersion_parameter},min_training_prediction_or_response{other.min_training_prediction_or_response},
     max_training_prediction_or_response{other.max_training_prediction_or_response},validation_tuning_metric{other.validation_tuning_metric},
     validation_indexes{other.validation_indexes}, quantile{other.quantile}, m_optimal{other.m_optimal},
-    calculate_custom_validation_error_function{other.calculate_custom_validation_error_function}
+    calculate_custom_validation_error_function{other.calculate_custom_validation_error_function},
+    calculate_custom_loss_function{other.calculate_custom_loss_function},calculate_custom_negative_gradient_function{other.calculate_custom_negative_gradient_function}
 {
 }
 
@@ -258,6 +266,8 @@ void APLRRegressor::throw_error_if_loss_function_does_not_exist()
         loss_function_exists=true;
     else if(loss_function=="weibull")
         loss_function_exists=true;
+    else if(loss_function=="custom_function")
+        loss_function_exists=true;
     if(!loss_function_exists)
         throw std::runtime_error("Loss function "+loss_function+" is not available in APLR.");   
 }
@@ -691,16 +701,28 @@ VectorXd APLRRegressor::calculate_neg_gradient_current(const VectorXd &sample_we
     }
     else if(loss_function=="weibull")
     {
-        output= dispersion_parameter / predictions_current.array() * ( (y_train.array()/predictions_current.array()).pow(dispersion_parameter) - 1);
+        output=dispersion_parameter / predictions_current.array() * ( (y_train.array()/predictions_current.array()).pow(dispersion_parameter) - 1);
     }    
+    else if(loss_function=="custom_function")
+    {
+        try
+        {
+            output=calculate_custom_negative_gradient_function(y_train, predictions_current, group_train);
+        }
+        catch(const std::exception& e)
+        {
+            std::string error_msg{"Error when calculating custom negative gradient function: " + static_cast<std::string>(e.what())};
+            throw std::runtime_error(error_msg);
+        }    
+    }
 
     if(link_function!="identity")
-        output=output.array()*differentiate_predictions().array();
+        output=output.array()*differentiate_predictions_wrt_linear_predictor().array();
 
     return output;
 }
 
-VectorXd APLRRegressor::differentiate_predictions()
+VectorXd APLRRegressor::differentiate_predictions_wrt_linear_predictor()
 {
     if(link_function=="logit")
         return 1.0/4.0 * (linear_predictor_current.array()/2.0).cosh().array().pow(-2);
@@ -1145,7 +1167,22 @@ void APLRRegressor::calculate_and_validate_validation_error(size_t boosting_step
 void APLRRegressor::calculate_validation_error(size_t boosting_step, const VectorXd &predictions)
 {
     if(validation_tuning_metric=="default")
-        validation_error_steps[boosting_step]=calculate_mean_error(calculate_errors(y_validation,predictions,sample_weight_validation,loss_function,dispersion_parameter,group_validation,unique_groups_validation,quantile),sample_weight_validation);
+    {
+        if(loss_function=="custom_function")
+        {
+            try
+            {
+                validation_error_steps[boosting_step] = calculate_custom_loss_function(y_validation, predictions, sample_weight_validation, group_validation);
+            }
+            catch(const std::exception& e)
+            {
+                std::string error_msg{"Error when calculating custom loss function: " + static_cast<std::string>(e.what())};
+                throw std::runtime_error(error_msg);
+            }
+        }
+        else
+            validation_error_steps[boosting_step]=calculate_mean_error(calculate_errors(y_validation,predictions,sample_weight_validation,loss_function,dispersion_parameter,group_validation,unique_groups_validation,quantile),sample_weight_validation);
+    }
     else if(validation_tuning_metric=="mse")
         validation_error_steps[boosting_step]=calculate_mean_error(calculate_errors(y_validation,predictions,sample_weight_validation,MSE_LOSS_FUNCTION),sample_weight_validation);
     else if(validation_tuning_metric=="mae")
@@ -1169,7 +1206,7 @@ void APLRRegressor::calculate_validation_error(size_t boosting_step, const Vecto
         }
         catch(const std::exception& e)
         {
-            std::string error_msg{"Error when calculating custom validation error: " + static_cast<std::string>(e.what())};
+            std::string error_msg{"Error when calculating custom validation error function: " + static_cast<std::string>(e.what())};
             throw std::runtime_error(error_msg);
         }
     }