Yi-Xin (Eshin) Liu avatar Yi-Xin (Eshin) Liu committed b72b600

Refactor ZODB storage scheme. Save space 500+ times.

Comments (0)

Files changed (5)

ngpy/ngofflattice_kooi.py

         self.n_SM = eval(ini.get('Data','n_SM'))
         self.r_seed = eval(ini.get('Data','r_seed'))
         self.r_test = eval(ini.get('Data','r_test'))
-        self.interval_sample = eval(ini.get('Data','interval_sample'))
+        self.interval_save = int(eval(ini.get('Data','interval_save')))
         self.interval_show = eval(ini.get('Data','interval_show'))
         # [Program]
         self.exe_name = ini.get('Program','exe_name')
     return False
 
 
-def particle_SM_nucleation(t,dn,r_test,r0,k,nu,pm,ps,pa,pi):
+def particle_SM_nucleation(t, dn, r_test, r0, k, nu, 
+                           pm, ps, pa, pi, max_try):
     n = int(dn)
-    x_low_limit = pm.o.x-pm.r
-    x_high_limit = pm.o.x+pm.r
-    y_low_limit = pm.o.y-pm.r
-    y_high_limit = pm.o.y+pm.r
+    x_low_limit = pm.o.x - pm.r
+    x_high_limit = pm.o.x + pm.r
+    y_low_limit = pm.o.y - pm.r
+    y_high_limit = pm.o.y + pm.r
+    p_new = []
     for i in xrange(n):
+        itry = 0
         while True:
-            x = np.random.uniform(x_low_limit,x_high_limit)
-            y = np.random.uniform(y_low_limit,y_high_limit)
-            test_pt = Vector2D(x,y)
+            x = np.random.uniform(x_low_limit, x_high_limit)
+            y = np.random.uniform(y_low_limit, y_high_limit)
+            test_pt = Vector2D(x, y)
             test_particle = Particle(o=test_pt,r=r_test)
             if pm.is_inner_point(test_pt):
-                #if not ps.is_inner_point(test_pt) and not is_in_particle_list (test_pt,pa) and not is_in_particle_list(test_pt,pi):
                 if (not test_particle.is_touch(ps) and not
                     is_touch_particle(test_particle,pa) and not
                     is_touch_particle(test_particle,pi)):
-                    pa.append(Particle(test_pt,r0,r0,t,k,nu))
+                    p = Particle(test_pt,r0,r0,t,k,nu)
+                    pa.append(p)
+                    p_new.append(p)
                     break
+            itry += 1
+            if itry > max_try:
+                break;
+    return p_new
 
 
 def particle_SM_growth(t,pm,ps,pa,pi):
     pao = pa[:]
     pio = pi[:]
+    p_stop = []
     for p in pao:
         p.grow_by_function(t)
         if (p.is_touch(ps) or is_touch_particle(p,pao) or
             is_touch_particle(p,pio)):
             pa.remove(p)
             pi.append(p)
+            p_stop.append(p)
+    return p_stop
 
 
 def init_draw(lx,ly,dx,dy,ps):
 dy     = %(dx)s
 dt     = 0.01 ; time step, in minute
 Nx     = 1/%(dt)s
-max_t  = 30.01
+max_t  = 100.01
 ;k_MA   = 0.25e5
 k_MA   = 0.012*%(dx)s/%(dt)s ; 0.26
-nu_MA  = 0.5
+nu_MA  = 1.5
 r0_SM  = %(dx)s/2 ; the size of new born nucleus is exactly one pixel
 k_SM   = 6.0 ; 0.3*dx/(Nx*dt) # nm/min
 nu_SM  = 1.0
 n_SM   = 6.0e-7 ; 50/(lx*ly)
 r_seed = 1000.0
 r_test = 20*%(r0_SM)s
-interval_sample = 1/%(dt)s
+interval_save = %(Nx)s/10
 interval_show = 3/%(dt)s
 
 [Program]
 exe_name = ngrun.py
-exe_dir = /export/home/lyx/opt/lyx/web/
-database = zconfig:///export/home/lyx/opt/lyx/web/zeo.conf
+exe_dir = /export/home/lyx/opt/lyx/ngpy/ngpy/
+database = zconfig:///export/home/lyx/opt/lyx/ngpy/ngpy/zeo.conf
 base_dir = .
 data_dir = 
 data_dir_suffix = 
