Merge pull request #39 from lewisamarshall/development

Serialization Fixes

.gitignore

@@ -59,3 +59,6 @@ target/
 
 # Ignore pycharm
 .idea
+
+# Ignore 
+.imdone

.travis.yml

@@ -15,7 +15,7 @@ before_install:
   - ./miniconda.sh -b
   - export PATH=/home/travis/miniconda/bin:$PATH
   - conda update --yes conda
-  - pip install coveralls
+  - pip install coveralls simplejson
   - mkdir examples
 install:
   - conda install --yes python=$TRAVIS_PYTHON_VERSION numpy scipy nose h5py coverage

benchmark.py

@@ -0,0 +1,31 @@
+import emigrate
+import ionize
+import cProfile
+
+solutions = [ionize.Solution(['hepes', 'tris'],
+                             [.05, .105]),
+             ionize.Solution(['caproic acid', 'tris', 'fluorescein',
+                              'alexa fluor 488', 'mops', 'acetic acid'],
+                             [.01, .1, .01, .01, .01, .01]),
+             ionize.Solution(['hydrochloric acid', 'tris'],
+                             [.04, .08]),
+             ]
+
+system = emigrate.Frame(dict(
+                         n_nodes=137,
+                         lengths=[.005, .001, .02],
+                         interface_length=.0005,
+                         solutions=solutions,
+                         current=-500.,
+                         ))
+
+
+solver = emigrate.Solver(system,
+                         filename='examples/benchmark.hdf5',
+                         precondition=True,
+                         flux_mode='slip')
+tmax = 200
+dt = 1
+ode_solver = 'dopri5'
+profile = True
+cProfile.run("solver.solve(dt, tmax)")

emigrate/Frame.py

@@ -2,6 +2,11 @@
 from scipy.special import erf
 import ionize
 import warnings
+try:
+    import simplejson as json
+except:
+    import json
+import ionize
 
 
 class Frame(object):
@@ -27,12 +32,36 @@ class Frame(object):
     pressure = 0
     bulk_flow = 0
 
+    # I/O
+    def _encode(self, obj):
+        if isinstance(obj, ionize.Ion):
+            ion = obj.serialize()
+            ion.update({'__ion__': True})
+            return ion
+        elif isinstance(obj, np.ndarray):
+            return {'__ndarray__': True, 'data': obj.tolist()}
+        # Let the base class default method raise the TypeError
+        return json.JSONEncoder.default(self, obj)
+
+    def _object_hook(self, obj):
+        if '__ion__' in obj:
+            return ionize.Ion(1, 1, 1, 1).deserialize(obj)
+        elif '__ndarray__' in obj:
+            return np.array(obj['data'])
+        return obj
+
     def __init__(self, constructor):
         """Initialize a Frame object."""
 
         # Try loading from file first.
         if isinstance(constructor, basestring):
-            self.open_file(constructor)
+            self.deserialize(constructor)
+            try:
+
+                self.deserialize(constructor)
+            except Exception as e:
+                raise e
+                # self.open_file(constructor)
 
         # Next look for a construction dictionary
         elif 'solutions' in constructor.keys():
@@ -44,45 +73,11 @@ def __init__(self, constructor):
 
         self._resolve_current()
 
-    def deserialize(self, constructor):
-        """Use a dictionary to update electrolyte properties."""
-        # Initialize optional properties
-        self.ions = constructor['ions']
-        for key, value in constructor.items():
-            if key in self.__class__.__dict__.keys():
-                if key in ['nodes', 'pH']:
-                    value = np.array(value)
-                elif key == 'concentrations':
-                    if isinstance(value, dict):
-                        value = np.array([value[i] for i in self.ions])
-                    else:
-                        value = np.array(value)
-                # print key, value
-                self.__dict__[key] = value
-            else:
-                warnings.warn('Unused initialization key: {}'.format(key))
-
     def serialize(self):
-        serial = self.__dict__.copy()
-        try:
-            serial['ions'] = [ion.name for ion in self.ions]
-        except:
-            serial['ions'] = self.ions
-        serial['nodes'] = self.nodes.tolist()
-        serial['concentrations'] = {i: c for i, c in
-                                    zip(serial['ions'],
-                                        self.concentrations.tolist()
-                                        )
-                                    }
-        serial['concentrations']['x'] = serial['nodes']
-        try:
-            serial['pH'] = serial['pH'].tolist()
-        except:
-            pass
-        serial['properties'] = {'pH': serial['pH']}
-        serial['properties']['x'] = serial['nodes']
-
-        return serial
+        return json.dumps(self.__dict__, default=self._encode)
+
+    def deserialize(self, source):
+        self.__dict__ = json.loads(source, object_hook=self._object_hook)
 
     def construct(self, constructor_input):
         """Construct electrophoretic system based on a set of solutions."""
