added the possibility to use a custom validation_tuning_metric

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)
+## 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)
 
 ### Constructor parameters
 
@@ -50,11 +50,19 @@ Limits 1) the number of terms already in the model that can be considered as int
 Specifies the variance power when ***loss_function*** is "tweedie". Specifies a dispersion parameter when ***loss_function*** is "negative_binomial", "cauchy" or "weibull". 
 
 #### validation_tuning_metric (default = "default")
-Specifies which metric to use for validating the model and tuning ***m***. Available options are "default" (using the same methodology as when calculating the training error), "mse", "mae", "negative_gini", "rankability" and "group_mse". The default is often a choice that fits well with respect to the ***loss_function*** chosen. However, if you want to use ***loss_function*** or ***dispersion_parameter*** as tuning parameters then the default is not suitable. "rankability" uses a methodology similar to the one described in https://towardsdatascience.com/how-to-calculate-roc-auc-score-for-regression-models-c0be4fdf76bb except that the metric is inverted and can be weighted by sample weights. "group_mse" requires that the "group" argument in the ***fit*** method is provided.
+Specifies which metric to use for validating the model and tuning ***m***. Available options are "default" (using the same methodology as when calculating the training error), "mse", "mae", "negative_gini", "rankability", "group_mse" and "custom_function". The default is often a choice that fits well with respect to the ***loss_function*** chosen. However, if you want to use ***loss_function*** or ***dispersion_parameter*** as tuning parameters then the default is not suitable. "rankability" uses a methodology similar to the one described in https://towardsdatascience.com/how-to-calculate-roc-auc-score-for-regression-models-c0be4fdf76bb except that the metric is inverted and can be weighted by sample weights. "group_mse" requires that the "group" argument in the ***fit*** method is provided. For "custom_function" see ***calculate_custom_validation_error_function*** below.
 
 #### quantile (default = 0.5)
 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:
+
+```
+def custom_validation_error_function(y, predictions, sample_weight, group):
+    squared_errors = (y-predictions)**2
+    return squared_errors.mean()
+```
 
 ## 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]=[])
 

aplr/aplr.py

@@ -1,11 +1,15 @@
 import numpy as np
 import numpy.typing as npt
-from typing import List
+from typing import List, Callable, Optional
 import aplr_cpp
 
 
 class APLRRegressor():
-    def __init__(self, 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):
+    def __init__(self, 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):
         self.m=m
         self.v=v
         self.random_state=random_state
@@ -23,6 +27,7 @@ def __init__(self, m:int=1000, v:float=0.1, random_state:int=0, loss_function:st
         self.dispersion_parameter=dispersion_parameter
         self.validation_tuning_metric=validation_tuning_metric
         self.quantile=quantile
+        self.calculate_custom_validation_error_function=calculate_custom_validation_error_function
 
         #Creating aplr_cpp and setting parameters
         self.APLRRegressor=aplr_cpp.APLRRegressor()
@@ -47,6 +52,7 @@ def __set_params_cpp(self):
         self.APLRRegressor.dispersion_parameter=self.dispersion_parameter
         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
 
     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()
@@ -116,7 +122,8 @@ def get_params(self, deep=True):
             "max_eligible_terms":self.max_eligible_terms,
             "dispersion_parameter":self.dispersion_parameter,
             "validation_tuning_metric":self.validation_tuning_metric,
-            "quantile":self.quantile
+            "quantile":self.quantile,
+            "calculate_custom_validation_error_function":self.calculate_custom_validation_error_function
         }
 
     #For sklearn
@@ -128,7 +135,9 @@ def set_params(self, **parameters):
 
 
 class APLRClassifier():
-    def __init__(self, m:int=9000, v:float=0.1, random_state:int=0, n_jobs:int=0, validation_ratio:float=0.2, bins:int=300, verbosity:int=0, 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):
+    def __init__(self, m:int=9000, v:float=0.1, random_state:int=0, n_jobs:int=0, validation_ratio:float=0.2, bins:int=300, verbosity:int=0, 
+                 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):
         self.m=m
         self.v=v
         self.random_state=random_state

cpp/APLRRegressor.h

