test_sisso_log_regressor.cc

// Copyright 2021 Thomas A. R. Purcell
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include <descriptor_identifier/solver/SISSOLogRegressor.hpp>
#include <boost/filesystem.hpp>
#include "gtest/gtest.h"
#include <random>

namespace
{
    class SISSOLogRegressorTests : public ::testing::Test
    {
    protected:
        void SetUp() override
        {
            allowed_op_maps::set_node_maps();
            node_value_arrs::initialize_d_matrix_arr();
            mpi_setup::init_mpi_env();

            node_value_arrs::initialize_values_arr(90, 10, 3, 2, true, false);

            _task_sizes_train = {90};
            _task_sizes_test = {10};
            _leave_out_inds = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};

            std::vector<double> value_1(90, 0.0);
            std::vector<double> value_2(90, 0.0);
            std::vector<double> value_3(90, 0.0);

            std::vector<double> test_value_1(10, 0.0);
            std::vector<double> test_value_2(10, 0.0);
            std::vector<double> test_value_3(10, 0.0);

            std::default_random_engine generator;
            std::uniform_real_distribution<double> distribution_feats(0.01, 100.0);
            std::uniform_real_distribution<double> distribution_params(0.9, 1.1);

            for(int ii = 0; ii < 90; ++ii)
            {
                value_1[ii] = distribution_feats(generator);
                value_2[ii] = distribution_feats(generator);
                value_3[ii] = distribution_feats(generator);
            }

            for(int ii = 0; ii < 10; ++ii)
            {
                test_value_1[ii] = distribution_feats(generator);
                test_value_2[ii] = distribution_feats(generator);
                test_value_3[ii] = distribution_feats(generator);
            }

            node_ptr feat_1 = std::make_shared<FeatureNode>(0, "A", value_1, test_value_1, Unit("m"));
            node_ptr feat_2 = std::make_shared<FeatureNode>(1, "B", value_2, test_value_2, Unit("m"));
            node_ptr feat_3 = std::make_shared<FeatureNode>(2, "C", value_3, test_value_3, Unit("s"));

            _phi_0 ={feat_1, feat_2, feat_3};

            double a00 = distribution_params(generator);
            double a10 = distribution_params(generator);
            double c00 = distribution_feats(generator);

            _prop = std::vector<double>(90, 0.0);
            std::transform(value_1.begin(), value_1.end(), value_2.begin(), _prop.begin(), [&c00, &a00, &a10](double v1, double v2){return c00 * std::pow(v1 * v1, a00) * std::pow(v2, a10);});

            _prop_test = std::vector<double>(10, 0.0);
            std::transform(test_value_1.begin(), test_value_1.end(), test_value_2.begin(), _prop_test.begin(), [&c00, &a00, &a10](double v1, double v2){return c00 * std::pow(v1 * v1, a00) * std::pow(v2, a10);});

            _prop_zero_int = std::vector<double>(90, 0.0);
            std::transform(value_1.begin(), value_1.end(), value_2.begin(), _prop_zero_int.begin(), [&a00, &a10](double v1, double v2){return std::pow(v1 * v1, a00) * std::pow(v2, a10);});

            _prop_test_zero_int = std::vector<double>(10, 0.0);
            std::transform(test_value_1.begin(), test_value_1.end(), test_value_2.begin(), _prop_test_zero_int.begin(), [&a00, &a10](double v1, double v2){return std::pow(v1 * v1, a00) * std::pow(v2, a10);});

            _allowed_ops = {"div", "add", "mult", "sub"};
            _allowed_param_ops = {};

        }
        std::vector<std::string> _allowed_param_ops;
        std::vector<std::string> _allowed_ops;
        std::vector<node_ptr> _phi_0;

        std::vector<double> _prop;
        std::vector<double> _prop_test;

        std::vector<double> _prop_zero_int;
        std::vector<double> _prop_test_zero_int;

