consolidate functions for interfacing with soar (8906e8eb) · Commits · Blais, Chris / SOAR_integration

tools/soar_utilities.py

0 → 100644

+510 −0

Original line number	Diff line number	Diff line
		from context import *
		import numpy as np
		import time
		import pickle
		import matplotlib.pyplot as plt
		import matplotlib

		# soar_libraries
		import soar.optimization.fEvaluate as feval
		from soar.graph_labeling.fDefineSystem import ProcessGraph
		from soar.graph_labeling.fPlotGraph import plot_graph
		from soar.optimization.fParetoDominance import pareto_dominance
		from soar.optimization.fFigurePareto import figure_pareto
		from soar.optimization.fParetoSearch import pareto_search


		def check_incidence(iIncidence):
		"""
		simple check that in nodes equal to out nodes (graph is closed)
		"""
		for row in iIncidence:
		assert sum(row) == 0, f"Incidence matrix row sums to {sum(row)}"


		def define_system(iIncidence=[], coordinates = [], verbose=5,
		arcSplitter =[[5 - 1, 6 - 1, 7 - 1]],
		arcEqualComposition = [[1-1, 2-1], [5 - 1, 6 - 1, 7 - 1, 8-1]],
		numberOfComponents = 2,
		):
		if verbose >= 0:
		print("***** DEFINE SYSTEM *****")
		pass

		strSystemID = "membrane1"
		print(strSystemID)
		# ==============================================================================================
		if verbose >= 1:
		print(" GRAPH STRUCTURE ")
		pass
		# ==============================================================================================

		# ---------------------------
		# Basic graph_labeling structure
		# ---------------------------

		iArcEnergy = []
		iArcReaction = []
		iArcSplitter = arcSplitter
		iArcEqualComposition = arcEqualComposition
		iNumberOfComponent = numberOfComponents

		# E C R S SPL
		if len(iIncidence) == 0:
		iIncidence = np.array(
		[
		[-1, +1, 0, 0, 0], # 1
		[0, -1, +1, 0, 0], # 2
		[0, 0, -1, +1, 0], # 3
		[+1, 0, 0, -1, 0], # 4
		[0, 0, 0, -1, +1], # 5
		[0, +1, 0, 0, -1], # 6
		[+1, 0, 0, 0, -1],
		]
		) # 7
		if len(coordinates) == 0:
		coordinates = np.array([[2, 1], [0, 0], [0.75, 0], [1.5, 0], [2, -1],
		[3,0]])

		iIncidence = iIncidence.T

		# ==============================================================================================
		if verbose >= 1:
		print(" STREAM COMPOSITION ")
		pass
		# ==============================================================================================

		[iNumberOfNode, iNumberOfArc] = iIncidence.shape
		ii = np.arange(iNumberOfComponent)
		t = iNumberOfComponent
		nn = np.arange(iNumberOfComponent) + 1 + t
		h = 2 * t + 1
		m = 2 * (iNumberOfComponent + 1)
		t = iNumberOfComponent

		Present = np.zeros([iNumberOfArc, m + 1])
		Present[:, m] = True # flow is always present
		Present[:, t] = True # energy is always present

		match strSystemID:
		case _:
		Present[:, 0] = True # component 1
		Present[:, 1] = True # component 2
		xArcReaction = np.zeros([iNumberOfArc, 1]).astype("bool")
		xArcReaction[iArcReaction] = True
		xArcEnergy = np.zeros([iNumberOfArc, 1]).astype("bool")
		xArcEnergy[iArcEnergy] = True
		xArcPhysical = ~np.logical_or(xArcReaction, xArcEnergy).squeeze()
		Present[xArcPhysical, 2] = True # component 3
		pass

		Present[:, nn] = Present[:, ii]
		Present[:, h] = Present[:, t]

		if verbose >= 1:
		print(" MEASUREMENTS ")
		pass

		xMeasured = np.zeros([iNumberOfArc, m + 1])
		xMeasured[Present == 0] = np.nan
		xKnown = np.zeros([iNumberOfArc, m + 1])
		xKnown[Present == 0] = np.nan


		one_component = True

		S = [[]]

		# ========================================================================
		# GROUND TRUTH
		# ========================================================================

		if verbose >= 1:
		print(" GROUND TRUTH ")
		pass


		xRefObservable = []
		xRefRedundant = []


		# ========================================================================
		# SPECIFY POSSIBLE SENSOR LOCATIONS AND WEIGHTS FOR OBSERVABILITY
		# ========================================================================

		xSensorCandidate = Present.copy()
		xSensorCandidate[:, :] = False
		xSensorCandidate[:, m] = np.all(
		np.vstack([Present[:, m], xKnown[:, m] != 1, xArcPhysical]),
		axis=0,
		)
		xSensorCandidate[:, 0] = np.all(
		np.vstack([Present[:, 0], xKnown[:, 1] != 1, xArcPhysical]),
		axis=0,
		)

		rWeightObservability = Present.copy()
		rWeightObservability[:, :] = 0
		rWeightObservability[:, m] = 1
		rWeightObservability[:, ii[0]] = 1
		rWeightObservability[:, nn[0]] = 1

		rWeightCost = xSensorCandidate.copy()
		rWeightCost = rWeightCost * 1.0
		rWeightCost[:, :] = 1.0
		rWeightRedundancy = xSensorCandidate.copy()
		rWeightRedundancy = rWeightRedundancy * 1.0

		process_graph = ProcessGraph(
		strSystemID=strSystemID,
		incidence_matrix=iIncidence,
		presence_matrix=Present,
		iNumberOfComponent=iNumberOfComponent,
		coordinates=coordinates,
		iArcReaction=iArcReaction,
		iArcEnergy=iArcEnergy,
		iArcSplitter=iArcSplitter,
		iArcEqualComposition=iArcEqualComposition,
		xKnown=xKnown,
		xMeasured=xMeasured,
		xSensorCandidate=xSensorCandidate,
		rWeightCost=rWeightCost,
		rWeightObservability=rWeightObservability,
		rWeightRedundancy=rWeightRedundancy,
		one_component=one_component,
		xRefObservable=xRefObservable,
		xRefRedundant=xRefRedundant,
		stoichiometry=S,
		)

		return process_graph

		def get_pareto_front(process_graph, plot_all=False, objective="linear",
		save_results=False, verbose=5):
		# Visualization setup
		if plot_all:
		matplotlib.use("Qt5Agg")
		plt.ion()
		plt.close("all")
		iSystemID = 99

		# for iSearchMethod in [1,2]:
		# match iSearchMethod:
		# case 0: # search by exhaustive enumeration
		# xExhaustive = True
		# feval.RedundancyDeficit = False
		# pass
		# case 1: # search with surplus redundanccy
		# xExhaustive = False
		# feval.RedundancyDeficit = False
		# pass
		# case 2: # search without surplus redundanccy
		# xExhaustive = False
		# feval.RedundancyDeficit = True
		# pass
		# case _:
		# print("Unknown option")
		# pass

		match objective:
		case "linear":
		iObj = 1
		# we want to enumerate all hydraulic layouts
		strObjective = "FlowOnly"
		xExhaustive = True
		feval.RedundancyDeficit = False
		iSearchMethod = 0
		case "bilinear":
		iObj = 2
		# too many layouts for flow + component, just do pareto
		strObjective = "FlowAndComponent1"
		xExhaustive = False
		feval.RedundancyDeficit = False
		iSearchMethod = 1
		pass
		case _:
		print("Unknown option")
		pass

		# ncriteria is if we are looking at red and obs, ncrit=2 is just obs, nCriteria 3 is all 3
		nCriteria = 3
		strProblem = (
		"Front_"
		+ str(iSystemID)
		+ "_"
		+ str(iObj)
		+ "_"
		+ str(nCriteria)
		)
		print(strProblem)

		tElapsedOuter = np.ones(5) * np.nan

		tStartOuter = time.perf_counter()
		tElapsedOuter[0] = time.perf_counter() - tStartOuter

		if plot_all:
		fig = plot_graph(process_graph)
		fig.canvas.draw()
		fig.canvas.flush_events()
		plt.show(block=False)

		process_graph.add_graph_info(verbose)
		tElapsedOuter[1] = time.perf_counter() - tStartOuter
		# iNumberOfArc = process_graph.iNumberOfArc

		match strObjective:
		case "FlowOnly":
		process_graph.xSensorCandidate[:, :-1] = 0
		process_graph.rWeightObservability[:, :-1] = 0
		process_graph.rWeightObservability[
		process_graph.xArcPhysical.squeeze(), :
		] = 1
		pass
		case _:
		pass

		X = process_graph.xSensorCandidate
		nSensor = np.sum(X).astype(int)
		nLayout = int(2**nSensor)
		# nRedundantMax = nSensor
		# nObservableMax=np.sum(process_graph.rWeightObservability)

		[I, J] = np.where(process_graph.xSensorCandidate)
		process_graph.ijSensorCandidate = np.where(
		process_graph.xSensorCandidate.T
		)[::-1]
		process_graph.nSensorCandidate = nSensor

		feval.criteria = np.ones([nLayout, nCriteria]) * np.nan
		feval.observableconfig = (
		np.ones(
		[
		nLayout,
		(process_graph.iNumberOfComponent + 1) * 2 + 1,
		]
		)
		* np.nan
		)

		# nLayout = 2
		iLayoutBnd = [0, nLayout]

		# ====================================
		tElapsedOuter[2] = time.perf_counter() - tStartOuter

		feval.lapse = np.ones([nLayout, nCriteria]) * np.nan

		xMeasured_default = process_graph.xMeasured.copy()

		if xExhaustive:
		for iLayout in range(iLayoutBnd[0], iLayoutBnd[-1]):
		print("layout {}".format(iLayout))
		process_graph.xMeasured = xMeasured_default.copy()
		feval.evaluate_layout(
		iLayout,
		process_graph,
		nCriteria,
		)
		pass

		crituniq = np.unique(feval.criteria, axis=0)
		crituniq[:, 1:] = -crituniq[:, 1:]
		xIncludeUniq = pareto_dominance(crituniq, crituniq)
		critfront = crituniq[xIncludeUniq, :]
		critfront[:, 1:] = -critfront[:, 1:]
		xLayoutOnFront = np.array(
		[
		np.any(np.all(x[None, :] == critfront, axis=1))
		for x in feval.criteria
		]
		)
		Solutions = np.where(xLayoutOnFront)[0]
		iObjectives = feval.criteria[Solutions, :]

		tElapsedOuter[3] = time.perf_counter() - tStartOuter
		pass

		else:
		# iLayoutBnd = [16, 31]
		def fbnd(iLayoutInterval):
		bounds = feval.fBound(
		iLayoutInterval,
		process_graph,
		nCriteria,
		)
		return bounds

		[LB, UB, xLayoutL, xLayoutU] = fbnd(iLayoutBnd)

		iTree = pareto_search(
		fbnd,
		iLayoutBnd,
		100,
		100000,
		strProblem,
		)

		tElapsedOuter[3] = time.perf_counter() - tStartOuter
		iObjectives = iTree[:, 2 + np.arange(nCriteria)]
		iObjectives[:, 1] = -iObjectives[:, 1]
		if nCriteria >= 3:
		iObjectives[:, 2] = -iObjectives[:, 2]
		if feval.RedundancyDeficit:
		iObjectives[:, 2] = (
		iObjectives[:, 2] + iObjectives[:, 0]
		)
		pass
		pass
		Solutions = iTree[:, 0]
		pass


		if plot_all:
		[_, _, _, _, _, fig, _] = figure_pareto(
		iObjectives,
		strProblem,
		)
		fig.canvas.draw()
		fig.canvas.flush_events()
		plt.show(block=False)

		criteria = feval.criteria.copy()
		observableconfig = feval.observableconfig.copy()

		# save a pickle file
		if save_results:
		folder = "./result/"
		if not os.path.exists(folder):
		print("Create result folder")
		os.mkdir(folder)
		pass

		filename = "".join(
		[
		folder,
		"res_Sys",
		str(iSystemID),
		"_Crit",
		str(nCriteria),
		"_Obj",
		strObjective,
		"_Search",
		str(iSearchMethod),
		".p",
		]
		)
		with open(filename, "wb") as filehandle:
		print("Store results")
		pickle.dump(
		[
		criteria,
		Solutions,
		process_graph,
		observableconfig,
		],
		filehandle,
		)
		print("Store results - done")
		pass


		return Solutions

		def get_one_result(process_graph, layout, objective="linear", verbose=5):

		iSystemID = 99

		match objective:
		case "linear":
		iObj = 1
		# we want to enumerate all hydraulic layouts
		strObjective = "FlowOnly"
		xExhaustive = True
		feval.RedundancyDeficit = False
		case "bilinear":
		iObj = 2
		# too many layouts for flow + component, just do pareto
		strObjective = "FlowAndComponent1"
		xExhaustive = False
		feval.RedundancyDeficit = False
		pass
		case _:
		print("Unknown option")
		pass

		# ncriteria is if we are looking at red and obs, ncrit=2 is just obs
		nCriteria = 2
		strProblem = (
		"Front_"
		+ str(iSystemID)
		+ "_"
		+ str(iObj)
		+ "_"
		+ str(nCriteria)
		)
		print(strProblem)

		tElapsedOuter = np.ones(5) * np.nan

		tStartOuter = time.perf_counter()
		tElapsedOuter[0] = time.perf_counter() - tStartOuter

		# this is where cutset enumeration happens
		process_graph.add_graph_info(verbose)

		tElapsedOuter[1] = time.perf_counter() - tStartOuter
		# iNumberOfArc = process_graph.iNumberOfArc

		match strObjective:
		case "FlowOnly":
		process_graph.xSensorCandidate[:, :-1] = 0
		process_graph.rWeightObservability[:, :-1] = 0
		process_graph.rWeightObservability[
		process_graph.xArcPhysical.squeeze(), :
		] = 1
		pass
		case _:
		pass

		X = process_graph.xSensorCandidate
		nSensor = np.sum(X).astype(int)
		nLayout = int(2**nSensor)
		# nRedundantMax = nSensor
		# nObservableMax=np.sum(process_graph.rWeightObservability)

		[I, J] = np.where(process_graph.xSensorCandidate)
		process_graph.ijSensorCandidate = np.where(
		process_graph.xSensorCandidate.T
		)[::-1]
		process_graph.nSensorCandidate = nSensor

		feval.criteria = np.ones([nLayout, nCriteria]) * np.nan
		feval.observableconfig = (
		np.ones(
		[
		nLayout,
		(process_graph.iNumberOfComponent + 1) * 2 + 1,
		]
		)
		* np.nan
		)

		iLayoutBnd = [0, nLayout]

		# ====================================
		tElapsedOuter[2] = time.perf_counter() - tStartOuter

		feval.lapse = np.ones([nLayout, nCriteria]) * np.nan

		xMeasured_default = process_graph.xMeasured.copy()

		# evaluate single layout by number
		process_graph.xMeasured = xMeasured_default.copy()
		sens, obs, red = feval.evaluate_obs_red(
		layout,
		process_graph,
		verbose = verbose,
		)

		return sens, obs, red