update the ccro nonlinear test. 5 test cases (c02a7a4e) · Commits · Blais, Chris / SOAR_integration

documenting_failures/build_test_cases_nl.py

0 → 100644

+132 −0

Original line number	Diff line number	Diff line

		import os
		import sys
		import pickle
		import random
		import traceback
		import yaml
		import numpy as np
		import copy

		# soar_path = os.path.join(r"C:\Users\tjf\Documents\01_gitlab_repos\soar\python")
		# if soar_path not in sys.path:
		# sys.path.insert(0, soar_path)
		# else:
		# print("soar toolbox already in path")

		# add
		module_path = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
		if module_path not in sys.path:
		sys.path.insert(0, module_path)
		else:
		print("path already in sys.path")

		from setup_example_soar import define_system, get_pareto_front, get_one_result

		# for plotting
		import matplotlib.pyplot as plt

		# compare to soar
		from test_case_nonlinear import constr_mem_nl, labelConstraint
		from numerical_labelling.labelling import jacsym

		def build_test_case_nl(ptype=1):
		"""
		ptype: problem type
		"""
		constr_nl, varbs_nl, coeff_nl = constr_mem_nl()
		# by default, we have all of the columns for the jacobian.
		match ptype:
		# all columns included, all constraints
		case 1:
		# use all columns, all rows.
		varbs_active = {varbs_nl, coeff_nl}
		constr_active = constr_nl
		case 2:
		# use all columns
		varbs_active = {varbs_nl, coeff_nl}

		constr_active = {}
		cindx = 0
		# exclude any nonlinear constraints
		for iconstr, constr in constr_nl.items():
		if not constr.type.startswith("nl_"):
		constr_active[cindx] = constr
		cindx += 1

		case 3:
		# use all columns
		varbs_active = {varbs_nl, coeff_nl}

		constr_active = {}
		cindx = 0
		# exclude any nonlinear constraints
		for iconstr, constr in constr_nl.items():
		if not constr.type == "nl_pressure":
		constr_active[cindx] = constr
		cindx += 1
		case 4:
		# use all columns
		varbs_active = {varbs_nl, coeff_nl}

		constr_active = {}
		cindx = 0
		# exclude any nonlinear constraints
		for iconstr, constr in constr_nl.items():
		if not constr.type == "nl_component":
		constr_active[cindx] = constr
		cindx += 1
		case 5:
		# use all columns
		varbs_active = {varbs_nl, coeff_nl}

		# use all constraints
		constr_active = constr_nl

		case _:
		# use all columns, all rows.
		varbs_active = {varbs_nl, coeff_nl}
		constr_active = constr_nl

		constr_active_sy = {}
		for iconstr, constr in constr_active.items():
		if isinstance(constr, labelConstraint):
		constr_active_sy[iconstr] = constr.add_expr
		else:
		constr_active_sy[iconstr] = constr

		jac_nl = jacsym(constr_active_sy.values(), varbs_active.values())

		return constr_nl, varbs_nl, coeff_nl, jac_nl

		def make_layout(bl_layout, varbs_nl, ptype=1):
		"""
		bl_layout: bilinear layout
		varbs_active: active variables (e.g. mass, pressure, coefficients, etc)
		type: test to be performed
		"""
		# get number of pressure variables
		nmass = sum(1 for varb in varbs_nl if varb.lower().startswith('m'))
		nconc = sum(1 for varb in varbs_nl if varb.lower().startswith('x'))
		npress = sum(1 for varb in varbs_nl if varb.lower().startswith('p'))
		ncoeff = len(varbs_nl) - nmass - nconc - npress
		match ptype:
		case 1:
		layout = bl_layout + [False]npress + [False]ncoeff
		case 2:
		layout = bl_layout + [False]npress + [False]ncoeff
		case 3:
		layout = bl_layout + [False]npress + [True]ncoeff
		case 4:
		layout = bl_layout + [False]npress + [True]ncoeff
		case 5:
		layout = bl_layout + [False]npress + [True]ncoeff
		case _:
		layout = bl_layout + [False]npress + [False]ncoeff

		return layout

documenting_failures/soar_compare_nonlinear.py

+211 −102