        std::vector<int> _task_sizes_train;
        std::vector<int> _task_sizes_test;
        std::vector<int> _leave_out_inds;
    };

    TEST_F(SISSOLogRegressorTests, FixInterceptFalseTest)
    {
#ifdef PARAMETERIZE
        std::shared_ptr<FeatureSpace> feat_space = std::make_shared<FeatureSpace>(
            mpi_setup::comm,
            _phi_0,
            _allowed_ops,
            _allowed_param_ops,
            _prop,
            _task_sizes_train,
            std::string("log_regression"),
            2,
            10,
            1,
            0,
            1.0,
            1e-50,
            1e50
        );
#else
        std::shared_ptr<FeatureSpace> feat_space = std::make_shared<FeatureSpace>(
            mpi_setup::comm,
            _phi_0,
            _allowed_ops,
            _prop,
            _task_sizes_train,
            std::string("log_regression"),
            2,
            10,
            1,
            0,
            1.0,
            1e-50,
            1e50
        );
#endif
        SISSOLogRegressor sisso(
            feat_space,
            "Property",
            Unit("m"),
            _prop,
            _prop_test,
            _task_sizes_train,
            _task_sizes_test,
            _leave_out_inds,
            2,
            2,
            3,
            false
        );
        std::vector<double> prop_comp(90, 0.0);
        std::transform(
            _prop.begin(),
            _prop.end(),
            sisso.prop().begin(),
            prop_comp.begin(),
            [](double p1, double p2){return std::abs(std::log(p1) - p2);}
        );
        EXPECT_FALSE(std::any_of(prop_comp.begin(), prop_comp.end(), [](double p){return p > 1e-10;}));

        std::transform(
            _prop_test.begin(),
            _prop_test.end(),
            sisso.prop_test().begin(),
            prop_comp.begin(),
            [](double p1, double p2){return std::abs(std::log(p1) - p2);}
        );
        EXPECT_FALSE(std::any_of(prop_comp.begin(), prop_comp.begin() + 10, [](double p){return p > 1e-10;}));

        EXPECT_EQ(sisso.n_samp(), 90);
        EXPECT_EQ(sisso.n_dim(), 2);
        EXPECT_EQ(sisso.n_residual(), 2);
        EXPECT_EQ(sisso.n_models_store(), 3);

        sisso.fit();

        EXPECT_EQ(sisso.models().size(), 2);
        EXPECT_EQ(sisso.models()[0].size(), 3);

        EXPECT_LT(sisso.models().back()[0].rmse(), 1e-7);
        EXPECT_LT(sisso.models().back()[0].test_rmse(), 1e-7);

        boost::filesystem::remove_all("feature_space/");
        boost::filesystem::remove_all("models/");
    }

    TEST_F(SISSOLogRegressorTests, FixInterceptTrueTest)
    {
#ifdef PARAMETERIZE
        std::shared_ptr<FeatureSpace> feat_space = std::make_shared<FeatureSpace>(
            mpi_setup::comm,
            _phi_0,
            _allowed_ops,
            _allowed_param_ops,
            _prop_zero_int,
            _task_sizes_train,
            "log_regression",
            2,
            10,
            1,
            0,
            1.0,
            1e-50,
            1e50
        );
#else
        std::shared_ptr<FeatureSpace> feat_space = std::make_shared<FeatureSpace>(
            mpi_setup::comm,
            _phi_0,
            _allowed_ops,
            _prop_zero_int,
            _task_sizes_train,
            "log_regression",
            2,
            10,
            1,
            0,
            1.0,
            1e-50,
            1e50
        );
#endif
        SISSOLogRegressor sisso(
            feat_space,
            "Property",
            Unit("m"),
            _prop_zero_int,
            _prop_test_zero_int,
            _task_sizes_train,
            _task_sizes_test,
            _leave_out_inds,
            2,
            2,
            3,
            true
        );

        std::vector<double> prop_comp(90, 0.0);
        std::transform(
            _prop_zero_int.begin(),
            _prop_zero_int.end(),
            sisso.prop().begin(),
            prop_comp.begin(),
            [](double p1, double p2){return std::abs(std::log(p1) - p2);}
        );
        EXPECT_FALSE(std::any_of(prop_comp.begin(), prop_comp.end(), [](double p){return p > 1e-10;}));

        std::transform(
            _prop_test_zero_int.begin(),
            _prop_test_zero_int.end(),
            sisso.prop_test().begin(),
            prop_comp.begin(),
            [](double p1, double p2){return std::abs(std::log(p1) - p2);}
        );
        EXPECT_FALSE(std::any_of(prop_comp.begin(), prop_comp.begin() + 10, [](double p){return p > 1e-10;}));

        EXPECT_EQ(sisso.n_samp(), 90);
        EXPECT_EQ(sisso.n_dim(), 2);
        EXPECT_EQ(sisso.n_residual(), 2);
        EXPECT_EQ(sisso.n_models_store(), 3);

        sisso.fit();

        EXPECT_EQ(sisso.models().size(), 2);
        EXPECT_EQ(sisso.models()[0].size(), 3);

        EXPECT_LT(sisso.models().back()[0].rmse(), 1e-7);
        EXPECT_LT(sisso.models().back()[0].test_rmse(), 1e-7);

        boost::filesystem::remove_all("feature_space/");
        boost::filesystem::remove_all("models/");
    }
}