diff --git a/MParT/Sigmoid.h b/MParT/Sigmoid.h index 17251e97..b91b5a6f 100644 --- a/MParT/Sigmoid.h +++ b/MParT/Sigmoid.h @@ -1,245 +1,324 @@ #ifndef MPART_SIGMOID_H #define MPART_SIGMOID_H +#include +#include "MParT/PositiveBijectors.h" +#include "MParT/Utilities/MathFunctions.h" #include "MParT/Utilities/KokkosSpaceMappings.h" #include "MParT/Utilities/Miscellaneous.h" -#include "MParT/Utilities/MathFunctions.h" - -#include - -namespace mpart{ +namespace mpart { /** * @brief A small namespace to store univariate functions used in @ref Sigmoid1d - * + * * A "sigmoid" class is expected to have at least three functions: - * + * * - \c Evaluate * - \c Derivative * - \c SecondDerivative */ 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 /** - * @brief Class to represent univariate function space spanned by sigmoids - * - * @details This class stores functions of the form - * \f[f_k(x;\mathbf{c},\mathbf{b},\mathbf{w})=b_0+b_1(x-c_1)+\sum_{j=2}^kw_js(b_j(x-c_j))\f] - * + * @brief Class to represent univariate function space spanned by sigmoids. + * Order should generally be >= 4. + * + * @details This class stores functions of the form + * \f[f_{k}(x;\mathbf{c},\mathbf{b},\mathbf{w})=b_0+b_1(x-c_1)+w_2g(-b_2(x-c_2))+w_3g(b_3(x-c_3))+\sum_{j=4}^kw_js(b_j(x-c_j))\f] + * * where \f$w_j,b_j\geq 0\f$ and \f$s\f$ is a monotone increasing function. * While it is expected that \f$s\f$ is a sigmoid (e.g., Logistic or Erf * function), any monotone function will allow this class to work as expected. * Without loss of generality we set \f$b_0,b_1\equiv 1\f$ and * \f$c_1\equiv0\f$-- for a linear basis, this will be equivalent to the class * of functions presented above. If \f$k\in\{0,1\}\f$, then we revert only to - * (increasing) linear functions. - * + * (increasing) linear functions. The function \f$g\f$ should be an "edge" term + * where it is a monotone increasing function with + * \f$\lim_{x\to-\infty}g(x)=0\f$ and \f$\lim_{x\to\infty}g(x)=\infty\f$. The + * difference here between \f$g\f$ and \f$s\f$ is that \f$s\f$ is generally + * expected to approach a constant value as \f$x\to\infty\f$. + * + * * In the above expression, we denote \f$\mathbf{c}\f$ as the "centers", * \f$\mathbf{b}\f$ as the "width" or "bandwidth", and \f$\mathbf{w}\f$ as the - * "weights". If we want a function class that works with up to "order n" - * sigmoids, then we need to store one "center+width+weight" for the first - * order, then two "center+width+weight"s for the second order, then three - * "center+width+weight"s for the third order and so on. Generally, the weights - * will be uniform, but the centers and widths of these sigmoids can affect - * behavior dramatically. If performing approximation in \f$L^2(\mu)\f$, a good - * heuristic is to place the weights at the quantiles of \f$\mu\f$ and create - * widths that are half the distance to the nearest neighbor. - * - * @todo Currently, \f$b_0,b_1,c_1\f$ are not accounted for in the calculation. - * + * "weights". If we want a function class that works with up to \f$n\f$ + * sigmoids, then we need to store two "center+width+weight" for the edge + * terms, then one "center+width+weight" for the first order, then two + * "center+width+weight"s for the second order (representing two sigmoids), + * then three "center+width+weight"s for the third order and so on. Generally, + * the weights will be uniform, but the centers and widths of these sigmoids can + * affect behavior dramatically. If performing approximation in \f$L^2(\mu)\f$, + * a good heuristic is to place the weights at the quantiles of \f$\mu\f$ and + * create widths that are half the distance to the nearest neighbor. + * + * Currently the minimum order one can construct is 1, i.e., an affine map. + * * @tparam MemorySpace Where the (nonlinear) parameters are stored * @tparam SigmoidType Class defining eval, @ref SigmoidTypes */ -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) - { - 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()); - } - // Arithmetic sum length calculation - double order_double = (MathSpace::sqrt(1+8*centers.extent(0))-1)/2; - int order = order_double; - if(order <= 0 || MathSpace::abs((double)order - order_double) > 1e-15) { - std::stringstream ss; - ss << "Incorrect length of centers/widths/weights."; - ss << "Length should be of form 1+2+3+...+n for some order n"; - ProcAgnosticError::error(ss.str().c_str()); - } - order_ = order; - } - - /** - * @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)); - if(centers.extent(0) != widths.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); - ProcAgnosticError::error(ss.str().c_str()); - } - // Arithmetic sum length calculation - double order_double = (MathSpace::sqrt(1+8*centers.extent(0))-1)/2; - int order = order_double; - if(order <= 0 || MathSpace::abs((double)order - order_double) > 1e-15) { - std::stringstream ss; - ss << "Incorrect length of centers/widths/weights, " << centers.extent(0) << "."; - ss << "Length should be of form 1+2+3+...+n for some order n>0"; - ProcAgnosticError::error(ss.str().c_str()); - } - Kokkos::parallel_for(centers.extent(0), KOKKOS_LAMBDA(unsigned int i){weights(i) = 1.;}); - order_ = order; - } - - /** - * @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()); - } - - int param_idx = 0; - for(int curr_order = 0; curr_order <= max_order; curr_order++) { - output[curr_order] = 0.; - for(int basis_idx = 0; basis_idx < curr_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()); - } - int param_idx = 0; - for(int curr_order = 0; curr_order <= max_order; curr_order++) { - output[curr_order] = 0.; - output_diff[curr_order] = 0.; - for(int basis_idx = 0; basis_idx < curr_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()); - } - output[0] = 0.; - output_diff[0] = 0.; - output_diff2[0] = 0.; - int param_idx = 0; - for(int curr_order = 0; 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; 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_;} - - private: - int order_; - Kokkos::View centers_, widths_, weights_; +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 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.; }); + 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()); + } + + 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++; + } + } + } + + /** + * @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; + + 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++; + } + } + } + + /** + * @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[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++; + } + } + } + + /** + * @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; + } + int order_; + static int constexpr START_SIGMOIDS_ORDER = 4; + static int constexpr START_SIGMOIDS_IDX = 2; + Kokkos::View centers_, widths_, weights_; }; -} +} // namespace mpart -#endif //MPART_SIGMOID_H \ No newline at end of file +#endif // MPART_SIGMOID_H \ No newline at end of file diff --git a/src/AffineMap.cpp b/src/AffineMap.cpp index 84a6bc09..1e4b7baa 100644 --- a/src/AffineMap.cpp +++ b/src/AffineMap.cpp @@ -42,7 +42,6 @@ AffineMap::AffineMap(StridedMatrix A, template void AffineMap::Factorize(){ - std::cout << "Factorizing AffineMap" << std::endl; if(A_.extent(0)!=A_.extent(1)){ StridedMatrix Asub = Kokkos::subview(A_, Kokkos::ALL(), std::make_pair(A_.extent(1)-A_.extent(0),A_.extent(1))); luSolver_.compute(Asub); diff --git a/tests/MultiIndices/Test_MultiIndexSet.cpp b/tests/MultiIndices/Test_MultiIndexSet.cpp index 78534302..a57c8c2f 100644 --- a/tests/MultiIndices/Test_MultiIndexSet.cpp +++ b/tests/MultiIndices/Test_MultiIndexSet.cpp @@ -500,7 +500,6 @@ TEST_CASE("MultiIndexSet ReducedMargin Test", "[MultiIndexSet_RM]") { REQUIRE(multi.Sum() == 1); REQUIRE(multi.Get(dim_idx) == 1); } - std::cout << std::endl; unsigned int P = 2; SECTION("OrderP_RM") { MultiIndexSet mset = MultiIndexSet::CreateTotalOrder(dim, P); diff --git a/tests/Test_MultivariateExpansionWorker.cpp b/tests/Test_MultivariateExpansionWorker.cpp index c448a9f7..c19b00cf 100644 --- a/tests/Test_MultivariateExpansionWorker.cpp +++ b/tests/Test_MultivariateExpansionWorker.cpp @@ -26,13 +26,15 @@ HomogeneousEval_T CreateEvaluator(int) { template<> OffdiagHomogeneousEval_T CreateEvaluator(int dim) { ProbabilistHermite offdiag; - const int order = 3; - const int params_size = order*(order+1)/2; + const int num_sigmoids = 3; + const int params_size = 2+num_sigmoids*(num_sigmoids+1)/2; Kokkos::View centers("Sigmoid centers", params_size); Kokkos::View widths("Sigmoid widths", params_size); Kokkos::View weights("Sigmoid weights", params_size); - int basis_idx = 0; - for(int curr_order = 1; curr_order <= order; curr_order++) { + centers(0) = -3; widths(0) = 1.5; weights(0) = 1.; + centers(1) = 3; widths(1) = 1.5; weights(1) = 1.; + int basis_idx = 2; + for(int curr_order = 1; curr_order <= num_sigmoids; curr_order++) { for(int j = 0; j void TestSigmoidGradients(Function Sigmoid, unsigned int N_grad_points, double fd_delta) { Kokkos::View gradPts("Gradient points", N_grad_points); @@ -35,12 +40,12 @@ void TestSigmoidGradients(Function Sigmoid, unsigned int N_grad_points, double f for(int j = 0; j <= max_order; j++) { double fd_diff = (output_pos_fd[j]-output[j])/fd_delta; double fd_diff2 = (output_diff_pos_fd[j] - output_diff[j])/fd_delta; - REQUIRE_THAT(output_deriv[j], WithinRel(output[j], 1e-12)); - REQUIRE_THAT(output_2deriv[j], WithinRel(output[j], 1e-12)); - REQUIRE_THAT(output_deriv_pos_fd[j], WithinRel(output_pos_fd[j], 1e-12)); - REQUIRE_THAT(output_diff[j], WithinRel(fd_diff, 10*fd_delta)); - REQUIRE_THAT(output_diff_2deriv[j], WithinRel(fd_diff, 20*fd_delta)); - REQUIRE_THAT(output_diff2[j], WithinAbs(fd_diff2, fd_delta)); + CheckNearZero(output_deriv[j], output[j], 1e-12); + CheckNearZero(output_2deriv[j], output[j], 1e-12); + CheckNearZero(output_deriv_pos_fd[j], output_pos_fd[j], 1e-12); + CheckNearZero(output_diff[j], fd_diff, sqrt(fd_delta)); + CheckNearZero(output_diff_2deriv[j], fd_diff, sqrt(fd_delta)); + CheckNearZero(output_diff2[j], fd_diff2, sqrt(fd_delta)); } } } @@ -59,6 +64,22 @@ TEMPLATE_TEST_CASE("Sigmoid1d","[sigmoid1d]", SigmoidTypes::Logistic) { weights = Kokkos::View("Sigmoid weights", N_wrong); CHECK_THROWS_AS((Sigmoid1d(centers, widths, weights)), std::invalid_argument); CHECK_THROWS_AS((Sigmoid1d(centers, widths)), std::invalid_argument); + int N_wrong_arr[4] = {0, 1, 2+(1+3)}; + for(int N_wrong : N_wrong_arr) { + centers = Kokkos::View("Sigmoid Centers", N_wrong); + widths = Kokkos::View("Sigmoid widths", N_wrong); + weights = Kokkos::View("Sigmoid weights", N_wrong); + CHECK_THROWS_AS((Sigmoid1d(centers, widths, weights)), std::invalid_argument); + CHECK_THROWS_AS((Sigmoid1d(centers, widths)), std::invalid_argument); + } + centers = Kokkos::View("Sigmoid Centers", 2); + widths = Kokkos::View("Sigmoid widths", 2); + Sigmoid1d Sigmoid = Sigmoid1d(centers, widths); + CHECK(Sigmoid.GetOrder() == 1 + 2); // Affine + 2 edge terms + centers = Kokkos::View("Sigmoid Centers", 3); + widths = Kokkos::View("Sigmoid widths", 3); + Sigmoid = Sigmoid1d(centers, widths); + CHECK(Sigmoid.GetOrder() == 1 + 2 + 1); // Affine + 2 edge terms + 1 sigmoid } unsigned int N_grad_points = 100; @@ -66,8 +87,9 @@ TEMPLATE_TEST_CASE("Sigmoid1d","[sigmoid1d]", SigmoidTypes::Logistic) { const double support_bound = 100.; SECTION("Single Sigmoid") { - const int order = 1; - const int param_length = order*(order+1)/2; + const int num_sigmoid = 1; + const int order = num_sigmoid+1+2; + const int param_length = 2 + num_sigmoid*(num_sigmoid+1)/2; Kokkos::View center("Sigmoid Center", param_length); Kokkos::View width("Sigmoid Width", param_length); Kokkos::View weight("Sigmoid Weight", param_length); @@ -81,8 +103,11 @@ TEMPLATE_TEST_CASE("Sigmoid1d","[sigmoid1d]", SigmoidTypes::Logistic) { const double grid_bdry = 5.; double eval_pts[num_pts_grid+3]; eval_pts[0] = -support_bound; eval_pts[1] = 0.; eval_pts[2] = support_bound; - for(int p = 0; p < num_pts_grid; p++) eval_pts[p+3] = -grid_bdry + 2*(p-3)*grid_bdry/num_pts_grid; - double expect_output[(order+1)*3] = {0., 0., 0., 0.5, 0., 1.}; + for(int p = 0; p < num_pts_grid; p++) eval_pts[p+3] = -grid_bdry + p*2*grid_bdry/(num_pts_grid-1); + double expect_output[(order+1)*3] = { + 1., eval_pts[0], eval_pts[0], 0., 0., + 1., eval_pts[1], -std::log(2), std::log(2), 0.5, + 1., eval_pts[2], 0., eval_pts[2], 1.}; double output[(order+1)*(num_pts_grid+3)]; int j = 0; @@ -96,9 +121,9 @@ TEMPLATE_TEST_CASE("Sigmoid1d","[sigmoid1d]", SigmoidTypes::Logistic) { double prev = 0.; for(; j < num_pts_grid + 3; j++) { Sigmoid.EvaluateAll(output+j*(order+1), order, eval_pts[j]); - // SPECIFIC FOR order==1 - CHECK(output[j*2] == 0.); - double next = output[j*2+1]; + CHECK(output[j*(order+1) ] == 1.); + CHECK(output[j*(order+1)+1] == eval_pts[j]); + double next = output[j*(order+1)+3]; CHECK(next > prev); prev = next; } @@ -106,13 +131,17 @@ TEMPLATE_TEST_CASE("Sigmoid1d","[sigmoid1d]", SigmoidTypes::Logistic) { } SECTION("Multiple Sigmoids") { - const int order = 3; - const int param_length = order*(order+1)/2; + const int num_sigmoids = 3; + const int order = num_sigmoids+1+2; + const int param_length = 2 + num_sigmoids*(num_sigmoids+1)/2; Kokkos::View center("Sigmoid Center", param_length); Kokkos::View width("Sigmoid Width", param_length); Kokkos::View weight("Sigmoid Weight", param_length); - int param_idx = 0; - for(int curr_order = 1; curr_order <= order; curr_order++) { + double edge_bound = 3.; + center(0) = -edge_bound; width(0) = 2*edge_bound/10; weight(0) = 1.; + center(1) = edge_bound; width(1) = 2*edge_bound/10; weight(1) = 1.; + int param_idx = 2; + for(int curr_order = 1; curr_order <= num_sigmoids; curr_order++) { for(int i = 0; i < curr_order; i++) { center(param_idx) = 4*(-(curr_order-1)/2 + i); width(param_idx) = 1/((double)i+1); @@ -127,22 +156,33 @@ TEMPLATE_TEST_CASE("Sigmoid1d","[sigmoid1d]", SigmoidTypes::Logistic) { const double grid_bdry = 5.; double eval_pts[num_pts_grid]; for(int p = 0; p < num_pts_grid; p++) { - eval_pts[p] = -grid_bdry + 2*p*grid_bdry/num_pts_grid; + eval_pts[p] = -grid_bdry + 2*p*grid_bdry/(num_pts_grid-1); } double output[order+1]; Sigmoid.EvaluateAll(output, order, -support_bound); - CHECK(output[0] == 0.); - for(int i = 1; i <= order; i++) REQUIRE_THAT(output[i], WithinAbs(0., 1e-12)); + CHECK(output[0] == 1.); + CHECK(output[1] == -support_bound); + CHECK(output[2] == -(support_bound-edge_bound)*width(0)); + CHECK(output[3] == 0.); + for(int i = 4; i <= order; i++) { + REQUIRE_THAT(output[i], WithinAbs(0., 1e-12)); + } Sigmoid.EvaluateAll(output, order, support_bound); - CHECK(output[0] == 0.); - for(int i = 1; i <= order; i++) REQUIRE_THAT(output[i], WithinAbs(1., 1e-12)); - double prev[order+1] = {0., 0., 0., 0.}; + CHECK(output[0] == 1.); + CHECK(output[1] == support_bound); + CHECK(output[2] == 0.); + CHECK(output[3] == (support_bound-edge_bound)*width(1)); + for(int i = 4; i <= order; i++) REQUIRE_THAT(output[i], WithinAbs(1., 1e-12)); + double prev[order+1] = {0.}; + // set prev for left edge term to negative infty + prev[2] = -2*support_bound; for(int j = 0; j < num_pts_grid; j++) { Sigmoid.EvaluateAll(output, order, eval_pts[j]); - CHECK(output[0] == 0.); - for(int curr_order = 1; curr_order <= order; curr_order++) { - CHECK(output[curr_order] > prev[curr_order]); - prev[curr_order] = output[curr_order]; + CHECK(output[0] == 1.); + CHECK(output[1] == eval_pts[j]); + for(int curr_sigmoid = 2; curr_sigmoid <= order; curr_sigmoid++) { + CHECK(output[curr_sigmoid] > prev[curr_sigmoid]); + prev[curr_sigmoid] = output[curr_sigmoid]; } } TestSigmoidGradients(Sigmoid, 100, 1e-7);