-#!/usr/bin/env python
 # -*- coding: utf-8 -*-
 """
     ngrun.py
 from .vector2d import Vector2D
 from .ngofflattice_kooi import Param as FileParam
 from .ngofflattice_kooi import calc_num_nucleation
-from .ngofflattice_kooi import particle_SM_nucleation,particle_SM_growth
+from .ngofflattice_kooi import particle_SM_nucleation, particle_SM_growth
 
-from .ngzodb import connect_zodb,setup_simulation
+from .ngzodb import connect_zodb, setup_simulation
 from .ngutil import now2str
 
+def save_frame(frames, SM_list_g, index, t, pa, pi, pm, p_new, p_stop):
+    SM_list = [str(p.ID) for p in pa + pi]
+    frame = PersistentMapping({
+                        't':t,
+                        'r_MA':pm.r,
+                        'SM':PersistentList(SM_list),
+                            })
+    # Add new SM particles to the global list
+    for p in p_new:
+        SM_list_g[p.ID] = p
+    # Update SM particles in the global list which stops growing
+    for p in p_stop:
+        if SM_list_g.has_key(p.ID):
+            SM_list_g[p.ID].te = t
+    frames[index] = frame
+    transaction.commit()
+
+
 def ngrun(zodb_URI,sim_id):
     db = connect_zodb(zodb_URI)
     sim_uuid = uuid.UUID(sim_id)
     simulations = db['simulations']
     simulation = simulations[sim_uuid]
-    p = simulation['parameter']
+    param = simulation['parameter']
+    lx = param.lx
+    ly = param.ly
+    dx = param.dx
+    dy = param.dy
+    k_MA = param.k_MA
+    I_SM = param.n_SM
+    k_SM = param.k_SM
+    r0_SM = param.r0_SM
+    nu_SM = param.nu_SM
+    dt = param.dt
+    max_t = param.max_t
+    r_seed = param.r_seed
+    r_test = param.r_test
+    interval_save = param.interval_save
 
     #UPDATE and ABORT are not handled Currently
     #Only NEW and UPDATE simulation can be run
     if simulation['status'] not in ('NEW','UPDATE'):
         return
 
-    o_M = Vector2D(p.lx/2,p.ly/2)
-    particle_MA = Particle(o_M,p.r_seed,p.r_seed,0,p.k_MA,p.nu_MA)
-    particle_seed = Particle(o_M,p.r_seed,p.r_seed,0,0,0)
+    o_M = Vector2D(param.lx/2, param.ly/2)
+    particle_MA = Particle(o_M, param.r_seed, param.r_seed, 
+                           0, param.k_MA, param.nu_MA)
+    particle_seed = Particle(o_M, param.r_seed, param.r_seed, 0, 0, 0)
     particle_SM_active = []
     particle_SM_inactive = []
     area_untransformed = []
     simulation['run_time'] = now2str()
     simulation['status'] = 'ACTIVE'
     simulation['particle_seed'] = particle_seed
+    simulation['particle_MA'] = particle_MA
+    simulation['particle_SM'] = PersistentMapping({})
+    ti = 2.0 * r_test / k_MA
+    particle_MA.r = r_seed + 2 * r_test
+    simulation['t_i'] = ti
     transaction.commit()
+    SM_list_global = simulation['particle_SM']
 
     if not simulation.has_key('frames'):
         simulation['frames'] = IOBTree.IOBTree()
         transaction.commit()
     frames = simulation['frames']
-    for i,t in np.ndenumerate(np.arange(p.dt,p.max_t+p.dt,p.dt)):
-        if 2 * particle_MA.r > p.lx:
+    for i,t in np.ndenumerate(np.arange(ti + param.dt, 
+                                        param.max_t + param.dt, param.dt)):
+        if 2 * particle_MA.r > lx or 2 * particle_MA.r > ly:
             break
         index, = i
 
-        I_SM = p.n_SM
-        dn = calc_num_nucleation(p.dt,I_SM,area_untransformed,
+        dn = calc_num_nucleation(dt, I_SM, area_untransformed,
                                  particle_MA,particle_seed,
                                  particle_SM_active,particle_SM_inactive)
+        N1 = len(particle_SM_active)
+        N2 = len(particle_SM_inactive)
         N = len(particle_SM_active) + len(particle_SM_inactive)
-        if dn > 0 and particle_MA.r - p.r_seed > 2 * p.r_test:
-            particle_SM_nucleation(t,dn,p.r_test,
-                                   p.r0_SM,p.k_SM,p.nu_SM,
-                                   particle_MA,particle_seed,
-                                   particle_SM_active,particle_SM_inactive)
+        print index,dn,N,N1,N2,particle_MA.r,particle_MA.r-r_seed - 2*r_test
+        p_new = []
+        if dn > 0 and particle_MA.r - r_seed > 2 * r_test:
+            max_try = int(lx*ly/(dx*dy))
+            p_new = particle_SM_nucleation(t, dn, r_test,
+                                    r0_SM, k_SM, nu_SM,
+                                    particle_MA, particle_seed,
+                                    particle_SM_active, 
+                                    particle_SM_inactive, max_try)
+            #print dn
+
+        dn = len(p_new) # the actual nucleated particle
         rr = 0.0
         for pc in particle_SM_active:
             rr += pc.r
-        rr *= p.k_SM
-        dr = (p.k_MA * p.dt - rr * p.dt / particle_MA.r -
-              dn * p.r0_SM * p.r0_SM / (2 * particle_MA.r))
+        rr *= k_SM
+        dr = k_MA * dt - rr * dt / particle_MA.r
+        if N > 0 and dn > 0:
+            dr -= dn * r0_SM * r0_SM / (2 * particle_MA.r)
         # The first occurence of nuclei may lead dr < 0,
         # which should be avoided
         if (N == 0 and dr > 0) or (N > 0):
             particle_MA.grow_by_dr(dr)
 
-        particle_SM_growth(t,particle_MA,particle_seed,
-                           particle_SM_active,particle_SM_inactive)
+        p_stop = particle_SM_growth(t, particle_MA, particle_seed,
+                                    particle_SM_active, 
+                                    particle_SM_inactive)
 
-        pa_list = IOBTree.IOBTree()
-        pi_list = IOBTree.IOBTree()
-        for i in range(len(particle_SM_active)):
-            pa_list[i] = copy.deepcopy(particle_SM_active[i])
-        for i in range(len(particle_SM_inactive)):
-            pi_list[i] = copy.deepcopy(particle_SM_inactive[i])
-        frame = PersistentMapping({
-            'particle_MA':copy.deepcopy(particle_MA),
-            'particle_SM_active':pa_list,
-            'particle_SM_inactive':pi_list
-            })
-        try:
-            frames[index] = frame
-            transaction.commit()
-        except KeyboardInterrupt:
-            transaction.abort()
-            simulation['status'] = 'ABORT'
-            simulation['abort_time'] = now2str()
-            transaction.commit()
-            sys.exit(0)
+        if index % interval_save == 0:
+            try:
+                save_frame(frames, SM_list_global, index, t,
+                           particle_SM_active, particle_SM_inactive,
+                           particle_MA, p_new, p_stop)
+            except KeyboardInterrupt:
+                transaction.abort()
+                simulation['status'] = 'ABORT'
+                simulation['abort_time'] = now2str()
+                transaction.commit()
+                sys.exit(0)
 
     simulation['status'] = 'FINISH'
     simulation['finish_time'] = now2str()
         simulation['abort_time'] = now2str()
         transaction.commit()
 
-
-if __name__ == '__main__':
-    params = FileParam('ngrc.ini')
-    db = connect_zodb(params.database)
-    sim_id = None
-    if not db.has_key('simulations'):
-        sim_id = setup_simulation(params)
-    simulations = db['simulations']
-    if not len(simulations):
-        sim_id = setup_simulation(params)
-    for k in simulations.keys():
-        simulation = simulations[k]
-        if simulation['status'] == 'NEW':
-            sim_id = k
-            break
-    if sim_id is None:
-        sim_id = setup_simulation(params)
-
-    ngrun(sim_id)
-
 
 
 def create_zodb(zodb_URI):
+    ''' zodb_URI example:
+            zeo://localhost:1234
+    '''
+
     db = connect_zodb(zodb_URI)
     if not db.has_key('simulations'):
         db['simulations'] = OOBTree.OOBTree()
 AUTHOUR:
     Yi-Xin Liu <liuyxpp@gmail.com>
     Fudan University
-REVISION:
-    2011.8.22
 """