@@ -158,9 +153,9 @@ def _create_concentrations(self, constructor):
                                              right_side +
                                              constructor['lengths'][idx])
                     ion_concentration += ((cs[idx]-cs[idx-1]) *
-                                          (erf((self.nodes-left_side)
-                                           / constructor['interface_length'])
-                                           / 2. + .5))
+                                          (erf((self.nodes-left_side) /
+                                           constructor['interface_length']) /
+                                           2. + .5))
 
             self.concentrations.append(ion_concentration)
         self.concentrations = np.array(self.concentrations)
@@ -170,7 +165,7 @@ def _resolve_current(self):
             cd = self.current/self.area
             if self.current_density and not self.current_density == cd:
                 warnings.warn('Current and current density do not match. '
-                              + 'Correcting current density.')
+                              'Correcting current density.')
             self.current_density = cd
         elif self.current_density:
             self.current = self.area * self.current_density
@@ -184,11 +179,11 @@ def open_file(self, filename):
     my_solutions = [ionize.Solution(['hydrochloric acid', 'tris'], [.05, .1]),
                     ionize.Solution(['caproic acid', 'tris'], [.05, .1])
                     ]
-    system = Electrolyte({'lengths': [0.01]*2,
-                         'solutions': my_solutions})
+    system = Frame({'lengths': [0.01]*2,
+                    'solutions': my_solutions})
     print system.concentrations
-    print Electrolyte(system.serialize()).serialize()['concentrations'].keys()
-    print Electrolyte(system.serialize()).concentrations
+    print Frame(system.serialize())
+    print Frame(system.serialize()).__dict__
     # # print system.domain
     # print system.ions
     # print system.concentrations

emigrate/Solver.py

@@ -50,7 +50,6 @@ class Solver(object):
     equilibrator = None
 
     # #TODO:50 Remove these
