Merge pull request #3 from scaroff/master

unit tests for two component constrained
llr-cta · Jul 17, 2018 · 0d6424e · 0d6424e
2 parents a3d99b5 + 5397ae9
commit 0d6424e
Show file tree

Hide file tree

Showing 2 changed files with 87 additions and 10 deletions.
diff --git a/src/math/pdf_1d_two_component_constrained.cpp b/src/math/pdf_1d_two_component_constrained.cpp
@@ -125,7 +125,7 @@ void TwoComponent1DConstraintPDF::set_parameter_values(ConstVecRef values)
 
   Eigen::VectorXd param_2G(5);
   param_2G[0] = prob_cpt1_;
-  param_2G[1] = 0;
+  param_2G[1] = mu1_;
   param_2G[2] = solNeg;
   param_2G[3] = mu2;
   param_2G[4] = sig2;
@@ -186,9 +186,8 @@ double TwoComponent1DConstraintPDF::value_and_parameter_gradient_1d(double x,  V
   double val = 0;
 
 
-  double dFdmu2Jacobian = 0, dFdppJacobian = 0, dFdresJacobian = 0, dFdnJacobian = 0; //dFdmu1Jacobian = 0,
+  double dFdmu1Jacobian = 0, dFdmu2Jacobian = 0, dFdppJacobian = 0, dFdresJacobian = 0, dFdnJacobian = 0;
 
-  gradient[1] = 0;
   const int OptParam = 5;
 
   if (fast_mode_){
@@ -197,28 +196,31 @@ double TwoComponent1DConstraintPDF::value_and_parameter_gradient_1d(double x,  V
 		  Eigen::MatrixXd JacobMatrix = Jacobian(x,mu1_,mu2_,pp_,res_,n_);
 		  for (int i=0; i < 5;i++){
 			  dFdppJacobian += gradient2G(i)*JacobMatrix(i,0);
-			  //~ dFdmu1Jacobian += gradient2G(i)*JacobMatrix(i,1);
+			  dFdmu1Jacobian += gradient2G(i)*JacobMatrix(i,1);
 			  dFdresJacobian += gradient2G(i)*JacobMatrix(i,2);
 			  dFdmu2Jacobian += gradient2G(i)*JacobMatrix(i,3);
 			  dFdnJacobian += gradient2G(i)*JacobMatrix(i,4);
 		  }
-		  gradient[0] = dFdppJacobian;
-		  gradient[2] = dFdresJacobian;
-		  gradient[3] = dFdmu2Jacobian;
-		  gradient[4] = dFdnJacobian;
-	 }
+	}
+	gradient[0] = dFdppJacobian;
+	gradient[1] = dFdmu1Jacobian;
+	gradient[2] = dFdresJacobian;
+	gradient[3] = dFdmu2Jacobian;
+	gradient[4] = dFdnJacobian;
+
   }
   else{
 	val = pdfTest_->value_and_parameter_gradient_1d(x,  gradient2G);
 	Eigen::MatrixXd JacobMatrix = Jacobian(x,mu1_,mu2_,pp_,res_,n_);
 	for (int i=0; i < 5;i++){
 		dFdppJacobian += gradient2G(i)*JacobMatrix(i,0);
-		//~ dFdmu1Jacobian += gradient2G(i)*JacobMatrix(i,1);
+		dFdmu1Jacobian += gradient2G(i)*JacobMatrix(i,1);
 		dFdresJacobian += gradient2G(i)*JacobMatrix(i,2);
 		dFdmu2Jacobian += gradient2G(i)*JacobMatrix(i,3);
 		dFdnJacobian += gradient2G(i)*JacobMatrix(i,4);
 		}
 	gradient[0] = dFdppJacobian;
+	gradient[1] = dFdmu1Jacobian;
 	gradient[2] = dFdresJacobian;
 	gradient[3] = dFdmu2Jacobian;
 	gradient[4] = dFdnJacobian;

diff --git a/unit_tests/math/test_pdf_1d.cpp b/unit_tests/math/test_pdf_1d.cpp
@@ -1080,6 +1080,81 @@ TEST(TestLogQuadSpline, BinnedParameterGradientCheck_ABZero) {
   }
 }
 