+
+import uuid
 import numpy as np
 from persistent import Persistent
 from .vector2d import Vector2D
 
 class Particle(Persistent):
-    def __init__(self,o=Vector2D(),r0=0.0,r=0.0,t0=0.0,k=1.0,nu=1.0):
-        self.o=o        # center of the circular particle
-        self.r0=float(r0)
-        self.r=float(r)
-        self.t0=float(t0)  # particle born time
-        self.k=float(k)
-        self.nu=float(nu)      # growth kinetics: r=k*(t-t0)^nu
+    def __init__(self, o=Vector2D(), 
+                 r0=0.0, r=0.0, t0=0.0, k=1.0, nu=1.0):
+        self.ID = uuid.uuid4()
+        self.o = o        # center of the circular particle
+        self.r0 = float(r0)
+        self.r = float(r)
+        self.t0 = float(t0)  # particle born time
+        self.te = 0.0 # The time when it becomes inactive
+        self.k = float(k)
+        self.nu = float(nu)      # growth kinetics: r=k*(t-t0)^nu
+
     def __str__(self):
         return 'Particle of radius '+str(self.r)+' locates at '+str(self.o)+' with growth kinetics r= '+str(self.k)+' * t^'+str(self.nu)
-    def __eq__(self,other):
+
+    def __eq__(self, other):
         return ((self.o==other.o) and (self.r==other.r))