-    area_variation = None
     ion_names = None
 
     def __init__(self, initial_condition, filename=False, precondition=False,
@@ -70,9 +69,7 @@ def __init__(self, initial_condition, filename=False, precondition=False,
 
         # Precondition if requested.
         if precondition:
-            preconditioner(self.state,
-                           self.fluxer,
-                           self.area_variation)
+            preconditioner(self.state, self.fluxer)
 
         # Create empty solution dictionaries
         self.ion_names = [ion.name for ion in self.initial_condition.ions]
@@ -115,7 +112,7 @@ def _write_solution(self):
         """Write the current state to solutions."""
         self.frame_series.add_frame(self.time, self.state)
 
-    def solve(self, interval=1., max_time=10, method='dopri5'):
+    def solve(self, interval=1., max_time=10):
         """Solve for a series of time points using an ODE solver."""
         if max_time is None:
             raise RuntimeError('Solving requires a finite maximum time.')
@@ -124,7 +121,7 @@ def solve(self, interval=1., max_time=10, method='dopri5'):
             pass
         return self.frame_series
 
-    def iterate(self, interval=1., max_time=None, method='dopri5'):
+    def iterate(self, interval=1., max_time=None):
         self._initialize_solver()
         while self.solver.successful():
             if self.solver.t >= max_time and max_time is not None:

emigrate/equilibration_schemes/VariablepH.py

@@ -157,6 +157,7 @@ def _calc_mobility(self):
 
     def _calc_diffusivity(self):
         """Calculate diffusivity."""
+        # TODO: Check if this works for low ionization fraction
         self.state.diffusivity = (self.state.absolute_mobility[:,
                                                                :, np.newaxis] *
                                   self.state.ionization_fraction /

emigrate/flux_schemes/Compact.py

@@ -30,7 +30,6 @@ def set_current(self):
 
     def set_E(self):
         """Calculate the electric field at each node."""
-        self.set_current()
         self.state.field = -self.state.current_density/self.conductivity()
 
     def diffusive_flux(self):

emigrate/flux_schemes/Differentiate.py

@@ -130,41 +130,3 @@ def construct_matrix(self, boundary_functions,
         construct = sp.csc_matrix(construct)
 
         return construct
-
-if __name__ == '__main__':
-    from scipy.special import erf
-    from matplotlib import pyplot as plot
-    Nt = 50
-    z = np.linspace(-1, 1, Nt)
-    test_functions = np.array([19*erf(z*3), 25*erf(z*2), 15*(erf(z*2) +
-                               .3*np.random.random(z.shape))/(10*z**2+1)])
-    my_diff = Differentiate(Nt, 1, method='6th-Order')
-    # my_diff = Differentiate(Nt, 1, method='dissipative')
-    # print my_diff.A1.todense(), '\n'
-    # print my_diff.B1.todense()
-
-    if False:
-        plot.plot(z, test_functions)
-        plot.show()
-
-    if False:
-        d1 = my_diff.first_derivative(test_functions.T)
-        d2 = my_diff.second_derivative(test_functions.T)
-
-    if True:
-        print my_diff.fA2
-        print my_diff.fM
-        print test_functions[2, :].shape
-        smoothed = my_diff.smooth(test_functions[2, :])
-        print smoothed.shape
-
-    if True:
-        for i in range(3):
-            plot.plot(z, np.ravel(test_functions[i, :]), label='test-function')
-        # plot.plot(z, np.ravel(d2[:, 0]))
-        plot.plot(z, np.ravel(smoothed), label='smoothed')
-        # plot.plot(z, np.ravel(d2[:, 1]))
-        # plot.plot(z, np.ravel(d1[:, 0]))
-        # plot.plot(z, np.ravel(d1[:, 1]))
-        plot.legend(loc="upper left")
-        plot.show()

emigrate/flux_schemes/SLIP.py

@@ -2,7 +2,6 @@
 import numpy as np
 from Fluxer import Fluxer
 from Flux_Limiter import Flux_Limiter
-# pylint: disable = W0232, E1101
 
 
 class SLIP(Fluxer):
@@ -40,8 +39,7 @@ def set_area_flux(self):
 
     def set_node_flux(self):
         """Calculate the flux of nodes."""
-        flux = self.first_derivative(self.node_cost() *
-                                     self.xz)
+        flux = self.first_derivative(self.node_cost() * self.xz)
         flux = self.differ.smooth(flux)
         flux *= self.pointwave
         flux[0, ] = flux[-1, ] = 0.
@@ -80,8 +78,8 @@ def electromigration_flux(self):
     def limit(self, concentrations, flux):
         self.set_alpha()
         limited_flux = 0.5*(flux[:, 3:] + flux[:, :-3]) +\
-            self.alpha[:, 1:-2] * (np.diff(concentrations, 1)[:, 1:-1]
-                                   - self.limiter(concentrations))
+            self.alpha[:, 1:-2] * (np.diff(concentrations, 1)[:, 1:-1] -
+                                   self.limiter(concentrations))
         return limited_flux
 
     def set_alpha(self):
@@ -108,18 +106,15 @@ def limiter(self, c):
 
     def node_cost(self):
         """Calculate the cost function of each node."""
-        self.set_Kag()
         deriv = np.fabs(self.cz)
-        if np.isnan(deriv).any():
-            print deriv
         cost = deriv / np.nanmax(deriv, 1)[:, np.newaxis]
-        cost = np.nanmax(cost, 0) + self.Kag
+        cost = np.nanmax(cost, 0) + self.Kag()
         return cost
 
-    def set_Kag(self):
+    def Kag(self):
         """Set the Kag parameter for spacing of low-gradient grid points."""
-        self.Kag = ((self.N-self.NI) / self.NI *
-                    self.Vthermal / abs(self.state.voltage))
+        return ((self.N-self.NI) / self.NI *
+                self.Vthermal / abs(self.state.voltage))
 
     # Helper Functions
     def set_derivatives(self):
@@ -156,8 +151,8 @@ def set_E(self):
 
     def conductivity(self):
         """Calculate the conductivty at each location."""
-        conductivity = np.sum(self.state.molar_conductivity
-                              * self.state.concentrations, 0)
+        conductivity = np.sum(self.state.molar_conductivity *
+                              self.state.concentrations, 0)
         conductivity += self.state.water_conductivity
         return conductivity
 

emigrate/preconditioner.py

@@ -2,16 +2,12 @@
 import numpy as np
 
 
-def preconditioner(state, fluxer, area_variation):
+def preconditioner(state, fluxer):
     """Precondition the system by spacing the grid points."""
     # set up the interpolator to get the new parameters
     concentration_interpolator = interp1d(state.nodes,
                                           state.concentrations,
                                           kind='cubic')
-    if area_variation:
-        area_interpolator = interp1d(state.nodes,
-                                     state.area,
-                                     kind='cubic')
 
     # Get the node cost from the flux calculator
     fluxer.update()
@@ -24,7 +20,10 @@ def preconditioner(state, fluxer, area_variation):
     # get the new grid parameters
     state.nodes = _precondition_x(state, cost)
     state.concentrations = concentration_interpolator(state.nodes)
-    if area_variation:
+    if np.size(state.area) > 1:
+        area_interpolator = interp1d(state.nodes,
+                                     state.area,
+                                     kind='cubic')
         state.area = area_interpolator(state.nodes)