+
+TEST(TestTwoComponent1DConstraintPDF, SetAndRecallParameters) {
+  pdf_1d::LimitedGaussianPDF gauss1_pdf(1.0, 9.0);
+  pdf_1d::LimitedGaussianPDF gauss2_pdf(1.0, 9.0);
+  pdf_1d::TwoComponent1DPDF TwoComponent(&gauss1_pdf, "gauss_low", &gauss2_pdf, "gauss_high");
+  pdf_1d::TwoComponent1DConstraintPDF pdf(&TwoComponent);
+
+  EXPECT_EQ(pdf.num_parameters(), 5U);
+  auto pvec = pdf.parameters();
+  EXPECT_EQ(pvec[0].name, "probability");
+  EXPECT_EQ(pvec[1].name, "mean_low_gauss");
+  EXPECT_EQ(pvec[2].name, "resolution");
+  EXPECT_EQ(pvec[3].name, "mean_high_gauss");
+  EXPECT_EQ(pvec[4].name, "Width_Ratio");
+  Eigen::VectorXd p(5);
+  p << 0.4, 0, 0.4, 5.0, 0.5;
+  pdf.set_parameter_values(p);
+  EXPECT_EQ(pdf.parameter_values(), p);
+}
+
+TEST(TestTwoComponent1DConstraintPDF, GradientCheck) {
+  pdf_1d::LimitedGaussianPDF gauss1_pdf(1.0, 9.0);
+  pdf_1d::LimitedGaussianPDF gauss2_pdf(1.0, 9.0);
+  pdf_1d::TwoComponent1DPDF TwoComponent(&gauss1_pdf, "gauss_low", &gauss2_pdf, "gauss_high");
+  pdf_1d::TwoComponent1DConstraintPDF pdf(&TwoComponent);
+  Eigen::VectorXd p(5);
+  p << 0.4, 0, 0.4, 5.0, 0.5;
+  pdf.set_parameter_values(p);
+  EXPECT_EQ(p, pdf.parameter_values());
+  for(double x = 0; x<10.0; x+=0.1)
+  {
+    double good;
+    EXPECT_TRUE(gradient_check(pdf, x, 1e-6, good));
+    EXPECT_LE(good, 0.5);
+  }
+}
+
+TEST(TestTwoComponent1DConstraintPDF, HessianCheck) {
+  pdf_1d::LimitedGaussianPDF gauss1_pdf(1.0, 9.0);
+  pdf_1d::LimitedGaussianPDF gauss2_pdf(1.0, 9.0);
+  pdf_1d::TwoComponent1DPDF TwoComponent(&gauss1_pdf, "gauss_low", &gauss2_pdf, "gauss_high");
+  pdf_1d::TwoComponent1DConstraintPDF pdf(&TwoComponent);
+  Eigen::VectorXd p(5);
+  p << 0.4, 0, 0.4, 5.0, 0.5;
+  pdf.set_parameter_values(p);
+  EXPECT_EQ(p, pdf.parameter_values());
+  for(double x = 0; x<10.0; x+=0.1)
+  {
+    double good;
+    EXPECT_TRUE(hessian_check(pdf, x, 1e-6, good));
+    EXPECT_LE(good, 0.5);
+  }
+}
+
+TEST(TestTwoComponent1DConstraintPDF, ParameterGradientCheck) {
+  pdf_1d::LimitedGaussianPDF gauss1_pdf(1.0, 9.0);
+  pdf_1d::LimitedGaussianPDF gauss2_pdf(1.0, 9.0);
+  pdf_1d::TwoComponent1DPDF TwoComponent(&gauss1_pdf, "gauss_low", &gauss2_pdf, "gauss_high");
+  pdf_1d::TwoComponent1DConstraintPDF pdf(&TwoComponent);
+  Eigen::VectorXd p(5);
+  p << 0.4, 0, 0.4, 5.0, 0.5;
+  Eigen::VectorXd dp(5);
+  dp << 1e-6, 1e-6, 1e-6, 1e-6, 1e-6;
+  for(double x = 0; x<10.0; x+=0.1)
+  {
+    Eigen::VectorXd good(5);
+    EXPECT_TRUE(gradient_check_par(pdf, x, p, dp, good));
+    EXPECT_LE(good(0), 0.5);
+    EXPECT_LE(good(1), 0.5);
+    EXPECT_LE(good(2), 0.5);
+    EXPECT_LE(good(3), 0.5);
+    EXPECT_LE(good(4), 0.5);
+  }
+}
+
 int main(int argc, char **argv) {
   ::testing::InitGoogleTest(&argc, argv);
   return RUN_ALL_TESTS();