From bfb63163119c12aaaae3a222f24ccb1f5a1b05d4 Mon Sep 17 00:00:00 2001 From: Daniel Sharp Date: Wed, 31 Jan 2024 11:37:34 -0500 Subject: [PATCH] Normalize sigmoid whitespace --- MParT/Sigmoid.h | 423 ++++++++++++++++++++++++------------------------ 1 file changed, 216 insertions(+), 207 deletions(-) diff --git a/MParT/Sigmoid.h b/MParT/Sigmoid.h index 2ded6d3e..b91b5a6f 100644 --- a/MParT/Sigmoid.h +++ b/MParT/Sigmoid.h @@ -20,21 +20,21 @@ namespace mpart { */ namespace SigmoidTypes { struct Logistic { - KOKKOS_INLINE_FUNCTION double static Evaluate(double x) { - return 0.5 + 0.5 * MathSpace::tanh(x / 2); - } - KOKKOS_INLINE_FUNCTION double static Inverse(double y) { - return y > 1 ? -MathSpace::log((1 - y) / y) : MathSpace::log(y / (1 - y)); - } - KOKKOS_INLINE_FUNCTION double static Derivative(double x) { - double fx = Evaluate(x); - return fx * (1 - fx); // Known expression for the derivative of this - } - KOKKOS_INLINE_FUNCTION double static SecondDerivative(double x) { - double fx = Evaluate(x); - return fx * (1 - fx) * - (1 - 2 * fx); // Known expression for the second derivative of this - } + KOKKOS_INLINE_FUNCTION double static Evaluate(double x) { + return 0.5 + 0.5 * MathSpace::tanh(x / 2); + } + KOKKOS_INLINE_FUNCTION double static Inverse(double y) { + return y > 1 ? -MathSpace::log((1 - y) / y) : MathSpace::log(y / (1 - y)); + } + KOKKOS_INLINE_FUNCTION double static Derivative(double x) { + double fx = Evaluate(x); + return fx * (1 - fx); // Known expression for the derivative of this + } + KOKKOS_INLINE_FUNCTION double static SecondDerivative(double x) { + double fx = Evaluate(x); + return fx * (1 - fx) * + (1 - 2 * fx); // Known expression for the second derivative of this + } }; } // namespace SigmoidTypes @@ -78,237 +78,246 @@ struct Logistic { template class Sigmoid1d { public: - /** - * @brief Construct a new Sigmoid 1d object - * - * Each input should be of length \f$n(n+1)/2\f$, where \f$n\f$ is the - * maximum order. - * - * @param centers Where to center the sigmoids - * @param widths How "wide" the sigmoids should be - * @param weights How much to weight the sigmoids linearly - */ - Sigmoid1d(Kokkos::View centers, - Kokkos::View widths, - Kokkos::View weights) - : centers_(centers), widths_(widths), weights_(weights) { - Validate(); - } + /** + * @brief Construct a new Sigmoid 1d object + * + * Each input should be of length \f$n(n+1)/2\f$, where \f$n\f$ is the + * maximum order. + * + * @param centers Where to center the sigmoids + * @param widths How "wide" the sigmoids should be + * @param weights How much to weight the sigmoids linearly + */ + Sigmoid1d(Kokkos::View centers, + Kokkos::View widths, + Kokkos::View weights) + : centers_(centers), widths_(widths), weights_(weights) { + Validate(); + } - /** - * @brief Construct a new Sigmoid 1d object from centers and widths - * - * Each input should be of length \f$n(n+1)/2\f$, where \f$n\f$ is the - * maximum order. - * - * @param centers - * @param widths - */ - Sigmoid1d(Kokkos::View centers, - Kokkos::View widths) - : centers_(centers), widths_(widths) { + /** + * @brief Construct a new Sigmoid 1d object from centers and widths + * + * Each input should be of length \f$n(n+1)/2\f$, where \f$n\f$ is the + * maximum order. + * + * @param centers + * @param widths + */ + Sigmoid1d(Kokkos::View centers, + Kokkos::View widths) + : centers_(centers), widths_(widths) { Kokkos::View weights ("Sigmoid weights", centers.extent(0)); - Kokkos::parallel_for(centers.extent(0), KOKKOS_LAMBDA(unsigned int i) { weights(i) = 1.; }); + Kokkos::parallel_for(centers.extent(0), KOKKOS_LAMBDA(unsigned int i) { weights(i) = 1.; }); weights_ = weights; Validate(); - } + } + + /** + * @brief Evaluate all sigmoids at one input + * + * @param output Place to put the output (size max_order+1) + * @param max_order Maximum order of basis function to evaluate + * @param input Point to evaluate function + */ + void EvaluateAll(double* output, int max_order, double input) const { + if (order_ < max_order) { + std::stringstream ss; + ss << "Sigmoid basis evaluation order too large.\n"; + ss << "Given order " << max_order << ", "; + ss << "can only evaluate up to order " << order_; + ProcAgnosticError::error( + ss.str().c_str()); + } - /** - * @brief Evaluate all sigmoids at one input - * - * @param output Place to put the output (size max_order+1) - * @param max_order Maximum order of basis function to evaluate - * @param input Point to evaluate function - */ - void EvaluateAll(double* output, int max_order, double input) const { - if (order_ < max_order) { - std::stringstream ss; - ss << "Sigmoid basis evaluation order too large.\n"; - ss << "Given order " << max_order << ", "; - ss << "can only evaluate up to order " << order_; - ProcAgnosticError::error( - ss.str().c_str()); - } + output[0] = 1.; + if (max_order == 0) return; + + output[1] = input; + if (max_order == 1) return; - output[0] = 1.; - if (max_order == 0) return; - output[1] = input; - if (max_order == 1) return; output[2] = -weights_(0)*EdgeType::Evaluate(-widths_(0)*(input-centers_(0))); if (max_order == 2) return; + output[3] = weights_(1)*EdgeType::Evaluate( widths_(1)*(input-centers_(1))); if (max_order == 3) return; - int param_idx = START_SIGMOIDS_IDX; - for (int curr_order = START_SIGMOIDS_ORDER; curr_order <= max_order; curr_order++) { - output[curr_order] = 0.; - for (int basis_idx = 0; basis_idx <= curr_order - START_SIGMOIDS_ORDER; basis_idx++) { - output[curr_order] += - weights_(param_idx) * - SigmoidType::Evaluate(widths_(param_idx) * - (input - centers_(param_idx))); - param_idx++; - } - } - } + int param_idx = START_SIGMOIDS_IDX; + for (int curr_order = START_SIGMOIDS_ORDER; curr_order <= max_order; curr_order++) { + output[curr_order] = 0.; + for (int basis_idx = 0; basis_idx <= curr_order - START_SIGMOIDS_ORDER; basis_idx++) { + output[curr_order] += + weights_(param_idx) * + SigmoidType::Evaluate(widths_(param_idx) * + (input - centers_(param_idx))); + param_idx++; + } + } + } + + /** + * @brief Evaluate all sigmoids up to given order and first derivatives + * + * @param output Storage for sigmoid evaluation, size max_order+1 + * @param output_diff Storage for sigmoid derivative, size max_order+1 + * @param max_order Number of sigmoids to evaluate + * @param input Where to evaluate sigmoids + */ + void EvaluateDerivatives(double* output, double* output_diff, int max_order, + double input) const { + if (order_ < max_order) { + std::stringstream ss; + ss << "Sigmoid basis evaluation order too large.\n"; + ss << "Given order " << max_order << ", "; + ss << "can only evaluate up to order " << order_; + ProcAgnosticError::error( + ss.str().c_str()); + } + + output[0] = 1.; + output_diff[0] = 0.; + if (max_order == 0) return; - /** - * @brief Evaluate all sigmoids up to given order and first derivatives - * - * @param output Storage for sigmoid evaluation, size max_order+1 - * @param output_diff Storage for sigmoid derivative, size max_order+1 - * @param max_order Number of sigmoids to evaluate - * @param input Where to evaluate sigmoids - */ - void EvaluateDerivatives(double* output, double* output_diff, int max_order, - double input) const { - if (order_ < max_order) { - std::stringstream ss; - ss << "Sigmoid basis evaluation order too large.\n"; - ss << "Given order " << max_order << ", "; - ss << "can only evaluate up to order " << order_; - ProcAgnosticError::error( - ss.str().c_str()); - } + output[1] = input; + output_diff[1] = 1.; + if (max_order == 1) return; - output[0] = 1.; - output_diff[0] = 0.; - if (max_order == 0) return; - output[1] = input; - output_diff[1] = 1.; - if (max_order == 1) return; output[2] = -weights_(0)*EdgeType::Evaluate(-widths_(0)*(input-centers_(0))); output_diff[2] = weights_(0)*widths_(0)*EdgeType::Derivative(-widths_(0)*(input-centers_(0))); if (max_order == 2) return; + output[3] = weights_(1)*EdgeType::Evaluate( widths_(1)*(input-centers_(1))); output_diff[3] = weights_(1)*widths_(1)*EdgeType::Derivative( widths_(1)*(input-centers_(1))); if (max_order == 3) return; - int param_idx = START_SIGMOIDS_IDX; - for (int curr_order = START_SIGMOIDS_ORDER; curr_order <= max_order; curr_order++) { - output[curr_order] = 0.; - output_diff[curr_order] = 0.; - for (int basis_idx = 0; basis_idx <= curr_order - START_SIGMOIDS_ORDER; basis_idx++) { - output[curr_order] += - weights_(param_idx) * - SigmoidType::Evaluate(widths_(param_idx) * - (input - centers_(param_idx))); - output_diff[curr_order] += - weights_(param_idx) * widths_(param_idx) * - SigmoidType::Derivative(widths_(param_idx) * - (input - centers_(param_idx))); - param_idx++; - } - } - } + int param_idx = START_SIGMOIDS_IDX; + for (int curr_order = START_SIGMOIDS_ORDER; curr_order <= max_order; curr_order++) { + output[curr_order] = 0.; + output_diff[curr_order] = 0.; + for (int basis_idx = 0; basis_idx <= curr_order - START_SIGMOIDS_ORDER; basis_idx++) { + output[curr_order] += + weights_(param_idx) * + SigmoidType::Evaluate(widths_(param_idx) * + (input - centers_(param_idx))); + output_diff[curr_order] += + weights_(param_idx) * widths_(param_idx) * + SigmoidType::Derivative(widths_(param_idx) * + (input - centers_(param_idx))); + param_idx++; + } + } + } - /** - * @brief Evaluate sigmoids up to given order and first+second derivatives - * - * @param output Storage for sigmoid evaluation, size max_order+1 - * @param output_diff Storage for sigmoid derivative, size max_order+1 - * @param output_diff2 Storage for sigmoid 2nd deriv, size max_order+1 - * @param max_order Maximum order of sigmoid to evaluate - * @param input Where to evaluate the sigmoids - */ - void EvaluateSecondDerivatives(double* output, double* output_diff, - double* output_diff2, int max_order, - double input) const { - if (order_ < max_order) { - std::stringstream ss; - ss << "Sigmoid basis evaluation order too large.\n"; - ss << "Given order " << max_order << ", "; - ss << "can only evaluate up to order " << order_; - ProcAgnosticError::error( - ss.str().c_str()); - } + /** + * @brief Evaluate sigmoids up to given order and first+second derivatives + * + * @param output Storage for sigmoid evaluation, size max_order+1 + * @param output_diff Storage for sigmoid derivative, size max_order+1 + * @param output_diff2 Storage for sigmoid 2nd deriv, size max_order+1 + * @param max_order Maximum order of sigmoid to evaluate + * @param input Where to evaluate the sigmoids + */ + void EvaluateSecondDerivatives(double* output, double* output_diff, + double* output_diff2, int max_order, + double input) const { + if (order_ < max_order) { + std::stringstream ss; + ss << "Sigmoid basis evaluation order too large.\n"; + ss << "Given order " << max_order << ", "; + ss << "can only evaluate up to order " << order_; + ProcAgnosticError::error( + ss.str().c_str()); + } + + output[0] = 1.; + output_diff[0] = 0.; + output_diff2[0] = 0.; + if (max_order == 0) return; + + output[1] = input; + output_diff[1] = 1.; + output_diff2[1] = 0.; + if (max_order == 1) return; - output[0] = 1.; - output_diff[0] = 0.; - output_diff2[0] = 0.; - if (max_order == 0) return; - output[1] = input; - output_diff[1] = 1.; - output_diff2[1] = 0.; - if (max_order == 1) return; output[2] = -weights_(0)*EdgeType::Evaluate(-widths_(0)*(input-centers_(0))); output_diff[2] = weights_(0)*widths_(0)*EdgeType::Derivative(-widths_(0)*(input-centers_(0))); output_diff2[2] = -weights_(0)*widths_(0)*widths_(0)*EdgeType::SecondDerivative(-widths_(0)*(input-centers_(0))); if (max_order == 2) return; + output[3] = weights_(1)*EdgeType::Evaluate( widths_(1)*(input-centers_(1))); output_diff[3] = weights_(1)*widths_(1)*EdgeType::Derivative( widths_(1)*(input-centers_(1))); output_diff2[3] = weights_(1)*widths_(1)*widths_(1)*EdgeType::SecondDerivative( widths_(1)*(input-centers_(1))); if (max_order == 3) return; - int param_idx = START_SIGMOIDS_IDX; - for (int curr_order = START_SIGMOIDS_ORDER; curr_order <= max_order; curr_order++) { - output[curr_order] = 0.; - output_diff[curr_order] = 0.; - output_diff2[curr_order] = 0.; - for (int basis_idx = 0; basis_idx <= curr_order - START_SIGMOIDS_ORDER; basis_idx++) { - output[curr_order] += - weights_(param_idx) * - SigmoidType::Evaluate(widths_(param_idx) * - (input - centers_(param_idx))); - output_diff[curr_order] += - weights_(param_idx) * widths_(param_idx) * - SigmoidType::Derivative(widths_(param_idx) * - (input - centers_(param_idx))); - output_diff2[curr_order] += - weights_(param_idx) * widths_(param_idx) * widths_(param_idx) * - SigmoidType::SecondDerivative(widths_(param_idx) * - (input - centers_(param_idx))); - param_idx++; - } - } - } + int param_idx = START_SIGMOIDS_IDX; + for (int curr_order = START_SIGMOIDS_ORDER; curr_order <= max_order; curr_order++) { + output[curr_order] = 0.; + output_diff[curr_order] = 0.; + output_diff2[curr_order] = 0.; + for (int basis_idx = 0; basis_idx <= curr_order - START_SIGMOIDS_ORDER; basis_idx++) { + output[curr_order] += + weights_(param_idx) * + SigmoidType::Evaluate(widths_(param_idx) * + (input - centers_(param_idx))); + output_diff[curr_order] += + weights_(param_idx) * widths_(param_idx) * + SigmoidType::Derivative(widths_(param_idx) * + (input - centers_(param_idx))); + output_diff2[curr_order] += + weights_(param_idx) * widths_(param_idx) * widths_(param_idx) * + SigmoidType::SecondDerivative(widths_(param_idx) * + (input - centers_(param_idx))); + param_idx++; + } + } + } - /** - * @brief Get the maximum order of this function class - * - * @return int - */ - int GetOrder() const { return order_; } + /** + * @brief Get the maximum order of this function class + * + * @return int + */ + int GetOrder() const { return order_; } private: void Validate() { if (centers_.extent(0) != widths_.extent(0) || - centers_.extent(0) != weights_.extent(0)) { - std::stringstream ss; - ss << "Sigmoid: incompatible dims of centers and widths.\n"; - ss << "centers: " << centers_.extent(0) << ", \n"; - ss << "widths: " << widths_.extent(0) << ",\n"; - ss << "weights: " << weights_.extent(0) << "\n"; - ProcAgnosticError::error( - ss.str().c_str()); - } - if (centers_.extent(0) < 2) { - std::stringstream ss; - ss << "Sigmoid: centers/widths/weights too short.\n"; - ss << "Length should be of form 2+(1+2+3+...+n) for some order n>=0"; - ProcAgnosticError::error( - ss.str().c_str()); - } - int n_sigmoid_centers = centers_.extent(0) - 2; - // Arithmetic sum length calculation - // Number of centers should be n_sigmoid_centers = num_sigmoids*(num_sigmoids+1)/2 - // Solve for num_sigmoids - double n_sig_double = (MathSpace::sqrt(1 + 8 * n_sigmoid_centers) - 1) / 2; - int n_sig = n_sig_double; - if (n_sig < 0 || MathSpace::abs((double)n_sig - n_sig_double) > 1e-15) { - std::stringstream ss; - ss << "Incorrect length of centers/widths/weights."; - ss << "Length should be of form 2+(1+2+3+...+n) for some order n"; - ProcAgnosticError::error( - ss.str().c_str()); - } - // one added for affine part of this, two added for edge terms - order_ = n_sig + 1 + 2; + centers_.extent(0) != weights_.extent(0)) { + std::stringstream ss; + ss << "Sigmoid: incompatible dims of centers and widths.\n"; + ss << "centers: " << centers_.extent(0) << ", \n"; + ss << "widths: " << widths_.extent(0) << ",\n"; + ss << "weights: " << weights_.extent(0) << "\n"; + ProcAgnosticError::error( + ss.str().c_str()); + } + if (centers_.extent(0) < 2) { + std::stringstream ss; + ss << "Sigmoid: centers/widths/weights too short.\n"; + ss << "Length should be of form 2+(1+2+3+...+n) for some order n>=0"; + ProcAgnosticError::error( + ss.str().c_str()); + } + int n_sigmoid_centers = centers_.extent(0) - 2; + // Arithmetic sum length calculation + // Number of centers should be n_sigmoid_centers = num_sigmoids*(num_sigmoids+1)/2 + // Solve for num_sigmoids + double n_sig_double = (MathSpace::sqrt(1 + 8 * n_sigmoid_centers) - 1) / 2; + int n_sig = n_sig_double; + if (n_sig < 0 || MathSpace::abs((double)n_sig - n_sig_double) > 1e-15) { + std::stringstream ss; + ss << "Incorrect length of centers/widths/weights."; + ss << "Length should be of form 2+(1+2+3+...+n) for some order n"; + ProcAgnosticError::error( + ss.str().c_str()); + } + // one added for affine part of this, two added for edge terms + order_ = n_sig + 1 + 2; } - int order_; + int order_; static int constexpr START_SIGMOIDS_ORDER = 4; static int constexpr START_SIGMOIDS_IDX = 2; - Kokkos::View centers_, widths_, weights_; + Kokkos::View centers_, widths_, weights_; }; } // namespace mpart