@@ -143,12 +143,14 @@ class APLRRegressor
     std::vector<size_t> validation_indexes;
     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;
 
     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);
+        std::string validation_tuning_metric="default", double quantile=0.5,
+        const std::function<double(VectorXd,VectorXd,VectorXd,VectorXi)> &calculate_custom_validation_error_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={},
@@ -174,15 +176,16 @@ class APLRRegressor
 APLRRegressor::APLRRegressor(size_t m,double v,uint_fast32_t random_state,std::string loss_function,std::string link_function,size_t n_jobs,
     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):
+    std::string validation_tuning_metric, double quantile, 
+    const std::function<double(VectorXd,VectorXd,VectorXd,VectorXi)> &calculate_custom_validation_error_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}
+        validation_indexes{std::vector<size_t>(0)}, quantile{quantile}, calculate_custom_validation_error_function{calculate_custom_validation_error_function}
 {
 }
 
@@ -197,7 +200,8 @@ APLRRegressor::APLRRegressor(const APLRRegressor &other):
     max_eligible_terms{other.max_eligible_terms},number_of_base_terms{other.number_of_base_terms},
     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}
+    validation_indexes{other.validation_indexes}, quantile{other.quantile}, m_optimal{other.m_optimal},
+    calculate_custom_validation_error_function{other.calculate_custom_validation_error_function}
 {
 }
 
@@ -1157,6 +1161,18 @@ void APLRRegressor::calculate_validation_error(size_t boosting_step, const Vecto
             throw std::runtime_error("When validation_tuning_metric is group_mse then the group argument in fit() must be provided.");    
         validation_error_steps[boosting_step]=calculate_mean_error(calculate_errors(y_validation,predictions,sample_weight_validation,"group_mse",dispersion_parameter,group_validation,unique_groups_validation,quantile),sample_weight_validation);
     }
+    else if(validation_tuning_metric=="custom_function")
+    {
+        try
+        {
+            validation_error_steps[boosting_step] = calculate_custom_validation_error_function(y_validation, predictions, sample_weight_validation, group_validation);
+        }
+        catch(const std::exception& e)
+        {
+            std::string error_msg{"Error when calculating custom validation error: " + static_cast<std::string>(e.what())};
+            throw std::runtime_error(error_msg);
+        }
+    }
     else
         throw std::runtime_error(validation_tuning_metric + " is an invalid validation_tuning_metric.");
 }

cpp/pythonbinding.cpp

@@ -5,24 +5,28 @@
 #include <pybind11/operators.h>
 #include <pybind11/eigen.h>
 #include <pybind11/iostream.h>
+#include <pybind11/functional.h>
 #include "APLRRegressor.h"
 #include "APLRClassifier.h"
 
 
 namespace py = pybind11;
 
