Merge pull request #5 from Jagirhussan/main

Parametric FTU construction

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"]

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

src/ftuutils/codegenerationutils.py

@@ -28,14 +28,21 @@ 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'
     variabledescription += "STATE_INFO = [\n"
     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