diff --git a/pyproject.toml b/pyproject.toml index 6676b62..894a07d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,8 @@ [build-system] requires = ["setuptools>=61.0", "setuptools_scm>=8.0"] build-backend = "setuptools.build_meta" +[tool.setuptools-git-versioning] +enabled = true [project] name = "ftuutils" dynamic = ["version"] diff --git a/src/ftuutils/base.py b/src/ftuutils/base.py index 6b8cd31..9a6b88f 100644 --- a/src/ftuutils/base.py +++ b/src/ftuutils/base.py @@ -345,6 +345,7 @@ def resolve(self,G,phsdefinitions,phsdata,generateCompositionJSON=True): Returns: bool,dict: True if a FTU was successfully constructed, FTU project as python dict """ nodeAttrib = nx.get_node_attributes(G,'phs') + nodeType = nx.get_node_attributes(G,'type') unresolved = False networks = dict() existingNets = dict() @@ -401,6 +402,8 @@ def resolve(self,G,phsdefinitions,phsdata,generateCompositionJSON=True): if not unresolved: for nd in curNodes: + if nodeType[nd] == 'out': + raise Exception(f"Unexpected workflow, boundary node up for resolution") #Get phstype and component connection information if nodeAttrib[nd] not in phsdefinitions: unresolved = True @@ -661,13 +664,14 @@ def resolve(self,G,phsdefinitions,phsdata,generateCompositionJSON=True): return (not unresolved),result - def composeCompositePHS(self,G,phsdefinitions=None,phsdata=None): + def composeCompositePHS(self,G,phsdefinitions=None,phsdata=None,substituteParameters=True): """Method to generate composite PHS from Args: G (networkx Graph or dict): FTU graph (networkx.Graph) or Dict description phsdefinitions (dict): Description of phs classes used in the FTU required if G is not dict phsdata (dict): Data about networks, external inputs etc used in the FTU required if G is not dict + substituteParameters (bool): Substitute phs parameters and simplify before composition, phs parameters are eliminated. Default True """ composition = None ftuGraph = None @@ -682,7 +686,7 @@ def composeCompositePHS(self,G,phsdefinitions=None,phsdata=None): composer = Composer() composer.loadComposition(composition) - composer.compose() + composer.compose(substituteParameters) if ftuGraph is not None: composer.setConnectivityGraph(ftuGraph) return composer @@ -786,15 +790,28 @@ def __init__(self,points,defaultPHS,conductivityTensor,edgeLengththeshold=1.5,de #Check edge lengths edge_len = np.array([np.linalg.norm(pdict[u]-pdict[v]) for u, v in G.edges()]) ecutoff = np.mean(edge_len)*self.edgeLengthThreshold - - for i,e in enumerate(G.edges()): - if edge_len[i]>ecutoff: - G.remove_edge(e[0],e[1]) - else: - evec = (pdict[e[1]] - pdict[e[0]])/edge_len[i] - wt = np.abs(conductivityTensor.dot(evec)) #Weight should be positive - projection can be negative due to direction - eweight[e] = wt/edge_len[i] #normalize by edge length - + if np.issubdtype(conductivityTensor.dtype,np.number): + for i,e in enumerate(G.edges()): + if edge_len[i]>ecutoff: + G.remove_edge(e[0],e[1]) + else: + evec = (pdict[e[1]] - pdict[e[0]]) + wt = np.abs(conductivityTensor.dot(evec)) #Weight should be positive - projection can be negative due to direction + eweight[e] = wt/edge_len[i] #normalize by edge length + else:#String or sympy + for i,e in enumerate(G.edges()): + if edge_len[i]>ecutoff: + G.remove_edge(e[0],e[1]) + else: + evec = (pdict[e[1]] - pdict[e[0]])/edge_len[i] + wt = [] + for j,x in enumerate(conductivityTensor): + if evec[j] != 0.0: + wt.append(f"{evec[j]/edge_len[i]}*"+x) + #Weight should be positive - projection can be negative due to direction + eweight[e] = "abs("+"+".join(wt)+")" #normalize by edge length + + nx.set_edge_attributes(G,eweight,self.getDefaultNetworkID()) self.graph = G diff --git a/src/ftuutils/codegenerationutils.py b/src/ftuutils/codegenerationutils.py index 25e5df1..104a498 100644 --- a/src/ftuutils/codegenerationutils.py +++ b/src/ftuutils/codegenerationutils.py @@ -28,7 +28,7 @@ def exportAsPython(composer): Args: composer (compositionutils.Composer): Composer instance that has the resolved FTU """ - numconstants,constantsubs,nonlinearrhsterms,inputs,arrayedinputs,arraymapping,uCapterms,ucapdescriptive,nonlineararrayedrhsterms,nonlinearrhstermsdescriptive,arrayedrhs,invarraymapping,rhs,ubaridxmap,cleaninputs = composer.generatePythonIntermediates() + numconstants,phsconstants,constantsubs,nonlinearrhsterms,inputs,arrayedinputs,arraymapping,uCapterms,ucapdescriptive,nonlineararrayedrhsterms,nonlinearrhstermsdescriptive,arrayedrhs,invarraymapping,rhs,ubaridxmap,ftuidmap,cleaninputs = composer.generatePythonIntermediates() stateVec = Matrix(composer.stateVec) # Generate metedata variabledescription = 'VOI_INFO = {"name": "t", "units": "second", "component": "main", "type": VariableType.VARIABLE_OF_INTEGRATION}\n' @@ -36,6 +36,13 @@ def exportAsPython(composer): for k, v in composer.statevalues.items(): variabledescription += f'\t{{"name": "{k}", "units": "{v["units"]}", "component": "main", "type": VariableType.STATE}},\n' variabledescription += "]\n\nVARIABLE_INFO = [\n" + for k, v in ubaridxmap.items(): + #TODO get the dimension + variabledescription += f'\t{{"name": "{k}", "units": "dimensionless", "component": "main", "type": VariableType.EXTERNAL_INPUT}},\n' + if len(ftuidmap)>0: + for k,v in ftuidmap.items(): + variabledescription += f'\t{{"name": "{k}", "units": "dimensionless", "component": "main", "type": VariableType.CONSTANT}},\n' + # Maintain this order when creating variables # Do constant subs, constant subs will have multiple values for the same constant due to precision definedNames = [] @@ -44,6 +51,11 @@ def exportAsPython(composer): if not v.name.startswith("-"): variabledescription += f'\t{{"name": "{v}", "units": "dimensionless", "component": "main", "type": VariableType.CONSTANT}},\n' definedNames.append(v) + for k,v in phsconstants.items(): + if k.name in arraymapping: + vunit = v['units'] + variabledescription += f'\t{{"name": "{k}", "units": "{vunit}", "component": "main", "type": VariableType.CONSTANT}},\n' + # TODO compute the units of calculated terms # Do uCap terms @@ -74,6 +86,13 @@ def exportAsPython(composer): STATE_COUNT = {len(composer.stateVec)} VARIABLE_COUNT = {numconstants} +def heaviside(x): + if x > 0: + return 1.0 + return 0.0 + +def Abs(x): + return np.fabs(x) class VariableType(Enum): VARIABLE_OF_INTEGRATION = 0 @@ -100,6 +119,12 @@ def initialise_variables(states, variables):\n""" except: stmt = f"\t{arraymapping[k.name]} = {v['value']} #{k}\n" pycode += stmt + for k, v in ubaridxmap.items(): + pycode += f"\t{v} = 0.0 #{k} External input\n" + + if len(ftuidmap)>0: + for k,v in ftuidmap.items(): + pycode += f"\t{v} = 1.0 #{k} This needs to be set for accurate simulation\n" # Do constant subs definedVariables = [] @@ -112,6 +137,16 @@ def initialise_variables(states, variables):\n""" stmt = f"\t{arraymapping[v.name]} = {k} #{v}\n" pycode += stmt definedVariables.append(v) + for v, k in phsconstants.items(): + if v not in definedVariables: + if v.name in arraymapping: + try: + stmt = f"\t{arraymapping[v.name]} = {float(k['value']):6f} #{v}\n" + except: + stmt = f"\t{arraymapping[v.name]} = {k['value']} #{v}\n" + pycode += stmt + definedVariables.append(v) + pycode += "\ndef compute_computed_constants(variables):\n\tpass\n\n" pycode += "def compute_variables(voi, states, rates, variables):\n\tt=voi #mapping to t\n" # Do uCap terms @@ -126,8 +161,8 @@ def initialise_variables(states, variables):\n""" pycode += f"\t#{nonlinearrhstermsdescriptive[k]}\n" pycode += f"\t{k} = {v}\n" for i, v in enumerate(arrayedrhs): - pycode += fr"\t#\dot{{{composer.stateVec[i]}}} = {_stringsubs(str(v),invarraymapping)} # {sympy.simplify(rhs[i,0])}\n" - pycode += fr"\trates[{i}] = {v}\n" + pycode += f"\t#\dot{{{composer.stateVec[i]}}} = {_stringsubs(str(v),invarraymapping)} # {sympy.simplify(rhs[i,0])}\n" + pycode += f"\trates[{i}] = {v}\n" # Do inputs pycode += "\ndef compute_inputs(voi,inputs,states,variables):\n" @@ -141,11 +176,6 @@ def initialise_variables(states, variables):\n""" pycode += f''' from math import exp -def Heaviside(x): - if x > 0: - return 1.0 - return 0.0 - def process_time_sensitive_events(voi, states, rates, variables): """Method to process events such as (re)setting inputs, updating switches etc Unline process_events, this method is called in rhs calculation @@ -204,10 +234,6 @@ def solve_model(starttime=0,stoptime=300,steps=300): initialise_variables(states,variables) # Set timespan to solve over voi = np.linspace(starttime, stoptime, steps) - print(f" c0 {{variables[44]}}") - print(f" ct {{variables[27]}}") - print(f" sigma {{variables[136]}}") - print(f" sqpi {{variables[135]}}") # Construct ODE object to solve r = ode(rhs) @@ -242,7 +268,7 @@ def solve_model(starttime=0,stoptime=300,steps=300): ax.plot(r[ix,:]) ax.title.set_text(f'{{ix//{(len(composer.stateVec)//len(composer.inodeIndexes))}+1}}') ix += {(len(composer.stateVec)//len(composer.inodeIndexes))} - if ix+{(len(composer.stateVec)//len(composer.inodeIndexes))} > result.shape[0]: + if ix+{(len(composer.stateVec)//len(composer.inodeIndexes))} > r.shape[0]: break plt.subplots_adjust(hspace=0.3) plt.subplots_adjust(wspace=0.3) @@ -259,7 +285,7 @@ def exportAsODEStepper(composer,modelName): Setup the code as a Python class, with the ability to step through time and set inputs """ - numconstants,constantsubs,nonlinearrhsterms,inputs,arrayedinputs,arraymapping,uCapterms,ucapdescriptive,nonlineararrayedrhsterms,nonlinearrhstermsdescriptive,arrayedrhs,invarraymapping,rhs_,ubaridxmap,cleaninputs = composer.generatePythonIntermediates() + numconstants,phsconstants,constantsubs,nonlinearrhsterms,inputs,arrayedinputs,arraymapping,uCapterms,ucapdescriptive,nonlineararrayedrhsterms,nonlinearrhstermsdescriptive,arrayedrhs,invarraymapping,rhs_,ubaridxmap,ftuidmap,cleaninputs = composer.generatePythonIntermediates() #Also handle constant subs revconstantsubs = dict() @@ -330,7 +356,9 @@ def exportAsODEStepper(composer,modelName): csyms = list(set(csyms)) csym_map = {} for kc in csyms: - csym_map[kc]=int(inputStateSymbolMap[kc].replace("states[","").replace("]","")) + if kc in inputStateSymbolMap: + csym_map[kc]=int(inputStateSymbolMap[kc].replace("states[","").replace("]","")) + #Map node label to input symbols and expressions - using subs as xreplace doesnt seem to handle ** subbedHam = sympy.simplify(cham.xreplace(parvals).xreplace(revconstantsubs)) supportEnergyCalculations[k] = {'cmap':csym_map,'expr':cinps,'ham':str(subbedHam),'rhsix':rhsidxmap,'cix':cix} @@ -346,7 +374,13 @@ def exportAsODEStepper(composer,modelName): #input name will be suffixed by the network id through which it communicates ik = '_'.join(k.split('_')[:-1]) inputhook[ik] = v - inputhookcode = f'class {modelName}Hooks():\n #PHS Input name and variable entry\n' + inputhookcode = f'class {modelName}Hooks():\n #FTU parameters\n' + inputhookcode += f" CELL_COUNT = {len(composer.inodeIndexes)}\n" + inputhookcode += " ftuparameterhooks={" + if len(ftuidmap)>0: + for k,v in ftuidmap.items(): + inputhookcode += f"\n '{k}':'{v}'," + inputhookcode += "\n },\n #PHS Input name and variable entry\n" inputhookcode += " inputhooks={" for k,v in inputhook.items(): inputhookcode += f"\n '{k}':'{v}'," @@ -378,7 +412,22 @@ def exportAsODEStepper(composer,modelName): for k1,v1 in stateNameToVectorIndexMap.items(): inputhookcode += f"\n '{k1}':[{','.join(map(str,v1))}]," inputhookcode += "\n }\n" - + if len(phsconstants)>0: + inputhookcode += " #Node, PHS parameter name and variable link\n phsparameterhooks={" + ngroup = {} + for k,v in phsconstants.items(): + nm = k.name.split('_') + grp = int(nm[-1]) + gname = '_'.join(nm[:-1]) + if not grp in ngroup: + ngroup[grp] = {} + ngroup[grp][gname] = arraymapping[k.name] + for n,v in ngroup.items(): + inputhookcode += f"\n '{n}':{{" + for g,a in v.items(): + inputhookcode += f"\n '{g}':'{a}'," + inputhookcode += f"\n }}," + inputhookcode += "\n }\n" phsIndexMap = dict() #Node ordering is in composer.inodeIndexes for si, k in enumerate(composer.inodeIndexes): @@ -420,6 +469,12 @@ def exportAsODEStepper(composer,modelName): for k1,v1 in svmap.items(): inputhookcode += f" '{k1}':{v1},\n" inputhookcode += f" }},\n" + if not composer.substituteParameters: + inputhookcode += f" 'parametervarmap':{{\n" + for k in composer.nodePHSData[k].parameters: + inputhookcode += f" '{k.name}':'{arraymapping[k.name]}',\n" + inputhookcode += f" }},\n" + inputhookcode += f" 'inputExpressions':{{\n" for k1,v1 in inexpr.items(): inputhookcode += f" '{k1}':{{'statevecindex':{v1[0]},'expr':'{v1[1]}'}},\n" @@ -433,12 +488,16 @@ def exportAsODEStepper(composer,modelName): from numpy import exp from scipy.integrate import ode -def Heaviside(x): +__version__ = "0.0.1" + +def heaviside(x): if x > 0: return 1.0 return 0.0 -__version__ = "0.0.1" +def Abs(x): + return np.fabs(x) + {inputhookcode} @@ -467,7 +526,12 @@ def __init__(self): except: stmt = f" {arraymapping[k.name]} = {v['value']} #{k}\n" pycode += stmt - + for k, v in ubaridxmap.items(): + pycode += f" {v} = 0.0 #{k} External input\n" + + if len(ftuidmap)>0: + for k,v in ftuidmap.items(): + pycode += f" {v} = 1.0 #{k} This needs to be set for accurate simulation\n" # Do constant subs definedVariables = [] for k, v in constantsubs.items(): @@ -479,6 +543,16 @@ def __init__(self): stmt = f" {arraymapping[v.name]} = {k} #{v}\n" pycode += stmt definedVariables.append(v) + for k, v in phsconstants.items(): + if k not in definedVariables: + if k.name in arraymapping: + try: + stmt = f" {arraymapping[k.name]} = {float(v['value']):6f} #{k}\n" + except: + stmt = f" {arraymapping[k.name]} = {v['value']} #{k}\n" + pycode += stmt + definedVariables.append(k) + pycode += "\n def compute_variables(self,voi):\n t=voi #mapping to t\n states, rates, variables = self.states,self.rates,self.variables\n" # Do uCap terms diff --git a/src/ftuutils/compositionutils.py b/src/ftuutils/compositionutils.py index 7a63699..ce52d44 100644 --- a/src/ftuutils/compositionutils.py +++ b/src/ftuutils/compositionutils.py @@ -490,6 +490,16 @@ def __str__(self): """ return res +def getNumerics(elem): + numericconstants = [] + if len(elem.free_symbols)>0: #numbers also get passed as expression, but will not have any free_symbols + for tm in elem.args: + numericconstants.extend(getNumerics(tm)) + elif isinstance(elem, sympy.Number): + if np.fabs(float(elem)) != 1.0: + numericconstants.append(elem) + return numericconstants + class Composer: """ @@ -525,7 +535,7 @@ def __init__(self) -> None: False # Will be set to true if it was computed by compose ) self.composition = None #Store composition information if used to construct the composer - + self.substituteParameters = True def setConnectivityGraph(self,graph:nx.Graph): """Set the base ftugraph from which the composition was created @@ -813,8 +823,28 @@ def compose(self, substituteParameters=True): bcsize = 0 networkLap = dict() networkAdj = dict() + self.substituteParameters = substituteParameters + self.requiredFTUparameters = False for n, g in self.networkGraphs.items(): - networkLap[n] = nx.laplacian_matrix(g.to_undirected()).todense() + try: + networkLap[n] = nx.laplacian_matrix(g.to_undirected()).todense() + except: + self.requiredFTUparameters = True + ug = g.to_undirected() + lap = sympy.zeros(ug.number_of_nodes(),ug.number_of_nodes()) + nodes = list(ug.nodes()) + for i,nd in enumerate(nodes): + nedges = list(ug.edges(nd,data='weight')) + ndeg = len(nedges) + lap[i,i] = ndeg + for ed in nedges: + if ed[0]==nd: + tar = nodes.index(ed[1]) + else: + tar = nodes.index(ed[0]) + lap[i,tar] = -sympy.sympify(ed[2]) + networkLap[n] = lap + # Construct node connectivity matrix, should be consistent with KCL nadj = np.zeros((g.number_of_nodes(), g.number_of_nodes())) for nd, nbrdict in g.adjacency(): @@ -1133,37 +1163,61 @@ def generatePythonIntermediates(self): rhsnumericconstants = [] self.raw_rhs = deepcopy(rhs) for elem in rhs: - if isinstance(elem, sympy.Expr): - for tm in elem.args: - rhsnumericconstants.extend( + nc = getNumerics(elem) + rhsnumericconstants.extend( [ term - for term in tm.args + for term in nc if term not in rhsfreesymbols and isinstance(term, sympy.Number) and np.fabs(float(term)) != 1.0 ] - ) - elif isinstance(elem, sympy.Number): - if np.fabs(float(elem)) != 1.0: - rhsnumericconstants.append(elem) + ) + + # Constants in uCapVec for i, s in enumerate(self.stateVec): elem = interioru[i] - if isinstance(elem, sympy.Expr): - for tm in elem.args: - rhsnumericconstants.extend( + nc = getNumerics(elem) + rhsnumericconstants.extend( [ term - for term in tm.args + for term in nc if term not in rhsfreesymbols and isinstance(term, sympy.Number) and np.fabs(float(term)) != 1.0 ] - ) - elif isinstance(elem, sympy.Number): - if np.fabs(float(elem)) != 1.0: - rhsnumericconstants.append(elem) + ) + + for v,ham in self.cellHamiltonians.items(): + nc = getNumerics(ham.hamiltonian) + rhsnumericconstants.extend( + [ + term + for term in nc + if term not in rhsfreesymbols + and isinstance(term, sympy.Number) + and np.fabs(float(term)) != 1.0 + ] + ) + #Find constants in composite parameters and get the list of nonlinear terms as well + nonlinearrhsterms = dict() + for c, t in self.compositeparameters.items(): + fs = t["value"].free_symbols + cvdict = dict() + if len(fs) > 0: #This is a nonlinearterm + nc = getNumerics(t["value"]) + rhsnumericconstants.extend( + [ + term + for term in nc + if term not in rhsfreesymbols + and isinstance(term, sympy.Number) + and np.fabs(float(term)) != 1.0 + ]) + nonlinearrhsterms[c] = t["value"] + else: + rhsnumericconstants.append(np.fabs(float(t["value"]))) constantsubs = dict() constCtr = 1 @@ -1171,48 +1225,63 @@ def generatePythonIntermediates(self): constantsubs[np.abs(c)] = f"c_{constCtr}" constCtr += 1 - # Remove constant entries that are same to given precision - constCtr = 1 - constantstoprecision = dict() - newkeys = dict() - for k, v in constantsubs.items(): - pk = f"{float(k):6f}" - if pk not in constantstoprecision: - constantstoprecision[pk] = sympy.Symbol(f"c_{constCtr}") - constCtr += 1 - newkeys[k] = constantstoprecision[pk] - - for k, v in newkeys.items(): - constantsubs[k] = v - + # # Remove constant entries that are same to given precision + # constCtr = 1 + # constantstoprecision = dict() + # newkeys = dict() + # for k, v in constantsubs.items(): + # pk = f"{float(k):6f}" + # if pk not in constantstoprecision: + # constantstoprecision[pk] = sympy.Symbol(f"c_{constCtr}") + # constCtr += 1 + # newkeys[k] = constantstoprecision[pk] + + # for k, v in newkeys.items(): + # constantsubs[k] = v + + #Convert to sympy Symbols for substitution + for k in constantsubs: + constantsubs[k] = sympy.Symbol(constantsubs[k]) + + for c in nonlinearrhsterms: + #Sympy handle Heaviside wierdly - so rename here and + vs = f"{nonlinearrhsterms[c]}" + if "Heaviside" in vs: + vn = sympy.sympify(vs.replace("Heaviside(","heaviside(")).xreplace(constantsubs) + vk = f"{vn}" + nonlinearrhsterms[c] = sympy.sympify(vk) + else: + v = nonlinearrhsterms[c].xreplace(constantsubs) + nonlinearrhsterms[c] = v + # Xreplace is faster than subs - no deep mathematical reasoning, ok for constant replacement cleanrhs = rhs.xreplace(constantsubs) cleaninputs = inputs.xreplace(constantsubs) # Generate python # Constants are contained in constantsubs, map and add all constants in compositeparameters that are used by functions in the composite parameters - nonlinearrhsterms = dict() - for c, t in self.compositeparameters.items(): - fs = t["value"].free_symbols - cvdict = dict() - if len(fs) > 0: - # Load the values of composite parameter parameters - for f in fs: - if f in self.compositeparameters: - if self.compositeparameters[f]["value"] in constantsubs: - cvdict[f] = constantsubs[ - self.compositeparameters[f]["value"] - ] - elif ( - np.fabs( - float(self.compositeparameters[f]["value"])) != 1.0 - ): - cvdict[f] = sympy.Symbol(f"c_{constCtr}") - constantsubs[self.compositeparameters[f]["value"]] = cvdict[ - f - ] - constCtr += 1 - nonlinearrhsterms[c] = t["value"].xreplace(cvdict) + # nonlinearrhsterms = dict() + # for c, t in self.compositeparameters.items(): + # fs = t["value"].free_symbols + # cvdict = dict() + # if len(fs) > 0: + # # Load the values of composite parameter parameters + # for f in fs: + # if f in self.compositeparameters: + # if self.compositeparameters[f]["value"] in constantsubs: + # cvdict[f] = constantsubs[ + # self.compositeparameters[f]["value"] + # ] + # elif ( + # np.fabs( + # float(self.compositeparameters[f]["value"])) != 1.0 + # ): + # cvdict[f] = sympy.Symbol(f"c_{constCtr}") + # constantsubs[self.compositeparameters[f]["value"]] = cvdict[ + # f + # ] + # constCtr += 1 + # nonlinearrhsterms[c] = t["value"].xreplace(cvdict) # Find all symbolic constants in nonlinearrhsterms for k, v in nonlinearrhsterms.items(): @@ -1223,41 +1292,54 @@ def generatePythonIntermediates(self): if self.compositeparameters[f]["value"] in constantsubs: cvdict[f] = constantsubs[self.compositeparameters[f]["value"]] elif np.fabs(float(self.compositeparameters[f]["value"])) != 1.0: - cvdict[f] = sympy.Symbol(f"c_{constCtr}") - constantsubs[self.compositeparameters[f] - ["value"]] = cvdict[f] + # cvdict[f] = sympy.Symbol(f"c_{constCtr}") + # constantsubs[self.compositeparameters[f] + # ["value"]] = cvdict[f] + #Maintain the same name + constantsubs[self.compositeparameters[f]["value"]] = f constCtr += 1 if len(cvdict) > 0: v = v.xreplace(cvdict) - # Store all numeric ones inside constantsubs - if isinstance(v, sympy.Expr): - for tm in v.args: - xtn = [ - term - for term in tm.args - if term not in fs and isinstance(term, sympy.Number) - ] - for c in xtn: - if np.abs(c) not in constantsubs and np.fabs(float(c)) != 1.0: - constantsubs[np.abs(c)] = sympy.Symbol( - f"c_{constCtr}") - constCtr += 1 - elif isinstance(v, sympy.Number): - if np.abs(v) not in constantsubs and np.fabs(float(v)) != 1.0: - constantsubs[np.abs(v)] = f"c_{constCtr}" - constCtr += 1 + # # Store all numeric ones inside constantsubs + # if isinstance(v, sympy.Expr): + # for tm in v.args: + # xtn = [ + # term + # for term in tm.args + # if term not in fs and isinstance(term, sympy.Number) + # ] + # for c in xtn: + # if np.abs(c) not in constantsubs and np.fabs(float(c)) != 1.0: + # constantsubs[np.abs(c)] = sympy.Symbol( + # f"c_{constCtr}") + # constCtr += 1 + # elif isinstance(v, sympy.Number): + # if np.abs(v) not in constantsubs and np.fabs(float(v)) != 1.0: + # constantsubs[np.abs(v)] = f"c_{constCtr}" + # constCtr += 1 # Remove constant entries that are same to given precision constCtr = 1 constantstoprecision = dict() newkeys = dict() skippedkeys = dict() + phsconstants = dict() + #Get all phs constants from composite parameters (all with numeric values) + #if not self.substituteParameters: + for k,v in self.compositeparameters.items(): + if len(v['value'].free_symbols)==0: + phsconstants[k] = v #float(v['value']) + for k, v in constantsubs.items(): pk = f"{float(k):6f}" if pk not in constantstoprecision: - constantstoprecision[pk] = sympy.Symbol(f"c_{constCtr}") + if v.name.startswith('c_'): + constantstoprecision[pk] = sympy.Symbol(f"c_{constCtr}") + else: + constantstoprecision[pk] = v constCtr += 1 - if v != constantstoprecision[pk]: + #Only for constant defined by the code and not phs constants + if v != constantstoprecision[pk] and not v in phsconstants: #v.name.startswith('c_'): skippedkeys[v] = constantstoprecision[pk] newkeys[k] = constantstoprecision[pk] @@ -1267,62 +1349,83 @@ def generatePythonIntermediates(self): # Handle elements that are functions, the multiplication operation with a state should be composition cleanedrhs = [] - for elem in cleanrhs: - newelem = 0 - # Logic works for Add expression - if not isinstance(elem, sympy.Add): - fs = elem.free_symbols - nterm = None - state = None - for f in fs: - if f in nonlinearrhsterms: - nterm = nonlinearrhsterms[f] - break - if nterm is not None: - # #Get the current state for the term - states = [] - for f in fs: + for relem in cleanrhs: + expandedelem = sympy.expand(relem) + #Look into each product term of a sum or product + # Sum need to be summed, product need to be multiplied + if isinstance(expandedelem,sympy.Add): + reducedelem = 0 + for elem in expandedelem.args: + #Each element is a product or free + estates = [] + enterms = [] + for f in elem.free_symbols: if f in stateVec.free_symbols: - states.append(f) - # If the state in term appears in nterm, divide - for st in states: - if st in nterm.free_symbols: - state = st - # handles just one! - break - if nterm is not None and state is not None: - newelem = elem / state - else: - newelem = term - cleanedrhs.append(sympy.simplify(newelem)) - else: - for term in elem.args: - fs = term.free_symbols - nterm = None - state = None - for f in fs: + estates.append(f) if f in nonlinearrhsterms: - nterm = nonlinearrhsterms[f] - break - if nterm is not None: - # #Get the current state for the term - states = [] - for f in fs: - if f in stateVec.free_symbols: - states.append(f) - # break - # If the state in term appears in nterm, divide - for st in states: - if st in nterm.free_symbols: - state = st - # handles just one! - break - if nterm is not None and state is not None: - nterm = term / state - newelem += nterm + enterms.append(f) + if len(estates)>0 and len(enterms)>0: + denom = 1 + for nt in enterms: + entf = nonlinearrhsterms[nt].free_symbols + for s in estates: + if s in entf: + denom *= s + reducedelem += sympy.simplify(elem/denom) else: - newelem += term - cleanedrhs.append(sympy.simplify(newelem)) + reducedelem += elem + # # Logic works for Multiplication expression + # newelem = 1 + # for term in elem.args: + # fs = term.free_symbols + # nterm = None + # state = None + # for f in fs: + # if f in nonlinearrhsterms: + # nterm = nonlinearrhsterms[f] + # break + # if nterm is not None: + # # #Get the current state for the term + # states = [] + # for f in fs: + # if f in stateVec.free_symbols: + # states.append(f) + # # break + # # If the state in term appears in nterm, divide + # for st in states: + # if st in nterm.free_symbols: + # state = st + # # handles just one! + # break + # if nterm is not None and state is not None: + # nterm = term / state + # newelem *= nterm + # else: + # newelem *= term + # reducedelem.append(sympy.simplify(newelem)) + cleanedrhs.append(reducedelem) + elif isinstance(expandedelem,sympy.Mul): # if its a product + reducedelem = 1 + estates = [] + enterms = [] + for f in expandedelem.free_symbols: + if f in stateVec.free_symbols: + estates.append(f) + if f in nonlinearrhsterms: + enterms.append(f) + if len(estates)>0 and len(enterms)>0: + denom = 1 + for nt in enterms: + entf = nonlinearrhsterms[nt].free_symbols + for s in estates: + if s in entf: + denom *= s + reducedelem *= sympy.simplify(expandedelem/denom) + else: + reducedelem *= expandedelem + cleanedrhs.append(reducedelem) + else: + cleanedrhs.append(relem) # Constants also contain u vector, however they are updated after each step # Do the update in compute_variables method, and initialise them in initialise_variables method @@ -1337,6 +1440,7 @@ def generatePythonIntermediates(self): arraymapping[s] = f"states[{i}]" invarraymapping[f"states[{i}]"] = s + numconstants = 0 # Do ubar first as numconstants change due to precision selection # ubar entries @@ -1347,6 +1451,16 @@ def generatePythonIntermediates(self): invarraymapping[f"variables[{numconstants}]"] = s.name ubaridxmap[s.name] = f"variables[{numconstants}]" numconstants += 1 + #Find any connectivity related symbols + ftuidmap = dict() + for s in interioru.free_symbols: + if not s.name.startswith("u_"): + arraysubs[s] = sympy.Symbol(f"variables[{numconstants}]") + arraymapping[s.name] = f"variables[{numconstants}]" + invarraymapping[f"variables[{numconstants}]"] = s.name + ftuidmap[s.name] = f"variables[{numconstants}]" + numconstants += 1 + # Multiple k's will have same v due to defined precision definedConstants = [] @@ -1357,7 +1471,32 @@ def generatePythonIntermediates(self): invarraymapping[f"variables[{numconstants}]"] = str(v) numconstants += 1 definedConstants.append(v) - + if not self.substituteParameters: + #Insert phs constants + for v,k in phsconstants.items(): + arraysubs[v] = sympy.Symbol(f"variables[{numconstants}]") + arraymapping[str(v)] = f"variables[{numconstants}]" + invarraymapping[f"variables[{numconstants}]"] = str(v) + numconstants += 1 + else: + #Reduce repeats + existingvalues = {} + newphsconstants = {} + for k,vdict in phsconstants.items(): + v = vdict['value'] + if float(v) in existingvalues: + arraysubs[k] = existingvalues[float(v)] + arraymapping[str(k)] = str(arraysubs[k]) + invarraymapping[arraymapping[str(k)] ] = str(k) + else: + arraysubs[k] = sympy.Symbol(f"variables[{numconstants}]") + arraymapping[str(k)] = f"variables[{numconstants}]" + invarraymapping[f"variables[{numconstants}]"] = str(k) + numconstants += 1 + existingvalues[v] = arraysubs[k] + newphsconstants[k] = vdict + phsconstants = newphsconstants + # uCap entries for s in self.stateVec: arraysubs[sympy.Symbol(f"u_{s}")] = sympy.Symbol( @@ -1401,11 +1540,11 @@ def generatePythonIntermediates(self): arrayedinputs = [] # Use cleanedrhs and not cleanrhs - nonlinear terms with functions are transformed to compositions and not multiplication for elem in cleanedrhs: - arrayedrhs.append(elem.xreplace(arraysubs)) + arrayedrhs.append(elem.xreplace(skippedkeys).xreplace(arraysubs)) for elem in cleaninputs: - arrayedinputs.append(elem.xreplace(arraysubs)) + arrayedinputs.append(elem.xreplace(skippedkeys).xreplace(arraysubs)) - return numconstants,constantsubs,nonlinearrhsterms,inputs,arrayedinputs,arraymapping,uCapterms,ucapdescriptive,nonlineararrayedrhsterms,nonlinearrhstermsdescriptive,arrayedrhs,invarraymapping,rhs,ubaridxmap,cleaninputs + return numconstants,phsconstants,constantsubs,nonlinearrhsterms,inputs,arrayedinputs,arraymapping,uCapterms,ucapdescriptive,nonlineararrayedrhsterms,nonlinearrhstermsdescriptive,arrayedrhs,invarraymapping,rhs,ubaridxmap,ftuidmap,cleaninputs def exportAsPython(self): """Export composed FTU as python code similar to OpenCOR export""" @@ -1513,9 +1652,17 @@ def ddecSetupOperators(self,weight='weight'): ew = v1[weight] break edgeWeights.append(ew) - sew = np.sqrt(ew) - adm[ei,nidx[ed[0]]] = -sew - adm[ei,nidx[ed[1]]] = +sew + if isinstance(ew, (int, float, complex)) and not isinstance(ew, bool): + sew = np.sqrt(ew) + adm[ei,nidx[ed[0]]] = -sew + adm[ei,nidx[ed[1]]] = +sew + else: + if adm.dtype=='float': + adm = np.zeros((numEdges,numNodes),dtype='str') + + adm[ei,nidx[ed[0]]] = f"-sqrt({ew})" + adm[ei,nidx[ed[1]]] = f"sqrt({ew})" + else: for ei,ed in enumerate(edges): ew = 1.0 @@ -1523,10 +1670,16 @@ def ddecSetupOperators(self,weight='weight'): if weight in ews: ew = ews[weight] edgeWeights.append(ew) - sew = np.sqrt(ew) + if isinstance(ew, (int, float, complex)) and not isinstance(ew, bool): + sew = np.sqrt(ew) + adm[ei,nidx[ed[0]]] = -sew + adm[ei,nidx[ed[1]]] = +sew + else: + if adm.dtype=='float': + adm = np.zeros((numEdges,numNodes),dtype='str') + adm[ei,nidx[ed[0]]] = f"-sqrt({ew})" + adm[ei,nidx[ed[1]]] = f"sqrt({ew})" - adm[ei,nidx[ed[0]]] = -sew - adm[ei,nidx[ed[1]]] = +sew return adm,edgeWeights diff --git a/src/ftuutils/simulationutils.py b/src/ftuutils/simulationutils.py index 695281d..62e903f 100644 --- a/src/ftuutils/simulationutils.py +++ b/src/ftuutils/simulationutils.py @@ -53,6 +53,24 @@ def visit_Name(self, node): return ast.Name(**{**node.__dict__, 'id':self.variableNameMap[node.id]}) else: return node + +def handleStars(code, numnodes): + """Replace code that has node[*]. with node[1..numnodes]. + + Args: + code (string): Python code + numnodes (int): Number of nodes + """ + cstrip = code.split('\n') + ncstrip = [] + for cs in cstrip: + if 'node[*]' in cs: + for ni in range(1,numnodes+1): + ncstrip.append(cs.replace('node[*]',f'node[{ni}]')) + else: + ncstrip.append(cs) + return '\n'.join(ncstrip) + def loadTemplateFile(filename): """Load template files stored as resources within ftuutils package""" @@ -95,22 +113,35 @@ def __init__(self,composer) -> None: code = f"{self.hamletcode}\nself.hamletNodes={self.modelname}Inputs.nodes\n" exec(compile(code,'','exec')) - + self.CELL_COUNT = self.inputhookInstance.CELL_COUNT #Number of nodes + self.variablemap = self.inputhookInstance.inputhooks for k,v in self.inputhookInstance.statehooks.items(): for sn,sm in v.items(): self.inputhookInstance.inputhooks[f"node[{k}].{sn}"] = sm - self.variablemap = self.inputhookInstance.inputhooks self.statenamehooks = self.inputhookInstance.statenamehooks self.phsnamehooks = self.inputhookInstance.phsnamehooks + self.phsparameterhooks = self.inputhookInstance.phsparameterhooks + for k,v in self.inputhookInstance.phsparameterhooks.items(): + for sn,sm in v.items(): + self.inputhookInstance.inputhooks[f"node[{k}].{sn}"] = sm + #Empty ftuparameterhooks def is returned as tuple ({},) + self.ftuparameterhooks = self.inputhookInstance.ftuparameterhooks + if isinstance(self.inputhookInstance.ftuparameterhooks,tuple): + self.ftuparameterhooks = self.inputhookInstance.ftuparameterhooks[0] + + for k,v in self.ftuparameterhooks.items(): + self.inputhookInstance.inputhooks[k] = v - def addExperiment(self,name,time,inputblock,preamble=""): + + def addExperiment(self,name,time,inputblock,parameterblock=None,preamble=""): """Add an experiment for simulation Args: name (string): Name of the experiment time (list): Simulation time block associated with the experiment [start,end, [numsteps]] inputblock (string): Python code that will be executed to generate input signals for FTU simulation + parameterblock (string): Python code to set parameter values. Default None preamble (str, optional): Any preamble python code that will be added to start of the generated code; ideally used to import packages that are used by the input block. Defaults to "". Raises: @@ -122,7 +153,7 @@ def addExperiment(self,name,time,inputblock,preamble=""): if "return " in isplit[-1]: raise Exception(f"inputblock code should not return, rather assign to input variables!") - self.experiments[name] = {'time':time,'process_time_sensitive_events':inputblock,'preamble':preamble} + self.experiments[name] = {'time':time,'process_time_sensitive_events':inputblock,'preamble':preamble,'parameters':parameterblock} def generate(self,targetDir,provenance={},defaultnetworkid=1): """Generate code for the experiments and store them in the target directory @@ -158,7 +189,8 @@ def generate(self,targetDir,provenance={},defaultnetworkid=1): numsteps = stepsize #Process the input code block to - eventCodex = ast.unparse(refactor.visit(ast.parse(v['process_time_sensitive_events']))).strip() + evcode = handleStars(v['process_time_sensitive_events'],self.CELL_COUNT) + eventCodex = ast.unparse(refactor.visit(ast.parse(evcode))).strip() cx = eventCodex.split('\n') #Indentation should match that of 'process_time_sensitive_events' function def indent = " " @@ -168,7 +200,27 @@ def generate(self,targetDir,provenance={},defaultnetworkid=1): eventCode = "" for c in cx: eventCode += f"{indent}{c}\n" - + + parameterupdates = '' + pblock = v['parameters'] + if not pblock is None: + pvcode = handleStars(pblock,self.CELL_COUNT) + pCodex = ast.unparse(refactor.visit(ast.parse(pvcode))).strip() + #This code will be in the __init_ block, variables, states and rates will be with respect to the instance + #so add self. prefix + pselfcx = pCodex.replace("variables","self.variables").replace("states","self.states").replace("rates","self.rates") + pcx = pselfcx.split('\n') + pvx = ast.unparse(ast.parse(pvcode)).strip().split('\n') #Get the statements for comments + #Indentation should match that of '__init__' function def + indent = " " + if pselfcx.startswith("def"): + pcx = pcx[1:] + pvx = pvx[1:] + indent = " " + parameterupdates = f"{indent}#Experiment specific parameters setting starts\n" + for ix,pc in enumerate(pcx): + parameterupdates += f"{indent}{pc}{indent}#{pvx[ix]}\n" + parameterupdates += f"{indent}#Experiment specific parameters setting ends" code = v['preamble'] if len(code)>0: code +='\n' @@ -184,6 +236,7 @@ class {self.modelname}_{k}({self.modelname}): """ def __init__(self) -> None: super().__init__() +{parameterupdates} self.cellHam = np.zeros(self.CELL_COUNT) self.energyInputs = np.zeros(self.CELL_COUNT) self.totalEnergyInputs = np.zeros(self.CELL_COUNT) diff --git a/tests/ftuworkflow.ipynb b/tests/ftuworkflow.ipynb index c0b95a0..58693e4 100644 --- a/tests/ftuworkflow.ipynb +++ b/tests/ftuworkflow.ipynb @@ -86,6 +86,9 @@ "from ftuutils.base import FTUDelaunayGraph \n", "#Set edge weights based on orientation with respect to conductivity tensor\n", "conductivity_tensor = np.array([1.2,0.9,0.5]) #Fibre sheet normal\n", + "#Parameterised conductivity tensor can also be provided. \n", + "#These values can be set at the time of creating an experiment\n", + "conductivity_tensor = np.array(['Df','Ds','Dn'],dtype=str) #Fibre sheet normal\n", "g = FTUDelaunayGraph(points,\"APN\",conductivity_tensor)" ] }, @@ -213,7 +216,15 @@ "G = g.getGraph()\n", "\n", "#Call the FTU composition logic to create a FTU with above information\n", - "composer = g.composeCompositePHS(G,phstypes,phsdata)" + "\n", + "#composer = g.composeCompositePHS(G,phstypes,phsdata)\n", + "\n", + "#The above call will create a composite phs, whose parameters are substituted in the final\n", + "#expression. Use this approach when the PHS parameters will not be changed to explored in the experiments\n", + "#When experiments with differing phs parameters need to created, the composite PHS\n", + "#can be created such that the parameters are not substituted at build time but resolved at runtime\n", + "#For such approaches use\n", + "composer = g.composeCompositePHS(G,phstypes,phsdata,substituteParameters=False)" ] }, { @@ -265,21 +276,28 @@ "For instance, in the above case\n", "\n", "input is `i_1`\n", + "\n", "the states are `Tai, Ta, u, v`\n", + "\n", "However since states are associated with cells/graph nodes. They need to be prefixed with their node number, for instance\n", "`node[1].Ta`\n", "\n", "The above refers to state value `Ta`, belonging to node[1], here `1` is the label assigned to the node in the graph generation process.\n", "\n", + "Similarly if the composite PHS is composed such that the individual phs parameters are not substituted at composition time, those parameters are also made available using the above syntax, for instance\n", + "`node[3].eps0`.\n", + "\n", "Where there are more than one celltypes, say `APN`, `FHN`...\n", - "The states names are prefixed by the celltype, as `APN_Ta,APN_u`, `FHN_u, FHN_v`\n", + "The states names and phs parameter names are prefixed by the celltype, as `APN_Ta,APN_u`, `FHN_u, FHN_v`.\n", + "\n", + "To apply a variation to all nodes, the special operator `*` can be used. For instance, `node[*].Ta = 5.0` will set the state value `Ta` for all nodes to `5.0`.\n", "\n", "Below we create a single experiment, where the boundary cells are stimulated with a current of `0.5` units between `100