optiphy/ana/sample.py

0001 #!/usr/bin/env python
0002 """
0003
0004 https://stackoverflow.com/questions/33976911/generate-a-random-sample-of-points-distributed-on-the-surface-of-a-unit-sphere
0005
0006 """
0007 import numpy as np
0008
0009
0010 def xy_grid_coordinates(nx=11, ny=11, sx=100., sy=100.):
0011     """
0012     :param nx:
0013     :param ny:
0014     :param sx:
0015     :param sy:
0016     :return xyz: (nx*ny,3) array of XY grid coordinates
0017     """
0018     x = np.linspace(-sx,sx,nx)
0019     y = np.linspace(-sy,sy,ny)
0020     xx, yy = np.meshgrid(x, y)
0021     zz = np.zeros( (nx, ny) )
0022     xyz = np.dstack( (xx,yy,zz) ).reshape(-1,3)
0023     return xyz
0024
0025
0026
0027 def sample_linear(n, lo=-1., hi=1.):
0028     """
0029     :param n: number of values
0030     :param lo: min value
0031     :param hi: max value
0032     :return ll: array of shape (n,) with random sample values uniformly between lo and hi
0033     """
0034     mi = (lo + hi)/2.
0035     ex = (hi - lo)/2.
0036     uu = 2.*np.random.rand(n) - 1.   # -1->1
0037     ll = mi + uu*ex
0038     return ll
0039
0040 def sample_linspace(n, lo=-1., hi=1. ):
0041     return np.linspace(lo, hi, n )
0042
0043
0044 def sample_disc(n, dtype=np.float64):
0045     """
0046     https://stackoverflow.com/questions/5837572/generate-a-random-point-within-a-circle-uniformly
0047
0048     NB the sqrt(u) for the radius avoids clumping in the middle
0049
0050     Circumference of circle 2.pi.r, area pi.r^2
0051
0052
0053     PDX(x) = 2x       Density for uniform radius sampling
0054
0055                  + 2
0056                 /
0057                /
0058               /
0059              /
0060             /    + 1
0061            /
0062           /
0063          /
0064         +--------+
0065         0        1
0066
0067
0068     CDF(x) = Integral 2x dx = x^2
0069     ICDF   sqrt(x)
0070
0071
0072     """
0073     a = np.zeros( (n, 3), dtype=dtype )
0074     r = np.sqrt(np.random.rand(n))
0075     t = 2.*np.pi*np.random.rand(n)
0076     a[:,0] = r*np.cos(t)
0077     a[:,1] = r*np.sin(t)
0078     return a
0079
0080
0081
0082
0083
0084
0085
0086 def sample_trig(n):
0087     """
0088     :param n: number of points
0089     """
0090     theta = 2*np.pi*np.random.rand(n)
0091     phi = np.arccos(2*np.random.rand(n)-1)
0092     x = np.cos(theta) * np.sin(phi)
0093     y = np.sin(theta) * np.sin(phi)
0094     z = np.cos(phi)
0095     return np.array([x,y,z])
0096
0097 def sample_normals(npoints):
0098     vec = np.random.randn(3, npoints)
0099     vec /= np.linalg.norm(vec, axis=0)
0100     return vec
0101
0102 def sample_reject(npoints):
0103     vec = np.zeros((3,npoints))
0104     abc = 2*np.random.rand(3,npoints)-1
0105     norms = np.linalg.norm(abc,axis=0)
0106     mymask = norms<=1
0107     abc = abc[:,mymask]/norms[mymask]
0108     k = abc.shape[1]
0109     vec[:,0:k] = abc
0110     while k<npoints:
0111        abc = 2*np.random.rand(3)-1
0112        norm = np.linalg.norm(abc)
0113        if 1e-5 <= norm <= 1:
0114            vec[:,k] = abc/norm
0115            k = k+1
0116     return vec
0117
0118
0119 if __name__ == '__main__':
0120
0121     n = 10
0122     a = sample_trig(n)
0123     b = sample_normals(n)
0124     c = sample_reject(n)
0125
0126