+std::function<double(VectorXd,VectorXd,VectorXd,VectorXi)> empty_calculate_custom_validation_error_function={};
+
 PYBIND11_MODULE(aplr_cpp, m) {
     py::class_<APLRRegressor>(m, "APLRRegressor",py::module_local())
         .def(py::init<int&,double&,int&,std::string&,std::string&,int&,double&,int&,int&,int&,int&,int&,int&,int&,int&,double&,std::string&,
-            double&>(),
+            double&,std::function<double(const VectorXd &y, const VectorXd &predictions, const VectorXd &sample_weight, const VectorXi &group)>&>(),
             py::arg("m")=1000,py::arg("v")=0.1,py::arg("random_state")=0,py::arg("loss_function")="mse",py::arg("link_function")="identity",
             py::arg("n_jobs")=0,py::arg("validation_ratio")=0.2,
             py::arg("reserved_terms_times_num_x")=100,py::arg("bins")=300,py::arg("verbosity")=0,
             py::arg("max_interaction_level")=1,py::arg("max_interactions")=100000,py::arg("min_observations_in_split")=20,
             py::arg("ineligible_boosting_steps_added")=10,py::arg("max_eligible_terms")=5,
             py::arg("dispersion_parameter")=1.5,
             py::arg("validation_tuning_metric")="default",
-            py::arg("quantile")=0.5
+            py::arg("quantile")=0.5,
+            py::arg("calculate_custom_validation_error_function")=empty_calculate_custom_validation_error_function
             )
         .def("fit", &APLRRegressor::fit,py::arg("X"),py::arg("y"),py::arg("sample_weight")=VectorXd(0),py::arg("X_names")=std::vector<std::string>(),
             py::arg("validation_set_indexes")=std::vector<size_t>(),py::arg("prioritized_predictors_indexes")=std::vector<size_t>(),
@@ -73,6 +77,7 @@ PYBIND11_MODULE(aplr_cpp, m) {
         .def_readwrite("validation_tuning_metric",&APLRRegressor::validation_tuning_metric)
         .def_readwrite("validation_indexes",&APLRRegressor::validation_indexes)
         .def_readwrite("quantile",&APLRRegressor::quantile)
+        .def_readwrite("calculate_custom_validation_error_function",&APLRRegressor::calculate_custom_validation_error_function)
         .def(py::pickle(
             [](const APLRRegressor &a) { // __getstate__
                 /* Return a tuple that fully encodes the state of the object */

cpp/test ALRRegressor gamma custom_validation.cpp

@@ -0,0 +1,67 @@
+#include <cmath>
+#include <iostream>
+#include <vector>
+#include <numeric>
+#include "../dependencies/eigen-master/Eigen/Dense"
+#include "APLRRegressor.h"
+#include "term.h"
+
+
+using namespace Eigen;
+
+double calculate_custom_validation_error(const VectorXd &y, const VectorXd &predictions, const VectorXd &sample_weight, const VectorXi &group)
+{
+    VectorXd error{(y.array()-predictions.array()).pow(2)};
+    return error.sum();
+}
+
+int main()
+{
+    std::vector<bool> tests;
+    tests.reserve(1000);
+
+    //Model
+    APLRRegressor model{APLRRegressor()};
+    model.m=100;
+    model.v=0.1;
+    model.bins=300;
+    model.n_jobs=0;
+    model.loss_function="gamma";
+    model.link_function="log";
+    model.verbosity=3;
+    model.max_interaction_level=0;
+    model.max_interactions=1000;
+    model.min_observations_in_split=20;
+    model.ineligible_boosting_steps_added=10;
+    model.max_eligible_terms=5;
+    model.validation_tuning_metric="custom_function";
+    model.calculate_custom_validation_error_function=calculate_custom_validation_error;
+
+    //Data    
+    MatrixXd X_train{load_csv_into_eigen_matrix<MatrixXd>("data/X_train.csv")};
+    MatrixXd X_test{load_csv_into_eigen_matrix<MatrixXd>("data/X_test.csv")}; 
+    VectorXd y_train{load_csv_into_eigen_matrix<MatrixXd>("data/y_train.csv")};    
+    VectorXd y_test{load_csv_into_eigen_matrix<MatrixXd>("data/y_test.csv")}; 
+
+    VectorXd sample_weight{VectorXd::Constant(y_train.size(),1.0)};
+
+    std::cout<<X_train;
+
+    //Fitting
+    //model.fit(X_train,y_train);
+    model.fit(X_train,y_train,sample_weight);
+    //model.fit(X_train,y_train,sample_weight,{},{0,1,2,3,4,5,10,static_cast<size_t>(y_train.size()-1)});
+    std::cout<<"feature importance\n"<<model.feature_importance<<"\n\n";
+
+    VectorXd predictions{model.predict(X_test)};
+    MatrixXd li{model.calculate_local_feature_importance(X_test)};
+
+    //Saving results
+    save_as_csv_file("data/output.csv",predictions);
+
+    std::cout<<predictions.mean()<<"\n\n";
+    tests.push_back(is_approximately_equal(predictions.mean(),23.6503,0.00001));
+
+    //Test summary
+    std::cout<<"\n\nTest summary\n"<<"Passed "<<std::accumulate(tests.begin(),tests.end(),0)<<" out of "<<tests.size()<<" tests.";
+}

setup.py

@@ -15,7 +15,7 @@
 
 setuptools.setup(
     name='aplr',
-    version='6.1.0',
+    version='6.2.0',
     description='Automatic Piecewise Linear Regression',
     ext_modules=[sfc_module],
     author="Mathias von Ottenbreit",