Merge pull request #42 from pyurdme/stochss

Stochss

build_docs/conf.py

@@ -56,7 +56,7 @@
 # The short X.Y version.
 version = '1.0'
 # The full version, including alpha/beta/rc tags.
-release = '1.0'
+release = '1.0.1'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

build_docs/index.rst

@@ -8,7 +8,7 @@ Welcome to PyURDME's documentation!
 
 PyURDME is a general software framework for modeling and simulation of stochastic reaction-diffusion processes on unstructured, tetrahedral (3D) and triangular (2D) meshes. Unstructured meshes allow for a more flexible handling of complex geometries compared to structured, Cartesian meshes. The current core simulation algorithm is based on the mesoscopic reaction-diffusion master equation (RDME) model.
 
-PyURDME was originally based on the URDME software package, it has been completely rewritten in python and updated.  It depends on the FEniCS libraries for mesh generation and assemply, see http://fenicsproject.org/.  The core simulation routines are implemetned in C, and requires GCC for compliation.  The default solver is an efficient implementation of the Next Subvolume Method (NSM).
+PyURDME was originally based on the URDME software package, it has been completely rewritten in python and updated.  It depends on the FEniCS libraries for mesh generation and assemply, see http://fenicsproject.org/.  The core simulation routines are implemented in C, and requires GCC for compliation.  The default solver is an efficient implementation of the Next Subvolume Method (NSM).
 
 Contents:
 

examples/simple_diffusion/circle.geo

@@ -0,0 +1,13 @@
+cl__1 = 1;
+Point(1) = {0, 0, 0, 1};
+Point(2) = {1, 0, 0, 1};
+Point(3) = {0, 1, 0, 1};
+Point(5) = {-1, 0, 0, 1};
+Point(6) = {0, -1, 0, 1};
+Circle(1) = {5, 1, 3};
+Circle(2) = {3, 1, 2};
+Circle(3) = {5, 1, 6};
+Circle(4) = {2, 1, 6};
+Line Loop(6) = {1, 2, 4, -3};
+Plane Surface(6) = {6};
+Coherence;

examples/simple_diffusion/simple_diffusion2.py

@@ -39,8 +39,7 @@ def __init__(self):
         self.add_species([A,B])
 
         # A circle
-        mesh = dolfin.UnitCircleMesh(20)
-        self.mesh = pyurdme.URDMEMesh(mesh)
+        self.mesh = pyurdme.URDMEMesh.read_dolfin_mesh("circle.xml")
 
         # A mesh function for the cells
         cell_function = dolfin.CellFunction("size_t",self.mesh)
@@ -77,28 +76,14 @@ def __init__(self):
 
     model = simple_diffusion2()
     result = model.run()
-    A = result.get_species("A")
-    #print numpy.sum(A,axis=1)
-    data = model.get_solver_datastructure()
-    u0 = model.u0
-    print numpy.sum(u0,axis=1)
-    ix = numpy.argmax(u0[0,:])
-
-    c = model.mesh.coordinates()
-    x = c[:,0]
-    print c[ix,:]
-    dof2vtx = dolfin.dof_to_vertex_map(model.mesh.get_function_space())
-    u0 = data["u0"]
-    ixdof = numpy.argmax(u0[0,:])
-    print ix, dof2vtx[ixdof]
-    print u0[0,ixdof]
-    print A[0,ix]
-    print A[1,ix]
-#print A0[ixdof]
 
-    #print numpy.max(x)
-    #print numpy.min(x)
-    # Dump timeseries in Paraview format
+    # Write mesh and subdomain files for the StochSS UI
+    sd = model.get_subdomain_vector()
+    with open("simple_diffusion_subdomains.txt",'w') as fd:
+        for ndx,val in enumerate(sd):
+            fd.write("{0},{1}\n".format(ndx,val))
+
+    # For visualization in Paraview
     result.export_to_vtk(species="B",folder_name="Bout")
     result.export_to_vtk(species="A",folder_name="Aout")
 

install_ubuntu.sh

@@ -10,5 +10,4 @@ add-apt-repository ppa:fenics-packages/fenics
 apt-get update
 apt-get -y install fenics
 apt-get -y install cython python-h5py
-apt-get -y install python-pip
-pip install https://www.github.com/ahellander/pyurdme/tarball/installer
+python setup.py install

pyurdme/model.py

@@ -326,7 +326,7 @@ def create_mass_action(self):
             raise ReactionError("Reaction: " +self.name + "A mass-action reaction cannot involve more than two species.")
 
         # Case EmptySet -> Y
-        propensity_function = 'return ' + self.marate.name;
+        propensity_function = self.marate.name;
 
         # There are only three ways to get 'total_stoch==2':
         for r in self.reactants:
@@ -345,7 +345,7 @@ def create_mass_action(self):
             propensity_function += "*vol"
 
 
-        self.propensity_function = propensity_function + ';'
+        self.propensity_function = "return " + propensity_function + ';'
 
     def set_type(self,type):
         if type not in {'mass-action','customized'}:

pyurdme/pyurdme.py

@@ -150,6 +150,7 @@ def __setstate__(self, state):
                     mesh.constrained_domain = compiled_object
                 if 'num_dof_voxels' in state['mesh']:
                     mesh.num_dof_voxels = state['mesh']['num_dof_voxels']
