Parametric expt design example and bug fixes

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
@@ -787,7 +790,7 @@ 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
-        if conductivityTensor.dtype!='<U1':
+        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])    

src/ftuutils/compositionutils.py

@@ -1652,20 +1652,34 @@ 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     
                 ews = self.ftugraph.get_edge_data(ed[0],ed[1])
                 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
 

src/ftuutils/simulationutils.py

@@ -205,17 +205,22 @@ def generate(self,targetDir,provenance={},defaultnetworkid=1):
             pblock = v['parameters']
             if not pblock is None:
                 pvcode = handleStars(pblock,self.CELL_COUNT)
-                pCodex = ast.unparse(refactor.visit(ast.parse(pvcode))).strip()            
-                pcx = pCodex.split('\n')
+                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 pCodex.startswith("def"):
+                if pselfcx.startswith("def"):
                     pcx = pcx[1:]
+                    pvx = pvx[1:]
                     indent = "    "
-                parameterupdates = "" 
-                for pc in pcx:
-                    parameterupdates += f"{indent}{pc}\n"
-
+                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'
@@ -231,7 +236,7 @@ class {self.modelname}_{k}({self.modelname}):
     """
     def __init__(self) -> None:
         super().__init__()
-        {parameterupdates}
+{parameterupdates}
         self.cellHam = np.zeros(self.CELL_COUNT)
         self.energyInputs = np.zeros(self.CELL_COUNT)
         self.totalEnergyInputs = np.zeros(self.CELL_COUNT)

tests/ftuworkflow.ipynb

@@ -88,7 +88,7 @@
     "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",
+    "conductivity_tensor = np.array(['Df','Ds','Dn'],dtype=str) #Fibre sheet normal\n",
     "g = FTUDelaunayGraph(points,\"APN\",conductivity_tensor)"
    ]
   },
@@ -216,7 +216,9 @@
     "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",
+    "\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",
@@ -295,7 +297,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -315,7 +317,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -329,8 +331,9 @@
     "Df = 1.5\n",
     "Ds = 0.9\n",
     "Dn = 0.7\n",
+    "node[1].eps0 = 1.0\n",
     "\"\"\"\n",
-    "#exptdesigned.addExperiment('test',[0,400,400],expt1,parameterblock=pblock)"
+    "exptdesigned.addExperiment('test',[0,400,400],expt1,parameterblock=pblock) #Note that using the same experiment name will overwrite previous record"
    ]
   },
   {
@@ -342,7 +345,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -364,7 +367,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [
     {