Original line number	Diff line number	Diff line
		@@ -3,7 +3,10 @@ import os
		import sys
		import pickle
		import random
		import traceback
		import yaml
		import numpy as np
		from time import time

		soar_path = os.path.join(r"C:\Users\tjf\Documents\01_gitlab_repos\soar\python")
		if soar_path not in sys.path:
		@@ -29,7 +32,7 @@ import matplotlib.pyplot as plt
		from numerical_labelling.labelling import *

		# compare to soar
		from test_case_nonlinear import constr_mem_nl, get_x0_nl
		from build_test_cases_nl import build_test_case_nl, make_layout
		import soar.optimization.fEvaluate as feval
		import matplotlib.pyplot as plt

		@@ -51,7 +54,51 @@ def check_incidence(iIncidence):
		for row in iIncidence:
		assert sum(row) == 0, f"Incidence matrix row sums to {sum(row)}"

		def test_nonlinear_soar(plot=False, random_sampling=False):
		def log_err(error_log, e, result, method=""):
		err_trace = traceback.format_exc()
		if result not in error_log.keys():
		error_log[result] = {}

		error_log[result][method] = {
		"error":e, "method":method, "trace":err_trace,
		}

		return error_log

		def log_fail(result, failures, error_log, method, or_method, or_soar, layout):
		"""
		method: the value we've gotten, and the method used (eg obs_rref, red_lu)
		or_method: observability or redundancy for method
		or_soar: obs or red from soar
		"""
		if result not in failures.keys():
		failures[result] = {}
		if method not in error_log[result].keys():
		failures[result][method] = {
		"method":or_method, "soar":or_soar, "layout":layout,
		}
		else:
		failures[result][method] = "error"

		return failures

		def log_success(result, successes, method, soar, lu, layout):
		"""
		method: results for method (could be tuple of obs and red)
		soar: results for soar
		lu: results for lu
		layout: bool array of measured vars
		"""
		successes[result]= {
		"method":method,
		"soar":soar,
		"lu":lu,
		"layout":layout,
		}

		return successes

		def test_nonlinear_soar(plot=False, random_sampling=False, ptype=1, verbose=1):
		# establish simple membrane system for testing
		iIncidence = np.array(
		[ # e P1 Mx Me1 Sp1 P2
		@@ -76,7 +123,7 @@ def test_nonlinear_soar(plot=False, random_sampling=False):
		[1.5,-1]] # p2
		)
		pgraph = define_system(
		iIncidence=iIncidence, coordinates = coordinates, verbose=5,
		iIncidence=iIncidence, coordinates = coordinates, verbose=verbose,
		arcSplitter =[[5 - 1, 6 - 1, 7 - 1]],
		arcEqualComposition = [[1-1, 2-1], [5 - 1, 6 - 1, 7 - 1, 8 - 1]],
		numberOfComponents = 3,
		@@ -85,17 +132,25 @@ def test_nonlinear_soar(plot=False, random_sampling=False):
		xMeasured_default = pgraph.xMeasured

		# load pareto optimal layouts if they exist, run and save if they don't
		results_path = os.path.join(os.path.dirname(__file__), "bilinear_pareto_front.pkl")
		results_path = os.path.join(module_path, "documenting_failures", "pareto_fronts", "bilinear_pareto_front.pkl")
		if os.path.exists(results_path):
		with open(results_path, "rb") as f:
		results = pickle.load(f)
		else:
		results = get_pareto_front(pgraph, objective="bilinear", verbose=5)
		results = get_pareto_front(pgraph, objective="nonlinear", verbose=5)
		results = [int(res) for res in results]
		with open(results_path, "wb") as f:
		pickle.dump(results, f)

		# use the bilinear pareto front for comparison
		# file the results for the pareto front.
		soar_or_results_path = os.path.join(module_path, "documenting_failures", "pareto_fronts", "bilinear_pareto_front_labels.pkl")
		if os.path.exists(soar_or_results_path):
		with open(soar_or_results_path, "rb") as f:
		soar_or_results = pickle.load(f)
		else:
		soar_or_results = {}

		# use the nonlinear pareto front for comparison
		# generate cases where we can tell what is observable in the nonlinear system
		# picking cases where the pressures are known.

		@@ -107,31 +162,74 @@ def test_nonlinear_soar(plot=False, random_sampling=False):
		plt.show(block=False)

		# x0 = get_x0_nl()
		x0_file = os.path.join(os.path.dirname(__file__), "x0_solved.pkl")
		with open(x0_file, "rb") as f:
		x0 = pickle.load(f)

		constr_nl, varbs_nl = constr_mem_nl()
		jac_nl = jacsym(constr_nl.values(), varbs_nl.values())
		# # x0_file = os.path.join(os.path.dirname(__file__), "x0_solved.pkl")
		# with open(x0_file, "rb") as f:
		# x0 = pickle.load(f)
		x0_file = os.path.join(module_path, "operating_point_search", "saved_x0", "x0_ccro_final_ans_rev2.yaml")
		with open(x0_file, "r") as f:
		x0 = yaml.safe_load(f)
		# typecheck
		for key, val in x0.items():
		x0[key] = int(val)

		constr_nl, varbs_nl, coeff_nl, jac_nl = build_test_case_nl(ptype=ptype)
		all_varbs = {varbs_nl, coeff_nl}
		# constr_nl, varbs_nl, consts_nl = constr_mem_nl()
		# jac_nl = jacsym(constr_nl.values(), varbs_nl.values())

		# needed for removing values from obs, red for each method (soar doesn't have pressure)
		nstreams = iIncidence.shape[0]
		ncoeffs = len(coeff_nl)

		failures = {}
		error_log = {}
		success_log = {}

		# random.seed = 42
		# if random_sampling:
		# results = random.sample(results, 20)

		# add 100 random layouts with higher number of measured sensors, so we have
		# redundant sensor layouts tested
		random.seed(42)
		results += random.choices(range(30000, 65536), k=200)

		# results = [
		# # 8512,
		# # 4169,
		# # 17538,
		# # 16646,
		# # 3104,
		# 21800, # adding redundant layouts
		# 65535,
		# 31679,
		# 56134,
		# ]

		# for message at end
		errs_rref = 0
		errs_lu = 0
		errs_qr = 0
		fails_rref_obs = 0
		fails_lu_obs = 0
		fails_qr_obs = 0
		fails_rref_red = 0
		fails_lu_red = 0
		fails_qr_red = 0

		if random_sampling:
		results = random.sample(results, 20)

		results = [
		8512,
		4169,
		17538,
		16646,
		3104,
		]
		for result in results:
		pgraph.xMeasured = xMeasured_default.copy()
		# sens, obs, red = evaluate_obs_red(int(result), pgraph)
		sens, obs, red = get_one_result(pgraph, result, objective="bilinear", verbose=5)
		# get the result if we have saved it
		if result in soar_or_results.keys():
		sens, obs, red = soar_or_results[result]
		else:
		sens, obs, red = get_one_result(
		pgraph, result, objective="bilinear", verbose=verbose
		)
		soar_or_results[result] = sens, obs, red

		# for concentrations we have it backwards.
		rev_sens = np.concatenate([sens[8:],sens[0:8]])
		# convert to list. in soar -1 means undetermined, right?
		obs_soar_mass = np.array([True if s==1 else False for s in obs[:,-1]])
		red_soar_mass = np.array([True if s==1 else False for s in red[:,-1]])
		@@ -142,55 +240,49 @@ def test_nonlinear_soar(plot=False, random_sampling=False):

		obs_soar = np.concatenate([obs_soar_mass, obs_soar_conc])
		red_soar = np.concatenate([red_soar_mass, red_soar_conc])
		layout = list(np.array(sens).astype(bool))
		bl_layout = list(np.array(rev_sens).astype(bool))

		# append the pressure sensors (all unmeasured) to the layout
		# for concentrations we have it backwards.
		rev_sens = np.concatenate([sens[8:],sens[0:8]])
		layout = layout+[False]*8

		# make the layout including pressure, other variables
		# layout = layout+[False]*8
		layout = make_layout(bl_layout, all_varbs, ptype=ptype)

		error_log[result] = {}
		try:
		obs_rref, red_rref = observability_redundancy_labelling_rref(jac_nl, layout, x0=x0)
		obs_rref_p, red_rref_p = observability_redundancy_labelling_rref(jac_nl, layout, x0=x0)
		except Exception as e:
		if result not in error_log.keys():
		error_log[result] = {}
		error_log[result]['rref'] = {
		"error":e, "type":"obs", "method":"rref",
		}
		obs_rref = np.full(len(obs_soar), fill_value=np.nan)
		red_rref = np.full(len(red_soar), fill_value=np.nan)
		err_log = log_err(error_log, e, result, method="rref")
		obs_rref_p = np.full(len(obs_soar)+8, fill_value=np.nan)
		red_rref_p = np.full(len(red_soar)+8, fill_value=np.nan)
		errs_rref+=1
		try:
		obs_lu, red_lu = observability_redundancy_labelling_lu(jac_nl, layout, x0=x0)
		obs_lu_p, red_lu_p = observability_redundancy_labelling_lu(jac_nl, layout, x0=x0)
		except Exception as e:
		if result not in error_log.keys():
		error_log[result] = {}
		error_log[result]['lu'] = {
		"error":e, "type":"obs", "method":"lu",
		}
		err_log = log_err(error_log, e, result, method="lu")
		# add in a nan array so we can run the other sections
		obs_lu = np.full(len(obs_soar), fill_value=np.nan)
		red_lu = np.full(len(red_soar), fill_value=np.nan)
		obs_lu_p = np.full(len(obs_soar)+8, fill_value=np.nan)
		red_lu_p = np.full(len(red_soar)+8, fill_value=np.nan)
		errs_lu+=1

		try:
		obs_qr, red_qr = observability_redundancy_labelling_qr(jac_nl, layout, x0=x0)
		obs_qr_p, red_qr_p = observability_redundancy_labelling_qr(jac_nl, layout, x0=x0)
		except Exception as e:
		if result not in error_log.keys():
		error_log[result] = {}
		error_log[result]['qr'] = {
		"error":e, "type":"obs", "method":"qr",
		}
		obs_qr = np.full(len(obs_soar), fill_value=np.nan)
		err_log = log_err(error_log, e, result, method="qr")
		obs_qr_p = np.full(len(obs_soar)+8, fill_value=np.nan)
		red_qr_p = np.full(len(red_soar)+8, fill_value=np.nan)
		errs_qr+=1

		# for comparison, remove the values for pressure, coeff
		# TODO: structure this like the SOAR matrix, so we can refer to column
		ncut = nstreams + ncoeffs
		obs_rref = obs_rref_p[:-ncut]
		red_rref = red_rref_p[:-ncut]
		obs_lu = obs_lu_p[:-ncut]
		red_lu = red_lu_p[:-ncut]
		obs_qr = obs_qr_p[:-ncut]
		# red_qr = red_qr_p[:-8]
		red_qr = np.full(len(red_soar), fill_value=np.nan)

		# for comparison, remove last 8 values
		obs_rref = obs_rref[:-8]
		red_rref = red_rref[:-8]
		obs_lu = obs_lu[:-8]
		red_lu = red_lu[:-8]
		obs_qr = obs_qr[:-8]
		red_qr = obs_qr[:-8]

		obs_rref_equ = np.array_equal(obs_rref, obs_soar)
		obs_lu_equ = np.array_equal(obs_lu, obs_soar)
		obs_qr_equ = np.array_equal(obs_qr, obs_soar)
		@@ -201,71 +293,88 @@ def test_nonlinear_soar(plot=False, random_sampling=False):

		failures[result] = {}
		if not obs_rref_equ:
		if result not in failures.keys():
		failures[result] = {}
		if "rref" not in error_log[result].keys():
		failures[result]['rref_obs'] = {
		"obs_method":obs_rref, "obs_soar":obs_soar, "type":"obs", "method":"rref",
		}
		else:
		failures[result]["rref_obs"] = "error"
		failures = log_fail(result, failures, error_log, "obs_rref", obs_rref, obs_soar, layout)
		fails_rref_obs += 1
		if not red_rref_equ:
		failures = log_fail(result, failures, error_log, "red_rref", red_rref, red_soar, layout)
		fails_rref_red += 1

		if not obs_lu_equ:
		if result not in failures.keys():
		failures[result] = {}
		if "lu" not in error_log[result].keys():
		failures[result]['lu_obs'] = {
		"obs_method":obs_lu, "obs_soar":obs_soar, "type":"obs", "method":"lu",
		}
		else:
		failures[result]['lu_obs'] = "error"
		failures = log_fail(result, failures, error_log, "obs_lu", obs_lu, obs_soar, layout)
		fails_lu_obs += 1
		if not red_lu_equ:
		failures = log_fail(result, failures, error_log, "red_lu", red_lu, red_soar, layout)
		fails_lu_red += 1

		if not obs_qr_equ:
		if result not in failures.keys():
		failures[result] = {}
		if "qr" not in error_log[result].keys():
		failures[result]['qr_obs'] = {
		"obs_method":obs_qr, "obs_soar":obs_soar, "type":"obs", "method":"qr",
		}
		else:
		failures[result]['qr_obs'] = "error"
		failures = log_fail(result, failures, error_log, "obs_qr", obs_qr, obs_soar, layout)
		fails_qr_obs += 1
		if not red_qr_equ:
		failures = log_fail(result, failures, error_log, "red_qr", red_qr, red_soar, layout)
		fails_qr_red += 1

		if (obs_rref_equ and red_rref_equ) and (not obs_lu_equ or not red_lu_equ):
		success_log = log_success(
		result=result,
		successes=success_log,
		method=(obs_rref, red_rref),
		soar=(obs_soar, red_soar),
		lu=(obs_lu, red_lu),
		layout = layout,
		)

		if not red_rref_equ:
		if result not in failures.keys():
		failures[result] = {}
		if "rref_red" not in error_log[result].keys():
		failures[result]['rref_red'] = {
		"method":red_rref, "soar":red_soar, "type":"red", "method":"rref",
		}
		else:
		failures[result]["rref_obs"] = "error"

		# save the soar results to save time:
		with open(soar_or_results_path, "wb") as f:
		pickle.dump(soar_or_results, f)

		print("!"*80)
		print(f"!!! Method {str(ptype)} results:")
		print("!"*80)

		if len(error_log) > 0:
		err_log = os.path.join(r".\documenting_failures\errors\nonlinear_errors.pkl")
		err_log = os.path.join(rf".\documenting_failures\errors\nonlinear_errors_type{str(ptype)}.pkl")
		if not os.path.exists(os.path.dirname(err_log)):
		os.mkdir(os.path.dirname(err_log))

		with open(err_log, "wb") as f:
		pickle.dump(obj=error_log, file=f)
		print(f"errors in {len(error_log)} layouts out of {len(results)}")
		print(f"errors out of {len(results)} layouts:\nrref: {errs_rref}\nLU: {errs_lu}\nQR: {errs_qr}")
		else:
		print("no errors")

		if len(failures) > 0:
		fail_log = os.path.join(r".\documenting_failures\failures\nonlinear_failures.pkl")
		fail_log = os.path.join(rf".\documenting_failures\failures\nonlinear_failures_type{str(ptype)}.pkl")
		if not os.path.exists(os.path.dirname(fail_log)):
		os.mkdir(os.path.dirname(fail_log))

		with open(fail_log, "wb") as f:
		pickle.dump(obj=failures, file=f)
		print(f"failures in {len(fail_log)} layouts out of {len(results)}")

		print(f"failures out of {len(results)} layouts obs:\nrref: {fails_rref_obs}\nLU: {fails_lu_obs}\nQR: {fails_qr_obs}")
		print(f"failures out of {len(results)} layouts red:\nrref: {fails_rref_red}\nLU: {fails_lu_red}\nQR: {fails_qr_red}")
		else:
		print("all layouts pass for linear system")

		if len(success_log) > 0:
		success_path = os.path.join(rf".\documenting_failures\successes\nonlinear_successes_type{str(ptype)}.pkl")
		if not os.path.exists(os.path.dirname(success_path)):
		os.mkdir(os.path.dirname(success_path))

		with open(success_path, "wb") as f:
		pickle.dump(obj=success_log, file=f)
		print(f"RREF agree SOAR, LU disagree in {len(success_log)} layouts out of {len(results)}")


		if __name__ == "__main__":
		test_nonlinear_soar(random_sampling = True)
		t1 = time()
		test_nonlinear_soar(random_sampling = False, ptype=1, verbose=0)
		test_nonlinear_soar(random_sampling = False, ptype=2, verbose=0)
		test_nonlinear_soar(random_sampling = False, ptype=3, verbose=0)
		test_nonlinear_soar(random_sampling = False, ptype=4, verbose=0)
		test_nonlinear_soar(random_sampling = False, ptype=5, verbose=0)
		t2 = time()

		print(f"total test time: {t2 - t1} seconds")

documenting_failures/test_case_nonlinear.py

+64 −23

File changed.

Preview size limit exceeded, changes collapsed.