From 3391004e7747a3ae7a8eb9257610d60d9f54c086 Mon Sep 17 00:00:00 2001 From: jialuw96 Date: Sun, 21 Apr 2024 21:48:11 -0700 Subject: [PATCH 1/7] changes and debuggings for testing --- measure_optimize.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/measure_optimize.py b/measure_optimize.py index 42e09ca..511bc12 100644 --- a/measure_optimize.py +++ b/measure_optimize.py @@ -483,7 +483,7 @@ def _check_sigma(self, error_cov, error_option): elif error_option == CovarianceStructure.time_correlation: # check the first dimension (length of DCMs) - if len(error_cov) != self.n_total_measurements: + if len(error_cov) != self.sens_info.total_measure_idx: raise ValueError( "error_cov must have the same length as n_total_measurements. Expect length:" + str(self.n_total_measurements) @@ -491,15 +491,15 @@ def _check_sigma(self, error_cov, error_option): # check the time correlation matrice shape for each DCM # loop over the index of DCM to retrieve the number of time points for DCM - for i in range(self.n_total_measurements): + for i in range(self.sens_info.total_measure_idx): # check row number - if len(error_cov[0]) != self.sens_info.Nt: + if len(error_cov[i]) != self.sens_info.Nt: raise ValueError( "error_cov[i] must have the shape Nt*Nt. Expect number of rows:" + str(self.sens_info.Nt) ) # check column number - if len(error_cov[0][0]) != self.sens_info.Nt: + if len(error_cov[i][i]) != self.sens_info.Nt: raise ValueError( "error_cov[i] must have the shape Nt*Nt. Expect number of columns:" + str(self.sens_info.Nt) @@ -697,7 +697,7 @@ def _build_sigma(self, error_cov, error_option): # different time correlation matrix for each measurement # no covariance between measurements elif error_option == CovarianceStructure.time_correlation: - for i in range(self.n_total_measurements): + for i in range(self.sens_info.total_measure_idx): # give the error covariance to Sigma # each measurement has a different time-correlation structure # that is why this is a 3D matrix @@ -707,9 +707,7 @@ def _build_sigma(self, error_cov, error_option): for t1 in range(self.sens_info.Nt): for t2 in range(self.sens_info.Nt): # for the ith measurement, the error matrix is error_cov[i] - Sigma[sigma_i_start + t1, sigma_i_start + t2] = error_cov[i][ - t1 - ][t2] + Sigma[sigma_i_start + t1, sigma_i_start + t2] = error_cov[i][t1][t2] if self.precompute_print_level >= 2: print("Error covariance matrix option:", error_option) @@ -1985,6 +1983,9 @@ def extract_store_sol(self, budget_opt, store_name): pickle.dump(fim_result, file2) file2.close() + # return the answer y, and the FIM + return ans_y, fim_result + def _locate_initial_file(self, budget): """ Given the budget, select which initial solution it should use by locating the closest budget that has stored solutions From 5e7b3f9fcf4bbe1465ea9ed4d44f11b013b45923 Mon Sep 17 00:00:00 2001 From: jialuw96 Date: Sun, 21 Apr 2024 21:48:32 -0700 Subject: [PATCH 2/7] add tests; needs comments --- test_cvxpy.py | 172 +++++++++++++++++++ test_example.py | 428 +++++++++++++++++++++++++++++++++++++++++++++++ test_function.py | 427 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 1027 insertions(+) create mode 100644 test_cvxpy.py create mode 100644 test_example.py create mode 100644 test_function.py diff --git a/test_cvxpy.py b/test_cvxpy.py new file mode 100644 index 0000000..781c633 --- /dev/null +++ b/test_cvxpy.py @@ -0,0 +1,172 @@ +import cvxpy as cp +import numpy as np +import pandas as pd +import pyomo.environ as pyo +from measure_optimize import ( + MeasurementOptimizer, + SensitivityData, + MeasurementData, + CovarianceStructure, + ObjectiveLib, +) +import pickle +import time +import matplotlib.pyplot as plt +#from idaes.core.util.model_diagnostics import DiagnosticsToolbox + + +class TestCvxpy(unittest.TestCase): + + def test_cv(self): + # set up problem formulation + + # number of time points for DCM + Nt = 8 + + # maximum manual measurement number for each measurement + max_manual_num = 10 + # minimal measurement interval + min_interval_num = 10.0 + # maximum manual measurement number for all measurements + total_max_manual_num = 10 + + # index of columns of SCM and DCM in Q + static_ind = [0, 1, 2] + dynamic_ind = [3, 4, 5] + # this index is the number of SCM + nubmer of DCM, not number of DCM timepoints + all_ind = static_ind + dynamic_ind + num_total_measure = len(all_ind) + + # meausrement names + all_names_strategy3 = [ + "CA.static", + "CB.static", + "CC.static", + "CA.dynamic", + "CB.dynamic", + "CC.dynamic", + ] + + # define static costs + static_cost = [2000, 2000, 2000, 200, 200, 200] # CA # CB # CC # CA # CB # CC + + # each static-cost measure has no per-sample cost + dynamic_cost = [0] * len(static_ind) + # each dynamic-cost measure costs $ 400 per sample + dynamic_cost.extend([400] * len(dynamic_ind)) + + # error covariance matrix + error_cov = [ + [1, 0.1, 0.1, 1, 0.1, 0.1], + [0.1, 4, 0.5, 0.1, 4, 0.5], + [0.1, 0.5, 8, 0.1, 0.5, 8], + [1, 0.1, 0.1, 1, 0.1, 0.1], + [0.1, 4, 0.5, 0.1, 4, 0.5], + [0.1, 0.5, 8, 0.1, 0.5, 8], + ] + + ## change the correlation of DCM VS SCM to half the original value + # variance + var_list = [1, 4, 8, 1, 4, 8] + # standard deviation + std_list = [np.sqrt(var_list[i]) for i in range(len(var_list))] + + # Generate correlation matrix + corr_original = [[0] * len(var_list) for i in range(len(var_list))] + + # compute correlation, corr[i,j] = cov[i,j]/std[i]/std[j] + for i in range(len(var_list)): + # loop over Nm + for j in range(len(var_list)): + corr_original[i][j] = error_cov[i][j] / std_list[i] / std_list[j] + + # DCM-SCM correlations are half of DCM-DCM or SCM-SCM correlation + corr_num = 0.5 + + # loop over SCM rows + for i in range(len(static_ind)): + # loop over DCM columns + for j in range(len(static_ind), len(var_list)): + # correlation matrix is symmetric + # we make DCM-SCM correlation half of its original value + corr_original[i][j] *= corr_num + corr_original[j][i] *= corr_num + + # update covariance matrix + # change back from correlation matrix + for i in range(len(var_list)): + for j in range(len(var_list)): + error_cov[i][j] = corr_original[i][j] * std_list[i] * std_list[j] + + ## define MeasurementData object + measure_info = MeasurementData( + all_names_strategy3, # name string + all_ind, # jac_index: measurement index in Q + static_cost, # static costs + dynamic_cost, # dynamic costs + min_interval_num, # minimal time interval between two timepoints + max_manual_num, # maximum number of timepoints for each measurement + total_max_manual_num, # maximum number of timepoints for all measurement + ) + + + # create data object to pre-compute Qs + # read jacobian from the source csv + # Nt is the number of time points for each measurement + jac_info = SensitivityData("./kinetics_source_data/Q_drop0.csv", Nt) + static_measurement_index = [ + 0, + 1, + 2, + ] # the index of CA, CB, CC in the jacobian array, considered as SCM + dynamic_measurement_index = [ + 0, + 1, + 2, + ] # the index of CA, CB, CC in the jacobian array, also considered as DCM + jac_info.get_jac_list( + static_measurement_index, # the index of SCMs in the jacobian array + dynamic_measurement_index, + ) # the index of DCMs in the jacobian array + + # use MeasurementOptimizer to pre-compute the unit FIMs + calculator = MeasurementOptimizer( + jac_info, # SensitivityData object + measure_info, # MeasurementData object + error_cov=error_cov, # error covariance matrix + error_opt=CovarianceStructure.measure_correlation, # error covariance options + print_level = 3 + ) + + # calculate a list of unit FIMs + calculator.assemble_unit_fims() + + # optimization options + objective = ObjectiveLib.A + + budget_opt = 5000 + + mixed_integer_opt = True + + file_store_name = "./cvxpy_results/cvxpy_A_" + #file_store_name = None + + num_dynamic_time = np.linspace(0, 60, 9) + + static_dynamic = [[0, 3], [1, 4], [2, 5]] # These pairs cannot be chosen simultaneously + time_interval_for_all = True + + # map the timepoint index to its real time + dynamic_time_dict = {} + for i, tim in enumerate(num_dynamic_time[1:]): + dynamic_time_dict[i] = np.round(tim, decimals=2) + + calculator.continuous_optimization_cvxpy(objective=objective, + mixed_integer = mixed_integer_opt, + budget=budget_opt, + static_dynamic_pair = static_dynamic, + time_interval_all_dynamic = time_interval_for_all, + num_dynamic_t_name=num_dynamic_time, # number of time points of DCMs + solver="MOSEK", + store_name = file_store_name) + diff --git a/test_example.py b/test_example.py new file mode 100644 index 0000000..4090e02 --- /dev/null +++ b/test_example.py @@ -0,0 +1,428 @@ +import numpy as np +import pandas as pd +import pyomo.environ as pyo +from measure_optimize import ( + MeasurementOptimizer, + SensitivityData, + MeasurementData, + CovarianceStructure, + ObjectiveLib, +) +import pickle +import time +import unittest + + + +class TestExample(unittest.TestCase): + + def test_eg(self): + ### STEP 1: set up measurement cost strategy + + # number of time points for DCM + Nt = 8 + + # maximum manual measurement number for each measurement + max_manual_num = 10 + # minimal measurement interval + min_interval_num = 10.0 + # maximum manual measurement number for all measurements + total_max_manual_num = 10 + + # index of columns of SCM and DCM in Q + # SCM measurements: # CA # CB # CC. This is the index of measurement column in the reformulated Jacobian + static_ind = [0, 1, 2] + # DCM measurements: # CA # CB # CC. This is the index of measurement column in the reformulated Jacobian + dynamic_ind = [3, 4, 5] + # this index is the number of SCM + nubmer of DCM, not number of DCM timepoints + all_ind = static_ind + dynamic_ind + num_total_measure = len(all_ind) + + # meausrement names + all_names_strategy3 = [ + "CA.static", + "CB.static", + "CC.static", + "CA.dynamic", + "CB.dynamic", + "CC.dynamic", + ] + + # define static costs in $ + # CA # CB # CC # CA # CB # CC + static_cost = [2000, 2000, 2000, 200, 200, 200] + + # each static-cost measure has no per-sample cost + dynamic_cost = [0] * len(static_ind) + # each dynamic-cost measure costs $ 400 per sample + dynamic_cost.extend([400] * len(dynamic_ind)) + + ## define MeasurementData object + measure_info = MeasurementData( + all_names_strategy3, # name string + all_ind, # jac_index: measurement index in Q + static_cost, # static costs + dynamic_cost, # dynamic costs + min_interval_num, # minimal time interval between two timepoints + max_manual_num, # maximum number of timepoints for each measurement + total_max_manual_num, # maximum number of timepoints for all measurement + ) + + + ### STEP 2: Read and create Jacobian object + + # error covariance matrix + error_cov = [ + [1, 0.1, 0.1, 1, 0.1, 0.1], + [0.1, 4, 0.5, 0.1, 4, 0.5], + [0.1, 0.5, 8, 0.1, 0.5, 8], + [1, 0.1, 0.1, 1, 0.1, 0.1], + [0.1, 4, 0.5, 0.1, 4, 0.5], + [0.1, 0.5, 8, 0.1, 0.5, 8], + ] + + ## change the correlation of DCM VS SCM to half the original value + # variance + var_list = [error_cov[i][i] for i in range(len(error_cov))] + # standard deviation + std_list = [np.sqrt(var_list[i]) for i in range(len(var_list))] + + # Generate correlation matrix + corr_original = [[0] * len(var_list) for i in range(len(var_list))] + + # compute correlation, corr[i,j] = cov[i,j]/std[i]/std[j] + for i in range(len(var_list)): + # loop over Nm + for j in range(len(var_list)): + corr_original[i][j] = error_cov[i][j] / std_list[i] / std_list[j] + + # DCM-SCM correlations are half of DCM-DCM or SCM-SCM correlation + corr_num = 0.5 + + # loop over SCM rows + for i in range(len(static_ind)): + # loop over DCM columns + for j in range(len(static_ind), len(var_list)): + # correlation matrix is symmetric + # we make DCM-SCM correlation half of its original value + corr_original[i][j] *= corr_num + corr_original[j][i] *= corr_num + + # update covariance matrix + # change back from correlation matrix + for i in range(len(var_list)): + for j in range(len(var_list)): + error_cov[i][j] = corr_original[i][j] * std_list[i] * std_list[j] + + + # create data object to pre-compute Qs + # read jacobian from the source csv + # Nt is the number of time points for each measurement + # This csv file looks like this: + # A1 A2 E1 E2 + #Unnamed: 0 + #1 -1.346695 -1.665335e-13 2.223380 -1.665335e-13 + # .... + + jac_info = SensitivityData("./kinetics_source_data/Q_drop0.csv", Nt) + static_measurement_index = [ + 0, + 1, + 2, + ] # the index of CA, CB, CC in the jacobian array, considered as SCM + dynamic_measurement_index = [ + 0, + 1, + 2, + ] # the index of CA, CB, CC in the jacobian array, also considered as DCM + jac_info.get_jac_list( + static_measurement_index, # the index of SCMs in the jacobian array + dynamic_measurement_index, + ) # the index of DCMs in the jacobian array + + + ### STEP 3: Create MeasurementOptimizer object, precomputation atom FIMs + + + # use MeasurementOptimizer to pre-compute the unit FIMs + calculator = MeasurementOptimizer( + jac_info, # SensitivityData object + measure_info, # MeasurementData object + error_cov=error_cov, # error covariance matrix + error_opt=CovarianceStructure.measure_correlation, # error covariance options + print_level=3, # I use highest here to see more information + ) + + # calculate a list of unit FIMs + calculator.assemble_unit_fims() + + + + + def kinetics_experiment(mip_option, + objective, + small_element=0.0001, + file_store_name="test", + initializer_option="lp_A", + rerun_all_paper_results=False, + linear_solver_opt="ma57"): + """ + This function runs the MO experiments. + + Arguments + --------- + mip_option: boolean, mixed integer problem or not + objective: objective function, chosen from ObjectiveLib.A or ObjectiveLib.D + small_element: scaler, the small element added to the diagonal of FIM + file_store_name: string, save result names with this string + initializer_option: string, choose what solutions to initialize from, lp_A, nlp_D, milp_A, minlp_D are options + rerun_all_paper_results: boolean, if run all results or just sensitivity test + linear_solver_opt: choose linear solver here. Directly comment this line or give None if default linear solver is used. + + Return + ------ + None. + """ + # optimization options + fixed_nlp_opt = False + mix_obj_option = False + alpha_opt = 0.9 + + sparse_opt = True + fix_opt = False + + + num_dynamic_time = np.linspace(0, 60, 9) + + static_dynamic = [[0, 3], [1, 4], [2, 5]] # These pairs cannot be chosen simultaneously + time_interval_for_all = True + + # map the timepoint index to its real time + dynamic_time_dict = {} + for i, tim in enumerate(num_dynamic_time[1:]): + dynamic_time_dict[i] = np.round(tim, decimals=2) + + if rerun_all_paper_results: + # give range of budgets for this case + budget_ranges = np.linspace(1000, 5000, 11) + else: + # give a trial ranges for a test; we use the first 3 budgets in budget_ranges + budget_ranges = [3800, 4600] + + # ==== initialization strategy ==== + # according to the initializer option, we provide different sets of initialization files + if initializer_option == "milp_A": + # all budgets + curr_results = np.linspace(1000, 5000, 11) + # initial solution file path and name + file_name_pre, file_name_end = "./kinetics_results/MILP_", "_a" + + elif initializer_option == "minlp_D": + # all budgets + curr_results = np.linspace(1000, 5000, 11) + # initial solution file path and name + file_name_pre, file_name_end = "./kinetics_results/MINLP_", "_d" + + elif initializer_option == "lp_A": + # all budgets + curr_results = np.linspace(1000, 5000, 41) + # initial solution file path and name + file_name_pre, file_name_end = "./kinetics_results/LP_", "_a" + + elif initializer_option == "nlp_D": + # all budgets + curr_results = np.linspace(1000, 5000, 41) + # initial solution file path and name + file_name_pre, file_name_end = "./kinetics_results/NLP_", "_d" + + # initialize the initial solution dict. key: budget. value: initial solution file name + # this initialization dictionary provides a more organized input format for initialization + initial_solution = {} + # loop over budget + for b in curr_results: + initial_solution[b] = file_name_pre + str(int(b)) + file_name_end + + # store objective and solutions + sol_list = [] + obj_list = [] + + # ===== run a test for a few budgets ===== + + # use a starting budget to create the model + start_budget = budget_ranges[0] + # timestamp for creating pyomo model + t1 = time.time() + # call the optimizer function to formulate the model and solve for the first time + # optimizer method will 1) create the model and save as self.mod 2) initialize the model + calculator.optimizer( + start_budget, # budget + initial_solution, # a collection of initializations + mixed_integer=mip_option, # if relaxing integer decisions + obj=objective, # objective function options, A or D + mix_obj=mix_obj_option, # if mixing A- and D-optimality to be the OF + alpha=alpha_opt, # the weight of A-optimality if using mixed obj + fixed_nlp=fixed_nlp_opt, # if it is a fixed NLP problem + fix=fix_opt, # if it is a squared problem + upper_diagonal_only=sparse_opt, # if only defining upper triangle part for symmetric matrix + num_dynamic_t_name=num_dynamic_time, # number of time points of DCMs + static_dynamic_pair=static_dynamic, # if one measurement can be both SCM and DCM + time_interval_all_dynamic=time_interval_for_all, # time interval for time points of DCMs + fim_diagonal_small_element=small_element, # a small element added for FIM diagonals to avoid ill-conditioning + print_level=1, + ) # print level for optimization part + # timestamp for solving pyomo model + t2 = time.time() + calculator.solve(mip_option=mip_option, objective=objective, linear_solver=linear_solver_opt) + # timestamp for finishing + t3 = time.time() + print("model and solver wall clock time:", t3 - t1) + print("solver wall clock time:", t3 - t2) + res, fim = calculator.extract_store_sol(start_budget, file_store_name) + sol_list.append(res) + obj_list.append(fim) + + + # continue to run the rest of budgets if the test goes well + for b in budget_ranges[1:]: + print("====Solving with budget:", b, "====") + # open the update toggle every time so no need to create model every time + calculator.update_budget(b) + # solve the model + calculator.solve(mip_option=mip_option, objective=objective, linear_solver=linear_solver_opt) + # extract and select solutions + res, fim = calculator.extract_store_sol(b, file_store_name) + sol_list.append(res) + obj_list.append(fim) + + return sol_list, obj_list + + ### Test A-optimality, LP + + ## The most often used options are listed as arguments here + mip_option_value = False # mixed integer problem or not + objective_value = ObjectiveLib.A # objective, ObjectiveLib.A or ObjectiveLib.D + small_element_value = 0.0001 # the small element added to the diagonal of FIM + # initialize with A-opt. MILP solutions + # choose what solutions to initialize from: + # minlp_D: initialize with minlp_D solutions + # milp_A: initialize with milp_A solutions + # lp_A: iniitalize with lp_A solution + # nlp_D: initialize with nlp_D solution + initializer_option_value = "lp_A" + # if run all results or just sensitivity test + rerun_all_paper_results_value = False + + + sol_set, obj_set = kinetics_experiment(mip_option_value, + objective_value, + small_element=small_element_value, + file_store_name=file_store_name_value, + initializer_option=initializer_option_value, + rerun_all_paper_results=rerun_all_paper_results_value, + linear_solver_opt=linear_solver_opt_value) + + + self.assertAlmostEqual(np.trace(obj_set[0]), 153.05, places=1) + self.assertAlmostEqual(np.linalg.det(obj_set[0]), 0.383, places=1) + self.assertAlmostEqual(np.trace(obj_set[1]), 168.659, places=1) + self.assertAlmostEqual(np.linalg.det(obj_set[1]), 5.3288, places=1) + + + ### Test A-optimality, MILP + + ## The most often used options are listed as arguments here + mip_option_value = True # mixed integer problem or not + objective_value = ObjectiveLib.A # objective, ObjectiveLib.A or ObjectiveLib.D + small_element_value = 0.0001 # the small element added to the diagonal of FIM + # initialize with A-opt. MILP solutions + # choose what solutions to initialize from: + # minlp_D: initialize with minlp_D solutions + # milp_A: initialize with milp_A solutions + # lp_A: iniitalize with lp_A solution + # nlp_D: initialize with nlp_D solution + initializer_option_value = "lp_A" + # if run all results or just sensitivity test + rerun_all_paper_results_value = False + + + sol_set, obj_set = kinetics_experiment(mip_option_value, + objective_value, + small_element=small_element_value, + file_store_name=file_store_name_value, + initializer_option=initializer_option_value, + rerun_all_paper_results=rerun_all_paper_results_value, + linear_solver_opt=linear_solver_opt_value) + + + self.assertAlmostEqual(np.trace(obj_set[0]), 118.92, places=1) + self.assertAlmostEqual(np.linalg.det(obj_set[0]), 0.1275, places=1) + self.assertAlmostEqual(np.trace(obj_set[1]), 168.459, places=1) + self.assertAlmostEqual(np.linalg.det(obj_set[1]), 3.539, places=1) + + + ### Test D-optimality, NLP + + ## The most often used options are listed as arguments here + mip_option_value = False # mixed integer problem or not + objective_value = ObjectiveLib.D # objective, ObjectiveLib.A or ObjectiveLib.D + small_element_value = 0.0001 # the small element added to the diagonal of FIM + # initialize with A-opt. MILP solutions + # choose what solutions to initialize from: + # minlp_D: initialize with minlp_D solutions + # milp_A: initialize with milp_A solutions + # lp_A: iniitalize with lp_A solution + # nlp_D: initialize with nlp_D solution + initializer_option_value = "milp_A" + # if run all results or just sensitivity test + rerun_all_paper_results_value = False + + + ### Test D-optimality + sol_set, obj_set = kinetics_experiment(mip_option_value, + objective_value, + small_element=small_element_value, + file_store_name=file_store_name_value, + initializer_option=initializer_option_value, + rerun_all_paper_results=rerun_all_paper_results_value, + linear_solver_opt=linear_solver_opt_value) + + + self.assertAlmostEqual(np.trace(obj_set[0]), 117.62, places=1) + self.assertAlmostEqual(np.linalg.det(obj_set[0]), 12.11, places=1) + self.assertAlmostEqual(np.trace(obj_set[1]), 139.544, places=1) + self.assertAlmostEqual(np.linalg.det(obj_set[1]), 21.5379, places=1) + + ### Test D-optimality, MINLP + + ## The most often used options are listed as arguments here + mip_option_value = True # mixed integer problem or not + objective_value = ObjectiveLib.D # objective, ObjectiveLib.A or ObjectiveLib.D + small_element_value = 0.0001 # the small element added to the diagonal of FIM + # initialize with A-opt. MILP solutions + # choose what solutions to initialize from: + # minlp_D: initialize with minlp_D solutions + # milp_A: initialize with milp_A solutions + # lp_A: iniitalize with lp_A solution + # nlp_D: initialize with nlp_D solution + initializer_option_value = "milp_A" + # if run all results or just sensitivity test + rerun_all_paper_results_value = False + + + ### Test D-optimality + sol_set, obj_set = kinetics_experiment(mip_option_value, + objective_value, + small_element=small_element_value, + file_store_name=file_store_name_value, + initializer_option=initializer_option_value, + rerun_all_paper_results=rerun_all_paper_results_value, + linear_solver_opt=linear_solver_opt_value) + + + self.assertAlmostEqual(np.trace(obj_set[0]), 110.1356, places=1) + self.assertAlmostEqual(np.linalg.det(obj_set[0]), 6.049, places=1) + self.assertAlmostEqual(np.trace(obj_set[1]), 124.0202, places=1) + self.assertAlmostEqual(np.linalg.det(obj_set[1]), 18.2956, places=1) + + if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/test_function.py b/test_function.py new file mode 100644 index 0000000..4d6e99f --- /dev/null +++ b/test_function.py @@ -0,0 +1,427 @@ +import numpy as np +import pandas as pd +import pyomo.environ as pyo +from measure_optimize import ( + MeasurementOptimizer, + SensitivityData, + MeasurementData, + CovarianceStructure, + ObjectiveLib, +) +import pickle +import time +import unittest + +class TestSensitivity(unittest.TestCase): + + def test_sens_data(self): + + ### Read and create Jacobian object + + + # create data object to pre-compute Qs + # read jacobian from the source csv + # Nt is the number of time points for each measurement + # This csv file looks like this: + # A1 A2 E1 E2 + #Unnamed: 0 + #1 -1.346695 -1.665335e-13 2.223380 -1.665335e-13 + # .... + + Nt = 8 + + jac_info = SensitivityData("./kinetics_source_data/Q_drop0.csv", Nt) + static_measurement_index = [ + 0, + 1, + 2, + ] # the index of CA, CB, CC in the jacobian array, considered as SCM + dynamic_measurement_index = [ + 0, + 1, + 2, + ] # the index of CA, CB, CC in the jacobian array, also considered as DCM + jac_info.get_jac_list( + static_measurement_index, # the index of SCMs in the jacobian array + dynamic_measurement_index, + ) # the index of DCMs in the jacobian array + + # test the shape + # 0,0 are both SCMs, shape should be Nt*Nt which is 8*8. same to 2,2 + self.assertTrue(np.shape(jac_info.jac)==(48,4)) + self.assertTrue(jac_info.total_measure_idx-6==0) + + # test if values are right + self.assertAlmostEqual(jac_info.jac[0,0], -1.346, places=2) + self.assertAlmostEqual(jac_info.jac[1,0], -1.039, places=2) + self.assertAlmostEqual(jac_info.jac[26,2], 1.6379, places=2) + self.assertAlmostEqual(jac_info.jac[47,2], -1.568, places=2) + self.assertAlmostEqual(jac_info.jac[47,3], -6.308, places=2) + + +class TestMeasurementError(unittest.TestCase): + + def test_measurement_error(self): + ### STEP 1: set up measurement cost strategy + + # number of time points for DCM + Nt = 8 + + # maximum manual measurement number for each measurement + max_manual_num = 10 + # minimal measurement interval + min_interval_num = 10.0 + # maximum manual measurement number for all measurements + total_max_manual_num = 10 + + # index of columns of SCM and DCM in Q + # SCM measurements: # CA # CB # CC. This is the index of measurement column in the reformulated Jacobian + static_ind = [0, 1, 2] + # DCM measurements: # CA # CB # CC. This is the index of measurement column in the reformulated Jacobian + dynamic_ind = [3, 4, 5] + # this index is the number of SCM + nubmer of DCM, not number of DCM timepoints + all_ind = static_ind + dynamic_ind + num_total_measure = len(all_ind) + + # meausrement names + all_names_strategy3 = [ + "CA.static", + "CB.static", + "CC.static", + "CA.dynamic", + "CB.dynamic", + "CC.dynamic", + ] + + # define static costs in $ + # CA # CB # CC # CA # CB # CC + static_cost = [2000, 2000, 2000, 200, 200, 200] + + # each static-cost measure has no per-sample cost + dynamic_cost = [0] * len(static_ind) + # each dynamic-cost measure costs $ 400 per sample + dynamic_cost.extend([400] * len(dynamic_ind)) + + ## define MeasurementData object + with self.assertRaises(ValueError): + measure_info = MeasurementData( + all_names_strategy3, # name string + num_total_measure, # jac_index: measurement index in Q + static_cost, # static costs + dynamic_cost, # dynamic costs + min_interval_num, # minimal time interval between two timepoints + max_manual_num, # maximum number of timepoints for each measurement + total_max_manual_num, # maximum number of timepoints for all measurement + ) + + ## define MeasurementData object + with self.assertRaises(ValueError): + measure_info = MeasurementData( + all_names_strategy3, # name string + num_total_measure, # jac_index: measurement index in Q + static_cost[0], # static costs + dynamic_cost, # dynamic costs + min_interval_num, # minimal time interval between two timepoints + max_manual_num, # maximum number of timepoints for each measurement + total_max_manual_num, # maximum number of timepoints for all measurement + ) + + ## define MeasurementData object + with self.assertRaises(ValueError): + measure_info = MeasurementData( + all_names_strategy3, # name string + num_total_measure, # jac_index: measurement index in Q + static_cost, # static costs + dynamic_cost[0], # dynamic costs + min_interval_num, # minimal time interval between two timepoints + max_manual_num, # maximum number of timepoints for each measurement + total_max_manual_num, # maximum number of timepoints for all measurement + ) + + ## define MeasurementData object + with self.assertRaises(ValueError): + measure_info = MeasurementData( + all_names_strategy3, # name string + num_total_measure, # jac_index: measurement index in Q + static_cost, # static costs + dynamic_cost, # dynamic costs + "10", # minimal time interval between two timepoints + max_manual_num, # maximum number of timepoints for each measurement + total_max_manual_num, # maximum number of timepoints for all measurement + ) + + ## define MeasurementData object + with self.assertRaises(ValueError): + measure_info = MeasurementData( + all_names_strategy3, # name string + num_total_measure, # jac_index: measurement index in Q + static_cost, # static costs + dynamic_cost[0], # dynamic costs + min_interval_num, # minimal time interval between two timepoints + "10", # maximum number of timepoints for each measurement + total_max_manual_num, # maximum number of timepoints for all measurement + ) + + ## define MeasurementData object + with self.assertRaises(ValueError): + measure_info = MeasurementData( + all_names_strategy3, # name string + num_total_measure, # jac_index: measurement index in Q + static_cost, # static costs + dynamic_cost[0], # dynamic costs + min_interval_num, # minimal time interval between two timepoints + max_manual_num, # maximum number of timepoints for each measurement + "20", # maximum number of timepoints for all measurement + ) + + +class TestCovariance(unittest.TestCase): + + def test_covariance(self): + + ### STEP 1: set up measurement cost strategy + + # number of time points for DCM + Nt = 8 + + # maximum manual measurement number for each measurement + max_manual_num = 10 + # minimal measurement interval + min_interval_num = 10.0 + # maximum manual measurement number for all measurements + total_max_manual_num = 10 + + # index of columns of SCM and DCM in Q + # SCM measurements: # CA # CB # CC. This is the index of measurement column in the reformulated Jacobian + static_ind = [0, 1, 2] + # DCM measurements: # CA # CB # CC. This is the index of measurement column in the reformulated Jacobian + dynamic_ind = [3, 4, 5] + # this index is the number of SCM + nubmer of DCM, not number of DCM timepoints + all_ind = static_ind + dynamic_ind + num_total_measure = len(all_ind) + + # meausrement names + all_names_strategy3 = [ + "CA.static", + "CB.static", + "CC.static", + "CA.dynamic", + "CB.dynamic", + "CC.dynamic", + ] + + # define static costs in $ + # CA # CB # CC # CA # CB # CC + static_cost = [2000, 2000, 2000, 200, 200, 200] + + # each static-cost measure has no per-sample cost + dynamic_cost = [0] * len(static_ind) + # each dynamic-cost measure costs $ 400 per sample + dynamic_cost.extend([400] * len(dynamic_ind)) + + ## define MeasurementData object + measure_info = MeasurementData( + all_names_strategy3, # name string + all_ind, # jac_index: measurement index in Q + static_cost, # static costs + dynamic_cost, # dynamic costs + min_interval_num, # minimal time interval between two timepoints + max_manual_num, # maximum number of timepoints for each measurement + total_max_manual_num, # maximum number of timepoints for all measurement + ) + + + ### STEP 2: Read and create Jacobian object + + + # create data object to pre-compute Qs + # read jacobian from the source csv + # Nt is the number of time points for each measurement + # This csv file looks like this: + # A1 A2 E1 E2 + #Unnamed: 0 + #1 -1.346695 -1.665335e-13 2.223380 -1.665335e-13 + # .... + + jac_info = SensitivityData("./kinetics_source_data/Q_drop0.csv", Nt) + static_measurement_index = [ + 0, + 1, + 2, + ] # the index of CA, CB, CC in the jacobian array, considered as SCM + dynamic_measurement_index = [ + 0, + 1, + 2, + ] # the index of CA, CB, CC in the jacobian array, also considered as DCM + jac_info.get_jac_list( + static_measurement_index, # the index of SCMs in the jacobian array + dynamic_measurement_index, + ) # the index of DCMs in the jacobian array + + + ### Test CovarianceStructure.identity (default) + + # use MeasurementOptimizer to pre-compute the unit FIMs + # do not define any error covariance structure, it should defaultly be identity matrix + calculator = MeasurementOptimizer( + jac_info, # SensitivityData object + measure_info, # MeasurementData object + print_level=0, # I use highest here to see more information + ) + + # test the shape + # 0,0 are both SCMs, shape should be Nt*Nt which is 8*8. same to 2,2 + self.assertTrue(np.shape(calculator.Sigma_inv[(0,0)])==(8,8)) + self.assertTrue(np.shape(calculator.Sigma_inv[(2,2)])==(8,8)) + # 2,3 are SCM - DCM combination, shape should be Nt*1 which is 8*1. same to 2, 26 + self.assertTrue(np.shape(calculator.Sigma_inv[(2,3)])==(8,1)) + self.assertTrue(np.shape(calculator.Sigma_inv[(2,26)])==(8,1)) + + # test if values are right + self.assertAlmostEqual(calculator.Sigma_inv[(2,2)][0,0], 1.0, places=2) + self.assertAlmostEqual(calculator.Sigma_inv[(3,3)], 1.0, places=1) + self.assertAlmostEqual(calculator.Sigma_inv[(26,26)], 1.0, places=2) + self.assertAlmostEqual(calculator.Sigma_inv[(2,3)][0,0], 0.0, places=3) + self.assertAlmostEqual(calculator.Sigma_inv[(2,26)][7,0], 0.0, places=3) + + + ### Test CovarianceStructure.variance + + # initialize the variance of each measurement CA, CB, CC + error_cov_initial = [1, 4, 8, 1, 4, 8] + error_cov = [] + # for each time point in each measurement, use the variance + for i in error_cov_initial: + error_cov.extend([i]*Nt) + + # use MeasurementOptimizer to pre-compute the unit FIMs + calculator = MeasurementOptimizer( + jac_info, # SensitivityData object + measure_info, # MeasurementData object + error_cov=error_cov, # error covariance matrix + error_opt=CovarianceStructure.variance, # error covariance options + print_level=0, # I use highest here to see more information + ) + + # test the shape + # 0,0 are both SCMs, shape should be Nt*Nt which is 8*8. same to 2,2 + self.assertTrue(np.shape(calculator.Sigma_inv[(0,0)])==(8,8)) + self.assertTrue(np.shape(calculator.Sigma_inv[(2,2)])==(8,8)) + # 2,3 are SCM - DCM combination, shape should be Nt*1 which is 8*1. same to 2, 26 + self.assertTrue(np.shape(calculator.Sigma_inv[(2,3)])==(8,1)) + self.assertTrue(np.shape(calculator.Sigma_inv[(2,26)])==(8,1)) + + # test if values are right + self.assertAlmostEqual(calculator.Sigma_inv[(2,2)][0,0], 0.125, places=2) + self.assertAlmostEqual(calculator.Sigma_inv[(3,3)], 1.0, places=1) + self.assertAlmostEqual(calculator.Sigma_inv[(26,26)], 0.125, places=2) + self.assertAlmostEqual(calculator.Sigma_inv[(2,3)][0,0], 0.0, places=3) + self.assertAlmostEqual(calculator.Sigma_inv[(2,26)][7,0], 0.0, places=3) + + ### Test CovarianceStructure.time_correlation + + # create the time-correlation matrix for one measurement, shape Nt*Nt + error_cov_initial = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8] + # create such matrix to be a diagonal matrix + error_cov_one_measure = [[0]*Nt for i in range(Nt)] + for i in range(Nt): + error_cov_one_measure[i][i] = error_cov_initial[i] + # for each measurement, use the same time-correlation matrix + error_cov = [] + for i in range(len(all_ind)): + error_cov.append(error_cov_one_measure) + + # use MeasurementOptimizer to pre-compute the unit FIMs + calculator = MeasurementOptimizer( + jac_info, # SensitivityData object + measure_info, # MeasurementData object + error_cov=error_cov, # error covariance matrix + error_opt=CovarianceStructure.time_correlation, # error covariance options + print_level=0, # I use highest here to see more information + ) + + # test the shape + # 0,0 are both SCMs, shape should be Nt*Nt which is 8*8. same to 2,2 + self.assertTrue(np.shape(calculator.Sigma_inv[(0,0)])==(8,8)) + self.assertTrue(np.shape(calculator.Sigma_inv[(2,2)])==(8,8)) + # 2,3 are SCM - DCM combination, shape should be Nt*1 which is 8*1. same to 2, 26 + self.assertTrue(np.shape(calculator.Sigma_inv[(2,3)])==(8,1)) + self.assertTrue(np.shape(calculator.Sigma_inv[(2,26)])==(8,1)) + + # test if values are right + self.assertAlmostEqual(calculator.Sigma_inv[(2,2)][0,0], 10.0, places=2) + self.assertAlmostEqual(calculator.Sigma_inv[(3,3)], 10.0, places=1) + self.assertAlmostEqual(calculator.Sigma_inv[(26,26)], 1.25, places=2) + self.assertAlmostEqual(calculator.Sigma_inv[(2,3)][0,0], 0.0, places=3) + self.assertAlmostEqual(calculator.Sigma_inv[(2,26)][7,0], 0.0, places=3) + + + ### Test CovarianceStructure.measure_correlation + + # error covariance matrix + error_cov = [[1.0, 0.1, 0.1, 0.5, 0.05, 0.05], + [0.1, 4.0, 0.5, 0.05, 2.0, 0.25], + [0.1, 0.5, 8.0, 0.05, 0.25, 4.0], + [0.5, 0.05, 0.05, 1.0, 0.1, 0.1], + [0.05, 2.0, 0.25, 0.1, 4.0, 0.5], + [0.05, 0.25, 4.0, 0.1, 0.5, 8.0]] + + + # use MeasurementOptimizer to pre-compute the unit FIMs + calculator = MeasurementOptimizer( + jac_info, # SensitivityData object + measure_info, # MeasurementData object + error_cov=error_cov, # error covariance matrix + error_opt=CovarianceStructure.measure_correlation, # error covariance options + print_level=0, # I use highest here to see more information + ) + + # test the shape + # 0,0 are both SCMs, shape should be Nt*Nt which is 8*8. same to 2,2 + self.assertTrue(np.shape(calculator.Sigma_inv[(0,0)])==(8,8)) + self.assertTrue(np.shape(calculator.Sigma_inv[(2,2)])==(8,8)) + # 2,3 are SCM - DCM combination, shape should be Nt*1 which is 8*1. same to 2, 26 + self.assertTrue(np.shape(calculator.Sigma_inv[(2,3)])==(8,1)) + self.assertTrue(np.shape(calculator.Sigma_inv[(2,26)])==(8,1)) + + # test if values are right + self.assertAlmostEqual(calculator.Sigma_inv[(2,2)][0,0], 0.168, places=2) + self.assertAlmostEqual(calculator.Sigma_inv[(3,3)], 1.3379, places=1) + self.assertAlmostEqual(calculator.Sigma_inv[(26,26)], 0.168, places=2) + self.assertAlmostEqual(calculator.Sigma_inv[(2,3)][0,0], 0.00737, places=3) + self.assertAlmostEqual(calculator.Sigma_inv[(2,26)][7,0], -0.08407, places=3) + + + ### Test CovarianceStructure.time_measure_correlation + + # error covariance matrix should be a sum(Nt)*sum(Nt) matrix + error_cov = [[1]*calculator.total_num_time for i in range(calculator.total_num_time)] + + # use MeasurementOptimizer to pre-compute the unit FIMs + calculator = MeasurementOptimizer( + jac_info, # SensitivityData object + measure_info, # MeasurementData object + error_cov=error_cov, # error covariance matrix + error_opt=CovarianceStructure.time_measure_correlation, # error covariance options + print_level=0, # I use highest here to see more information + ) + + # test the shape + # 0,0 are both SCMs, shape should be Nt*Nt which is 8*8. same to 2,2 + self.assertTrue(np.shape(calculator.Sigma_inv[(0,0)])==(8,8)) + self.assertTrue(np.shape(calculator.Sigma_inv[(2,2)])==(8,8)) + # 2,3 are SCM - DCM combination, shape should be Nt*1 which is 8*1. same to 2, 26 + self.assertTrue(np.shape(calculator.Sigma_inv[(2,3)])==(8,1)) + self.assertTrue(np.shape(calculator.Sigma_inv[(2,26)])==(8,1)) + + # test if values are right + self.assertAlmostEqual(calculator.Sigma_inv[(2,2)][0,0], 0.000434, places=4) + self.assertAlmostEqual(calculator.Sigma_inv[(3,3)], 0.000434, places=4) + self.assertAlmostEqual(calculator.Sigma_inv[(26,26)], 0.000434, places=4) + self.assertAlmostEqual(calculator.Sigma_inv[(2,3)][0,0], 0.000434, places=4) + self.assertAlmostEqual(calculator.Sigma_inv[(2,26)][7,0], 0.000434, places=4) + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file From 57e654d69f7e3804749b6217d147d1a3299aca3d Mon Sep 17 00:00:00 2001 From: jialuw96 Date: Sun, 21 Apr 2024 21:56:51 -0700 Subject: [PATCH 3/7] add test --- test_cvxpy.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test_cvxpy.py b/test_cvxpy.py index 781c633..967ca28 100644 --- a/test_cvxpy.py +++ b/test_cvxpy.py @@ -169,4 +169,6 @@ def test_cv(self): num_dynamic_t_name=num_dynamic_time, # number of time points of DCMs solver="MOSEK", store_name = file_store_name) + + From 78bbd425eea666902b46ffd20b052dbbaa6dca92 Mon Sep 17 00:00:00 2001 From: jialuw96 Date: Sun, 21 Apr 2024 22:15:45 -0700 Subject: [PATCH 4/7] add cvxpy tests --- measure_optimize.py | 255 ++++++++++++++++++++++++++++++++++++++------ test_cvxpy.py | 39 ++++++- 2 files changed, 257 insertions(+), 37 deletions(-) diff --git a/measure_optimize.py b/measure_optimize.py index 511bc12..832f75a 100644 --- a/measure_optimize.py +++ b/measure_optimize.py @@ -16,6 +16,7 @@ from dataclasses import dataclass import pickle from scipy.spatial import cKDTree +import cvxpy as cp class CovarianceStructure(Enum): @@ -54,6 +55,7 @@ class ObjectiveLib(Enum): A = 0 D = 1 + E = 2 @classmethod def has_value(cls, value): @@ -931,18 +933,18 @@ def _measure_matrix(self, measurement_vector): measurement_matrix: a full measurement matrix, construct the weights for covariances """ # check if measurement vector legal - if len(measurement_vector) != self.total_no_measure: + if len(measurement_vector) != self.num_measure_dynamic_flatten: raise ValueError( "Measurement vector is of wrong shape, expecting length of" - + str(self.total_no_measure) + + str(self.num_measure_dynamic_flatten) ) # initialize measurement matrix as a 2D array - measurement_matrix = np.zeros((self.total_no_measure, self.total_no_measure)) + measurement_matrix = np.zeros((self.num_measure_dynamic_flatten, self.num_measure_dynamic_flatten)) # loop over total measurement index - for i in range(self.total_no_measure): - for j in range(self.total_no_measure): + for i in range(self.num_measure_dynamic_flatten): + for j in range(self.num_measure_dynamic_flatten): measurement_matrix[i, j] = min( measurement_vector[i], measurement_vector[j] ) @@ -1387,7 +1389,7 @@ def dynamic_install_init(m, j): def dynamic_manual_num(m): """the timepoints for each dynamical measurement should be smaller than a given limit - Eq. 11f in paper + Eq. 11g in paper """ start = ( self.n_static_measurements + i * self.dynamic_Nt @@ -1608,12 +1610,12 @@ def compute_fim(self, measurement_vector): fim = np.zeros((self.no_param, self.no_param)) # go over all measurement index - for m1_rank in range(self.total_no_measure): + for m1_rank in range(self.num_measure_dynamic_flatten): # get the corresponding gradient vector for this measurement # use np.matrix to keep the dimensions for the vector, otherwise numpy makes the vector 1D instead of Np*1 jac_m1 = np.matrix(self.jac[m1_rank]) # go over all measurement index - for m2_rank in range(self.total_no_measure): + for m2_rank in range(self.num_measure_dynamic_flatten): # get the corresponding gradient vector for this measurement jac_m2 = np.matrix(self.jac[m2_rank]) # solution of if these two measurements are selected @@ -1983,9 +1985,41 @@ def extract_store_sol(self, budget_opt, store_name): pickle.dump(fim_result, file2) file2.close() - # return the answer y, and the FIM return ans_y, fim_result + def extract_store_sol_cvxpy(self, budget_opt, y_matrice, obj, fim, if_install_y, store_name): + """ + Extract the solution, FIM, objective function from Cvxpy model, and store in two pickle files + "store_name + budget_opt": pickle file storing the y solution + "store_name + fim + budget_opt": pickle file storing the FIM + + Arguments + --------- + :param budget_opt: budget of the problem + :param y_matrice: the Cvxpy y matrice solution + :param obj: the Cvxpy calculated objective function + :param fim: the FIM calculated by Cvxpy + :param if_install_y: DCM installaion binary variable solution + :param store_name: if not None, store the solution and FIM in pickle file with the given name + """ + + # extract solution + ans_y, sol_y = self.extract_solutions(y_matrice=y_matrice, if_install_y=if_install_y) + + # extract FIM + if self.precompute_print_level >= 2: + print_fim(fim) + print("CVXPY calculated objective value:", obj.value) + + if store_name: + file = open(store_name + str(budget_opt), "wb") + pickle.dump(ans_y, file) + file.close() + + file2 = open(store_name + "fim_" + str(budget_opt), "wb") + pickle.dump(fim, file2) + file2.close() + def _locate_initial_file(self, budget): """ Given the budget, select which initial solution it should use by locating the closest budget that has stored solutions @@ -2260,34 +2294,51 @@ def cost_computation(m): print("warmstart initialize total dynamic: ", total_dynamic_initial) print("warmstart initialize cost:", cost_init) - def continuous_optimization_cvxpy(self, objective="D", budget=100, solver=None): + def continuous_optimization_cvxpy(self, objective=ObjectiveLib.D, + budget=3000, + mixed_integer = False, + static_dynamic_pair = None, + time_interval_all_dynamic = False, + num_dynamic_t_name = None, + surrogate_obj = True, + solver=None, + store_name = None): """ This optimization problem can also be formulated and solved in the CVXPY framework. - This is a generalization code for CVXPY problems for reference, not currently used for the paper. + This is a generalization code for CVXPY problems for Eq. 11 in MO paper. Arguments --------- :param objective: can choose from 'D', 'A', 'E' for now. if defined others or None, use A-optimality. - :param cost_budget: give a total limit for costs. + :param budget: give a total limit for costs. + :param mixed_integer: if the problem is a mixed-integer problem, or relaxed + :param static_dynamic_pair: list of lists + a list of the name of measurements, that are selected as either dynamic or static measurements. + :param time_interval_all_dynamic: boolean + if True, the minimal time interval applies for all dynamical measurements + :param num_dynamic_t_name: list + a list of the exact time points for the dynamic-cost measurements time points + :param surrogate_obj: boolean + if True, using a surrogate variable to be the objective, the objective function < surrogate variable + if False, directly use the objective function to be the objective :param solver: default to be MOSEK. Look for CVXPY document for more solver information. + :param store_name: string + name used to store the solution and objective function Returns ------- None """ - # compute Atomic FIM - self.assemble_unit_fims() - - # evaluate fim in the CVXPY framework + # evaluate fim in the CVXPY framework, Eq.(11a) in MO paper def eval_fim(y): """Evaluate FIM from y solution FIM = sum(cov_y[i,j]*unit FIM[i,j]) for all i, j in n_responses """ fim = sum( - y[i, j] * self.unit_fims[i * self.total_no_measure + j] - for i in range(self.total_no_measure) - for j in range(self.total_no_measure) + y[i, j] * self.unit_fims[i * self.num_measure_dynamic_flatten + j] + for i in range(self.num_measure_dynamic_flatten) + for j in range(self.num_measure_dynamic_flatten) ) return fim @@ -2302,40 +2353,167 @@ def d_opt(y): return cp.log_det(fim) def e_opt(y): - """E-optimality as OBJ""" + """E-optimality as OBJ + This OBJ is not used for now (pending test) + """ fim = eval_fim(y) return -cp.lambda_min(fim) # construct variables - y_matrice = cp.Variable( - (self.total_no_measure, self.total_no_measure), nonneg=True - ) + if not mixed_integer: # continuous variables + y_matrice = cp.Variable( + (self.num_measure_dynamic_flatten, self.num_measure_dynamic_flatten), nonneg=True + ) + else: # discrete binary variables + y_matrice = cp.Variable( + (self.num_measure_dynamic_flatten, self.num_measure_dynamic_flatten), boolean=True + ) + + # boolean for if a DCM is installed + if_install_y = cp.Variable(self.n_dynamic_measurements, nonneg = True) + + # surrogate objective + surro_obj = cp.Variable(nonneg = True) # cost limit p_cons = [ - sum(y_matrice[i, i] * self.cost[i] for i in range(self.total_no_measure)) + sum(y_matrice[i, i] * self.cost_list[i] for i in range(self.num_measure_dynamic_flatten)) # static and dynamic timepoint cost + + sum(if_install_y[i]*self.dynamic_install_cost[i] for i in range(self.n_dynamic_measurements)) # dynamic installation costs <= budget ] # loop over all measurement index - for k in range(self.total_no_measure): + for k in range(self.num_measure_dynamic_flatten): # loop over all measurement index - for l in range(self.total_no_measure): - # y[k,l] = y[k]*y[l] relaxation + for l in range(self.num_measure_dynamic_flatten): + # y[k,l] = y[k]*y[l] McCormick relaxation p_cons += [y_matrice[k, l] <= y_matrice[k, k]] p_cons += [y_matrice[k, l] <= y_matrice[l, l]] p_cons += [y_matrice[k, k] + y_matrice[l, l] - 1 <= y_matrice[k, l]] + # symmetric matrix p_cons += [y_matrice.T == y_matrice] + #total number of manual dynamical measurements number + #Eq. 11f in paper + p_cons += [sum(y_matrice[i,i] for i in range(self.n_static_measurements, self.sens_info.total_measure_idx)) + <= self.manual_number] + + # loop over dynamical measurements + #the timepoints for each dynamical measurement should be smaller than a given limit + #Eq. 11g in paper + for i in range(self.n_dynamic_measurements): + start = ( + self.n_static_measurements + i * self.sens_info.Nt + ) # the start index of this dynamical measurement + end = ( + self.n_static_measurements + (i + 1) * self.sens_info.Nt + ) # the end index of this dynamical measurement + # all time points for this DCM should be smaller than L_dynamic + p_cons += [sum(y_matrice[j, j] for j in range(start, end)) <= self.each_manual_number] + + # if install a DCM + for i in range(self.n_dynamic_measurements): + start = ( + self.n_static_measurements + i * self.sens_info.Nt + ) # the start index of this dynamical measurement + end = ( + self.n_static_measurements + (i + 1) * self.sens_info.Nt + ) # the end index of this dynamical measurement + # all time points for this DCM should be smaller than L_dynamic + for j in range(start, end): + p_cons += [if_install_y[i] >= y_matrice[j,j]] + + # constrains for if installing DCM + p_cons += [if_install_y[i] <= sum(y_matrice[k,k] for k in range(start, end))] + p_cons += [if_install_y[i] <= 1 ] + + # if some measurements can only be dynamic or static + if static_dynamic_pair is not None: + # loop over the index of the static, and dynamic measurements + for i, pair in enumerate(static_dynamic_pair): + # if_install_y needs a padding of self.n_static_measurements + p_cons += [if_install_y[pair[1]-self.n_static_measurements] + y_matrice[pair[0], pair[0]] <= 1] + + # pair time index and real time + # for e.g., time 2h is the index 16, dynamic[16] = 2 + # this is for the time interval between two DCMs computation + dynamic_time = {} + # loop over time index + for i in range(self.sens_info.Nt): + # index: real time + dynamic_time[i] = num_dynamic_t_name[i] + #self.dynamic_time = dynamic_time + + # if there is minimal interval constraint + # Eq. 11h in the paper + if self.min_time_interval is not None: + # if this constraint applies to all dynamic measurements + if time_interval_all_dynamic: + for t in range(self.sens_info.Nt): + # end time is an open end of the region, so another constraint needs to be added to include end_time + # if dynamic_time[t]+discretize_time <= end_time+0.1*discretize_time: + sumi = 0 + count = 0 + # get the timepoints in this interval + while (count + t < self.sens_info.Nt) and ( + dynamic_time[count + t] - dynamic_time[t] + ) < self.min_time_interval: + for i in range(self.n_dynamic_measurements): + surro_idx = ( + self.n_static_measurements + + i + * self.sens_info.Nt + + t + + count + ) + sumi += y_matrice[surro_idx, surro_idx] + count += 1 + + p_cons += [sumi <= 1] + + # if this constraint applies to each dynamic measurements, in a local way + else: + for i in range(self.n_dynamic_measurements): + for t in range(self.sens_info.Nt): + # end time is an open end of the region, so another constraint needs to be added to include end_time + # if dynamic_time[t]+discretize_time <= end_time+0.1*discretize_time: + + # sumi is the summation of all measurements selected during this time interval + sumi = 0 + # count helps us go through each time points in this time interval + count = 0 + # get timepoints in this interval + while (count + t < self.sens_info.Nt) and ( + dynamic_time[count + t] - dynamic_time[t] + ) < self.min_time_interval: + # surro_idx gets the index of the current time point + surro_idx = ( + self.n_static_measurements + + i + * self.sens_info.Nt + + t + + count + ) + # sum up all timepoints selections + sumi += y_matrice[surro_idx, surro_idx] + count += 1 + + p_cons += [ sumi <= 1 ] + + # D-optimality - if objective == "D": - obj = cp.Maximize(d_opt(y_matrice)) + if objective == ObjectiveLib.D: + if surrogate_obj: + p_cons += [surro_obj <= d_opt(y_matrice)] + obj = cp.Maximize(surro_obj) + else: + obj = cp.Maximize(d_opt(y_matrice)) # E-optimality - elif objective == "E": + elif objective == ObjectiveLib.E: obj = cp.Maximize(e_opt(y_matrice)) # A-optimality else: - if self.verbose: + if self.precompute_print_level >= 2: print("Use A-optimality (Trace).") obj = cp.Maximize(a_opt(y_matrice)) @@ -2346,9 +2524,11 @@ def e_opt(y): else: problem.solve(solver=solver, verbose=True) - self.__solution_analysis(y_matrice, obj.value) + self.extract_store_sol_cvxpy(budget,y_matrice, obj, eval_fim(y_matrice).value, if_install_y, store_name) + + return obj - def extract_solutions(self): + def extract_solutions(self, y_matrice=None, if_install_y=None): """ Extract and show solutions from a solved Pyomo model mod is an argument because we @@ -2371,9 +2551,14 @@ def extract_solutions(self): for i in range(self.num_measure_dynamic_flatten): # loop over the measurement choice index for j in range(i, self.num_measure_dynamic_flatten): - cov = pyo.value(self.mod.cov_y[i, j]) + if not y_matrice: + cov = pyo.value(self.mod.cov_y[i, j]) + else: + cov = y_matrice[i,j].value # give value to its symmetric part ans_y[i, j] = ans_y[j, i] = cov + + # round small errors to integers # loop over all measurement choice @@ -2404,7 +2589,7 @@ def extract_solutions(self): ] ) # group DCM time points - sol_y = np.reshape(sol_y, (self.n_dynamic_measurements, self.dynamic_Nt)) + sol_y = np.reshape(sol_y, (self.n_dynamic_measurements, self.sens_info.Nt)) np.around(sol_y) # loop over each DCM for r in range(len(sol_y)): diff --git a/test_cvxpy.py b/test_cvxpy.py index 967ca28..e66fa30 100644 --- a/test_cvxpy.py +++ b/test_cvxpy.py @@ -12,6 +12,7 @@ import pickle import time import matplotlib.pyplot as plt +import unittest #from idaes.core.util.model_diagnostics import DiagnosticsToolbox @@ -144,7 +145,7 @@ def test_cv(self): # optimization options objective = ObjectiveLib.A - budget_opt = 5000 + budget_opt = 3000 mixed_integer_opt = True @@ -161,7 +162,7 @@ def test_cv(self): for i, tim in enumerate(num_dynamic_time[1:]): dynamic_time_dict[i] = np.round(tim, decimals=2) - calculator.continuous_optimization_cvxpy(objective=objective, + obj = calculator.continuous_optimization_cvxpy(objective=objective, mixed_integer = mixed_integer_opt, budget=budget_opt, static_dynamic_pair = static_dynamic, @@ -170,5 +171,39 @@ def test_cv(self): solver="MOSEK", store_name = file_store_name) + self.assertAlmostEqual(obj.value, 108.3066, places=1) + + + # optimization options + objective = ObjectiveLib.D + + budget_opt = 3000 + + mixed_integer_opt = False + + file_store_name = "./cvxpy_results/cvxpy_A_" + #file_store_name = None + + num_dynamic_time = np.linspace(0, 60, 9) + + static_dynamic = [[0, 3], [1, 4], [2, 5]] # These pairs cannot be chosen simultaneously + time_interval_for_all = True + + # map the timepoint index to its real time + dynamic_time_dict = {} + for i, tim in enumerate(num_dynamic_time[1:]): + dynamic_time_dict[i] = np.round(tim, decimals=2) + + obj = calculator.continuous_optimization_cvxpy(objective=objective, + mixed_integer = mixed_integer_opt, + budget=budget_opt, + static_dynamic_pair = static_dynamic, + time_interval_all_dynamic = time_interval_for_all, + num_dynamic_t_name=num_dynamic_time, # number of time points of DCMs + solver="MOSEK", + store_name = file_store_name) + + self.assertAlmostEqual(obj.value, 1.591, places=1) + From 87e72a253566d48b0c97be7873ea2d1ad204130d Mon Sep 17 00:00:00 2001 From: jialuw96 Date: Sun, 21 Apr 2024 22:26:16 -0700 Subject: [PATCH 5/7] test are uploaded --- test_example.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/test_example.py b/test_example.py index 4090e02..d34f4ed 100644 --- a/test_example.py +++ b/test_example.py @@ -311,12 +311,12 @@ def kinetics_experiment(mip_option, initializer_option_value = "lp_A" # if run all results or just sensitivity test rerun_all_paper_results_value = False + linear_solver_opt_value = "ma57" sol_set, obj_set = kinetics_experiment(mip_option_value, objective_value, small_element=small_element_value, - file_store_name=file_store_name_value, initializer_option=initializer_option_value, rerun_all_paper_results=rerun_all_paper_results_value, linear_solver_opt=linear_solver_opt_value) @@ -348,7 +348,6 @@ def kinetics_experiment(mip_option, sol_set, obj_set = kinetics_experiment(mip_option_value, objective_value, small_element=small_element_value, - file_store_name=file_store_name_value, initializer_option=initializer_option_value, rerun_all_paper_results=rerun_all_paper_results_value, linear_solver_opt=linear_solver_opt_value) @@ -381,7 +380,6 @@ def kinetics_experiment(mip_option, sol_set, obj_set = kinetics_experiment(mip_option_value, objective_value, small_element=small_element_value, - file_store_name=file_store_name_value, initializer_option=initializer_option_value, rerun_all_paper_results=rerun_all_paper_results_value, linear_solver_opt=linear_solver_opt_value) @@ -413,7 +411,6 @@ def kinetics_experiment(mip_option, sol_set, obj_set = kinetics_experiment(mip_option_value, objective_value, small_element=small_element_value, - file_store_name=file_store_name_value, initializer_option=initializer_option_value, rerun_all_paper_results=rerun_all_paper_results_value, linear_solver_opt=linear_solver_opt_value) From 4de24bc3c80d3c6f897855b3f2c995b3e9ae16c9 Mon Sep 17 00:00:00 2001 From: jialuw96 Date: Sun, 21 Apr 2024 22:36:31 -0700 Subject: [PATCH 6/7] finish commenting --- test_cvxpy.py | 14 ++++++++++---- test_example.py | 15 ++++++++------- test_function.py | 26 +++++++++++++------------- 3 files changed, 31 insertions(+), 24 deletions(-) diff --git a/test_cvxpy.py b/test_cvxpy.py index e66fa30..64f9479 100644 --- a/test_cvxpy.py +++ b/test_cvxpy.py @@ -17,7 +17,9 @@ class TestCvxpy(unittest.TestCase): - + """ + Test CVXPY A-optimality MILP problem, and D-optimality NLP problem + """ def test_cv(self): # set up problem formulation @@ -142,13 +144,15 @@ def test_cv(self): # calculate a list of unit FIMs calculator.assemble_unit_fims() + + ### Test A-optimality MILP problem # optimization options objective = ObjectiveLib.A budget_opt = 3000 mixed_integer_opt = True - + file_store_name = "./cvxpy_results/cvxpy_A_" #file_store_name = None @@ -170,10 +174,11 @@ def test_cv(self): num_dynamic_t_name=num_dynamic_time, # number of time points of DCMs solver="MOSEK", store_name = file_store_name) - + # test objective function value self.assertAlmostEqual(obj.value, 108.3066, places=1) + ### Test D-optimality NLP problem # optimization options objective = ObjectiveLib.D @@ -181,7 +186,7 @@ def test_cv(self): mixed_integer_opt = False - file_store_name = "./cvxpy_results/cvxpy_A_" + file_store_name = "./cvxpy_results/cvxpy_D_" #file_store_name = None num_dynamic_time = np.linspace(0, 60, 9) @@ -203,6 +208,7 @@ def test_cv(self): solver="MOSEK", store_name = file_store_name) + # test objective function value self.assertAlmostEqual(obj.value, 1.591, places=1) diff --git a/test_example.py b/test_example.py index d34f4ed..3712a2d 100644 --- a/test_example.py +++ b/test_example.py @@ -15,7 +15,8 @@ class TestExample(unittest.TestCase): - + """Test kinetics example, A-optimality LP and MILP, D-optimality NLP and MINLP + """ def test_eg(self): ### STEP 1: set up measurement cost strategy @@ -313,7 +314,7 @@ def kinetics_experiment(mip_option, rerun_all_paper_results_value = False linear_solver_opt_value = "ma57" - + # create object sol_set, obj_set = kinetics_experiment(mip_option_value, objective_value, small_element=small_element_value, @@ -321,7 +322,7 @@ def kinetics_experiment(mip_option, rerun_all_paper_results=rerun_all_paper_results_value, linear_solver_opt=linear_solver_opt_value) - + # test objective functions in A- and D- for both budgets self.assertAlmostEqual(np.trace(obj_set[0]), 153.05, places=1) self.assertAlmostEqual(np.linalg.det(obj_set[0]), 0.383, places=1) self.assertAlmostEqual(np.trace(obj_set[1]), 168.659, places=1) @@ -344,7 +345,7 @@ def kinetics_experiment(mip_option, # if run all results or just sensitivity test rerun_all_paper_results_value = False - + # create example object sol_set, obj_set = kinetics_experiment(mip_option_value, objective_value, small_element=small_element_value, @@ -352,7 +353,7 @@ def kinetics_experiment(mip_option, rerun_all_paper_results=rerun_all_paper_results_value, linear_solver_opt=linear_solver_opt_value) - + # test objective functions in A- and D- for both budgets self.assertAlmostEqual(np.trace(obj_set[0]), 118.92, places=1) self.assertAlmostEqual(np.linalg.det(obj_set[0]), 0.1275, places=1) self.assertAlmostEqual(np.trace(obj_set[1]), 168.459, places=1) @@ -384,7 +385,7 @@ def kinetics_experiment(mip_option, rerun_all_paper_results=rerun_all_paper_results_value, linear_solver_opt=linear_solver_opt_value) - + # test objective functions in A- and D- for both budgets self.assertAlmostEqual(np.trace(obj_set[0]), 117.62, places=1) self.assertAlmostEqual(np.linalg.det(obj_set[0]), 12.11, places=1) self.assertAlmostEqual(np.trace(obj_set[1]), 139.544, places=1) @@ -415,7 +416,7 @@ def kinetics_experiment(mip_option, rerun_all_paper_results=rerun_all_paper_results_value, linear_solver_opt=linear_solver_opt_value) - + # test objective functions in A- and D- for both budgets self.assertAlmostEqual(np.trace(obj_set[0]), 110.1356, places=1) self.assertAlmostEqual(np.linalg.det(obj_set[0]), 6.049, places=1) self.assertAlmostEqual(np.trace(obj_set[1]), 124.0202, places=1) diff --git a/test_function.py b/test_function.py index 4d6e99f..e7083fa 100644 --- a/test_function.py +++ b/test_function.py @@ -13,12 +13,11 @@ import unittest class TestSensitivity(unittest.TestCase): - + """Test sensitivity object by checking the shape and values of Jacobian matrix generated + """ def test_sens_data(self): - ### Read and create Jacobian object - # create data object to pre-compute Qs # read jacobian from the source csv # Nt is the number of time points for each measurement @@ -60,10 +59,10 @@ def test_sens_data(self): class TestMeasurementError(unittest.TestCase): - + """Test the measurement object by checking if the code throws + error as expected when given a wrong type of input + """ def test_measurement_error(self): - ### STEP 1: set up measurement cost strategy - # number of time points for DCM Nt = 8 @@ -102,7 +101,7 @@ def test_measurement_error(self): # each dynamic-cost measure costs $ 400 per sample dynamic_cost.extend([400] * len(dynamic_ind)) - ## define MeasurementData object + ## define MeasurementData object with wrong jac_index number with self.assertRaises(ValueError): measure_info = MeasurementData( all_names_strategy3, # name string @@ -114,7 +113,7 @@ def test_measurement_error(self): total_max_manual_num, # maximum number of timepoints for all measurement ) - ## define MeasurementData object + ## define MeasurementData object with wrong type of static costs with self.assertRaises(ValueError): measure_info = MeasurementData( all_names_strategy3, # name string @@ -126,7 +125,7 @@ def test_measurement_error(self): total_max_manual_num, # maximum number of timepoints for all measurement ) - ## define MeasurementData object + ## define MeasurementData object with wrong type of dynamic costs with self.assertRaises(ValueError): measure_info = MeasurementData( all_names_strategy3, # name string @@ -138,7 +137,7 @@ def test_measurement_error(self): total_max_manual_num, # maximum number of timepoints for all measurement ) - ## define MeasurementData object + ## define MeasurementData object with wrong type of minimal time interval with self.assertRaises(ValueError): measure_info = MeasurementData( all_names_strategy3, # name string @@ -150,7 +149,7 @@ def test_measurement_error(self): total_max_manual_num, # maximum number of timepoints for all measurement ) - ## define MeasurementData object + ## define MeasurementData object with wrong type of max # of timepoints with self.assertRaises(ValueError): measure_info = MeasurementData( all_names_strategy3, # name string @@ -162,7 +161,7 @@ def test_measurement_error(self): total_max_manual_num, # maximum number of timepoints for all measurement ) - ## define MeasurementData object + ## define MeasurementData object with wrong type of max # of total timepoints with self.assertRaises(ValueError): measure_info = MeasurementData( all_names_strategy3, # name string @@ -176,7 +175,8 @@ def test_measurement_error(self): class TestCovariance(unittest.TestCase): - + """Test if all five types of covariances generate the error covariance matrix correctly + """ def test_covariance(self): ### STEP 1: set up measurement cost strategy From 8cf0772d53fd66753b0ffefcc925a4d4ebd4981a Mon Sep 17 00:00:00 2001 From: Jialu <51490104+jialuw96@users.noreply.github.com> Date: Sun, 21 Apr 2024 22:42:07 -0700 Subject: [PATCH 7/7] Update README.md --- README.md | 61 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/README.md b/README.md index 6ac7623..babcd35 100644 --- a/README.md +++ b/README.md @@ -58,6 +58,20 @@ The following instructions assume you have anaconda installed. We suggest create `conda install jupyter notebook` + +### (Optional) `Cvxpy` environment setting +- this is needed only for constructing and solving the problem with `cvxpy`. + + Step 1: same as step 1 above to create the environment + + Step 2: same as step 3 above to install dependencies + + Step 3: install cvxpy with: + + `conda install -c conda-forge cvxpy` + + Step 4: install `Mosek`. You need to validate a liscence for using this solver, see the link: https://docs.mosek.com/latest/install/installation.html + ### Software versions we use for the results `Python`: 3.8 @@ -70,6 +84,8 @@ The following instructions assume you have anaconda installed. We suggest create `CyIpopt`: 1.3.0 +`Cvxpy`: 1.2.1 + ## Code - `measure_optimize.py`: Measurement optimization optimization framework @@ -82,6 +98,21 @@ The following instructions assume you have anaconda installed. We suggest create - `draw_figure.ipynb`: Generates all results figures in the manuscript +- `cvxpy_problem.py`: Kinetics case study implemented with `Cvxpy` + +## Unit tests contents and instructions + +- `test_example.py`: test the kinetics example in `Pyomo` in A- and D-optimality, mixed integer and relaxed +- `test_cvxpy.py`: (`Cvxpy` and `Mosek` needed, see above for instructions) test the kinetics example in `Cvxpy` in A-optimality MILP and D-optimality NLP problems +- `test_function.py`: test the `SensitivityData`, `MeasurementData` objects, and all five types of covariance matrix + +To run unittests, run the following command in the environment created above: + +`python -m unittest test_example` + +`python -m unittest test_cvxpy` + +`python -m unittest test_function` ## Example to run code and reproduce figures for case studies @@ -129,7 +160,37 @@ Setup the scripts to reproduce result files and figures from the paper: - `plot_one_solution` receives and draws the solution of one measurement under four strategies. To reproduce result figure like Fig. S-2 in paper, call it 6 times to draw all 6 figures and combine to a panel figure. + +### Kinetics case study with Cvxpy + +- Step 1: run `cvxpy_problem.py` +- Step 2: with `mip_option` and `objective`, choose to run the A-optimality or D-optimality, mixed-integer or relaxed problem +- Step 3: with `test`, set up the budget ranges as you want to try. + + If `test` is False: we use the budget range [1000, 5000] with a 400 discretization, + i.e. [1000, 1400, 1800, ..., 5000] for mixed-integer problems. + + Otherwise, we use three budget [3000, 5000] to do a test run. + +- Step 4: store results for drawing figures + + To do this, define the param `file_store_name` with a string you given, for e.g., "MINLP_result_". + + Then both the solutions and the FIM of the results are stored separately. + + For e.g., if running in the range [1000, 5000], the stored files will be: + + MINLP_result_1000, MINLP_result_fim_1000, + + ... + + MINLP_result_5000, MINLP_result_fim_5000, +- Step 7: use draw_figure.ipynb to read stored FIM and solutions + + - `read_fim` receives the string name, for.e.g. `MINLP_result_`, and budget ranges, returns a list of A- and D-optimality values of the given FIMs + + - `plot_data` receives the Cvxpy solution and Pyomo solution, and draws them on the same figure ### Rotary bed case study