diff --git a/CellModeller/Biophysics/BacterialModels/CLBacterium.cl b/CellModeller/Biophysics/BacterialModels/CLBacterium.cl index 36e243b9..f10e096d 100644 --- a/CellModeller/Biophysics/BacterialModels/CLBacterium.cl +++ b/CellModeller/Biophysics/BacterialModels/CLBacterium.cl @@ -122,6 +122,45 @@ void print_matrix(float4* m) { */ } +void cyl_inertia_tensor(float muA, float l, float4 axis, float4 res[]) { + // first find the inv inertia tensor for a capsule along the x axis + float diag = (muA*l*l*l) / 12.f; + + // now find the matrix that transforms from the x-axis to the axis + float4 x_axis = {1.f, 0.f, 0.f, 0.f}; + float4 y_axis = {0.f, 1.f, 0.f, 0.f}; + float4 z_axis = {0.f, 0.f, 1.f, 0.f}; + float rot_ang = acos(dot(x_axis, axis)); + + float4 y_prime = y_axis; + float4 z_prime = z_axis; + + if (rot_ang > EPSILON) { + y_prime = rot(z_prime, rot_ang, y_axis); + z_prime = normalize(cross(x_axis, axis)); + } + + float4 M[4]; + M[0] = axis; + M[1] = y_prime; + M[2] = z_prime; + M[3] = 0.f; + + float4 MT[4]; + transpose(M, MT); + + float4 I[4]; + I[0].s0 = 0.f; + I[1].s1 = diag; + I[2].s2 = diag; + I[3].s3 = 0.f; + + // M^T(I)M + float4 IM[4]; + matmulmat(I, M, IM); + matmulmat(MT, IM, res); +} + void cyl_inv_inertia_tensor(float muA, float l, float4 axis, float4 res[]) { // first find the inv inertia tensor for a capsule along the x axis float diag = 12.f/(muA*l*l*l); @@ -298,6 +337,94 @@ float pt_to_plane_dist(float4 v, float4 n, float4 p) return dot((p-v), n); } +// intest distance from a point p to a sphere at point v +// with radius r +float pt_to_sphere_dist(float4 v, float rad, float4 p) +{ + return length(p - v) - rad; +} + +__kernel void find_sphere_contacts(const int max_cells, + const int max_contacts, + const int n_spheres, + __global const float4* sphere_pts, + __global const float* sphere_coeffs, + __global const float* sphere_rads, + __global const float* sphere_norms, + __global const float4* centers, + __global const float4* dirs, + __global const float* lens, + __global const float* rads, + __global int* n_cts, + __global int* frs, + __global int* tos, + __global float* dists, + __global float4* pts, + __global float4* norms, + __global float* reldists, + __global float* stiff) +{ + int i = get_global_id(0); + + // collision count + int k = n_cts[i]; //keep existing contacts + + float4 end1 = centers[i] - 0.5f*lens[i]*dirs[i]; // 'left' end of the cell + float4 end2 = centers[i] + 0.5f*lens[i]*dirs[i]; // 'right' end of the cell + + for (int n = 0; n < n_spheres; n++) { // loop through all spheres + int to1 = -2*n - 1; // 'to' if left end has contact with sphere n + int to2 = to1 - 1; // 'to' if right end has contact with sphere n + + float dist1 = sphere_norms[n] * pt_to_sphere_dist(sphere_pts[n], sphere_rads[n], end1)-rads[i]; + float dist2 = sphere_norms[n] * pt_to_sphere_dist(sphere_pts[n], sphere_rads[n], end2)-rads[i]; + + // check for old contacts with this sphere + int cti1 = -1; + int cti2 = -1; + for (int m = i*max_contacts; m < i*max_contacts+n_cts[i]; m++) { + if (tos[m] == to1) cti1 = m; + else if (tos[m] == to2) cti2 = m; + } + + bool two_pts = ((cti1 >= 0) || (dist1<0.f) ) && ( (cti2 >= 0) || (dist2<0.f) ); + float stiffness = two_pts*ISQRT2 + (!two_pts)*1.0; + + // if we're in contact, or we were in contact, recompute + if ((cti1 >= 0) || (dist1<0.f) ){ + // need to make a new contact + if (cti1 < 0) { + cti1 = i*max_contacts+k; + k++; + } + + frs[cti1] = i; + tos[cti1] = to1; + dists[cti1] = dist1; + pts[cti1] = end1; // FIXME: not really the right point + norms[cti1] = sphere_norms[n] * normalize(-end1 + sphere_pts[n]); + reldists[cti1] = stiffness*sphere_coeffs[n]*dist1; + stiff[cti1] = stiffness*sphere_coeffs[n]; + } + + if ( (cti2 >= 0) || (dist2<0.f) ){ + if (cti2 < 0) { + cti2 = i*max_contacts+k; + k++; + } + + frs[cti2] = i; + tos[cti2] = to2; + dists[cti2] = dist2; + pts[cti2] = end2; + norms[cti2] = sphere_norms[n] * normalize(-end2 + sphere_pts[n]); + reldists[cti2] = stiffness*sphere_coeffs[n]*dist2; + stiff[cti2] = stiffness*sphere_coeffs[n]; + } + } + n_cts[i] = k; +} + __kernel void find_plane_contacts(const int max_cells, const int max_contacts, @@ -636,8 +763,6 @@ __kernel void collect_tos(const int max_cells, __kernel void build_matrix(const int max_contacts, - const float muA, - const float gamma, __global const float4* centers, __global const float4* dirs, __global const float* lens, @@ -662,21 +787,20 @@ __kernel void build_matrix(const int max_contacts, float4 r_a = pts[i]-centers[a]; float8 fr_ent = 0.f; - float4 Ia[4]; - cyl_inv_inertia_tensor(muA, lens[a]+2.f*rads[a], dirs[a], Ia); float4 nxr_a = 0.f; nxr_a = cross(norms[i], r_a); - fr_ent.s012 = norms[i].s012/(muA*(lens[a]+2.f*rads[a])); - fr_ent.s3 = -dot(nxr_a, Ia[0]); - fr_ent.s4 = -dot(nxr_a, Ia[1]); - fr_ent.s5 = -dot(nxr_a, Ia[2]); - fr_ent.s6 = (1.f/gamma) * dot(dirs[a], r_a) * dot(dirs[a], norms[i])/(lens[a]+2.f*rads[a]); + fr_ent.s012 = norms[i].s012; + fr_ent.s345 = -nxr_a.s012; + //fr_ent.s3 = -dot(nxr_a, Ia[0]); + //fr_ent.s4 = -dot(nxr_a, Ia[1]); + //fr_ent.s5 = -dot(nxr_a, Ia[2]); + fr_ent.s6 = dot(dirs[a], r_a) * dot(dirs[a], norms[i])/(lens[a]+2.f*rads[a]); fr_ents[i] = fr_ent * stiff[i]; int b = tos[i]; - // plane contacts have no to_ent, and have negative indices + // plane and sphere contacts have no to_ent, and have negative indices if (b < 0) { to_ents[i] = 0.f; return; @@ -685,26 +809,25 @@ __kernel void build_matrix(const int max_contacts, float4 r_b = pts[i]-centers[b]; float8 to_ent = 0.f; - float4 Ib[4]; - cyl_inv_inertia_tensor(muA, lens[b]+2.f*rads[b], dirs[b], Ib); float4 nxr_b = 0.f; nxr_b = cross(norms[i], r_b); - to_ent.s012 = norms[i].s012/(muA*(lens[b]+2.f*rads[b])); - to_ent.s3 = -dot(nxr_b, Ib[0]); - to_ent.s4 = -dot(nxr_b, Ib[1]); - to_ent.s5 = -dot(nxr_b, Ib[2]); - to_ent.s6 = (1.f/gamma) * dot(dirs[b], r_b) * dot(dirs[b], norms[i])/(lens[b]+2.f*rads[b]); + to_ent.s012 = norms[i].s012; + to_ent.s345 = -nxr_b.s012; + //to_ent.s3 = -dot(nxr_b, Ib[0]); + //to_ent.s4 = -dot(nxr_b, Ib[1]); + //to_ent.s5 = -dot(nxr_b, Ib[2]); + to_ent.s6 = dot(dirs[b], r_b) * dot(dirs[b], norms[i])/(lens[b]+2.f*rads[b]); to_ents[i] = to_ent * stiff[i]; } -__kernel void calculate_Mx(const int max_contacts, +__kernel void calculate_Bx(const int max_contacts, __global const int* frs, __global const int* tos, __global const float8* fr_ents, __global const float8* to_ents, __global const float8* deltap, - __global float* Mx) + __global float* Bx) { int id = get_global_id(0); int ct = get_global_id(1); @@ -716,34 +839,68 @@ __kernel void calculate_Mx(const int max_contacts, //my machine can't dot float8s... float res0123 = dot(fr_ents[i].s0123, deltap[a].s0123) - dot(to_ents_i.s0123, deltap[b].s0123); float res4567 = dot(fr_ents[i].s4567, deltap[a].s4567) - dot(to_ents_i.s4567, deltap[b].s4567); - Mx[i] = res0123 + res4567; + Bx[i] = res0123 + res4567; } -__kernel void calculate_MTMx(const int max_contacts, +__kernel void calculate_BTBx(const int max_contacts, __global const int* n_cts, __global const int* n_cell_tos, __global const int* cell_tos, __global const float8* fr_ents, __global const float8* to_ents, - __global const float* Mx, - __global float8* MTMx) + __global const float* Bx, + __global float8* BTBx) { int i = get_global_id(0); int base = i*max_contacts; float8 res = 0.f; for (int k = base; k < base+n_cts[i]; k++) { float8 oldres = res; - res += fr_ents[k]*Mx[k]; + res += fr_ents[k]*Bx[k]; } for (int k = base; k < base+n_cell_tos[i]; k++) { int n = cell_tos[k]; if (n < 0) continue; - res -= to_ents[n]*Mx[n]; + res -= to_ents[n]*Bx[n]; } - MTMx[i] = res; + BTBx[i] = res; +} + +__kernel void calculate_Mx(const float muA, + const float gamma, + __global const float4* dirs, + __global const float* lens, + __global const float* rads, + __global const float8* x, + __global float8* Mx) +{ + int i = get_global_id(0); + float l = lens[i] + 2.f * rads[i]; + + float8 xi = x[i]; + float8 v = 0.f; + v.s012 = xi.s012 * muA * l; + + float4 I[4]; + cyl_inertia_tensor(muA, l, dirs[i], I); + float4 L = 0.f; + L.s012 = xi.s345; + float4 w = matmul(I, L); + v.s345 = w.s012; + + //float4 w = 0.f; + //w.s012 = xi.s345; + //float4 a = dirs[i]; + //float4 vv = (w - a*dot(a,w)) * muA * l*l*l / 12.f; + //v.s345 = vv.s012; + + v.s6 = xi.s6 * gamma; + + Mx[i] = v; } + __kernel void calculate_Minv_x(const float muA, const float gamma, __global const float4* dirs, @@ -856,14 +1013,17 @@ __kernel void add_impulse(const float muA, float4 dplin = 0.f; dplin.s012 = deltap_i.s012; - dcenters[i] += dplin/(muA*(lens[i]+2.f*rads[i])); + dcenters[i] += dplin; // FIXME: should probably store these so we don't recompute them - float4 Iinv[4]; - cyl_inv_inertia_tensor(muA, len_i+2.f*rad_i, dir_i, Iinv); - float4 dL = 0.f; - dL.s012 = deltap_i.s345; - float4 dpang = matmul(Iinv, dL); + //float4 Iinv[4]; + //cyl_inv_inertia_tensor(muA, len_i+2.f*rad_i, dir_i, Iinv); + //float4 dL = 0.f; + //dL.s012 = deltap_i.s345; + //float4 dpang = matmul(Iinv, dL); + + float4 dpang = 0.f; + dpang.s012 = deltap_i.s345; float dpangmag = length(dpang); if (dpangmag>ANG_LIMIT) { @@ -872,8 +1032,8 @@ __kernel void add_impulse(const float muA, dangs[i] += dpang; float dplen = deltap_i.s6; - //dlens[i] += max(0.f, target_dlens[i] + dplen/gamma); - dlens[i] = max(0.f, dlens[i]+dplen/gamma); + //dlens[i] += max(0.f, target_dlens[i] + dplen); + dlens[i] = max(0.f, dlens[i]+dplen); } diff --git a/CellModeller/Biophysics/BacterialModels/CLBacterium.py b/CellModeller/Biophysics/BacterialModels/CLBacterium.py index 384e3fc7..85d8180d 100644 --- a/CellModeller/Biophysics/BacterialModels/CLBacterium.py +++ b/CellModeller/Biophysics/BacterialModels/CLBacterium.py @@ -4,6 +4,7 @@ import pyopencl as cl import pyopencl.array as cl_array from pyopencl.array import vec +from pyopencl.array import max as device_max from pyopencl.elementwise import ElementwiseKernel from pyopencl.reduction import ReductionKernel import random @@ -21,14 +22,14 @@ def __init__(self, simulator, max_substeps=8, max_cells=10000, max_contacts=24, - max_planes=4, + max_planes=1, + max_spheres=1, max_sqs=192**2, grid_spacing=5.0, muA=1.0, gamma=10.0, dt=None, cgs_tol=5e-3, - reg_param=0.1, jitter_z=True, alternate_divisions=False, printing=True, @@ -49,19 +50,20 @@ def __init__(self, simulator, self.max_cells = max_cells self.max_contacts = max_contacts self.max_planes = max_planes + self.max_spheres = max_spheres self.max_sqs = max_sqs self.grid_spacing = grid_spacing self.muA = muA self.gamma = gamma self.dt = dt self.cgs_tol = cgs_tol - self.reg_param = numpy.float32(reg_param) self.max_substeps = max_substeps self.n_cells = 0 self.n_cts = 0 self.n_planes = 0 + self.n_spheres = 0 self.next_id = 0 @@ -88,6 +90,7 @@ def reset(self): self.n_cells=0 self.n_cts=0 self.n_planes=0 + self.n_spheres=0 def setRegulator(self, regulator): self.regulator = regulator @@ -111,13 +114,13 @@ def addCell(self, cellState, pos=(0,0,0), dir=(1,0,0), rad=0.5, **kwargs): # Eventually prob better to have a generic editCell() that deals with this stuff # def moveCell(self, cellState, delta_pos): - print "cell idx = %d"%cellState.idx + print("cell idx = %d"%cellState.idx) i = cellState.idx cid = cellState.id - print "cell center = " - print self.cell_centers[i] - print "delta_pos" - print delta_pos + print("cell center = ") + print(self.cell_centers[i]) + print("delta_pos") + print(delta_pos) pos = numpy.array(tuple(self.cell_centers[i])) pos[0:3] += numpy.array(tuple(delta_pos)) self.cell_centers[i] = pos @@ -133,6 +136,15 @@ def addPlane(self, pt, norm, coeff): self.plane_coeffs[pidx] = coeff self.set_planes() + def addSphere(self, pt, rad, coeff, norm): + sidx = self.n_spheres + self.n_spheres += 1 + self.sphere_pts[sidx] = tuple(pt)+(0,) + self.sphere_rads[sidx] = rad + self.sphere_coeffs[sidx] = coeff + self.sphere_norms[sidx] = norm + self.set_spheres() + def hasNeighbours(self): return False @@ -149,7 +161,7 @@ def init_cl(self): def init_kernels(self): """Set up the OpenCL kernels.""" from pkg_resources import resource_string - kernel_src = resource_string(__name__, 'CLBacterium.cl') + kernel_src = resource_string(__name__, 'CLBacterium.cl').decode() self.program = cl.Program(self.context, kernel_src).build(cache_dir=False) # Some kernels that seem like they should be built into pyopencl... @@ -165,6 +177,13 @@ def init_kernels(self): self.vsubkx = ElementwiseKernel(self.context, "float8 *res, const float k, const float8 *in1, const float8 *in2", "res[i] = in1[i] - k*in2[i]", "vecsubkx") + self.vmulk = ElementwiseKernel(self.context, + "float8 *res, const float k, const float8 *in1", + "res[i] = k*in1[i]", "vecmulk") + self.vnorm = ElementwiseKernel(self.context, + "float8 *res, const float8 *in1", + "res[i] = dot(in1[i], in1[i]", "vecnorm") + # cell geometry kernels self.calc_cell_area = ElementwiseKernel(self.context, "float* res, float* r, float* l", @@ -233,6 +252,17 @@ def init_data(self): self.plane_coeffs = numpy.zeros(plane_geom, numpy.float32) self.plane_coeffs_dev = cl_array.zeros(self.queue, plane_geom, numpy.float32) + # constraint spheres + sphere_geom = (self.max_spheres,) + self.sphere_pts = numpy.zeros(sphere_geom, vec.float4) + self.sphere_pts_dev = cl_array.zeros(self.queue, sphere_geom, vec.float4) + self.sphere_rads = numpy.zeros(sphere_geom, numpy.float32) + self.sphere_rads_dev = cl_array.zeros(self.queue, sphere_geom, numpy.float32) + self.sphere_coeffs = numpy.zeros(sphere_geom, numpy.float32) + self.sphere_coeffs_dev = cl_array.zeros(self.queue, sphere_geom, numpy.float32) + self.sphere_norms = numpy.zeros(sphere_geom, numpy.float32) + self.sphere_norms_dev = cl_array.zeros(self.queue, sphere_geom, numpy.float32) + # contact data ct_geom = (self.max_cells, self.max_contacts) self.ct_frs = numpy.zeros(ct_geom, numpy.int32) @@ -275,8 +305,8 @@ def init_data(self): self.deltap_dev = cl_array.zeros(self.queue, cell_geom, vec.float8) self.Mx = numpy.zeros(mat_geom, numpy.float32) self.Mx_dev = cl_array.zeros(self.queue, mat_geom, numpy.float32) - self.MTMx = numpy.zeros(cell_geom, vec.float8) - self.MTMx_dev = cl_array.zeros(self.queue, cell_geom, vec.float8) + self.BTBx = numpy.zeros(cell_geom, vec.float8) + self.BTBx_dev = cl_array.zeros(self.queue, cell_geom, vec.float8) self.Minvx_dev = cl_array.zeros(self.queue, cell_geom, vec.float8) # CGS intermediates @@ -287,7 +317,7 @@ def init_data(self): def load_from_cellstates(self, cell_states): - for (cid,cs) in cell_states.items(): + for (cid,cs) in list(cell_states.items()): i = cs.idx self.cell_centers[i] = tuple(cs.pos)+(0,) self.cell_dirs[i] = tuple(cs.dir)+(0,) @@ -301,10 +331,10 @@ def load_from_cellstates(self, cell_states): def load_test_data(self): import CellModeller.Biophysics.BacterialModels.CLData as data - self.cell_centers.put(range(len(data.pos)), data.pos) - self.cell_dirs.put(range(len(data.dirs)), data.dirs) - self.cell_lens.put(range(len(data.lens)), data.lens) - self.cell_rads.put(range(len(data.rads)), data.rads) + self.cell_centers.put(list(range(len(data.pos))), data.pos) + self.cell_dirs.put(list(range(len(data.dirs))), data.dirs) + self.cell_lens.put(list(range(len(data.lens))), data.lens) + self.cell_rads.put(list(range(len(data.rads))), data.rads) self.n_cells = data.n_cells self.set_cells() @@ -427,6 +457,14 @@ def set_planes(self): self.plane_coeffs_dev[0:self.n_planes].set(self.plane_coeffs[0:self.n_planes]) + def set_spheres(self): + """Copy sphere pts and coeffs to the device from local.""" + self.sphere_pts_dev[0:self.n_spheres].set(self.sphere_pts[0:self.n_spheres]) + self.sphere_rads_dev[0:self.n_spheres].set(self.sphere_rads[0:self.n_spheres]) + self.sphere_coeffs_dev[0:self.n_spheres].set(self.sphere_coeffs[0:self.n_spheres]) + self.sphere_norms_dev[0:self.n_spheres].set(self.sphere_norms[0:self.n_spheres]) + + def get_cts(self): """Copy contact froms, tos, dists, pts, and norms from the device.""" self.ct_frs[0:self.n_cts] = self.ct_frs_dev[0:self.n_cts].get() @@ -455,13 +493,13 @@ def matrixTest(self): if ii!=self.n_cells-1 or jj!=6: opstring = opstring + '\t' opstring = opstring + '\n' - print "MTM" - print opstring + print("MTM") + print(opstring) open('CellModeller/Biophysics/BacterialModels/matrix.mat', 'w').write(opstring) def dump_cell_data(self, n): - import cPickle + import pickle filename = 'data/data-%04i.pickle'%n outfile = open(filename, 'wb') data = (self.n_cells, @@ -470,7 +508,7 @@ def dump_cell_data(self, n): self.cell_lens_dev.get(), self.cell_rads_dev.get(), self.parents), - cPickle.dump(data, outfile, protocol=-1) + pickle.dump(data, outfile, protocol=-1) def dydt(self): self.set_cells() @@ -479,7 +517,7 @@ def finish(self): # pull cells from the device and update simulator if self.simulator: self.get_cells() - for state in self.simulator.cellStates.values(): + for state in list(self.simulator.cellStates.values()): self.updateCellState(state) def progress_init(self, dt): @@ -489,7 +527,7 @@ def progress_init(self, dt): self.n_ticks = int(math.ceil(dt/self.dt)) else: self.n_ticks = 1 - #print "n_ticks = %d"%(self.n_ticks) + # print("n_ticks = %d"%(self.n_ticks)) self.actual_dt = dt / float(self.n_ticks) self.progress_initialised = True @@ -508,7 +546,7 @@ def progress_finalise(self): self.minutes_elapsed = (numpy.float32(self.seconds_elapsed) / 60.0) self.hours_elapsed = (numpy.float32(self.minutes_elapsed) / 60.0) if self.frame_no % 10 == 0: - print '% 8i % 8i cells % 8i contacts %f hour(s) or %f minute(s) or %f second(s)' % (self.frame_no, self.n_cells, self.n_cts, self.hours_elapsed, self.minutes_elapsed, self.seconds_elapsed) + print('% 8i % 8i cells % 8i contacts %f hour(s) or %f minute(s) or %f second(s)' % (self.frame_no, self.n_cells, self.n_cts, self.hours_elapsed, self.minutes_elapsed, self.seconds_elapsed)) # pull cells from the device and update simulator if self.simulator: self.get_cells() @@ -517,7 +555,7 @@ def progress_finalise(self): # TJR: add flag for this cos a bit time consuming if self.computeNeighbours: self.updateCellNeighbours(self.simulator.idxToId) - for state in self.simulator.cellStates.values(): + for state in list(self.simulator.cellStates.values()): self.updateCellState(state) def step(self, dt): @@ -576,11 +614,12 @@ def sub_tick(self, dt): self.collect_tos() self.sub_tick_i += 1 + alpha = 10**(self.sub_tick_i) new_cts = self.n_cts - old_n_cts if (new_cts>0 or self.sub_tick_i==0) and self.sub_tick_i:\t' + if type(val) in [float, np.float32, np.float64]: + txt += '%g'%val + elif type(val) in [list, tuple, np.array]: + txt += ', '.join(['%g'%v for v in val]) + else: + txt += str(val) + txt += '
' self.selectedCell.emit(txt) self.updateGL() @@ -228,7 +267,7 @@ def translate(self, _trans): # currently the translation is not great due to the way PyGLWidget works if False: #self.sim and self.sim.cellStates.has_key(cid): self.sim.moveCell(cid, _trans) - print "Called self.sim.moveCell" + print("Called self.sim.moveCell") self.updateSelectedCell() else: # Translate the object by _trans @@ -254,16 +293,17 @@ def paintGL(self): #glScalef(s,s,s) # Draw a grid in xy plane + glEnable(GL_DEPTH_TEST) glDisable(GL_LIGHTING) glColor3f(1.0, 1.0, 1.0) glEnable(GL_LINE_SMOOTH) glLineWidth(1.0) glBegin(GL_LINES) - for i in range(5): - glVertex(-20, (i-2)*10) - glVertex(20, (i-2)*10) - glVertex((i-2)*10, -20) - glVertex((i-2)*10, 20) + for i in range(25): + glVertex(-120, (i-12)*10, 0) + glVertex(120, (i-12)*10, 0) + glVertex((i-12)*10, -120, 0) + glVertex((i-12)*10, 120, 0) glEnd() # Draw x,y,z axes @@ -271,13 +311,13 @@ def paintGL(self): glBegin(GL_LINES) glColor3f(1.0,0.0,0.0) glVertex(0,0,0) - glVertex(5,0,0) + glVertex(25,0,0) glColor3f(0.0,1.0,0.0) glVertex(0,0,0) - glVertex(0,5,0) + glVertex(0,25,0) glColor3f(0.0,0.0,1.0) glVertex(0,0,0) - glVertex(0,0,5) + glVertex(0,0,25) glEnd() glEnable(GL_DEPTH_TEST) glEnable(GL_LIGHTING) diff --git a/CellModeller/GUI/PyGLGUI.ui b/CellModeller/GUI/PyGLGUI.ui index 6a6dd7eb..3fb41184 100644 --- a/CellModeller/GUI/PyGLGUI.ui +++ b/CellModeller/GUI/PyGLGUI.ui @@ -75,6 +75,7 @@ + @@ -116,6 +117,17 @@ Ctrl+O + + + Load Geometry + + + Load geometry from pickle file + + + Ctrl+l + + Load Pickle @@ -224,6 +236,22 @@ + + actionLoadGeometry + triggered() + PyGLCMViewer + loadGeometry() + + + -1 + -1 + + + 589 + 357 + + + actionLoadPickle triggered() diff --git a/CellModeller/GUI/PyGLWidget.py b/CellModeller/GUI/PyGLWidget.py index 726e939b..e10c9d11 100644 --- a/CellModeller/GUI/PyGLWidget.py +++ b/CellModeller/GUI/PyGLWidget.py @@ -32,7 +32,7 @@ # #=============================================================================== -from PyQt4 import QtCore, QtGui, QtOpenGL +from PyQt5 import QtCore, QtGui, QtOpenGL import math import numpy import numpy.linalg as linalg @@ -67,14 +67,15 @@ def __init__(self, parent = None): self.last_point_ok_ = False self.last_point_3D_ = [1.0, 0.0, 0.0] self.isInRotation_ = False - self.pickSize = 3 + self.pickSize = 18 + self.pix_ratio = 1. # connections #self.signalGLMatrixChanged.connect(self.printModelViewMatrix) @QtCore.pyqtSlot() def printModelViewMatrix(self): - print self.modelview_matrix_ + print(self.modelview_matrix_) def initializeGL(self): # OpenGL state @@ -110,8 +111,9 @@ def set_pick_projection(self, x, y, _near, _far, _fovy): self.makeCurrent() glMatrixMode( GL_PROJECTION ) glLoadIdentity() - viewport =glGetIntegerv(GL_VIEWPORT) - gluPickMatrix(x, viewport[3]-y, self.pickSize, self.pickSize, viewport); + viewport = glGetIntegerv(GL_VIEWPORT) + #gluPickMatrix(x, viewport[3]-y, self.pickSize, self.pickSize, viewport); + gluPickMatrix(x*self.pix_ratio, viewport[3]-y*self.pix_ratio, self.pickSize, self.pickSize, viewport); gluPerspective( self.fovy_, float(self.width()) / float(self.height()), self.near_, self.far_ ) def set_center(self, _cog): @@ -219,7 +221,7 @@ def map_to_sphere(self, _v2D): def wheelEvent(self, _event): # Use the mouse wheel to zoom in/out - d = - float(_event.delta()) / 200.0 * self.radius_ + d = - float(_event.angleDelta().y()) / 200.0 * self.radius_ self.translate([0.0, 0.0, d]) self.updateGL() _event.accept() @@ -238,13 +240,16 @@ def selectName(self, point): buf = glRenderMode(GL_RENDER) selectedName = -1 closest_z = 1.0 + print('buf ', buf) for hit_record in buf: min_depth, max_depth, names = hit_record + print('Names ', names) if min_depth < closest_z: closest_z = min_depth for name in names: if name: selectedName = name + print('Selected name ', selectedName) glMatrixMode( GL_PROJECTION ) glPopMatrix() return selectedName @@ -282,7 +287,7 @@ def mouseMoveEvent(self, _event): # move in z direction if (((_event.buttons() & QtCore.Qt.LeftButton) and (_event.buttons() & QtCore.Qt.MidButton)) or (_event.buttons() & QtCore.Qt.LeftButton and _event.modifiers() & QtCore.Qt.ControlModifier)): - print "translating in Z" + print("translating in Z") value_y = self.radius_ * dy * 2.0 / h self.translate([0.0, 0.0, value_y]) # move in x,y direction diff --git a/CellModeller/GUI/Renderers.py b/CellModeller/GUI/Renderers.py index 74878cf1..f5dc3128 100644 --- a/CellModeller/GUI/Renderers.py +++ b/CellModeller/GUI/Renderers.py @@ -3,16 +3,123 @@ from OpenGL.arrays import vbo import math import numpy -import numpy.linalg as la import random +#import pygame # helper module for text rendering - turned off to see if it suppresses segfaults +#from pyopencl.array import vec + +class GLSphereRenderer: + def __init__(self, sim, properties=None, scales = None): + self.ncells_list = 0 + self.ncells_names_list = 0 + self.dlist = None + self.dlist_names = None + self.cellcol = [1, 1, 1] + self.sim = sim + self.quad = gluNewQuadric() + self.properties = properties + self.scales = scales + + def init_gl(self): + pass + + def build_list(self, cells): + if self.dlist: + glDeleteLists(self.dlist, 1) + index = glGenLists(1) + glNewList(index, GL_COMPILE) + self.render_cells() + glEndList() + self.dlist = index + + def build_list_names(self, cells): + if self.dlist_names: + glDeleteLists(self.dlist_names, 1) + index = glGenLists(1) + glNewList(index, GL_COMPILE) + self.render_cell_names() + glEndList() + self.dlist_names = index + + def render_gl(self, selection=None): + cells = self.sim.cellStates.values() + states = self.sim.cellStates.items() + # FIXED ============================================================================= + # Before, the renderer would only draw cells when the number of cells changed. + # Now it draws them whenever render_gl is called (by paintGL in PyGLCMViewer.py) + #if len(cells)!=self.ncells_list or len(cells)<500: + # self.build_list(cells) + # self.ncells_list = len(cells) + + #if len(cells)!=self.ncells_list: + # self.build_list(cells) + # self.ncells_list = len(cells) + #==================================================================================== + #glCallList(self.dlist) + self.render_cells(selection=selection) + + + def renderNames_gl(self, selection=None): + #cells = self.sim.cellStates.values() + #if len(cells)!=self.ncells_names_list or len(cells)<500: + # self.build_list_names(cells) + # self.ncells_names_list = len(cells) + #glCallList(self.dlist_names) + #for cell in cells: self.render_cell_name(cell, selection) + self.render_cell_names() + + def render_cell_names(self): + # glDisable(GL_DEPTH_TEST) + glDisable(GL_LIGHTING) + for cell in self.sim.cellStates.values(): + l = cell.length + p = cell.pos + + cid = cell.id + glPushName(cid) + + glMatrixMode(GL_MODELVIEW) + glPushMatrix() + glTranslatef(p[0],p[1],p[2]) + gluSphere(self.quad, r, 8, 8) + glPopMatrix() + + glPopName() + + glEnable(GL_LIGHTING) + + def render_cell(self, cell, selection): + r = cell.radius + p = cell.pos + cid = cell.id + if selection==cid: + cellcol = [1,0,0] + else: + cellcol = cell.color + + glColor3fv(cellcol) + + glMatrixMode(GL_MODELVIEW) + glPushMatrix() + glTranslatef(p[0],p[1],p[2]) + gluSphere(self.quad, r, 8, 8) + glPopMatrix() + + def render_cells(self, selection=None): + + # PLACEHOLDER + + #glDisable(GL_DEPTH_TEST) + glDisable(GL_LIGHTING) + cells = self.sim.cellStates.values() + for cell in cells: + self.render_cell(cell, selection) class GLGridRenderer: - def __init__(self, sig, integ, rng=None, alpha=0.5): + def __init__(self, sig, integ, rng=None, alpha=1.): self.sig = sig self.integ = integ self.rng = rng self.alpha = alpha - self.size = self.sig.gridSize self.dim = self.sig.gridDim[1:] self.len = self.dim[0]*self.dim[1] @@ -45,19 +152,18 @@ def render_gl(self, selection=None): mx = numpy.max(self.imageData) mn = numpy.min(self.imageData) if mx>mn: - scale = 255.0/(mx-mn) + scale = 255/(mx-mn) else: - scale = 1.0 - #print "mx = %f, mn = %f, scale = %f"%(mx,mn,scale) - self.imageData = numpy.clip((self.imageData - mn)*scale, 0.0, 255.0) + scale = 1 + self.imageData = (self.imageData - mn)*scale #print "Signal grid range = %f to %f"%(mn,mx) for s in range(self.sig.nSignals): self.byteImageData[0:self.dim[0],0:self.dim[1],s] = self.imageData[s,:,:].astype(numpy.uint8) - glDisable(GL_CULL_FACE) + glEnable(GL_TEXTURE_2D) glDisable(GL_LIGHTING) - glDisable(GL_DEPTH_TEST) + #glDisable(GL_DEPTH_TEST) glBindTexture(GL_TEXTURE_2D, self.texture) #Load the Data into the Texture glTexImage2D( GL_TEXTURE_2D, 0, GL_RGB, self.texDim, self.texDim, 0, GL_RGB, GL_UNSIGNED_BYTE, self.byteImageData ) @@ -80,7 +186,7 @@ def render_gl(self, selection=None): glVertex3f(self.orig[0] , self.orig[1]+self.dim[1]*self.size[1], 0.0) glEnd() glEnable(GL_LIGHTING) - glEnable(GL_DEPTH_TEST) + #glEnable(GL_DEPTH_TEST) glDisable(GL_BLEND) glDisable(GL_TEXTURE_2D) @@ -91,36 +197,36 @@ def __init__(self, sim, chanIdx): self.wallcol = [0, 0, 0] self.nodecol = [0, 0, 0] self.cellcol = [1, 1, 1] - self.chanIdx = chanIdx + self.chanIdx = chanIdx def init_gl(self): pass def renderNames_gl(self): - pass + pass def render_gl(self, selection=None): - # Render signals[chanIdx] as rgb - # - maxSig = [0,0,0] + # Render signals[chanIdx] as rgb + # + maxSig = [0,0,0] for i in range(len(self.chanIdx)): for (cid,cellState) in self.sim.cellStates.items(): - s = cellState.signals[self.chanIdx[i]] - if s > maxSig[i]: - maxSig[i] = s - if maxSig[i]==0: - maxSig[i]=1.0 - - for (cid,cellState) in self.sim.cellStates.items(): - col = [0,0,0] - for i in range(len(self.chanIdx)): - col[i] = cellState.signals[self.chanIdx[i]]/maxSig[i] - glColor3fv(col) - glBegin(GL_POLYGON) - for p in cellState.nodep: - glVertex2d(p[0], p[1]) - glEnd() - + s = cellState.signals[self.chanIdx[i]] + if s > maxSig[i]: + maxSig[i] = s + if maxSig[i]==0: + maxSig[i]=1.0 + + for (cid,cellState) in self.sim.cellStates.items(): + col = [0,0,0] + for i in range(len(self.chanIdx)): + col[i] = cellState.signals[self.chanIdx[i]]/maxSig[i] + glColor3fv(col) + glBegin(GL_POLYGON) + for p in cellState.nodep: + glVertex2d(p[0], p[1]) + glEnd() + class GLPlantRenderer: def __init__(self, sim): self.sim = sim @@ -132,8 +238,8 @@ def init_gl(self): pass def renderNames_gl(self): - # Render cell_ids for selection - # + # Render cell_ids for selection + # for (cid,cellState) in self.sim.cellStates.items(): glColor3fv(self.cellcol) glPushName(cid) @@ -141,11 +247,11 @@ def renderNames_gl(self): for p in cellState.nodep: glVertex2d(p[0], p[1]) glEnd() - glPopName() + glPopName() def render_gl(self, selection=None): # Render model using PyOpenGL - # + # # principle axes #scale=10.0 #c=cell.get_centre() @@ -158,40 +264,39 @@ def render_gl(self, selection=None): #glVertex2d(c[0],c[1]) #glVertex2d(c[0]+cell.pa2[0]*scale,c[1]+cell.pa2[1]*scale) #glEnd() - glDisable(GL_LIGHTING) - glDisable(GL_DEPTH_TEST) - cellStates = self.sim.cellStates + glDisable(GL_LIGHTING) + glDisable(GL_DEPTH_TEST) + cellStates = self.sim.cellStates - #for (cid,cellState) in cellStates.items(): - # cellcol = cellState.color - # glColor3fv(cellcol) + #for (cid,cellState) in cellStates.items(): + # cellcol = cellState.color + # glColor3fv(cellcol) # glBegin(GL_POLYGON) # for p in cellState.nodep: # glVertex2d(p[0], p[1]) # glEnd() - - for (cid,cellState) in cellStates.items(): - col = self.wallcol - for (wp1,wp2) in cellState.wallp: - glColor3fv(col) - glLineWidth(2) - glBegin(GL_LINES) - glVertex2d(wp1[0], wp1[1]) - glVertex2d(wp2[0], wp2[1]) - glEnd() - - if selection>0 and cellStates.has_key(selection): - scell = cellStates[selection] - glColor3fv([1,0,0]) - for (wp1,wp2) in scell.wallp: - glLineWidth(2) - glBegin(GL_LINES) - glVertex2d(wp1[0], wp1[1]) - glVertex2d(wp2[0], wp2[1]) - glEnd() - - # Draw nodes as points - #glColor3fv(self.nodecol) + for (cid,cellState) in cellStates.items(): + col = self.wallcol + for (wp1,wp2) in cellState.wallp: + glColor3fv(col) + glLineWidth(2) + glBegin(GL_LINES) + glVertex2d(wp1[0], wp1[1]) + glVertex2d(wp2[0], wp2[1]) + glEnd() + + if selection>0 and cellStates.has_key(selection): + scell = cellStates[selection] + glColor3fv([1,0,0]) + for (wp1,wp2) in scell.wallp: + glLineWidth(2) + glBegin(GL_LINES) + glVertex2d(wp1[0], wp1[1]) + glVertex2d(wp2[0], wp2[1]) + glEnd() + + # Draw nodes as points + #glColor3fv(self.nodecol) #glPointSize(3) #for node in self.model.nodes: @@ -211,280 +316,953 @@ def render_gl(self, selection=None): class GLBacteriumRenderer: - def __init__(self, sim, properties=None, scales = None): - self.ncells_list = 0 - self.ncells_names_list = 0 - self.names_list_step = -1 - self.list_step = -1 - self.dlist = None - self.dlist_names = None - self.cellcol = [1, 1, 1] - self.sim = sim - self.quad = gluNewQuadric() - self.properties = properties - self.scales = scales - - def init_gl(self): - pass - - def build_list(self, cells): - if self.dlist: - glDeleteLists(self.dlist, 1) - index = glGenLists(1) - glNewList(index, GL_COMPILE) - self.render_cells() - glEndList() - self.dlist = index - - def build_list_names(self, cells): - if self.dlist_names: - glDeleteLists(self.dlist_names, 1) - index = glGenLists(1) - glNewList(index, GL_COMPILE) - self.render_cell_names() - glEndList() - self.dlist_names = index - - def render_gl(self, selection=None): - # Uncomment this to directly draw all cells each frame - #self.render_cells() - - cells = self.sim.cellStates.values() - if self.sim.stepNum > self.list_step: - self.build_list(cells) - #self.ncells_list = len(cells) - self.list_step = self.sim.stepNum - - glCallList(self.dlist) + def __init__(self, sim, properties=None, scales = None): + self.ncells_list = 0 + self.ncells_names_list = 0 + self.dlist = None + self.dlist_names = None + self.cellcol = [1, 1, 1] + self.sim = sim + self.quad = gluNewQuadric() + self.properties = properties + self.scales = scales + self.last_rendered_step = -1 + + def init_gl(self): + pass - if self.sim.cellStates.has_key(selection): - self.render_selected_cell(selection) - self.render_constraints() - - def render_constraints(self): - phys = self.sim.phys - glDisable(GL_LIGHTING) - glEnable(GL_DEPTH_TEST) - glEnable(GL_LINE_SMOOTH) - glLineWidth(1.0) - #glBegin(GL_LINES) - glColor3f(1,0,0) - for i in range(phys.n_planes): - p = phys.plane_pts[i] - nn = phys.plane_norms[i] - n = [nn[i] for i in range(3)] - - # Make an orthonormal basis in the plane - v1 = numpy.cross([0,0,1],n) - if la.norm(v1)<1e-6: - v1 = numpy.cross([0,1,0],n) - v2 = numpy.cross(v1,n) - - glBegin(GL_LINE_STRIP) - glVertex3f(p[0], p[1], p[2]) - #glVertex3f(p[0]+n[0]*5, p[1]+n[1]*5, p[2]+n[2]*5) - glVertex3f(p[0]+v1[0]*5, p[1]+v1[1]*5, p[2]+v1[2]*5) - glVertex3f(p[0]+(v1[0]+v2[0])*5, p[1]+(v1[1]+v2[1])*5, p[2]+(v1[2]+v2[2])*5) - glVertex3f(p[0]+v2[0]*5, p[1]+v2[1]*5, p[2]+v2[2]*5) - glVertex3f(p[0], p[1], p[2]) - glEnd() - #glEnd() - glEnable(GL_LIGHTING) - glDisable(GL_DEPTH_TEST) - - - def render_selected_cell(self, selection): - glEnable(GL_DEPTH_TEST) - cellcol = [1,0,0] - cell = self.sim.cellStates[selection] - l = cell.length - r = cell.radius - - (e1,e2) = cell.ends - ae1 = numpy.array(e1) - ae2 = numpy.array(e2) - zaxis = numpy.array([0,0,1]) - caxis = numpy.array(cell.dir) #(ae2-ae1)/l - rotaxis = numpy.cross(caxis, zaxis) - rotangle = numpy.arccos(numpy.dot(caxis,zaxis)) - - # draw the outlines antialiased in RED - glColor3f(1.0, 0.0, 0.0) - glEnable(GL_BLEND) - glBlendFunc(GL_SRC_ALPHA ,GL_ONE_MINUS_SRC_ALPHA) - glEnable(GL_LINE_SMOOTH) - glLineWidth(8.0) - # draw wireframe for back facing polygons and cull front-facing ones - glPolygonMode(GL_BACK, GL_FILL) - glEnable(GL_CULL_FACE) - glCullFace(GL_FRONT) - glDepthFunc(GL_LEQUAL) - - glMatrixMode(GL_MODELVIEW) - glPushMatrix() - glTranslatef(e1[0],e1[1],e1[2]) - gluSphere(self.quad, r, 8, 8) - glRotatef(-rotangle*180.0/numpy.pi, rotaxis[0], rotaxis[1], rotaxis[2]) - gluCylinder(self.quad, r, r , l, 8, 1) - glPopMatrix() - glPushMatrix() - glTranslatef(e2[0],e2[1],e2[2]) - gluSphere(self.quad, r, 8, 8) - glPopMatrix() - - glDisable(GL_DEPTH_TEST) - - - - def renderNames_gl(self, selection=None): - #self.render_cell_names() - cells = self.sim.cellStates.values() - if self.sim.stepNum > self.names_list_step: - self.build_list_names(cells) - self.names_list_step = self.sim.stepNum - glCallList(self.dlist_names) - - def render_cell_names(self): - glEnable(GL_DEPTH_TEST) - glDisable(GL_LIGHTING) - glDepthFunc(GL_LESS) - glDisable(GL_CULL_FACE) - glPolygonMode(GL_FRONT, GL_FILL) - glDisable(GL_LINE_SMOOTH) - glDisable(GL_BLEND) - - for cell in self.sim.cellStates.values(): - l = cell.length - r = cell.radius - - (e1,e2) = cell.ends - ae1 = numpy.array(e1) - ae2 = numpy.array(e2) - zaxis = numpy.array([0,0,1]) - caxis = numpy.array(cell.dir) #(ae2-ae1)/l - rotaxis = numpy.cross(caxis, zaxis) - rotangle = numpy.arccos(numpy.dot(caxis,zaxis)) - - cid = cell.id - glPushName(cid) - - glMatrixMode(GL_MODELVIEW) - glPushMatrix() - glTranslatef(e1[0],e1[1],e1[2]) - glRotatef(-rotangle*180.0/numpy.pi, rotaxis[0], rotaxis[1], rotaxis[2]) - gluCylinder(self.quad, r, r , l, 8, 1) - gluSphere(self.quad, r, 8, 8) - glPopMatrix() - glPushMatrix() - glTranslatef(e2[0],e2[1],e2[2]) - gluSphere(self.quad, r, 8, 8) - glPopMatrix() - - glPopName() - - glEnable(GL_LIGHTING) - glDisable(GL_DEPTH_TEST) - - - def render_cells(self, selection=None): - glEnable(GL_DEPTH_TEST) - glDisable(GL_LIGHTING) - - cells = self.sim.cellStates.values() - for cell in cells: - l = cell.length - #r = cell.radius*2.0 - r = cell.radius - - (e1,e2) = cell.ends - ae1 = numpy.array(e1) - ae2 = numpy.array(e2) - zaxis = numpy.array([0,0,1]) - caxis = numpy.array(cell.dir) #(ae2-ae1)/l - rotaxis = numpy.cross(caxis, zaxis) - rotangle = numpy.arccos(numpy.dot(caxis,zaxis)) + #def render_text(self, position, textString, fontSize): + # pygame.font.init() + # font = pygame.font.Font (None, fontSize) + # textSurface = font.render(textString, True, (0,0,0,225), (255,255,255,255)) + # textData = pygame.image.tostring(textSurface, "RGBA", True) + # glRasterPos3d(*position) + # glDrawPixels(textSurface.get_width(), textSurface.get_height(), GL_RGBA, GL_UNSIGNED_BYTE, textData) + + def build_list(self, cells): + if self.dlist: + glDeleteLists(self.dlist, 1) + index = glGenLists(1) + glNewList(index, GL_COMPILE) + self.render_cells() + glEndList() + self.dlist = index + + def build_list_names(self, cells): + if self.dlist_names: + glDeleteLists(self.dlist_names, 1) + index = glGenLists(1) + glNewList(index, GL_COMPILE) + self.render_cell_names() + glEndList() + self.dlist_names = index + + def render_gl(self, selection=None): + cells = self.sim.cellStates.values() + states = self.sim.cellStates.items() + # FIXED ============================================================================= + # Before, the renderer would only draw cells when the number of cells changed. + # Now it draws them whenever render_gl is called (by paintGL in PyGLCMViewer.py) + #if len(cells)!=self.ncells_list or len(cells)<500: + # self.build_list(cells) + # self.ncells_list = len(cells) + + #if len(cells)!=self.ncells_list: + if self.last_rendered_step0: - glColor4f(s.signals[0],0,0,1) + if len(s.signals)>0: + glColor4f(s.signals[0],0,0,1) else: glColor4fv(c.color) #c.color[0:3]) else: @@ -786,7 +1564,7 @@ def render_cell_name(self, cell): def render_gl(self, selection=None): - cells = self.sim.cellStates.values() + cells = self.sim.cellStates.values() for cell in cells: self.render_cell(cell, selection) def renderNames_gl(self, selection=None): @@ -895,9 +1673,9 @@ def render_gl(self, selection=None): glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA) cells = self.sim.cellStates.values() - if len(cells)!=self.ncells_list: - self.build_list(cells) - self.ncells_list = len(cells) + #if len(cells)!=self.ncells_list: + self.build_list(cells) + #self.ncells_list = len(cells) glCallList(self.dlist) #for cell in cells: self.render_cell(cell, selection) diff --git a/CellModeller/Integration/CLCrankNicIntegrator.py b/CellModeller/Integration/CLCrankNicIntegrator.py index 3b673aa8..39ee8655 100644 --- a/CellModeller/Integration/CLCrankNicIntegrator.py +++ b/CellModeller/Integration/CLCrankNicIntegrator.py @@ -7,6 +7,7 @@ import pyopencl.array as cl_array from pyopencl.array import vec import math +from functools import reduce def unique_stable(ar, return_index=False, return_inverse=False): @@ -102,7 +103,7 @@ def __init__(self, sim, nSignals, nSpecies, maxCells, sig, greensThreshold=1e-12 # set the species for existing states to views of the levels array cs = self.cellStates - for id,c in cs.items(): + for id,c in list(cs.items()): c.species = self.specLevel[c.idx,:] @@ -135,7 +136,7 @@ def computeGreensFunc(self): self.greensFunc = self.greensFunc[:, min(inds[:,1]):max(inds[:,1])+1, \ min(inds[:,2]):max(inds[:,2])+1, \ min(inds[:,3]):max(inds[:,3])+1] - print "Truncated Green's function size is " + str(self.greensFunc.shape) + print("Truncated Green's function size is " + str(self.greensFunc.shape)) def addCell(self, cellState): @@ -203,7 +204,7 @@ def initKernels(self): sigRateKernel = self.regul.sigRateCL() #kernel_src = open('CellModeller/Integration/CLCrankNicIntegrator.cl', 'r').read() from pkg_resources import resource_string - kernel_src = resource_string(__name__, 'CLCrankNicIntegrator.cl') + kernel_src = resource_string(__name__, 'CLCrankNicIntegrator.cl').decode() # substitute user defined kernel code, and number of signals kernel_src = kernel_src % {'sigKernel': sigRateKernel, 'specKernel': specRateKernel, @@ -292,7 +293,7 @@ def dydt(self): def step(self, dt): if dt!=self.dt: - print "I can only integrate at fixed dt!" + print("I can only integrate at fixed dt!") return self.nCells = len(self.cellStates) @@ -301,9 +302,9 @@ def step(self, dt): s = self.specLevel[self.nCells-1] except IndexError: # Could resize here, then would have to rebuild views - print "Number of cells exceeded " \ + print("Number of cells exceeded " \ + self.__class__.__name__ \ - + "::maxCells (" + self.maxCells + ")" + + "::maxCells (" + self.maxCells + ")") self.dataLen = self.signalDataLen + self.nCells*self.nSpecies @@ -343,7 +344,7 @@ def step(self, dt): # c.signals = self.signalling.signals(c, self.signalLevel) # Update cellType array - for (id,c) in self.cellStates.items(): + for (id,c) in list(self.cellStates.items()): self.celltype[c.idx] = numpy.int32(c.cellType) self.celltype_dev.set(self.celltype) @@ -366,7 +367,7 @@ def setLevels(self, SSLevel, cellSigData): self.specLevel_dev.set(self.specLevel) self.cellSigLevels_dev.set(self.cellSigLevels) cs = self.cellStates - for id,c in cs.items(): #make sure everything is correct here + for id,c in list(cs.items()): #make sure everything is correct here c.species = self.specLevel[c.idx,:] c.signals = self.cellSigLevels[c.idx,:] self.celltype[c.idx] = numpy.int32(c.cellType) diff --git a/CellModeller/Integration/CLEulerIntegrator.py b/CellModeller/Integration/CLEulerIntegrator.py index 9b630d2b..3db4d153 100644 --- a/CellModeller/Integration/CLEulerIntegrator.py +++ b/CellModeller/Integration/CLEulerIntegrator.py @@ -36,7 +36,7 @@ def __init__(self, sim, nSpecies, maxCells, regul=None): # set the species for existing states to views of the levels array cs = self.cellStates - for id,c in cs.items(): + for id,c in list(cs.items()): c.species = self.specLevel[c.idx,:] @@ -95,8 +95,7 @@ def initKernels(self): # Get user defined kernel source specRateKernel = self.regul.specRateCL() from pkg_resources import resource_string - kernel_src = resource_string(__name__, 'CLEulerIntegrator.cl') - #kernel_src = open('CellModeller/Integration/CLEulerIntegrator.cl', 'r').read() + kernel_src = resource_string(__name__, 'CLEulerIntegrator.cl').decode() # substitute user defined kernel code, and number of signals kernel_src = kernel_src%(specRateKernel) self.program = cl.Program(self.context, kernel_src).build(cache_dir=False) @@ -119,7 +118,7 @@ def dydt(self): def step(self, dt): if dt!=self.dt: - print "I can only integrate at fixed dt!" + print("I can only integrate at fixed dt!") return self.nCells = len(self.cellStates) @@ -128,15 +127,15 @@ def step(self, dt): s = self.specLevel[self.nCells-1] except IndexError: # Could resize here, then would have to rebuild views - print "Number of cells exceeded " \ + print("Number of cells exceeded " \ + self.__class__.__name__ \ - + "::maxCells (" + self.maxCells + ")" + + "::maxCells (" + self.maxCells + ")") self.dataLen = self.nCells*self.nSpecies self.cellStates = self.sim.cellStates cs = self.cellStates - for id,c in cs.items(): + for id,c in list(cs.items()): self.effgrow[c.idx] = numpy.float32(c.effGrowth) self.effgrow_dev.set(self.effgrow) @@ -160,7 +159,7 @@ def setLevels(self, specLevel): self.makeViews() self.specLevel_dev.set(self.specLevel) cs = self.cellStates - for id,c in cs.items(): + for id,c in list(cs.items()): c.species = self.specLevel[c.idx,:] self.celltype[c.idx] = numpy.int32(c.cellType) self.celltype_dev.set(self.celltype) diff --git a/CellModeller/Integration/CLEulerSigIntegrator.py b/CellModeller/Integration/CLEulerSigIntegrator.py index 8986e388..39a972b9 100644 --- a/CellModeller/Integration/CLEulerSigIntegrator.py +++ b/CellModeller/Integration/CLEulerSigIntegrator.py @@ -7,6 +7,7 @@ import pyopencl.array as cl_array from pyopencl.array import vec import math +from functools import reduce def unique_stable(ar, return_index=False, return_inverse=False): @@ -62,8 +63,6 @@ def __init__(self, sim, nSignals, nSpecies, maxCells, sig, regul=None, boundcond self.regul = regul self.boundcond = boundcond - self.cellStates = sim.cellStates - self.nCells = len(self.cellStates) self.nSpecies = nSpecies self.nSignals = nSignals @@ -98,10 +97,18 @@ def __init__(self, sim, nSignals, nSpecies, maxCells, sig, regul=None, boundcond self.initArrays() #self.initKernels() + self.setCellStates(sim.cellStates) + + + def setCellStates(self, cs): # set the species for existing states to views of the levels array - cs = self.cellStates - for id,c in cs.items(): + self.cellStates = cs + self.nCells = len(self.cellStates) + self.cellSigLevels[:,:] = 0 + for id,c in list(cs.items()): c.species = self.specLevel[c.idx,:] + c.signals = self.cellSigLevels[c.idx,:] + self.celltype[c.idx] = numpy.int32(c.cellType) def makeViews(self): @@ -177,7 +184,7 @@ def initKernels(self): sigRateKernel = self.regul.sigRateCL() #kernel_src = open('CellModeller/Integration/CLCrankNicIntegrator.cl', 'r').read() from pkg_resources import resource_string - kernel_src = resource_string(__name__, 'CLCrankNicIntegrator.cl') + kernel_src = resource_string(__name__, 'CLCrankNicIntegrator.cl').decode() # substitute user defined kernel code, and number of signals kernel_src = kernel_src % {'sigKernel': sigRateKernel, 'specKernel': specRateKernel, @@ -266,7 +273,7 @@ def dydt(self): def step(self, dt): if dt!=self.dt: - print "I can only integrate at fixed dt!" + print("I can only integrate at fixed dt!") return self.nCells = len(self.cellStates) @@ -275,9 +282,9 @@ def step(self, dt): s = self.specLevel[self.nCells-1] except IndexError: # Could resize here, then would have to rebuild views - print "Number of cells exceeded " \ + print("Number of cells exceeded " \ + self.__class__.__name__ \ - + "::maxCells (" + self.maxCells + ")" + + "::maxCells (" + self.maxCells + ")") self.dataLen = self.signalDataLen + self.nCells*self.nSpecies @@ -312,7 +319,7 @@ def step(self, dt): # c.signals = self.signalling.signals(c, self.signalLevel) # Update cellType array - for (id,c) in self.cellStates.items(): + for (id,c) in list(self.cellStates.items()): self.celltype[c.idx] = numpy.int32(c.cellType) self.celltype_dev.set(self.celltype) @@ -335,7 +342,7 @@ def setLevels(self, SSLevel, cellSigData): self.specLevel_dev.set(self.specLevel) self.cellSigLevels_dev.set(self.cellSigLevels) cs = self.cellStates - for id,c in cs.items(): #make sure everything is correct here + for id,c in list(cs.items()): #make sure everything is correct here c.species = self.specLevel[c.idx,:] c.signals = self.cellSigLevels[c.idx,:] self.celltype[c.idx] = numpy.int32(c.cellType) diff --git a/CellModeller/Regulation/ModuleRegulator.py b/CellModeller/Regulation/ModuleRegulator.py index 009ec390..4ef2bf1c 100644 --- a/CellModeller/Regulation/ModuleRegulator.py +++ b/CellModeller/Regulation/ModuleRegulator.py @@ -40,7 +40,7 @@ def speciesRates(self, cstate, speciesLevels, signalLevels): return self.module.speciesRates(cstate, speciesLevels, signalLevels) def initSpeciesLevels(self, levels): - csv = self.cellStates.values() + csv = list(self.cellStates.values()) nCells = len(csv) for i in range(nCells): levels[i,:] = csv[i].species @@ -49,8 +49,8 @@ def step(self, dt=0): try: self.module.update(self.cellStates) except Exception as e: - print "Problem with regulation module " + self.modName - print e + print("Problem with regulation module " + self.modName) + print(e) def divide(self, pState, d1State, d2State): # Call the module's optional divide function diff --git a/CellModeller/Regulation/SBMLImport.py b/CellModeller/Regulation/SBMLImport.py index 4d7d230a..6bdcb3e6 100644 --- a/CellModeller/Regulation/SBMLImport.py +++ b/CellModeller/Regulation/SBMLImport.py @@ -1,7 +1,7 @@ from libsbml import * import math import new -import urllib2 +import urllib.request, urllib.error, urllib.parse # indicates whether we should perform SBML document consistency check. @@ -11,8 +11,8 @@ # get the array index of given species ID: def getSpecIndex(specIdToArrayIndexDict, specId): - if not specIdToArrayIndexDict.has_key(specId): - specIdToArrayIndexDict[specId] = len(specIdToArrayIndexDict.keys()) + if specId not in specIdToArrayIndexDict: + specIdToArrayIndexDict[specId] = len(list(specIdToArrayIndexDict.keys())) return specIdToArrayIndexDict[specId] @@ -31,7 +31,7 @@ def checkSBMLConsistency(document): if numFailures > 0: failureMsg = "" for failureNum in range(numFailures): - print "Failure " + str(failureNum) + ": " + document.getError(failureNum).getShortMessage() + "\n" + print("Failure " + str(failureNum) + ": " + document.getError(failureNum).getShortMessage() + "\n") raise Exception(failureMsg) @@ -40,11 +40,11 @@ def SBMLModelFromSBMLFile(sbmlFile): reader = SBMLReader() document = reader.readSBML(sbmlFile) if document.getNumErrors()>0: - print "Errors in reading SBML file..." + print("Errors in reading SBML file...") checkSBMLConsistency(document) model = document.getModel() if not model: - print "No model!" + print("No model!") return model @@ -120,7 +120,7 @@ def pythonMathFromASTNode(astNode, kineticLaw, model): # something weird -- return 1 and issue warning: else: - print "WARNING: unknown identifier " + nodeName + ", defaulting to 1. Globals: " + str(model.getNumParameters()) + ". Locals: " + str(kineticLaw.getNumParameters()) + print("WARNING: unknown identifier " + nodeName + ", defaulting to 1. Globals: " + str(model.getNumParameters()) + ". Locals: " + str(kineticLaw.getNumParameters())) return "1" else: @@ -184,7 +184,7 @@ def pythonStrFromSBMLModel(model): reactantSpecId = reactantSpecRef.getSpecies() # add key to hash table if it doesn't exist: - if not specIdToExpDict.has_key(reactantSpecId): + if reactantSpecId not in specIdToExpDict: specIdToExpDict[reactantSpecId] = "" if stoich == 1: @@ -204,7 +204,7 @@ def pythonStrFromSBMLModel(model): productSpecId = productSpecRef.getSpecies() # add key to hash table if it doesn't exist: - if not specIdToExpDict.has_key(productSpecId): + if productSpecId not in specIdToExpDict: specIdToExpDict[productSpecId] = "" if stoich == 1: @@ -217,7 +217,7 @@ def pythonStrFromSBMLModel(model): pythonVarDefs = "" pythonRateUpdateExps = "" specIdLst = [0]*len(specIdToExpDict) - for specId in specIdToExpDict.keys(): + for specId in list(specIdToExpDict.keys()): # get index into species array: index = getSpecIndex(specIdToArrayIndexDict, specId) @@ -302,13 +302,13 @@ def pythonStrFromSBMLModel(model): # create python module from code dynamically. def pythonModuleFromPythonStr(pythonStr): module = new.module("sbmlPythonEncoding") - exec pythonStr in module.__dict__ + exec(pythonStr, module.__dict__) return module # get text from given channel. def TextFromChannel(channelName): - f = urllib2.urlopen('http://async-message-passer.appspot.com/?content_only=1&channel_name=' + channelName) + f = urllib.request.urlopen('http://async-message-passer.appspot.com/?content_only=1&channel_name=' + channelName) return f.read() @@ -318,7 +318,7 @@ def pythonModuleFromFile(fileName): checkSBMLConsistency(document) model = document.getModel() pstr = pythonStrFromSBMLModel(model) - print "Got python str: " + pstr + print("Got python str: " + pstr) pythonModule = pythonModuleFromPythonStr(pstr) return pythonModule @@ -344,7 +344,7 @@ def pythonModuleFromChannel(channelName): pythonStr = pythonStrFromSBMLModel(model) - print "Got python str: " + pythonStr + print("Got python str: " + pythonStr) pythonModule = pythonModuleFromPythonStr(pythonStr) return pythonModule diff --git a/CellModeller/Signalling/GridDiffusion.py b/CellModeller/Signalling/GridDiffusion.py index f8292dd2..f863e897 100644 --- a/CellModeller/Signalling/GridDiffusion.py +++ b/CellModeller/Signalling/GridDiffusion.py @@ -2,6 +2,7 @@ import numpy from scipy.ndimage import laplace from scipy.ndimage.filters import convolve +from functools import reduce class GridDiffusion: @@ -39,8 +40,10 @@ def __init__(self, sim, nSignals, gridDim, gridSize, gridOrig, D, adv=None, init self.dV = reduce(lambda x, y: x * y, self.gridSize) self.regul = regul - self.cellStates = sim.cellStates + self.setCellStates(sim.cellStates) + def setCellStates(self, cs): + self.cellStates = cs def setBiophysics(self, biophysics): self.biophys = biophysics @@ -49,8 +52,8 @@ def setRegulator(self, regul): self.regul = regul def flattenIdx(self, idx): - print "idx = " + str(idx) - print "flat idx = " + str(idx[2] + idx[1]*self.gridDim[3] + idx[0]*self.gridDim[2]*self.gridDim[3]) + print("idx = " + str(idx)) + print("flat idx = " + str(idx[2] + idx[1]*self.gridDim[3] + idx[0]*self.gridDim[2]*self.gridDim[3])) return idx[2] + idx[1]*self.gridDim[3] + idx[0]*self.gridDim[2]*self.gridDim[3] def idxFromPos(self, p): @@ -76,7 +79,7 @@ def trilinearWeights(self, p): w[1,0,1] = dx*(1-dy)*dz if x>=-1 and x=0 and y=-1 and z=0 and x=-1 and y=-1 and z=-1 and x=-1 and y=-1 and zidmax: @@ -273,7 +267,7 @@ def reset(self): if not self.moduleStr: #This will take up any changes made in the model file - reload(self.module) + importlib.reload(self.module) else: # TJR: Module loaded from pickle, cannot reset? pass @@ -349,7 +343,7 @@ def addCell(self, cellType=0, cellAdh=0, length=3.5, **kwargs): # Eventually prob better to have a generic editCell() that deals with this stuff # def moveCell(self, cid, delta_pos): - if self.cellStates.has_key(cid): + if cid in self.cellStates: self.phys.moveCell(self.cellStates[cid], delta_pos) ## Proceed to the next simulation step @@ -357,7 +351,8 @@ def moveCell(self, cid, delta_pos): def step(self): self.reg.step(self.dt) states = dict(self.cellStates) - for (cid,state) in states.items(): + for (cid,state) in list(states.items()): + state.time = self.stepNum * self.dt if state.divideFlag: self.divide(state) #neighbours no longer current @@ -399,18 +394,40 @@ def writePickle(self, csv=False): data['moduleStr'] = self.moduleOutput data['moduleName'] = self.moduleName if self.integ: - print("Writing new pickle format") + # print("Writing new pickle format") data['specData'] = self.integ.levels - data['sigGrid'] = self.integ.signalLevel + if self.sig: data['sigGridOrig'] = self.sig.gridOrig data['sigGridDim'] = self.sig.gridDim data['sigGridSize'] = self.sig.gridSize - if self.sig: + if self.sig and self.integ: + data['sigGrid'] = self.integ.signalLevel data['sigData'] = self.integ.cellSigLevels data['sigGrid'] = self.integ.signalLevel - cPickle.dump(data, outfile, protocol=-1) + pickle.dump(data, outfile, protocol=-1) #output csv file with cell pos,dir,len - sig? + # Populate simulation from saved data pickle + def loadGeometryFromPickle(self, data): + self.setCellStates(data['cellStates']) + self.lineage = data['lineage'] + idx_map = {} + id_map = {} + idmax = 0 + for id,state in data['cellStates'].items(): + idx_map[state.id] = state.idx + id_map[state.idx] = state.id + if id>idmax: + idmax=id + self.idToIdx = idx_map + self.idxToId = id_map + self._next_id = idmax+1 + self._next_idx = len(data['cellStates']) + if self.integ: + self.integ.setCellStates(self.cellStates) + if self.sig: + self.integ.setCellStates(self.cellStates) + # Populate simulation from saved data pickle def loadFromPickle(self, data): self.setCellStates(data['cellStates']) @@ -419,7 +436,7 @@ def loadFromPickle(self, data): idx_map = {} id_map = {} idmax = 0 - for id,state in data['cellStates'].iteritems(): + for id,state in data['cellStates'].items(): idx_map[state.id] = state.idx id_map[state.idx] = state.id if id>idmax: @@ -429,9 +446,9 @@ def loadFromPickle(self, data): self._next_id = idmax+1 self._next_idx = len(data['cellStates']) if self.integ: - if data.has_key('sigData'): + if 'sigData' in data: self.integ.setLevels(data['specData'],data['sigData']) - elif data.has_key('specData'): + elif 'specData' in data: self.integ.setLevels(data['specData']) diff --git a/Examples/ACS2012/EdgeDetectorChamber.py b/Examples/ACS2012/EdgeDetectorChamber.py index 63048088..e8698d97 100644 --- a/Examples/ACS2012/EdgeDetectorChamber.py +++ b/Examples/ACS2012/EdgeDetectorChamber.py @@ -91,9 +91,9 @@ def sigRateCL(): def update(cells): if len(cells) > max_cells-16: - print 'reached cell limit' + print('reached cell limit') exit() - for (i,cell) in cells.iteritems(): + for (i,cell) in cells.items(): cell.color = [0.1, 0.1+cell.species[3]/10.0, 0.1+cell.species[4]*20.0] if cell.volume > getattr(cell, 'target_volume', 3.0): cell.divideFlag = True diff --git a/Examples/Conjugation.py b/Examples/Conjugation.py index a2896e42..b8b9a9a4 100644 --- a/Examples/Conjugation.py +++ b/Examples/Conjugation.py @@ -10,7 +10,7 @@ def setup(sim): # Set biophysics module - biophys = CLBacterium(sim, jitter_z=False, max_cells=50000,gamma=10.0,reg_param=0.1, cgs_tol=1E-5,compNeighbours=True) + biophys = CLBacterium(sim, jitter_z=False, max_cells=50000,gamma=10.0, cgs_tol=1E-5,compNeighbours=True) # Set up regulation module regul = ModuleRegulator(sim, sim.moduleName) @@ -40,7 +40,7 @@ def init(cell): def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/TimRudgeThesis/Meristem.py b/Examples/TimRudgeThesis/Meristem.py index a35645f5..6b116468 100644 --- a/Examples/TimRudgeThesis/Meristem.py +++ b/Examples/TimRudgeThesis/Meristem.py @@ -33,7 +33,7 @@ def init(cell): def max_x_coord(cells): mx = -1000 - for i,cell in cells.items(): + for i,cell in list(cells.items()): mx = max(mx, cell.pos[0]) return mx @@ -43,7 +43,7 @@ def update(cells): global growth_rad rightmost = -1 edge = max_x_coord(cells) - for (i,cell) in cells.items(): + for (i,cell) in list(cells.items()): dist = edge - cell.pos[0] if dist < growth_rad: cell.growthRate = 0.5 #(growth_rad-dist)/growth_rad * 0.5 diff --git a/Examples/Tutorial_1/Tutorial_1a.py b/Examples/Tutorial_1/Tutorial_1a.py index ecb2d2cd..081282b2 100644 --- a/Examples/Tutorial_1/Tutorial_1a.py +++ b/Examples/Tutorial_1/Tutorial_1a.py @@ -33,7 +33,7 @@ def init(cell): def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/Tutorial_1/Tutorial_1b.py b/Examples/Tutorial_1/Tutorial_1b.py index af3118ca..5f529c9d 100644 --- a/Examples/Tutorial_1/Tutorial_1b.py +++ b/Examples/Tutorial_1/Tutorial_1b.py @@ -32,7 +32,7 @@ def setup(sim): def max_y_coord(cells): #finds the largest y-coordinate in the colony my = 0.0 - for i,cell in cells.items(): + for i,cell in list(cells.items()): my = max(my, cell.pos[1]) return my @@ -46,7 +46,7 @@ def init(cell): def update(cells): #Iterate through each cell and flag cells that reach target size for division maxy = max_y_coord(cells) - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): dist = maxy - cell.pos[1] growthZone = 5.0 #width of the growth zone if dist < growthZone: diff --git a/Examples/Tutorial_1/Tutorial_1c.py b/Examples/Tutorial_1/Tutorial_1c.py index 26569e97..1eb30636 100644 --- a/Examples/Tutorial_1/Tutorial_1c.py +++ b/Examples/Tutorial_1/Tutorial_1c.py @@ -37,7 +37,7 @@ def init(cell): def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/Tutorial_2/Tutorial_2a.py b/Examples/Tutorial_2/Tutorial_2a.py index b92da170..8855a2f4 100644 --- a/Examples/Tutorial_2/Tutorial_2a.py +++ b/Examples/Tutorial_2/Tutorial_2a.py @@ -52,7 +52,7 @@ def specRateCL(): def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [numpy.clip(cell.species[0]/6.0,0.0,1.0), 1.0, 0.1] if cell.volume > cell.targetVol: a = 1 diff --git a/Examples/Tutorial_2/Tutorial_2b.py b/Examples/Tutorial_2/Tutorial_2b.py index 2eeadd21..b584e5f1 100644 --- a/Examples/Tutorial_2/Tutorial_2b.py +++ b/Examples/Tutorial_2/Tutorial_2b.py @@ -57,7 +57,7 @@ def specRateCL(): def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [numpy.clip((cell.species[0]-0.53)/0.05,0.0,1.0), 0.1, 0.1] if cell.volume > cell.targetVol: a = 1 diff --git a/Examples/Tutorial_3/Tutorial_3.py b/Examples/Tutorial_3/Tutorial_3.py index d724f996..61e13795 100644 --- a/Examples/Tutorial_3/Tutorial_3.py +++ b/Examples/Tutorial_3/Tutorial_3.py @@ -98,7 +98,7 @@ def update(cells): #Iterate through each cell and flag cells that reach target size for division v_max = 0.9 Km = 0.1 #This can be though of as the degree of mutualism - turns mutualism off when set to 0. - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [1,1,1] #[0.1+cell.species[0]/3.0, 0.1+cell.species[1]/3.0, 0.1] if cell.cellType==0: cell.growthRate = 0.1 + v_max * cell.species[1] / (Km + cell.species[1]) diff --git a/Examples/colorWalk_planes_3d.py b/Examples/colorWalk_planes_3d.py index dfb84911..358c934e 100644 --- a/Examples/colorWalk_planes_3d.py +++ b/Examples/colorWalk_planes_3d.py @@ -57,7 +57,7 @@ def numSpecies(): return 0 def update(cells): - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/ex1_simpleGrowth.py b/Examples/ex1_simpleGrowth.py index 9bc8cb28..1b845a75 100644 --- a/Examples/ex1_simpleGrowth.py +++ b/Examples/ex1_simpleGrowth.py @@ -9,7 +9,7 @@ def setup(sim): # Set biophysics, signalling, and regulation models - biophys = CLBacterium(sim, jitter_z=False) + biophys = CLBacterium(sim, jitter_z=False, max_cells=100000, gamma=100.) # use this file for reg too regul = ModuleRegulator(sim, sim.moduleName) @@ -20,14 +20,14 @@ def setup(sim): sim.addCell(cellType=0, pos=(0,0,0), dir=(1,0,0)) # Add some objects to draw the models - if sim.is_gui: - from CellModeller.GUI import Renderers - therenderer = Renderers.GLBacteriumRenderer(sim) - sim.addRenderer(therenderer) - else: - print "Running in batch mode: no display will be output" - - sim.pickleSteps = 10 + #if sim.is_gui: + from CellModeller.GUI import Renderers + therenderer = Renderers.GLBacteriumRenderer(sim) + sim.addRenderer(therenderer) + #else: + # print("Running in batch mode: no display will be output") + + sim.pickleSteps = 100 sim.saveOutput = True def init(cell): @@ -39,10 +39,19 @@ def init(cell): def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): if cell.volume > cell.targetVol: cell.divideFlag = True + gr = cell.strainRate/0.05 + cgr = gr - 0.5 + # Return value is tuple of colors, (fill, stroke) + #if cgr>0: + # cell.color = [1, 1-cgr*2, 1-cgr*2] + #else: + # cell.color = [1+cgr*2, 1+cgr*2, 1] + + def divide(parent, d1, d2): # Specify target cell size that triggers cell division d1.targetVol = 3.5 + random.uniform(0.0,0.5) diff --git a/Examples/ex1_simpleGrowth2D.py b/Examples/ex1_simpleGrowth2D.py index 251c9ea3..406b63a2 100644 --- a/Examples/ex1_simpleGrowth2D.py +++ b/Examples/ex1_simpleGrowth2D.py @@ -33,7 +33,7 @@ def init(cell): def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [cell.cellType*0.6+0.1, 1.0-cell.cellType*0.6, 0.3] if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/ex1a_simpleGrowth2D.py b/Examples/ex1a_simpleGrowth2D.py index f195e809..a8d22413 100644 --- a/Examples/ex1a_simpleGrowth2D.py +++ b/Examples/ex1a_simpleGrowth2D.py @@ -36,7 +36,7 @@ def init(cell): def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [cell.cellType*0.6+0.1, 1.0-cell.cellType*0.6, 0.3] if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/ex1a_simpleGrowth2Types.py b/Examples/ex1a_simpleGrowth2Types.py index 3fd04709..20e9d247 100644 --- a/Examples/ex1a_simpleGrowth2Types.py +++ b/Examples/ex1a_simpleGrowth2Types.py @@ -36,7 +36,7 @@ def init(cell): def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [cell.cellType*0.6+0.1, 1.0-cell.cellType*0.6, 0.3] if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/ex1b_simpleGrowth2D.py b/Examples/ex1b_simpleGrowth2D.py index 7f244927..285c60e7 100644 --- a/Examples/ex1b_simpleGrowth2D.py +++ b/Examples/ex1b_simpleGrowth2D.py @@ -36,7 +36,7 @@ def init(cell): def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [cell.cellType*0.6+0.1, 1.0-cell.cellType*0.6, 0.3] if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/ex1b_simpleGrowthRoundCell.py b/Examples/ex1b_simpleGrowthRoundCell.py index f96feea5..e65b7c8b 100644 --- a/Examples/ex1b_simpleGrowthRoundCell.py +++ b/Examples/ex1b_simpleGrowthRoundCell.py @@ -8,7 +8,7 @@ def setup(sim): # Set biophysics, signalling, and regulation models - biophys = CLBacterium(sim, jitter_z=False) + biophys = CLBacterium(sim, jitter_z=True) # use this file for reg too regul = ModuleRegulator(sim, sim.moduleName) @@ -34,7 +34,7 @@ def init(cell): def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [cell.cellType*0.6+0.1, 1.0-cell.cellType*0.6, 0.3] if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/ex2_constGene.py b/Examples/ex2_constGene.py index ff21cb45..baa98232 100644 --- a/Examples/ex2_constGene.py +++ b/Examples/ex2_constGene.py @@ -53,7 +53,7 @@ def specRateCL(): # Add def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [0.1+cell.species[0]/20.0, 0.1, 0.1] # Add/change if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/ex2a_dilution.py b/Examples/ex2a_dilution.py index 93acae3d..013ae89c 100644 --- a/Examples/ex2a_dilution.py +++ b/Examples/ex2a_dilution.py @@ -52,7 +52,7 @@ def specRateCL(): # Add def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [0.1+cell.species[0]*10.0, 0.1, 0.1] # Add/change if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/ex2b_diluteRepression.py b/Examples/ex2b_diluteRepression.py index e203bf91..1b698365 100644 --- a/Examples/ex2b_diluteRepression.py +++ b/Examples/ex2b_diluteRepression.py @@ -57,7 +57,7 @@ def specRateCL(): # Add def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [0.1+cell.species[0], 0.1+cell.species[1]*0.1, 0.1] # Add/change if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/ex3_simpleSignal.py b/Examples/ex3_simpleSignal.py index 6dc58d5f..be408c61 100644 --- a/Examples/ex3_simpleSignal.py +++ b/Examples/ex3_simpleSignal.py @@ -89,7 +89,7 @@ def sigRateCL(): #Add def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [0.1+cell.species[0]/20.0, 0.1, 0.1] if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/ex4_simpleCellCellSignaling.py b/Examples/ex4_simpleCellCellSignaling.py index b88285fb..3eed694f 100644 --- a/Examples/ex4_simpleCellCellSignaling.py +++ b/Examples/ex4_simpleCellCellSignaling.py @@ -99,7 +99,7 @@ def sigRateCL(): #Add def update(cells): #Iterate through each cell and flag cells that reach target size for division - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [0.1+cell.species[1]/20.0, 0.1+cell.species[2]/20.0, 0.1] if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/ex5_colonySector.py b/Examples/ex5_colonySector.py index 35fad734..97965b7f 100644 --- a/Examples/ex5_colonySector.py +++ b/Examples/ex5_colonySector.py @@ -1,36 +1,41 @@ import random from CellModeller.Regulation.ModuleRegulator import ModuleRegulator from CellModeller.Biophysics.BacterialModels.CLBacterium import CLBacterium +from CellModeller.GUI import Renderers import numpy import math +N0 = 10 def setup(sim): # Set biophysics, signalling, and regulation models - biophys = CLBacterium(sim, jitter_z=False) + biophys = CLBacterium(sim, jitter_z=False, gamma = 100, max_cells=100000, max_planes=1) regul = ModuleRegulator(sim, sim.moduleName) # use this file for reg too # Only biophys and regulation sim.init(biophys, regul, None, None) - sim.addCell(cellType=0, pos=(0,0,0)) + #biophys.addPlane((0,0,0),(0,0,1),1.0) #Base plane + #biophys.addPlane((10,0,0),(-1,0,0),1.0) + #biophys.addPlane((-10,0,0),(1,0,0),1.0) + #biophys.addPlane((0,10,0),(0,-1,0),1.0) + #biophys.addPlane((0,-10,0),(0,1,0),1.0) - if sim.is_gui: - # Add some objects to draw the models - from CellModeller.GUI import Renderers - therenderer = Renderers.GLBacteriumRenderer(sim) - sim.addRenderer(therenderer) + sim.addCell(cellType=0, pos=(0,0,0)) - sim.pickleSteps = 10 + # Add some objects to draw the models + therenderer = Renderers.GLBacteriumRenderer(sim) + sim.addRenderer(therenderer) + sim.pickleSteps = 1 def init(cell): cell.targetVol = 3.5 + random.uniform(0.0,0.5) cell.growthRate = 1.0 - cell.n_a = 3 - cell.n_b = 3 + cell.n_a = N0//2 + cell.n_b = N0 - cell.n_a def update(cells): - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [0.1, cell.n_a/3.0, cell.n_b/3.0] if cell.volume > cell.targetVol: cell.divideFlag = True @@ -44,15 +49,15 @@ def divide(parent, d1, d2): d1.n_b = 0 d2.n_a = 0 d2.n_b = 0 - for p in plasmids[:6]: + for p in plasmids[:N0]: if p == 0: d1.n_a +=1 else: d1.n_b +=1 - for p in plasmids[6:12]: + for p in plasmids[N0:2*N0]: if p == 0: d2.n_a +=1 else: d2.n_b +=1 - assert parent.n_a + parent.n_b == 6 - assert d1.n_a + d1.n_b == 6 - assert d2.n_a + d2.n_b == 6 + assert parent.n_a + parent.n_b == N0 + assert d1.n_a + d1.n_b == N0 + assert d2.n_a + d2.n_b == N0 assert parent.n_a*2 == d1.n_a+d2.n_a assert parent.n_b*2 == d1.n_b+d2.n_b assert parent.n_a > 0 or (d1.n_a == 0 and d2.n_a == 0) diff --git a/Examples/ex5_colonySector_3d.py b/Examples/ex5_colonySector_3d.py index 4dfecde4..9848c315 100644 --- a/Examples/ex5_colonySector_3d.py +++ b/Examples/ex5_colonySector_3d.py @@ -32,7 +32,7 @@ def init(cell): cell.n_b = 3 def update(cells): - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [0.1, cell.n_a/3.0, cell.n_b/3.0] if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/load.py b/Examples/load.py index e679f53c..01b7470b 100644 --- a/Examples/load.py +++ b/Examples/load.py @@ -19,7 +19,6 @@ def setup(sim): max_contacts=32, \ max_sqs=128**2, \ jitter_z=False, \ - reg_param=2, \ gamma=10) biophys.addPlane((0,0,-0.5), (0,0,1), 0.2) @@ -34,13 +33,13 @@ def setup(sim): # Or pop up a dialog to choose a pickle if dataFileName=='': - print "Please edit the model file to specify a pickle file." - print " -- No data was loaded and there are no cells in this simulation!" + print("Please edit the model file to specify a pickle file.") + print(" -- No data was loaded and there are no cells in this simulation!") else: # Import the data and load into Simulator - print "Loading data from pickle file: %s"%(dataFileName) - import cPickle - data = cPickle.load(open(dataFileName,'r')) + print("Loading data from pickle file: %s"%(dataFileName)) + import pickle + data = pickle.load(open(dataFileName,'r')) sim.loadFromPickle(data) if sim.is_gui: @@ -54,7 +53,7 @@ def init(cell): cell.growthRate = 0.5 def update(cells): - for (id, cell) in cells.iteritems(): + for (id, cell) in cells.items(): cell.color = [1,0,0] if cell.volume > cell.targetVol: cell.divideFlag = True diff --git a/Examples/sphere_constraints.py b/Examples/sphere_constraints.py new file mode 100644 index 00000000..7d7e2d5f --- /dev/null +++ b/Examples/sphere_constraints.py @@ -0,0 +1,53 @@ +import random +from CellModeller.Regulation.ModuleRegulator import ModuleRegulator +from CellModeller.Biophysics.BacterialModels.CLBacterium import CLBacterium +from CellModeller.GUI import Renderers +import numpy as np +import math + +max_cells = 50000 + + +def setup(sim): + # Set biophysics, signalling, and regulation models + + #max_sqr + #biophys = CLBacterium(sim, max_cells=max_cells, jitter_z=False, max_sqs=256**2) + biophys = CLBacterium(sim, max_cells=max_cells, jitter_z=False, max_spheres=1, gamma=100.) + + # use this file for reg too + regul = ModuleRegulator(sim) + # Only biophys and regulation + sim.init(biophys, regul, None, None) + + # Specify the initial cell and its location in the simulation + sim.addCell(cellType=0, pos=(0,0,0)) + + biophys.addSphere((0,0,0), 10, 1.0, -1) + + + # Add some objects to draw the models + therenderer = Renderers.GLBacteriumRenderer(sim) + sim.addRenderer(therenderer) + + sim.pickleSteps = 5 + +def init(cell): + # Specify mean and distribution of initial cell size + cell.targetVol = 3.0 + random.uniform(0.0,0.5) + # Specify growth rate of cells + #cell.growthRate = 0.6 + cell.growthRate = 1.0 + cell.color = [0, 1, 0] + +def update(cells): + #Iterate through each cell and flag cells that reach target size for division + for (id, cell) in cells.items(): + if cell.volume > cell.targetVol: + cell.divideFlag = True + +def divide(parent, d1, d2): + # Specify target cell size that triggers cell division + d1.targetVol = 3.0 + random.uniform(0.0,0.5) + d2.targetVol = 3.0 + random.uniform(0.0,0.5) + diff --git a/Scripts/CellModellerGUI.py b/Scripts/CellModellerGUI.py index 14f04325..23f3b699 100644 --- a/Scripts/CellModellerGUI.py +++ b/Scripts/CellModellerGUI.py @@ -6,8 +6,9 @@ # Tim Rudge # Jan 2011 -from PyQt4.QtGui import QApplication -from PyQt4 import uic +from PyQt5.QtWidgets import QApplication +from PyQt5 import uic +from PyQt5.QtCore import * import CellModeller.GUI.Renderers from CellModeller import Simulator @@ -26,6 +27,11 @@ ui.show() ui.raise_() cmv = ui.PyGLCMViewer +pix_ratio = qapp.devicePixelRatio() +cmv.setPixelRatio(pix_ratio) +label = ui.label +label.setTextFormat(Qt.RichText) +label.setAlignment(Qt.AlignJustify) # Load a model if specified if len(sys.argv) > 1: cmv.loadModelFile(sys.argv[1]) diff --git a/Scripts/Draw2DPDF.py b/Scripts/Draw2DPDF.py index 34a64488..af920463 100644 --- a/Scripts/Draw2DPDF.py +++ b/Scripts/Draw2DPDF.py @@ -8,7 +8,7 @@ from reportlab.lib.colors import Color import numpy -import cPickle +import pickle mxsig0 = 0 @@ -19,13 +19,13 @@ class CellModellerPDFGenerator(Canvas): def __init__(self, name, data, bg_color): self.name = name self.states = data.get('cellStates') - self.signals = False - if 'specData' in data: - self.signals = True - self.signal_levels = data.get('sigGrid') - self.signal_grid_orig = data.get('sigGridOrig') - self.signal_grid_dim = data.get('sigGridDim') - self.signal_grid_size = data.get('sigGridSize') + + self.signal_levels = data.get('sigGrid', None) + self.signal_grid_orig = data.get('sigGridOrig', None) + self.signal_grid_dim = data.get('sigGridDim', None) + self.signal_grid_size = data.get('sigGridSize', None) + self.signals = self.signal_levels is not None + self.parents = data.get('lineage') self.data = data self.bg_color = bg_color @@ -79,7 +79,7 @@ def draw_capsule(self, p, d, l, r, fill_color, stroke_color): self.restoreState() def draw_cells(self): - for id, state in self.states.items(): + for id, state in list(self.states.items()): p = state.pos d = state.dir l = state.length @@ -93,7 +93,7 @@ def draw_chamber(self): self.line(-100, -16, 100, -16) self.line(-100, 16, 100, 16) - def draw_signals(self, index=0, scale=0.0192, z=2): + def draw_signals(self, index=0, scale=0.0192, z=0): ''' Function for drawing signal grids, currently limited to 1 signal a plane at a fixed z-axis level through the grid @@ -109,12 +109,13 @@ def draw_signals(self, index=0, scale=0.0192, z=2): self.signal_grid_dim, \ self.signal_levels levels = levels.reshape(dim) - l = map(float,l) + mx = levels[index,:,:,z].max() + l = list(map(float,l)) for i in range(dim[1]): x = l[0]*i + orig[0] for j in range(dim[2]): y = l[1]*j + orig[1] - lvls = levels[index,i,j,z]/scale + lvls = levels[index,i,j,z]/mx mxsig0 = max(lvls, mxsig0) self.setFillColorRGB(lvls, 0, 0) self.rect(x-l[0]/2.0, y-l[1]/2.0, l[0], l[1], stroke=0, fill=1) @@ -140,7 +141,7 @@ def computeBox(self): mnx = -20 mxy = 20 mny = -20 - for (id,s) in self.states.iteritems(): + for (id,s) in list(self.states.items()): pos = s.pos l = s.length # add/sub length to keep cell in frame mxx = max(mxx,pos[0]+l) @@ -155,8 +156,8 @@ def computeBox(self): def importPickle(fname): if fname[-7:]=='.pickle': - print 'Importing CellModeller pickle file: %s'%fname - data = cPickle.load(open(fname, 'rb')) + print(('Importing CellModeller pickle file: %s'%fname)) + data = pickle.load(open(fname, 'rb')) # Check for old-style pickle that is tuple, # just extract cellStates from 1st element @@ -182,24 +183,24 @@ def main(): # e.g. outline color, page size, etc. # # For now, put these options into variables here: - bg_color = Color(1.0,1.0,1.0,alpha=1.0) + bg_color = Color(0,0,0,alpha=1.0) # For now just assume a list of files infns = sys.argv[1:] for infn in infns: # File names if infn[-7:]!='.pickle': - print 'Ignoring file %s, because its not a pickle...'%(infn) + print(('Ignoring file %s, because its not a pickle...'%(infn))) continue - outfn = string.replace(infn, '.pickle', '.pdf') + outfn = infn.replace('.pickle', '.pdf') outfn = os.path.basename(outfn) # Put output in this dir - print 'Processing %s to generate %s'%(infn,outfn) + print(('Processing %s to generate %s'%(infn,outfn))) # Import data data = importPickle(infn) if not data: - print "Problem importing data!" + print("Problem importing data!") return # Create a pdf canvas thing @@ -211,14 +212,14 @@ def main(): '''(w,h) = pdf.computeBox() sqrt2 = math.sqrt(2) world = (w/sqrt2,h/sqrt2)''' - world = (250,250) + world = (150,150) # Page setup page = (20,20) center = (0,0) # Render pdf - print 'Rendering PDF output to %s'%outfn + print(('Rendering PDF output to %s'%outfn)) pdf.draw_frame(outfn, world, page, center) if __name__ == "__main__": diff --git a/Scripts/LengthHistogram.py b/Scripts/LengthHistogram.py index 767cfbdd..7253bcfd 100644 --- a/Scripts/LengthHistogram.py +++ b/Scripts/LengthHistogram.py @@ -3,7 +3,7 @@ import math import numpy as np sys.path.append('.') -import cPickle +import pickle import CellModeller #import matplotlib.pyplot as plt @@ -15,12 +15,12 @@ def rad_pos(cellstate): return np.sqrt(cellstate.pos[0]*cellstate.pos[0]+cellstate.pos[1]*cellstate.pos[1]) def lengthHist(pickle, bins, file=False): - print('opening '+ pickle) - data = cPickle.load(open(pickle,'r')) + print(('opening '+ pickle)) + data = pickle.load(open(pickle,'r')) cs = data['cellStates'] it = iter(cs) n = len(cs) - print 'Number of cells = '+str(n) + print(('Number of cells = '+str(n))) lens = [] r = [] for it in cs: @@ -44,7 +44,7 @@ def dirLengthHist(dir,bins,file=False): number = f[5:-7] if file: fout = open('LengthData'+number+'.csv') - print 'step number = ', number + print(('step number = ', number)) n, lens = lengthHist(dir + f,bins,file) if file: fout.close() diff --git a/Scripts/batch.py b/Scripts/batch.py index f94eaf52..c63223a3 100644 --- a/Scripts/batch.py +++ b/Scripts/batch.py @@ -6,7 +6,7 @@ from CellModeller.Simulator import Simulator -max_cells = 10000 +max_cells = 50000 cell_buffer = 256 def simulate(modfilename, platform, device, steps=50): @@ -20,7 +20,7 @@ def simulate(modfilename, platform, device, steps=50): def main(): # Get module name to load if len(sys.argv)<2: - print "Please specify a model (.py) file" + print("Please specify a model (.py) file") exit(0) else: moduleName = sys.argv[1] @@ -31,17 +31,17 @@ def main(): import pyopencl as cl # Platform platforms = cl.get_platforms() - print "Select OpenCL platform:" + print("Select OpenCL platform:") for i in range(len(platforms)): - print 'press '+str(i)+' for '+str(platforms[i]) - platnum = int(input('Platform Number: ')) + print(('press '+str(i)+' for '+str(platforms[i]))) + platnum = int(eval(input('Platform Number: '))) # Device devices = platforms[platnum].get_devices() - print "Select OpenCL device:" + print("Select OpenCL device:") for i in range(len(devices)): - print 'press '+str(i)+' for '+str(devices[i]) - devnum = int(input('Device Number: ')) + print(('press '+str(i)+' for '+str(devices[i]))) + devnum = int(eval(input('Device Number: '))) else: platnum = int(sys.argv[2]) devnum = int(sys.argv[3]) diff --git a/Scripts/batchFile.py b/Scripts/batchFile.py index d62a4a34..50472678 100644 --- a/Scripts/batchFile.py +++ b/Scripts/batchFile.py @@ -10,7 +10,7 @@ sys.path.append('./') -print os.getcwd() +print((os.getcwd())) import CellModeller.AdaptiveSimulator @@ -39,14 +39,14 @@ try: os.mkdir(pickleDir) except OSError: - print pickleDir, 'exists' + print((pickleDir, 'exists')) pickleSetDir = os.path.join(pickleDir, mod_name) try: os.mkdir(pickleSetDir) except: - print pickleSetDir, 'exists' + print((pickleSetDir, 'exists')) folderName = time.strftime('%Y%m%d-%H%M%S', time.localtime()) pickleFileRoot = os.path.join(pickleSetDir, folderName) @@ -56,9 +56,9 @@ def simulate(mod_name, steps=50): - print 'simulate' + print('simulate') sim = Simulator(mod_name, 0.25, None, pickleSteps=50, pickleFileRoot=pickleFileRoot) - print 'start' + print('start') sim.phys.set_cells() sim.phys.calc_cell_geom() diff --git a/Scripts/contactGraph.py b/Scripts/contactGraph.py index d14a4296..2c38ccbd 100644 --- a/Scripts/contactGraph.py +++ b/Scripts/contactGraph.py @@ -2,7 +2,7 @@ import os import math import numpy as np -import cPickle +import pickle import CellModeller import subprocess import string @@ -56,8 +56,8 @@ def draw_graph(self): self.saveState() self.setLineWidth(0.003*units.cm) self.setStrokeColor((0.0,0.0,0.0)) - positions = networkx.get_node_attributes(G,'pos').values() - colors = networkx.get_node_attributes(G,'color').values() + positions = list(networkx.get_node_attributes(G,'pos').values()) + colors = list(networkx.get_node_attributes(G,'color').values()) for u,v,di in self.network.edges_iter(data=True): self.drawEdge(positions[u],positions[v]) for n in self.network.nodes_iter(): @@ -104,13 +104,13 @@ def get_current_contacts(G, data): G = networkx.Graph() fname = sys.argv[1] -data = cPickle.load(open(fname,'rb')) +data = pickle.load(open(fname,'rb')) cs = data['cellStates'] it = iter(cs) n = len(cs) oname = fname.replace('.pickle','_graph.pdf') -print "num_cells = "+str(n) +print(("num_cells = "+str(n))) cell_type={} pos_dict={} @@ -120,13 +120,13 @@ def get_current_contacts(G, data): get_current_contacts(G, data) -print "num_contacts = " + str(networkx.number_of_edges(G)) -degrees = G.degree().values() -print "mean degree = " + str(np.mean(degrees)) +print(("num_contacts = " + str(networkx.number_of_edges(G)))) +degrees = list(G.degree().values()) +print(("mean degree = " + str(np.mean(degrees)))) -if networkx.get_edge_attributes(G,'color').items(): - edges,ecolor = zip(*networkx.get_edge_attributes(G,'color').items()) -ncolor = networkx.get_node_attributes(G,'color').values() +if list(networkx.get_edge_attributes(G,'color').items()): + edges,ecolor = list(zip(*list(networkx.get_edge_attributes(G,'color').items()))) +ncolor = list(networkx.get_node_attributes(G,'color').values()) pdf = CellModellerPDFGenerator(oname, G, bg_color) world = (120,120) diff --git a/Scripts/multi_batch.py b/Scripts/multi_batch.py index 38ecf6fb..781ca5aa 100644 --- a/Scripts/multi_batch.py +++ b/Scripts/multi_batch.py @@ -20,7 +20,7 @@ def simulate(modfilename, platform, device, steps=50): def main(): # Get module name to load if len(sys.argv)<2: - print "Please specify a model (.py) file" + print("Please specify a model (.py) file") exit(0) else: moduleName = sys.argv[1] @@ -31,17 +31,17 @@ def main(): import pyopencl as cl # Platform platforms = cl.get_platforms() - print "Select OpenCL platform:" + print("Select OpenCL platform:") for i in range(len(platforms)): - print 'press '+str(i)+' for '+str(platforms[i]) - platnum = int(input('Platform Number: ')) + print(('press '+str(i)+' for '+str(platforms[i]))) + platnum = int(eval(input('Platform Number: '))) # Device devices = platforms[platnum].get_devices() - print "Select OpenCL device:" + print("Select OpenCL device:") for i in range(len(devices)): - print 'press '+str(i)+' for '+str(devices[i]) - devnum = int(input('Device Number: ')) + print(('press '+str(i)+' for '+str(devices[i]))) + devnum = int(eval(input('Device Number: '))) else: platnum = int(sys.argv[2]) devnum = int(sys.argv[3]) diff --git a/Scripts/printTiming.py b/Scripts/printTiming.py index 38953c48..f6f3996c 100644 --- a/Scripts/printTiming.py +++ b/Scripts/printTiming.py @@ -1,18 +1,18 @@ import sys import os sys.path.append('.') -import cPickle +import pickle import CellModeller dir = os.path.join('data', sys.argv[1]) -print "No. cells\ttime (s)" -print "-------------------" +print("No. cells\ttime (s)") +print("-------------------") for f in os.listdir(dir): if 'pickle' in f: ff = os.path.join(dir,f) - (cs,lin) = cPickle.load(open(ff,'r')) + (cs,lin) = pickle.load(open(ff,'r')) n = len(cs) t = os.path.getmtime(ff) - print "%i\t%f"%(n,t) + print(("%i\t%f"%(n,t))) diff --git a/Scripts/spatial_analyze.py b/Scripts/spatial_analyze.py index b178ea1b..15b86ecb 100644 --- a/Scripts/spatial_analyze.py +++ b/Scripts/spatial_analyze.py @@ -3,12 +3,12 @@ import math import numpy as np sys.path.append('.') -import cPickle +import pickle import CellModeller import matplotlib.pyplot as plt pname= sys.argv[1] -print('opening '+ pname) +print(('opening '+ pname)) fout=open('spatial.txt', "w") fout.write("r spec_level \n") @@ -18,12 +18,12 @@ narray=np.zeros(bin_num) specArray=np.zeros(bin_num) -data = cPickle.load(open(pname,'r')) +data = pickle.load(open(pname,'r')) cs = data['cellStates'] it = iter(cs) n = len(cs) -rArray = np.multiply(range(0,bin_num),rad_max/bin_num) +rArray = np.multiply(list(range(0,bin_num)),rad_max/bin_num) for it in cs: r = np.sqrt(cs[it].pos[0]*cs[it].pos[0]+cs[it].pos[1]*cs[it].pos[1]) diff --git a/Scripts/video.sh b/Scripts/video.sh index 2d7073ec..9e269fae 100755 --- a/Scripts/video.sh +++ b/Scripts/video.sh @@ -7,18 +7,19 @@ # Run Draw2DPDF to generate pdf files for f in $( ls *.pickle ); do echo Processing: $f - $CMPATH/bin/cmpython $CMPATH/Scripts/Draw2DPDF.py $f + python $HOME/Code/CellModeller/Scripts/Draw2DPDF.py $f done # Convert and resize etc. pdf files into jpegs for f in $( ls *.pdf ); do NAME=`basename $f .pdf` convert \ - -verbose \ + -colorspace RGB \ + -verbose \ -density 150 \ $NAME.pdf \ $NAME.png done # Run ffmpeg to generate video file -ffmpeg -framerate 7 -i step-%04d0.png -vf scale=1920:1080 -r 24 $1 +ffmpeg -framerate 7 -i %*.png -vf scale=1920:1080 -r 24 $1 diff --git a/setup.py b/setup.py index 17d59036..b0240950 100644 --- a/setup.py +++ b/setup.py @@ -22,10 +22,8 @@ setup(name='CellModeller', - install_requires=['numpy==1.9.2', 'scipy', 'pyopengl', 'pyopencl==2014.1', 'reportlab', 'matplotlib'], - setup_requires=['numpy==1.9.2', 'scipy', 'pyopengl', 'pyopencl==2014.1', 'reportlab', 'matplotlib'], - #setup_requires=['python=2'], - #setup_requires=['numpy', 'pyopengl', 'mako', 'pyopencl'], + install_requires=['numpy', 'scipy', 'pyopengl', 'mako', 'pyqt5', 'pyopencl', 'reportlab', 'matplotlib'], + setup_requires=['numpy', 'scipy', 'pyopengl', 'mako', 'pyqt5', 'pyopencl', 'reportlab', 'matplotlib'], packages=['CellModeller', 'CellModeller.Biophysics', 'CellModeller.Biophysics.BacterialModels', @@ -35,5 +33,6 @@ 'CellModeller.Signalling', 'CellModeller.GUI'], package_data={'':['*.cl','*.ui']}, + python_requires='>=3', version=str(version_git) )