+                mesh.init(2, 0)
                 self.__dict__['mesh'] = mesh
             except Exception as e:
                 print "Error unpickling model, could not recreate the mesh."
@@ -377,7 +378,6 @@ def get_subdomain_vector(self, subdomains={}):
 
         if subdomains == {}:
             self.sd = sd
-            print subdomains
         else:
             for dim, subdomain in subdomains.items():
                 if dim == 0:
@@ -1240,6 +1240,10 @@ def __getstate__(self):
 
         state = self.__dict__
         state["filecontents"] = filecontents
+        state["v2d"] = self.get_v2d()
+        state["d2v"] = self.get_d2v()
+
+
         for key, item in state.items():
             try:
                 pickle.dumps(item)
@@ -1266,11 +1270,26 @@ def __setstate__(self, state):
         for k,v in state.items():
             self.__dict__[k] = v
 
+    def get_v2d(self):
+        """ Return the vertex-to-dof mapping. """
+        if not hasattr(self, 'v2d'):
+            fs = self.model.mesh.get_function_space()
+            self.v2d = dolfin.vertex_to_dof_map(fs)
+
+        return self.v2d
+
+    def get_d2v(self):
+        """ Return the dof-to-vertex mapping. """
+        if not hasattr(self, 'd2v'):
+            fs = self.model.mesh.get_function_space()
+            self.d2v = dolfin.dof_to_vertex_map(fs)
+
+        return self.d2v
+
     def _reorder_dof_to_voxel(self, M, num_species=None):
         """ Reorder the colums of M from dof ordering to vertex ordering. """
 
-        fs = self.model.mesh.get_function_space()
-        v2d = dolfin.vertex_to_dof_map(fs)
+        v2d = self.get_v2d()
         if len(M.shape) == 1:
             num_timepoints = 1
         else:
@@ -1425,8 +1444,8 @@ def _initialize_sol(self):
             raise URDMEError("URDMEResult.model must be set before the sol attribute can be accessed.")
         numvox = self.model.mesh.num_vertices()
         fs = self.model.mesh.get_function_space()
-        vertex_to_dof_map = dolfin.vertex_to_dof_map(fs)
-        dof_to_vertex_map = dolfin.dof_to_vertex_map(fs)
+        vertex_to_dof_map = self.get_v2d()
+        dof_to_vertex_map = self.get_d2v()
 
         # The result is loaded in dolfin Functions, one for each species and time point
         for i, spec in enumerate(self.model.listOfSpecies):
@@ -1608,8 +1627,7 @@ def _copynumber_to_concentration(self,copy_number_data):
             where the volume unit is defined by the user input.
         """
 
-        fs = self.model.mesh.get_function_space()
-        v2d = dolfin.vertex_to_dof_map(fs)
+        v2d = self.get_v2d()
         shape = numpy.shape(copy_number_data)
         if len(shape) == 1:
             shape = (1,shape[0])
@@ -1639,7 +1657,8 @@ def _compute_solution_colors(self,species, time_index):
         # Convert RGB to HEX
         colors= []
         for row in crgba:
-            colors.append(self._rgb_to_hex(tuple(list(row[1:]))))
+            # get R,G,B of RGBA
+            colors.append(self._rgb_to_hex(tuple(list(row[0:3]))))
 
         # Convert Hex to Decimal
         for i,c in enumerate(colors):
@@ -1912,15 +1931,14 @@ def run(self, number_of_trajectories=1, seed=None, input_file=None, loaddata=Fal
             self.infile_name = input_file
             self.delete_infile = False
 
-        outfile = tempfile.NamedTemporaryFile(delete=False)
-        outfile.close()
-
         if not os.path.exists(self.infile_name):
             raise URDMEError("input file not found.")
 
         # Execute the solver
-        urdme_solver_cmd = [self.solver_dir + self.propfilename + '.' + self.NAME, self.infile_name, outfile.name]
         for run_ndx in range(number_of_trajectories):
+            outfile = tempfile.NamedTemporaryFile(delete=False)
+            outfile.close()
+            urdme_solver_cmd = [self.solver_dir + self.propfilename + '.' + self.NAME, self.infile_name, outfile.name]
 
             if seed is not None:
                 urdme_solver_cmd.append(str(seed+run_ndx))

pyurdme/urdme/build/Makefile.nsm

@@ -21,7 +21,15 @@ DMAT = -DURDME_LIBMAT
 CC	   = mpicc
 LINKER = mpicc 
 
-CFLAGS = -gdwarf-2 -m64 -c  -O3 -std=gnu99 $(MEMOPTS) -Wall -Wlong-long -Wformat -Wpointer-arith  $(INCLUDE) $(DMAT)
+LBITS := $(shell getconf LONG_BIT)
+ifeq ($(LBITS),64)
+    # 64 bit stuff
+    CFLAGS = -gdwarf-2 -m64 -c  -O3 -std=gnu99 $(MEMOPTS) -Wall -Wlong-long -Wformat -Wpointer-arith  $(INCLUDE) $(DMAT)
+else
+    # 32 bit stuff
+    CFLAGS = -gdwarf-2 -m32 -c  -O3 -std=gnu99 $(MEMOPTS) -Wall -Wlong-long -Wformat -Wpointer-arith  $(INCLUDE) $(DMAT)
+endif
+
 
 LFLAGS_BASE= -l hdf5 -l hdf5_hl 
 LFLAGS = $(LFLAGS_BASE) -lm