Parametric experimental design

src/ftuutils/base.py

@@ -668,6 +668,7 @@ def composeCompositePHS(self,G,phsdefinitions=None,phsdata=None,substituteParame
             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
@@ -786,15 +787,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 conductivityTensor.dtype!='<U1':
+            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,phsconstants,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,phsconstants,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()
@@ -348,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 += "    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}',"
@@ -380,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):
@@ -441,12 +488,16 @@ def exportAsODEStepper(composer,modelName):
 from numpy import exp
 from scipy.integrate import ode
 
+__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}
 
@@ -475,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():
@@ -487,15 +543,15 @@ def __init__(self):
                     stmt = f"        {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:
+    for k, v in phsconstants.items():
+        if k not in definedVariables:
+            if k.name in arraymapping:
                 try:
-                    stmt = f"        {arraymapping[v.name]} = {float(k):6f}  #{v}\n"
+                    stmt = f"        {arraymapping[k.name]} = {float(v['value']):6f}  #{k}\n"
                 except:
-                    stmt = f"        {arraymapping[v.name]} = {k}  #{v}\n"
+                    stmt = f"        {arraymapping[k.name]} = {v['value']}  #{k}\n"
                 pycode += stmt
-                definedVariables.append(v)
+                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"

src/ftuutils/compositionutils.py

@@ -824,8 +824,27 @@ def compose(self, substituteParameters=True):
         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():
@@ -1220,7 +1239,7 @@ def generatePythonIntermediates(self):
         # for k, v in newkeys.items():
         #     constantsubs[k] = v
 
-        #Goto sympy 
+        #Convert to sympy Symbols for substitution
         for k in constantsubs:
             constantsubs[k] = sympy.Symbol(constantsubs[k])
 
@@ -1309,13 +1328,9 @@ def generatePythonIntermediates(self):
         #if not self.substituteParameters:
         for k,v in self.compositeparameters.items():
             if len(v['value'].free_symbols)==0:
-                phsconstants[k] = float(v['value'])
+                phsconstants[k] = v #float(v['value'])
 
         for k, v in constantsubs.items():
-            #Phs named constants may get elimimated if they share similar values, so store
-            # if not v.name.startswith('c_'):
-            #    phsconstants[v] = k
-
             pk = f"{float(k):6f}"
             if pk not in constantstoprecision:
                 if v.name.startswith('c_'):
@@ -1425,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
@@ -1435,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 = []
@@ -1456,9 +1482,10 @@ def generatePythonIntermediates(self):
             #Reduce repeats
             existingvalues = {}
             newphsconstants = {}
-            for k,v in phsconstants.items():
-                if v in existingvalues:
-                    arraysubs[k] = existingvalues[v]
+            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:
@@ -1467,7 +1494,7 @@ def generatePythonIntermediates(self):
                     invarraymapping[f"variables[{numconstants}]"] = str(k)
                     numconstants += 1    
                     existingvalues[v] = arraysubs[k]  
-                    newphsconstants[k] = v            
+                    newphsconstants[k] = vdict            
             phsconstants = newphsconstants
 
         # uCap entries
@@ -1517,7 +1544,7 @@ def generatePythonIntermediates(self):
         for elem in cleaninputs:
             arrayedinputs.append(elem.xreplace(skippedkeys).xreplace(arraysubs))
 
-        return numconstants,phsconstants,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"""