Multimor optimization (#75)

* add unit tests for mscombo eq and ineq constraints * changes to fix optimization. But now calcs have errors. Need to reconcile. * Fixed mscombo bmdl and bmdu calcs * Code cleanup and consolidation * Fixed handling of no selected models in multitumor * Updated selected index value for no selection to -9999 * removed commented code * fix python to match interface --------- Co-authored-by: Andy Shapiro <[email protected]> Co-authored-by: Andy Shapiro <[email protected]>
USEPA · Nov 1, 2024 · 8646354 · 8646354
1 parent b0be96d
commit 8646354
Show file tree

Hide file tree

Showing 9 changed files with 361 additions and 248 deletions.
diff --git a/src/bmdscore/bmds_helper.cpp b/src/bmdscore/bmds_helper.cpp
diff --git a/src/bmdscore/bmds_helper.h b/src/bmdscore/bmds_helper.h
@@ -304,7 +304,7 @@ struct python_multitumor_analysis{
 
 struct python_multitumor_result{
   int ndatasets; //number of models for each
-  std::vector<bool> validResult;
+  bool validResult;  //true if BMD and slope factor were both calculated for multitumor model
   std::vector<int> nmodels; //# of models per dataset (size ndatasets)
   std::vector<std::vector<python_dichotomous_model_result>> models;  //Individual model fits for each dataset nmodels[i]*ndatasets
   std::vector<int> selectedModelIndex;
@@ -482,9 +482,8 @@ double zeroin_nested(double ax,double bx, double tol,
 	      double (*f)(std::vector<double> &, double, double, struct nestedObjData*),
 	      std::vector<double> &Parms, double ck, struct nestedObjData *objData);	
 double BMD_func(int n, double p[], double x, double ck);
-double getclmt(python_multitumor_analysis *pyAnal, python_multitumor_result *pyRes, double Dose, double target, double maxDose, std::vector<double> xParms, bool isBMDL);
-double BMDL_combofunc(struct python_multitumor_analysis *pyAnal, struct python_multitumor_result *pyRes, double Dose, double D, double LR, double gtol, int *is_zero);
-double BMDU_combofunc(struct python_multitumor_analysis *pyAnal, struct python_multitumor_result *pyRes, double Dose, double D, double LR, double gtol, int *is_zero);
+double getclmt(python_multitumor_analysis *pyAnal, python_multitumor_result *pyRes, double Dose, double target, double maxDose, std::vector<double> xParms, std::vector<double> fixedParm, bool isBMDL);
+void multitumorCLs(struct python_multitumor_analysis *pyAnal, struct python_multitumor_result *pyRes, double Dose, double D, double LR, double gtol, int *is_zero);
 void Multistage_ComboBMD (struct python_multitumor_analysis *pyAnal, struct python_multitumor_result *pyRes);
 double objfunc_bmdl(const std::vector<double> &x, std::vector<double> &grad, void *my_func_data);
 double objfunc_bmdu(const std::vector<double> &x, std::vector<double> &grad, void *my_func_data);

diff --git a/src/pybmds/bmdscore.pyi b/src/pybmds/bmdscore.pyi
@@ -373,7 +373,7 @@ class python_multitumor_result:
     nmodels: list[int]
     selectedModelIndex: list[int]
     slopeFactor: float
-    validResult: list[bool]
+    validResult: bool
     def __init__(self) -> None: ...
     def setSlopeFactor(self, arg0: float) -> None: ...
 

diff --git a/src/pybmds/models/multi_tumor.py b/src/pybmds/models/multi_tumor.py
@@ -231,7 +231,7 @@ def execute(self):
         self.structs.execute()
         for i, models in enumerate(self.structs.result.models):
             for j, model in enumerate(models):
-                if model.bmdsRes.BMD > 0:
+                if model.bmdsRes.validResult:
                     bmr = self.structs.analysis.models[i][j].BMR
                     model.bmdsRes.setSlopeFactor(bmr)
         self.results = MultitumorResult.from_model(self)

diff --git a/src/pybmds/types/multi_tumor.py b/src/pybmds/types/multi_tumor.py
@@ -45,7 +45,7 @@ class MultitumorResult(BaseModel):
     models: list[list[BmdModelDichotomousSchema]]  # all degrees for all datasets
     selected_model_indexes: list[int | None]
     slope_factor: float
-    valid_result: list[bool]
+    valid_result: bool
 
     @classmethod
     def from_model(cls, model) -> Self:

diff --git a/src/tests/src/test_cpp.cpp b/src/tests/src/test_cpp.cpp
@@ -12,6 +12,7 @@
 //#include "test_cpp.h"
 #include <string>
 #include <iostream>
+#include <iomanip>
 
 void runOldDichoAnalysis();
 void runOldContAnalysis();
@@ -29,6 +30,8 @@ std::vector<double> getMultitumorPrior(int degree, int prior_cols);
 std::vector<double> getNestedPrior(int ngrp, int prior_cols, bool restricted);
 void Nlogist_probs_test();
 void Nlogist_lk_test();
+void runTestMTInequalityConstraint();
+void runTestMTEqualityConstraint();
 
 bool showResultsOverride = true;
 
@@ -42,12 +45,12 @@ int main(void){
 //  runPythonDichoAnalysis();
 //  runPythonDichoMA();
 //  runPythonContAnalysis();
-//  runPythonMultitumorAnalysis();
-  runPythonNestedAnalysis();
+  runPythonMultitumorAnalysis();
+//  runPythonNestedAnalysis();
 //  Nlogist_probs_test();
 //  Nlogist_lk_test();
 ////  runTestMultitumorModel();
-
+//  runTestMTEqualityConstraint();
   return 0;
 
 }
@@ -4042,6 +4045,9 @@ void runPythonMultitumorAnalysis(){
   int BMD_type = 1;        // 1 = extra ; added otherwise
   double BMR = 0.1;
   double alpha = 0.05;
+  std::vector<std::vector<double>> doses; 
+  std::vector<std::vector<double>> Y; 
+  std::vector<std::vector<double>> n_group;
 
   //data 1st test
 //  std::vector<double> doses1 = {0,50,100,150,200};
@@ -4053,6 +4059,41 @@ void runPythonMultitumorAnalysis(){
 //  std::vector<double> doses3 = {0,50,100,150,200};
 //  std::vector<double> Y3 = {1,68,78,88,98};
 //  std::vector<double> n_group3 = {100,100,100,100,100};
+//  std::vector<int> n = {5,5,5};
+//  std::vector<int> degree = {0,0,0};
+//  doses.push_back(doses1);
+//  doses.push_back(doses2);
+//  doses.push_back(doses3);
+//  Y.push_back(Y1);
+//  Y.push_back(Y2);
+//  Y.push_back(Y3);
+//  n_group.push_back(n_group1);
+//  n_group.push_back(n_group2);
+//  n_group.push_back(n_group3); 
+
+  //online test
+  std::vector<double> doses1 = {0,50,100,200,400};
+  std::vector<double> Y1 = {0,1,2,10,19};
+  std::vector<double> n_group1 = {20,20,20,20,20};
+  std::vector<double> doses2 = {0,50,100,200,400};
+  std::vector<double> Y2 = {0,1,2,4,11};
+  std::vector<double> n_group2 = {20,20,20,20,20};
+  std::vector<double> doses3 = {0,50,100,200,400};
+  std::vector<double> Y3 = {0,2,2,6,9};
+  std::vector<double> n_group3 = {20,20,20,20,20};
+  std::vector<int> n = {5,5,5};
+  std::vector<int> degree = {0,0,0};
+  //std::vector<int> degree = {2,2,2};
+  doses.push_back(doses1);
+  doses.push_back(doses2);
+  doses.push_back(doses3);
+  Y.push_back(Y1);
+  Y.push_back(Y2);
+  Y.push_back(Y3);
+  n_group.push_back(n_group1);
+  n_group.push_back(n_group2);
+  n_group.push_back(n_group3);
+
 
 
   //data test one recommended model 
@@ -4062,32 +4103,54 @@ void runPythonMultitumorAnalysis(){
 //  std::vector<double> doses2 = {0,50,100,200,400};
 //  std::vector<double> Y2 = {1,68,78,88,98};
 //  std::vector<double> n_group2 = {100,100,100,100,100};
+//  std::vector<int> n = {5,5};
+//  std::vector<int> degree = {0,0};
+//  doses.push_back(doses1);
+//  doses.push_back(doses2);
+//  Y.push_back(Y1);
+//  Y.push_back(Y2);
+//  n_group.push_back(n_group1);
+//  n_group.push_back(n_group2);
 
 
-  //data test no recommended model 
-  std::vector<double> doses1 = {0,50,100,200,400};
-  std::vector<double> Y1 = {1,68,78,88,98};
-  std::vector<double> n_group1 = {100,100,100,100,100};
+  //data test one dataset
+//  std::vector<double> doses1 = {0,50,100,150,200};
+//  std::vector<double> Y1 = {0,5,30,65,90};
+//  std::vector<double> n_group1 = {100,100,100,100,100};
+//  std::vector<int> n = {5};
+//  std::vector<int> degree = {0};
+//  doses.push_back(doses1);
+//  Y.push_back(Y1);
+//  n_group.push_back(n_group1);
 
-  std::vector<std::vector<double>> doses; 
-  std::vector<std::vector<double>> Y; 
-  std::vector<std::vector<double>> n_group;
-  doses.push_back(doses1);
+  //data test one recommended model 
+//  std::vector<double> doses1 = {0,50,100,150,200};
+//  std::vector<double> Y1 = {0,5,30,65,90};
+//  std::vector<double> n_group1 = {100,100,100,100,100};
+//  std::vector<double> doses2 = {0,50,100,200,400};
+//  std::vector<double> Y2 = {1,68,78,88,98};
+//  std::vector<double> n_group2 = {100,100,100,100,100};
+//  std::vector<int> n = {5,5};
+//  std::vector<int> degree = {0,0};
+//  doses.push_back(doses1);
 //  doses.push_back(doses2);
-//  doses.push_back(doses3);
-  Y.push_back(Y1);
+//  Y.push_back(Y1);
 //  Y.push_back(Y2);
-//  Y.push_back(Y3);
-  n_group.push_back(n_group1);
+//  n_group.push_back(n_group1);
 //  n_group.push_back(n_group2);
-//  n_group.push_back(n_group3); 
 
-//  std::vector<int> n = {5,5,5};
-//  std::vector<int> degree = {0,0,0};
-//  std::vector<int> n = {5,5};
-//  std::vector<int> degree = {0,0};
-  std::vector<int> n(doses.size(), 5);
-  std::vector<int> degree(doses.size(), 0);
+//  //data test no recommended model 
+//  std::vector<double> doses1 = {0,50,100,200,400};
+//  std::vector<double> Y1 = {1,68,78,88,98};
+//  std::vector<double> n_group1 = {100,100,100,100,100};
+//  std::vector<int> n = {5};
+//  std::vector<int> degree = {0};
+//  doses.push_back(doses1);
+//  Y.push_back(Y1);
+//  n_group.push_back(n_group1);
+
+
+
 /////////////////////////////////////////////////
 ////END USER INPUT
 ////////////////////////////////////////////////////
@@ -4192,16 +4255,16 @@ void runPythonMultitumorAnalysis(){
   }
 
   std::cout<<std::endl;
-  //if (res.valid){
+  if (res.validResult){
     std::cout<<"BMD:  "<<res.BMD<<std::endl;
     std::cout<<"BMDL: "<<res.BMDL<<std::endl;
     std::cout<<"BMDU: "<<res.BMDU<<std::endl;
     std::cout<<"slope factor: "<<res.slopeFactor<<std::endl;
     std::cout<<"combined LL: "<<res.combined_LL<<std::endl;
     std::cout<<"combined LL constant: "<<res.combined_LL_const<<std::endl;
-  //} else {
-  //  std::cout<<"Multitumor analysis failed"<<std::endl;
-  //}
+  } else {
+    std::cout<<"Multitumor analysis failed"<<std::endl;
+  }
 
 }
 
@@ -4815,3 +4878,113 @@ void Nlogist_lk_test(){
 //        
 //}
 
+
+void runTestMTInequalityConstraint(){
+
+  //USER INPUT
+  //double target = 565.17482044575650;
+  //double target = -104.341181352210455;
+  double target = -103.61387669479440;
+
+  const std::vector<double> x = {-2.8269246586783039,1.3081031558978955E-021,0.69271407326031031 ,1.1269537527262939E-021,-0.0000000000000000,0.54665428796084770,0.16803414281323573,-0.0000000000000000,0.40754846974911302,2.0785118843572268};
+  //const std::vector<double> x = {-2.8269246587415990,6.5404535862508677E-022,0.69271208814523166,5.6347646693637056E-022,-0.0000000000000000,0.54665480790790499,0.16803682933352726,1.9124921929208785E-021,0.40755047235337671,2.0785001207417775};
+
+  std::vector<int> degree = {2,2,2};
+  //data
+  std::vector<double> doses1 = {0,50,100,200,400};
+  std::vector<double> Y1 = {0,1,2,10,19};
+  std::vector<double> n_group1 = {20,20,20,20,20};
+  std::vector<double> doses2 = {0,50,100,200,400};
+  std::vector<double> Y2 = {0,1,2,4,11};
+  std::vector<double> n_group2 = {20,20,20,20,20};
+  std::vector<double> doses3 = {0,50,100,200,400};
+  std::vector<double> Y3 = {0,2,2,6,9};
+  std::vector<double> n_group3 = {20,20,20,20,20};
+
+  std::vector<std::vector<double>> doses;
+  std::vector<std::vector<double>> Y;
+  std::vector<std::vector<double>> n_group;
+  doses.push_back(doses1);
+  doses.push_back(doses2);
+  doses.push_back(doses3);
+  Y.push_back(Y1);
+  Y.push_back(Y2);
+  Y.push_back(Y3);
+  n_group.push_back(n_group1);
+  n_group.push_back(n_group2);
+  n_group.push_back(n_group3);
+  int nT = doses.size();
+
+  //END USER INPUT
+
+  double maxDose = 0;
+  for (int i=0; i<nT; i++){
+    double tmpMax = *std::max_element(doses[i].begin(), doses[i].end());
+    if (tmpMax > maxDose) maxDose = tmpMax;
+  }
+  std::vector<double> grad(x.size());
+   std::vector<int> nObs;
+   struct msComboInEq ineq1;
+   ineq1.nT = nT;
+   ineq1.target = target;
+
+   for(int i=0; i<nT; i++){
+     std::vector<double> scaledDose = doses[i];
+     nObs.push_back(scaledDose.size());
+     std::cout<<"dataset "<<i<<std::endl;
+     for (int j=0; j<scaledDose.size(); j++){
+       scaledDose[j] /= maxDose;
+       std::cout<<"j:"<<j<<", dose:"<<scaledDose[j]<<std::endl;
+     } 
+     ineq1.doses.push_back(scaledDose);
+     ineq1.Y.push_back(Y[i]);
+     ineq1.n_group.push_back(n_group[i]);
+   }
+   ineq1.nObs = nObs;
+   ineq1.degree = degree;
+   for (int i=0; i<nT; i++){
+     std::cout<<"i:"<<i<<", nObs:"<<nObs[i]<<std::endl;
+   }
+   std::cout<<"target:"<<target<<std::endl;
+
+
+   std::cout << std::setprecision(15);
+   double ineqVal = myInequalityConstraint1(x, grad, &ineq1);
+   std::cout<<"ineqVal = "<<ineqVal<<std::endl;
+
+   std::cout<<"ineq grad:"<<std::endl;
+   for (int i=0; i<grad.size(); i++){
+      std::cout<<"i:"<<i<<", grad:"<<grad[i]<<std::endl;
+   }
+}
+
+
+void runTestMTEqualityConstraint(){
+
+//  //USER INPUT
+
+  const std::vector<double> x = {-2.8269246586783039,1.3081031558978955E-021,0.69271407326031031 ,1.1269537527262939E-021,-0.0000000000000000,0.54665428796084770,0.16803414281323573,-0.0000000000000000,0.40754846974911302,2.0785118843572268};
+  //const std::vector<double> x = {-2.8269246587415990,6.5404535862508677E-022,0.69271208814523166,5.6347646693637056E-022,-0.0000000000000000,0.54665480790790499,0.16803682933352726,1.9124921929208785E-021,0.40755047235337671,2.0785001207417775};
+
+  std::vector<int> degree = {2,2,2};
+  int nT = degree.size(); 
+  double bmr = 0.1;
+  std::vector<double> grad(x.size());
+
+  //END USER INPUT
+
+
+   struct msComboEq eq1;
+   eq1.bmr = bmr;
+   eq1.nT = nT;
+   eq1.degree = degree;
+
+   std::cout << std::setprecision(18);
+   double ineqVal = myEqualityConstraint(x, grad, &eq1);
+   std::cout<<"eqVal = "<<ineqVal<<std::endl;
+
+   std::cout<<"eq grad:"<<std::endl;
+   for (int i=0; i<grad.size(); i++){
+      std::cout<<"i:"<<i<<", grad:"<<grad[i]<<std::endl;
+   }
+}