Changeset 032646d in sasmodels
 Timestamp:
 Feb 2, 2018 8:03:01 AM (6 years ago)
 Branches:
 master, core_shell_microgels, magnetic_model, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 e526a9d
 Parents:
 8bd379a
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

explore/realspace.py
r8bd379a r032646d 9 9 import numpy as np 10 10 from numpy import pi, radians, sin, cos, sqrt 11 from numpy.random import poisson, uniform 11 from numpy.random import poisson, uniform, randn, rand 12 12 from numpy.polynomial.legendre import leggauss 13 13 from scipy.integrate import simps … … 91 91 raise NotImplementedError() 92 92 93 def dims(self): 94 # type: () > float, float, float 95 raise NotImplementedError() 96 93 97 def rotate(self, theta, phi, psi): 94 98 self.rotation = rotation(theta, phi, psi) * self.rotation … … 125 129 for s2 in shapes] 126 130 self.r_max = max(distances + [s.r_max for s in shapes]) 127 128 def volume(self): 129 return sum(shape.volume() for shape in self.shapes) 131 self.volume = sum(shape.volume for shape in self.shapes) 130 132 131 133 def sample(self, density): … … 142 144 self._scale = np.array([a/2, b/2, c/2])[None, :] 143 145 self.r_max = sqrt(a**2 + b**2 + c**2) 144 145 def volume(self): 146 return self.a*self.b*self.c 146 self.dims = a, b, c 147 self.volume = a*b*c 147 148 148 149 def sample(self, density): 149 num_points = poisson(density*self. a*self.b*self.c)150 num_points = poisson(density*self.volume) 150 151 points = self._scale*uniform(1, 1, size=(num_points, 3)) 151 152 values = self.value.repeat(points.shape[0]) … … 161 162 self._scale = np.array([ra, rb, length/2])[None, :] 162 163 self.r_max = sqrt(4*max(ra, rb)**2 + length**2) 163 164 def volume(self): 165 return pi*self.ra*self.rb*self.length 164 self.dims = 2*ra, 2*rb, length 165 self.volume = pi*ra*rb*length 166 166 167 167 def sample(self, density): 168 # density of the bounding box 168 # randomly sample from a box of side length 2*r, excluding anything 169 # not in the cylinder 169 170 num_points = poisson(density*4*self.ra*self.rb*self.length) 170 171 points = uniform(1, 1, size=(num_points, 3)) … … 183 184 self._scale = np.array([ra, rb, rc])[None, :] 184 185 self.r_max = 2*max(ra, rb, rc) 185 186 def volume(self): 187 return 4*pi/3 * self.ra * self.rb * self.rc 186 self.dims = 2*ra, 2*rb, 2*rc 187 self.volume = 4*pi/3 * ra * rb * rc 188 188 189 189 def sample(self, density): … … 203 203 self.rotate(*orientation) 204 204 self.shift(*center) 205 helix_length = helix_pitch * tube_length/sqrt(helix_radius**2 + helix_pitch**2) 206 total_radius = self.helix_radius + self.tube_radius 205 207 self.helix_radius, self.helix_pitch = helix_radius, helix_pitch 206 208 self.tube_radius, self.tube_length = tube_radius, tube_length 207 helix_length = helix_pitch * tube_length/sqrt(helix_radius**2 + helix_pitch**2) 208 self.r_max = sqrt((2*helix_radius + 2*tube_radius)*2 209 + (helix_length + 2*tube_radius)**2) 210 211 def volume(self): 209 self.r_max = sqrt(4*total_radius + (helix_length + 2*tube_radius)**2) 210 self.dims = 2*total_radius, 2*total_radius, helix_length 212 211 # small tube radius approximation; for larger tubes need to account 213 212 # for the fact that the inner length is much shorter than the outer 214 213 # length 215 returnpi*self.tube_radius**2*self.tube_length214 self.volume = pi*self.tube_radius**2*self.tube_length 216 215 217 216 def points(self, density): … … 262 261 side_c2 = copy(side_c).shift(0, 0, cdc) 263 262 shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 263 shape.dims = 2*da+a, 2*db+b, 2*dc+c 264 264 return shape 265 265 … … 497 497 Fq = sas_2J1x_x(qab*radius) * sas_sinx_x(qc*length/2) 498 498 Iq[k] = np.sum(w*Fq**2) 499 Iq = Iq /Iq[0]499 Iq = Iq 500 500 return Iq 501 501 … … 509 509 def sphere_Iq(q, radius): 510 510 Iq = sas_3j1x_x(q*radius)**2 511 return Iq /Iq[0]511 return Iq 512 512 513 513 def box_Iq(q, a, b, c): … … 527 527 outer_sum += outer_w * inner_sum 528 528 Iq = outer_sum / 4 # = outer*um*zm*8.0/(4.0*M_PI) 529 return Iq /Iq[0]529 return Iq 530 530 531 531 def box_Iqxy(qx, qy, a, b, c, view=(0, 0, 0)): … … 598 598 #  Test cases  599 599 600 def check_cylinder(radius=25, length=125, rho=2.):600 def build_cylinder(radius=25, length=125, rho=2.): 601 601 shape = EllipticalCylinder(radius, radius, length, rho) 602 602 fn = lambda q: cylinder_Iq(q, radius, length)*rho**2 … … 604 604 return shape, fn, fn_xy 605 605 606 def check_cylinder_lattice(radius=25, length=125, rho=2.): 607 nx, dx = 1, 2*radius 608 ny, dy = 30, 2*radius 609 nz, dz = 30, length 610 dx, dy, dz = 2*dx, 2*dy, 2*dz 611 def center(*args): 612 sigma = 0.333 613 space = 2 614 return [(space*n+np.random.randn()*sigma)*x for n, x in args] 615 t0 = time.time() 616 shapes = [EllipticalCylinder(radius, radius, length, rho, 617 #center=(ix*dx, iy*dy, iz*dz) 618 orientation=np.random.randn(3)*0, 619 center=center((ix, dx), (iy, dy), (iz, dz)) 620 ) 621 for ix in range(nx) 622 for iy in range(ny) 623 for iz in range(nz)] 624 shape = Composite(shapes) 625 print("generate points time", time.time()  t0) 626 fn = None 627 fn_xy = lambda qx, qy, view: cylinder_Iqxy(qx, qy, radius, length, view=view) 628 return shape, fn, fn_xy 629 630 def check_sphere(radius=125, rho=2): 606 def build_sphere(radius=125, rho=2): 631 607 shape = TriaxialEllipsoid(radius, radius, radius, rho) 632 fn = lambda q: cylinder_Iq(q, radius, length)*rho**2608 fn = lambda q: sphere_Iq(q, radius)*rho**2 633 609 fn_xy = lambda qx, qy, view: sphere_Iq(np.sqrt(qx**2+qy**2), radius)*rho**2 634 610 return shape, fn, fn_xy 635 611 636 def check_box(a=10, b=20, c=30, rho=2.):612 def build_box(a=10, b=20, c=30, rho=2.): 637 613 shape = Box(a, b, c, rho) 638 614 fn = lambda q: box_Iq(q, a, b, c)*rho**2 … … 640 616 return shape, fn, fn_xy 641 617 642 def check_box_lattice(a=10, b=20, c=30, rho=2.): 643 nx, dx = 3, a 644 ny, dy = 5, b 645 nz, dz = 5, c 646 dx, dy, dz = 2*dx, 2*dy, 2*dz 647 def center(*args): 648 sigma = 0.333 649 space = 2 650 return [(space*n+np.random.randn()*sigma)*x for n, x in args] 651 t0 = time.time() 652 shapes = [Box(a, b, c, rho, 653 #center=(ix*dx, iy*dy, iz*dz) 654 orientation=np.random.randn(3)*10, 655 center=center((ix, dx), (iy, dy), (iz, dz)) 656 ) 657 for ix in range(nx) 658 for iy in range(ny) 659 for iz in range(nz)] 660 shape = Composite(shapes) 661 print("generate points time", time.time()  t0) 662 fn = None 663 fn_xy = lambda qx, qy, view: box_Iqxy(qx, qy, a, b, c, view=view) 664 return shape, fn, fn_xy 665 666 667 def check_csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 618 def build_csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 668 619 shape = csbox(a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 669 620 fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) … … 672 623 return shape, fn, fn_xy 673 624 625 def build_cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 626 shuffle=0, rotate=0): 627 a, b, c = shape.dims 628 shapes = [copy(shape) 629 .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 630 (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 631 (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 632 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 633 for ix in range(nx) 634 for iy in range(ny) 635 for iz in range(nz)] 636 lattice = Composite(shapes) 637 return lattice 638 674 639 675 640 SHAPE_FUNCTIONS = OrderedDict([ 676 ("cylinder", check_cylinder), 677 ("sphere", check_sphere), 678 ("box", check_box), 679 ("csbox", check_csbox), 680 ("multicyl", check_cylinder_lattice), 681 ("multibox", check_box_lattice), 641 ("cylinder", build_cylinder), 642 ("sphere", build_sphere), 643 ("box", build_box), 644 ("csbox", build_csbox), 682 645 ]) 683 646 SHAPES = list(SHAPE_FUNCTIONS.keys()) … … 686 649 mesh=100, qmax=1.0, r_step=0.01, samples=5000): 687 650 rho_solvent = 0 688 qmin = qmax/100 0.651 qmin = qmax/100. 689 652 q = np.logspace(np.log10(qmin), np.log10(qmax), mesh) 690 653 r = shape.r_bins(q, r_step=r_step) 691 sampling_density = samples / shape.volume ()654 sampling_density = samples / shape.volume 692 655 rho, points = shape.sample(sampling_density) 693 656 t0 = time.time() … … 710 673 qy = np.linspace(0.0, qmax, mesh) 711 674 Qx, Qy = np.meshgrid(qx, qy) 712 sampling_density = samples / shape.volume ()713 #t0 = time.time()675 sampling_density = samples / shape.volume 676 t0 = time.time() 714 677 rho, points = shape.sample(sampling_density) 715 #print("sampletime", time.time()  t0)678 print("point generation time", time.time()  t0) 716 679 t0 = time.time() 717 680 Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) … … 734 697 formatter_class=argparse.ArgumentDefaultsHelpFormatter, 735 698 ) 736 parser.add_argument('d', 'dim', type=int, default=1, help='dimension 1 or 2') 737 parser.add_argument('m', 'mesh', type=int, default=100, help='number of mesh points') 738 parser.add_argument('s', 'samples', type=int, default=5000, help="number of sample points") 739 parser.add_argument('q', 'qmax', type=float, default=0.5, help='max q') 740 parser.add_argument('v', 'view', type=str, default='0,0,0', help='theta,phi,psi angles') 741 parser.add_argument('p', 'plot', action='store_true', help='plot points') 742 parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], help='oriented shape') 699 parser.add_argument('d', 'dim', type=int, default=1, 700 help='dimension 1 or 2') 701 parser.add_argument('m', 'mesh', type=int, default=100, 702 help='number of mesh points') 703 parser.add_argument('s', 'samples', type=int, default=5000, 704 help="number of sample points") 705 parser.add_argument('q', 'qmax', type=float, default=0.5, 706 help='max q') 707 parser.add_argument('v', 'view', type=str, default='0,0,0', 708 help='theta,phi,psi angles') 709 parser.add_argument('n', 'lattice', type=str, default='1,1,1', 710 help='lattice size') 711 parser.add_argument('z', 'spacing', type=str, default='2,2,2', 712 help='lattice spacing') 713 parser.add_argument('r', 'rotate', type=float, default=0., 714 help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise") 715 parser.add_argument('w', 'shuffle', type=float, default=0., 716 help="position relative to lattice, gaussian < 0.3, uniform otherwise") 717 parser.add_argument('p', 'plot', action='store_true', 718 help='plot points') 719 parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], 720 help='oriented shape') 743 721 parser.add_argument('pars', type=str, nargs='*', help='shape parameters') 744 722 opts = parser.parse_args() 745 723 pars = {key: float(value) for p in opts.pars for key, value in [p.split('=')]} 724 nx, ny, nz = [int(v) for v in opts.lattice.split(',')] 725 dx, dy, dz = [float(v) for v in opts.spacing.split(',')] 726 shuffle, rotate = opts.shuffle, opts.rotate 746 727 shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars) 728 if nx > 1 or ny > 1 or nz > 1: 729 shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 747 730 title = "%s(%s)" % (opts.shape, " ".join(opts.pars)) 748 731 if opts.dim == 1:
Note: See TracChangeset
for help on using the changeset viewer.