+
     def r2(self):
-        return self.r*self.r
+        return self.r * self.r
+
     def area(self):
-        return np.pi*self.r*self.r
+        return np.pi * self.r * self.r
+
     def perimeter(self):
-        return 2*np.pi*self.r
-    def grow_by_function(self,t):
-        if self.nu==0.0:
+        return 2 * np.pi * self.r
+
+    def grow_by_function(self, t):
+        if self.nu == 0.0:
             pass
-        elif self.nu==1.0:
-            self.r=self.r0+self.k*(t-self.t0)
+        elif self.nu == 1.0:
+            self.r = self.r0 + self.k * (t - self.t0)
         else:
-            self.r=self.r0+self.k*(t-self.t0)**self.nu
-    def grow_by_dr(self,dr):
+            self.r = self.r0 + self.k * (t - self.t0)**self.nu
+
+    def grow_by_dr(self, dr):
         self.r += dr
-    def is_inner_point(self,pt):
+
+    def is_inner_point(self, pt):
         '''Boundary points are also inner points'''
         return (self.o.distance2(pt) <= self.r2())
-    def is_boundary_point(self,pt):
-        return (abs(self.o.distance2(pt)-self.r2())<1e-12)
-    def is_touch(self,other):
+
+    def is_boundary_point(self, pt):
+        return (abs(self.o.distance2(pt) - self.r2()) < 1e-12)
+
+    def is_touch(self, other):
         '''is self touches other'''
-        d2=self.o.distance2(other.o)
-        rr=self.r+other.r
+        d2 = self.o.distance2(other.o)
+        rr = self.r + other.r
         rr *= rr
         return (d2 <= rr)
-    def draw_on_lattice(self,color,lx,ly,lattice):
-        Lx,Ly=np.shape(lattice)
-        dx=lx/Lx
-        dy=ly/Ly
-        Ox=int(self.o.x/dx)
-        Oy=int(self.o.y/dy)
-        R=int(self.r/dx) # dx=dy is assumed. isotropic square cell
-        grids=np.indices((2*R+1,2*R+1)).reshape(2,-1).T - R
+
+    def draw_on_lattice(self, color, lx, ly, lattice):
+        Lx,Ly = np.shape(lattice)
+        dx = lx / Lx
+        dy =ly / Ly
+        Ox = int(self.o.x / dx)
+        Oy = int(self.o.y / dy)
+        R = int(self.r / dx) # dx=dy is assumed. isotropic square cell
+        grids = np.indices((2*R+1, 2*R+1)).reshape(2,-1).T - R
         for i in grids:
-            ix,iy=i
-            if (ix*ix+iy*iy<=R*R) and (0<=ix+Ox<Lx) and (0<=iy+Oy<Ly):
+            ix, iy = i
+            if (ix * ix + iy * iy <= R*R) and (0 <= ix + Ox < Lx) and (0 <= iy + Oy < Ly):
                 lattice[ix+Ox,iy+Oy]=color
 
+
 def test():
     p0=Particle() # particle at (0,0) with zero radius
     print p0
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.