Skip to content

Commit

Permalink
Working sigmoids
Browse files Browse the repository at this point in the history
  • Loading branch information
@dannys4
dannys4 committed Sep 12, 2023
1 parent fd09255 commit b77ca66
Show file tree
Hide file tree
Showing 2 changed files with 88 additions and 40 deletions.
14 changes: 8 additions & 6 deletions MParT/Sigmoid.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,8 @@ template<typename MemorySpace, typename SigmoidType>
class Sigmoid1d: public ParameterizedFunctionBase<MemorySpace>
{
public:
Sigmoid1d(StridedVector<const double, MemorySpace> centers, StridedVector<const double, MemorySpace> widths): ParameterizedFunctionBase<MemorySpace>(1, 1, widths.extent(0))
Sigmoid1d(StridedVector<const double, MemorySpace> centers, StridedVector<const double, MemorySpace> widths):
centers_(centers), widths_(widths), ParameterizedFunctionBase<MemorySpace>(1, 1, widths.extent(0))
{
if(centers.extent(0) != widths.extent(0)) {
std::stringstream ss;
Expand All @@ -44,12 +45,12 @@ class Sigmoid1d: public ParameterizedFunctionBase<MemorySpace>
}

void EvaluateImpl(StridedMatrix<const double, MemorySpace> const& pts, StridedMatrix<double, MemorySpace> out) override {
Kokkos::parallel_for(pts.extent(1), KOKKOS_CLASS_LAMBDA (unsigned int pointInd) {
Kokkos::parallel_for(pts.extent(1), KOKKOS_CLASS_LAMBDA (unsigned int sample_index) {
double eval_pt = 0.;
for(int coeff_index = 0; coeff_index < this->numCoeffs; coeff_index++){
eval_pt += this->savedCoeffs(coeff_index)*SigmoidType::Evaluate(pts(0,pointInd));
eval_pt += this->savedCoeffs(coeff_index)*SigmoidType::Evaluate(widths_(coeff_index)*(pts(0,sample_index)-centers_(coeff_index)));
}
out(0,pointInd) = eval_pt;
out(0,sample_index) = eval_pt;
});
}

Expand All @@ -58,7 +59,7 @@ class Sigmoid1d: public ParameterizedFunctionBase<MemorySpace>
Kokkos::parallel_for(pts.extent(1), KOKKOS_CLASS_LAMBDA(unsigned int sample_index) {
double grad_pt = 0.;
for(int coeff_index = 0; coeff_index < this->numCoeffs; coeff_index++){
grad_pt += this->savedCoeffs(coeff_index)*SigmoidType::Derivative(pts(0,sample_index));
grad_pt += this->savedCoeffs(coeff_index)*widths_(coeff_index)*SigmoidType::Derivative(widths_(coeff_index)*(pts(0,sample_index)-centers_(coeff_index)));
}
out(0,sample_index) = sens(0,sample_index)*grad_pt;
});
Expand All @@ -67,12 +68,13 @@ class Sigmoid1d: public ParameterizedFunctionBase<MemorySpace>
void CoeffGradImpl(StridedMatrix<const double, MemorySpace> const& sens,StridedMatrix<const double, MemorySpace> const& pts, StridedMatrix<double, MemorySpace> out) override {
Kokkos::MDRangePolicy<Kokkos::Rank<2>, ExecutionSpace> policy ({0,0}, {this->numCoeffs, (int)sens.extent(1)});
Kokkos::parallel_for(policy, KOKKOS_CLASS_LAMBDA(unsigned int coeff_index, unsigned int sample_index) {
out(coeff_index,sample_index) = sens(0,sample_index)*SigmoidType::Evaluate(pts(0,sample_index));
out(coeff_index,sample_index) = sens(0,sample_index)*SigmoidType::Evaluate(widths_(coeff_index)*(pts(0,sample_index)-centers_(coeff_index)));
});
}

private:
using ExecutionSpace = typename MemoryToExecution<MemorySpace>::Space;
Kokkos::View<const double*, MemorySpace> centers_, widths_;
};
}

Expand Down
114 changes: 80 additions & 34 deletions tests/Test_Sigmoid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,43 @@ using namespace mpart;
using namespace Catch;
using MemorySpace = Kokkos::HostSpace;

template<typename Function>
void TestSigmoidGradients(Function Sigmoid, unsigned int N_grad_points, double fd_delta, Kokkos::View<double*, MemorySpace> coeff) {
Kokkos::View<double**, MemorySpace> gradPts("Gradient points", 1, N_grad_points);
Kokkos::View<double**, MemorySpace> gradPts_plus_delta("Gradient points plus delta", 1, N_grad_points);
Kokkos::View<double**, MemorySpace> sens("Sensitivities", 1, N_grad_points);
Kokkos::parallel_for(N_grad_points, KOKKOS_LAMBDA(unsigned int point_index) {
double gradPt = 3.0*(-1.0 + 2*((double) point_index)/((double) N_grad_points-1));
double sensPt = 2.0 + ((double) point_index)/((double) N_grad_points-1);
gradPts(0,point_index) = gradPt;
gradPts_plus_delta(0,point_index) = gradPt + fd_delta;
sens(0,point_index) = sensPt;
});
Kokkos::View<double**, MemorySpace> input_grad = Sigmoid.Gradient(sens, gradPts);
Kokkos::View<double**, MemorySpace> coeff_grad = Sigmoid.CoeffGrad(sens, gradPts);
Kokkos::View<double**, MemorySpace> gradPts_eval = Sigmoid.Evaluate(gradPts);
Kokkos::View<double**, MemorySpace> gradPts_plus_delta_eval = Sigmoid.Evaluate(gradPts_plus_delta);
Kokkos::View<double**, MemorySpace> gradPts_coeff_plus_delta_eval ("Evals for each directional add", coeff.extent(0), N_grad_points);
for(int j = 0; j < coeff.extent(0); j++) {
coeff(j) += fd_delta;
Kokkos::View<double**, MemorySpace> coeff_eval_perturb_j = Sigmoid.Evaluate(gradPts);
Kokkos::parallel_for(N_grad_points, KOKKOS_LAMBDA(unsigned int i){
gradPts_coeff_plus_delta_eval(j,i) = coeff_eval_perturb_j(0,i);
});
coeff(j) -= fd_delta;
}
double input_grad_error_accumulator = 0., coeff_grad_error_accumulator = 0.;
for(int i = 0; i < N_grad_points; i++) {
double input_grad_i_fd = (gradPts_plus_delta_eval(0,i) - gradPts_eval(0,i))/fd_delta;
input_grad_error_accumulator += fabs(input_grad_i_fd*sens(0,i) - input_grad(0,i));
for(int j = 0; j < coeff.extent(0); j++) {
double coeff_grad_ij_fd = (gradPts_coeff_plus_delta_eval(j,i) - gradPts_eval(0,i))/fd_delta;
coeff_grad_error_accumulator += fabs(coeff_grad_ij_fd*sens(0,i) - coeff_grad(j,i));
}
}
CHECK(input_grad_error_accumulator < 3*N_grad_points*fd_delta);
CHECK(coeff_grad_error_accumulator < 3*N_grad_points*coeff.extent(0)*fd_delta);
}

TEMPLATE_TEST_CASE("Sigmoid1d","[sigmoid1d]", Logistic) {
SECTION("Initialization") {
Expand All @@ -13,6 +50,11 @@ TEMPLATE_TEST_CASE("Sigmoid1d","[sigmoid1d]", Logistic) {

CHECK_THROWS_AS((Sigmoid1d<MemorySpace,TestType>(centers, widths)), std::invalid_argument);
}

unsigned int N_grad_points = 100;
double fd_delta = 1e-6;
const double support_bound = 100.;

SECTION("Single Sigmoid") {
Kokkos::View<double*, MemorySpace> center("Sigmoid Centers", 1);
Kokkos::View<double*, MemorySpace> width("Sigmoid Centers", 1);
Expand All @@ -23,49 +65,53 @@ TEMPLATE_TEST_CASE("Sigmoid1d","[sigmoid1d]", Logistic) {
coeff(0) = (double) coeff_int;
Sigmoid.WrapCoeffs(coeff);
Kokkos::View<double**, MemorySpace> evalPts("Input point", 1, 3);
evalPts(0,0) = -100; evalPts(0,1) = 0.0; evalPts(0,2) = 100;
evalPts(0,0) = -support_bound; evalPts(0,1) = 0.0; evalPts(0,2) = support_bound;
StridedMatrix<double,MemorySpace> out = Sigmoid.Evaluate(evalPts);
double approxTol = 1e-5;
REQUIRE_THAT(out(0,0), Matchers::WithinAbs(coeff_int*0.0, approxTol));
REQUIRE_THAT(out(0,1), Matchers::WithinAbs(coeff_int*0.5, approxTol));
REQUIRE_THAT(out(0,2), Matchers::WithinAbs(coeff_int*1.0, approxTol));
}

unsigned int N_grad_points = 100;
double fd_delta = 1e-5;
Kokkos::View<double**, MemorySpace> gradPts("Gradient points", 1, N_grad_points);
Kokkos::View<double**, MemorySpace> gradPts_plus_delta("Gradient points plus delta", 1, N_grad_points);
Kokkos::View<double**, MemorySpace> sens("Sensitivities", 1, N_grad_points);
Kokkos::parallel_for(N_grad_points, KOKKOS_LAMBDA(unsigned int point_index) {
double gradPt = 3.0*(-1.0 + 2*((double) point_index)/((double) N_grad_points-1));
double sensPt = 2.0 + ((double) point_index)/((double) N_grad_points-1);
gradPts(0,point_index) = gradPt;
gradPts_plus_delta(0,point_index) = gradPt + fd_delta;
sens(0,point_index) = sensPt;
});
Kokkos::View<double**, MemorySpace> input_grad = Sigmoid.Gradient(sens, gradPts);
Kokkos::View<double**, MemorySpace> coeff_grad = Sigmoid.CoeffGrad(sens, gradPts);
Kokkos::View<double**, MemorySpace> gradPts_eval = Sigmoid.Evaluate(gradPts);
Kokkos::View<double**, MemorySpace> gradPts_plus_delta_eval = Sigmoid.Evaluate(gradPts_plus_delta);
coeff(0) += fd_delta;
Kokkos::View<double**, MemorySpace> gradPts_plus_coeff_delta_eval = Sigmoid.Evaluate(gradPts);
double input_grad_error_accumulator = 0., coeff_grad_error_accumulator = 0.;
for(int i = 0; i < N_grad_points; i++) {
double input_grad_i_fd = (gradPts_plus_delta_eval(0,i) - gradPts_eval(0,i))/fd_delta;
double coeff_grad_i_fd = (gradPts_plus_coeff_delta_eval(0,i) - gradPts_eval(0,i))/fd_delta;
input_grad_error_accumulator += fabs(input_grad_i_fd*sens(0,i) - input_grad(0,i));
coeff_grad_error_accumulator += fabs(coeff_grad_i_fd*sens(0,i) - coeff_grad(0,i));
}
CHECK(input_grad_error_accumulator < N_grad_points*2*fd_delta);
CHECK(coeff_grad_error_accumulator < N_grad_points*2*fd_delta);
TestSigmoidGradients(Sigmoid, N_grad_points, fd_delta, coeff);
}
SECTION("Multiple Sigmoids") {
int N_Sigmoid = 3;
Kokkos::View<double*, MemorySpace> centers("Sigmoid Centers", N_Sigmoid);
Kokkos::View<double*, MemorySpace> widths("Sigmoid Centers", N_Sigmoid);
int N_sigmoid = 3;
Kokkos::View<double*, MemorySpace> centers("Sigmoid Centers", N_sigmoid);
Kokkos::View<double*, MemorySpace> widths("Sigmoid Centers", N_sigmoid);
Sigmoid1d<MemorySpace,TestType> Sigmoid (centers, widths);
Kokkos::View<double*, MemorySpace> coeffs("Sigmoid Coefficients", N_Sigmoid);
for(int j = 0; j < N_Sigmoid; j++) coeffs(j) = 0.5*(j+1);
Kokkos::View<double*, MemorySpace> coeffs("Sigmoid Coefficients", N_sigmoid);
for(int j = 0; j < N_sigmoid; j++) {
centers(j) = 1-0.75*j;
widths(j) = 2*(j+1);
}
Sigmoid.WrapCoeffs(coeffs);
for(int coeff_int = 1; coeff_int <= 2; coeff_int++) {
double coeff_sum = 0.;
for(int j = 1; j < N_sigmoid; j++) {
double coeff_j = 0.5*((double) N_sigmoid*(j+1));
coeffs(j) = coeff_j;
coeff_sum += coeff_j;
}
Kokkos::View<double**, MemorySpace> evalPts("Input point", 1, N_sigmoid+2);
evalPts(0,0) = -support_bound; evalPts(0,1) = support_bound;
Kokkos::View<double*, MemorySpace> expectedEval("Expected output", N_sigmoid);
for(int j = 0; j < N_sigmoid; j++) {
double pt_j = centers(j);
evalPts(0,j+2) = pt_j;
double exp_eval_j = 0.;
for(int k = 0; k < N_sigmoid; k++) {
exp_eval_j += coeffs(k)*TestType::Evaluate(widths(k)*(pt_j - centers(k)));
}
expectedEval(j) = exp_eval_j;
}
StridedMatrix<double,MemorySpace> out = Sigmoid.Evaluate(evalPts);
double approxTol = 1e-5;
REQUIRE_THAT(out(0,0), Matchers::WithinAbs(0.0 , approxTol));
REQUIRE_THAT(out(0,1), Matchers::WithinAbs(coeff_sum, approxTol));
for(int j = 0; j < N_sigmoid; j++) {
REQUIRE_THAT(out(0,j+2), Matchers::WithinAbs(expectedEval(j), approxTol));
}
}
TestSigmoidGradients(Sigmoid, N_grad_points, fd_delta, coeffs);
}
}

0 comments on commit b77ca66

Please sign in to comment.