matscipy.electrochemistry.steric_correction
Enforces minimum distances on coordinates within discrete distribtution.
Copyright 2020 IMTEK Simulation University of Freiburg
Authors:
Johannes Hoermann <johannes.hoermann@imtek-uni-freiburg.de>
Examples
Benchmark different scipy optimizers for the steric correction problem:
>>> # measures of box
>>> xsize = ysize = 5e-9 # nm, SI units
>>> zsize = 10e-9 # nm, SI units
>>>
>>> # get continuum distribution, z direction
>>> x = np.linspace(0, zsize, 2000)
>>> c = [0.1,0.1]
>>> z = [1,-1]
>>> u = 0.05
>>>
>>> phi = potential(x, c, z, u)
>>> C = concentration(x, c, z, u)
>>> rho = charge_density(x, c, z, u)
>>>
>>> # create distribution functions
>>> distributions = [interpolate.interp1d(x,c) for c in C]
>>>
>>> # sample discrete coordinate set
>>> box = np.array([xsize, ysize, zsize])m
>>> sample_size = 100
>>>
>>> samples = [ continuous2discrete(
>>> distribution=d, box=box, count=sample_size) for d in distributions ]
>>>
>>> # apply penalty for steric overlap
>>> x = np.vstack(samples)
>>>
>>> box = np.array([[0.,0.,0],box]) # needs lower corner
>>>
>>> n = x.shape[0]
>>> dim = x.shape[1]
>>>
>>> # benchmark methods
>>> mindsq, (p1,p2) = scipy_distance_based_closest_pair(x)
>>> pmin = np.min(x,axis=0)
>>> pmax = np.max(x,axis=0)
>>> mind = np.sqrt(mindsq)
>>> logger.info("Minimum pair-wise distance in sample: {}".format(mind))
>>> logger.info("First sample point in pair: ({:8.4e},{:8.4e},{:8.4e})".format(*p1))
>>> logger.info("Second sample point in pair ({:8.4e},{:8.4e},{:8.4e})".format(*p2))
>>> logger.info("Box lower boundary: ({:8.4e},{:8.4e},{:8.4e})".format(*box[0]))
>>> logger.info("Minimum coordinates in sample: ({:8.4e},{:8.4e},{:8.4e})".format(*pmin))
>>> logger.info("Maximum coordinates in sample: ({:8.4e},{:8.4e},{:8.4e})".format(*pmax))
>>> logger.info("Box upper boundary: ({:8.4e},{:8.4e},{:8.4e})".format(*box[1]))
>>>
>>> # stats: method, x, res, dt, mind, p1, p2 , pmin, pmax
>>> stats = [('initial',x,None,0,mind,p1,p2,pmin,pmax)]
>>>
>>> r = 4e-10 # 4 Angstrom steric radius
>>> logger.info("Steric radius: {:8.4e}".format(r))
>>>
>>> methods = [
>>> 'Powell',
>>> 'CG',
>>> 'BFGS',
>>> 'L-BFGS-B'
>>> ]
>>>
>>> for m in methods:
>>> try:
>>> logger.info("### {} ###".format(m))
>>> t0 = time.perf_counter()
>>> x1, res = apply_steric_correction(x,box=box,r=r,method=m)
>>> t1 = time.perf_counter()
>>> dt = t1 - t0
>>> logger.info("{} s runtime".format(dt))
>>>
>>> mindsq, (p1,p2) = scipy_distance_based_closest_pair(x1)
>>> mind = np.sqrt(mindsq)
>>> pmin = np.min(x1,axis=0)
>>> pmax = np.max(x1,axis=0)
>>>
>>> stats.append([m,x1,res,dt,mind,p1,p2,pmin,pmax])
>>>
>>> logger.info("Minimum pair-wise distance in final configuration: {:8.4e}".format(mind))
>>> logger.info("First sample point in pair: ({:8.4e},{:8.4e},{:8.4e})".format(*p1))
>>> logger.info("Second sample point in pair ({:8.4e},{:8.4e},{:8.4e})".format(*p2))
>>> logger.info("Box lower boundary: ({:8.4e},{:8.4e},{:8.4e})".format(*box[0]))
>>> logger.info("Minimum coordinates in sample: ({:8.4e},{:8.4e},{:8.4e})".format(*pmin))
>>> logger.info("Maximum coordinates in sample: ({:8.4e},{:8.4e},{:8.4e})".format(*pmax))
>>> logger.info("Box upper boundary: ({:8.4e},{:8.4e},{:8.4e})".format(*box[1]))
>>> except:
>>> logger.warn("{} failed.".format(m))
>>> continue
>>>
>>> stats_df = pd.DataFrame( [ {
>>> 'method': s[0],
>>> 'runtime': s[3],
>>> 'mind': s[4],
>>> **{'p1{:d}'.format(i): c for i,c in enumerate(s[5]) },
>>> **{'p2{:d}'.format(i): c for i,c in enumerate(s[6]) },
>>> **{'pmin{:d}'.format(i): c for i,c in enumerate(s[7]) },
>>> **{'pmax{:d}'.format(i): c for i,c in enumerate(s[8]) }
>>> } for s in stats] )
>>>
>>> print(stats_df.to_string(float_format='%8.6g'))
method runtime mind p10 p11 p12 p20 p21 p22 pmin0 pmin1 pmin2 pmax0 pmax1 pmax2
0 initial 0 1.15674e-10 2.02188e-09 4.87564e-10 5.21835e-09 2.03505e-09 3.72691e-10 5.22171e-09 1.17135e-12 1.49124e-10 6.34126e-12 4.98407e-09 4.99037e-09 9.86069e-09
1 Powell 75.2704 8.02318e-10 4.23954e-09 3.36242e-09 8.80092e-09 4.31183e-09 2.56345e-09 8.81278e-09 4.01789e-10 4.0081e-10 4.2045e-10 4.59284e-09 4.54413e-09 9.5924e-09
2 CG 27.0756 7.9992e-10 3.39218e-09 4.00079e-09 8.27255e-09 3.86337e-09 4.27807e-09 7.68863e-09 4.00018e-10 4.00146e-10 4.00565e-10 4.59941e-09 4.59989e-09 9.59931e-09
3 BFGS 19.0255 7.99527e-10 1.82802e-09 3.54397e-09 9.69736e-10 2.41411e-09 3.936e-09 1.34664e-09 4.00514e-10 4.01874e-10 4.0002e-10 4.59695e-09 4.59998e-09 9.58155e-09
4 L-BFGS-B 11.7869 7.99675e-10 4.34395e-09 3.94096e-09 1.28996e-09 4.44064e-09 3.15999e-09 1.14778e-09 4.12146e-10 4.01506e-10 4.03583e-10 4.6e-09 4.59898e-09 9.5982e-09
Functions
|
Enforce steric constraints on coordinate distribution within box. |
|
Constraint function. |
|
Constraint function. |
Find coordinate pair with minimum distance squared ||xi-xj||^2. |
|
|
Target function. |
|
Target function. |
|
Target function. |
Find coordinate pair with minimum distance ||xi-xj||. |
|
|
Find coordinate pair with minimum distance squared ||xi-xj||^2. |
Find coordinate pair with minimum distance ||xi-xj||. |
|
|
Target function. |
Classes
|
Lazy evaluation for log messages. |
- class matscipy.electrochemistry.steric_correction.DeferredMessage(msg, func, *args, **kwargs)
Bases:
objectLazy evaluation for log messages.
- __init__(msg, func, *args, **kwargs)
- matscipy.electrochemistry.steric_correction.brute_force_closest_pair(x)
Find coordinate pair with minimum distance squared ||xi-xj||^2.
- Parameters:
x ((N,dim) ndarray) – coordinates
- Returns:
float, (ndarray, ndarray)
- Return type:
minimum distance squared and coodinates pair
Examples
Compare the performance of closest pair algorithms:
>>> from matscipy.electrochemistry.steric_distribution import scipy_distance_based_target_function >>> from matscipy.electrochemistry.steric_distribution import numpy_only_target_function >>> from matscipy.electrochemistry.steric_distribution import brute_force_target_function >>> import itertools >>> import pandas as pd >>> import scipy.spatial.distance >>> import timeit >>> >>> funcs = [ >>> brute_force_closest_pair, >>> scipy_distance_based_closest_pair, >>> planar_closest_pair ] >>> func_names = ['brute','scipy','planar'] >>> stats = [] >>> N = 1000 >>> dim = 3 >>> for k in range(5): >>> x = np.random.rand(N,dim) >>> lambdas = [ (lambda x=x,f=f: f(x)) for f in funcs ] >>> rets = [ f() for f in lambdas ] >>> vals = [ v[0] for v in rets ] >>> coords = [ c for v in rets for p in v[1] for c in p ] >>> times = [ timeit.timeit(f,number=1) for f in lambdas ] >>> diffs = scipy.spatial.distance.pdist( >>> np.atleast_2d(vals).T,metric='euclidean') >>> stats.append((*vals,*diffs,*times,*coords)) >>> >>> func_name_tuples = list(itertools.combinations(func_names,2)) >>> diff_names = [ 'd_{:s}_{:s}'.format(f1,f2) for (f1,f2) in func_name_tuples ] >>> perf_names = [ 't_{:s}'.format(f) for f in func_names ] >>> coord_names = [ >>> 'p{:d}{:s}_{:s}'.format(i,a,f) for f in func_names for i in (1,2) for a in ('x','y','z') ] >>> float_fields = [*func_names,*diff_names,*perf_names,*coord_names] >>> dtypes = [ (field, 'f4') for field in float_fields ] >>> labeled_stats = np.array(stats,dtype=dtypes) >>> stats_df = pd.DataFrame(labeled_stats) >>> print(stats_df.T.to_string(float_format='%8.6g')) 0 1 2 3 4 brute 2.24089e-05 5.61002e-05 8.51047e-05 3.48424e-05 5.37235e-05 scipy 2.24089e-05 5.61002e-05 8.51047e-05 3.48424e-05 5.37235e-05 planar 2.24089e-05 5.61002e-05 8.51047e-05 3.48424e-05 5.37235e-05 d_brute_scipy 0 0 0 0 0 d_brute_planar 0 0 0 0 0 d_scipy_planar 0 0 0 0 0 t_brute 4.02697 3.85543 4.1414 3.90338 3.86993 t_scipy 0.00708364 0.00698962 0.00762594 0.00703242 0.00703579 t_planar 0.38302 0.39462 0.434342 0.407233 0.420773 p1x_brute 0.132014 0.331441 0.553405 0.534633 0.977582 p1y_brute 0.599688 0.186959 0.90897 0.575864 0.636278 p1z_brute 0.49631 0.993856 0.246418 0.853567 0.411793 p2x_brute 0.134631 0.333526 0.55322 0.534493 0.977561 p2y_brute 0.603598 0.179771 0.915063 0.576894 0.629313 p2z_brute 0.496833 0.994145 0.239493 0.859377 0.409509 p1x_scipy 0.132014 0.331441 0.553405 0.534633 0.977582 p1y_scipy 0.599688 0.186959 0.90897 0.575864 0.636278 p1z_scipy 0.49631 0.993856 0.246418 0.853567 0.411793 p2x_scipy 0.134631 0.333526 0.55322 0.534493 0.977561 p2y_scipy 0.603598 0.179771 0.915063 0.576894 0.629313 p2z_scipy 0.496833 0.994145 0.239493 0.859377 0.409509 p1x_planar 0.132014 0.331441 0.55322 0.534633 0.977561 p1y_planar 0.599688 0.186959 0.915063 0.575864 0.629313 p1z_planar 0.49631 0.993856 0.239493 0.853567 0.409509 p2x_planar 0.134631 0.333526 0.553405 0.534493 0.977582 p2y_planar 0.603598 0.179771 0.90897 0.576894 0.636278 p2z_planar 0.496833 0.994145 0.246418 0.859377 0.411793
- matscipy.electrochemistry.steric_correction.recursive_closest_pair(x, y)
Find coordinate pair with minimum distance squared ||xi-xj||^2.
Find coordinate pair with minimum distance squared ||xi-xj||^2 with one point from x and the other point from y
- Parameters:
x ((N,dim) ndarray) – coordinates
- Returns:
float, (ndarray, ndarray)
- Return type:
minimum distance squared and coordinate pair
- matscipy.electrochemistry.steric_correction.planar_closest_pair(x)
Find coordinate pair with minimum distance ||xi-xj||.
ATTENTION: this implementation tackles the planar problem!
- Parameters:
x ((N,dim) ndarray) – coordinates
- Returns:
float, (ndarray, ndarray)
- Return type:
minimum distance squared and coordinates pair
- matscipy.electrochemistry.steric_correction.scipy_distance_based_closest_pair(x)
Find coordinate pair with minimum distance ||xi-xj||.
- Parameters:
x ((N,dim) ndarray) – coordinates
- Returns:
float, (ndarray, ndarray)
- Return type:
minimum distance squared and coodinates pair
Examples
Handling condensed distance matrix indices:
>>> c = np.array([1, 2, 3, 4, 5, 6]) >>> print(c) [1 2 3 4 5 6]
>>> d = scipy.spatial.distance.squareform(c) >>> print(d) [[0 1 2 3] [1 0 4 5] [2 4 0 6] [3 5 6 0]]
>>> I,J = np.tril_indices(d.shape[0],-1) >>> print(I,J) [1 2 2 3 3 3] [0 0 1 0 1 2]
>>> I,J = np.triu_indices(d.shape[0],1) >>> print(I,J) [0 0 0 1 1 2] [1 2 3 2 3 3]
>>> print(d[I]) [1 2 4 3 5 6]
- matscipy.electrochemistry.steric_correction.brute_force_target_function(x, r=1.0, constraints=None)
Target function. Penalize dense packing for coordinates ||xi-xj||<ri+rj.
- Parameters:
x ((N,dim) ndarray) – particle coordinates
r (float or (N,) ndarray, optional (default=1.0)) – steric radii of particles
- Returns:
float
- Return type:
target function value
Examples
Compare performance of target functions:
>>> from matscipy.electrochemistry.steric_distribution import scipy_distance_based_target_function >>> from matscipy.electrochemistry.steric_distribution import numpy_only_target_function >>> from matscipy.electrochemistry.steric_distribution import brute_force_target_function >>> import itertools >>> import pandas as pd >>> import scipy.spatial.distance >>> import timeit >>> >>> funcs = [ >>> brute_force_target_function, >>> numpy_only_target_function, >>> scipy_distance_based_target_function ] >>> func_names = ['brute','numpy','scipy'] >>> >>> stats = [] >>> K = np.exp(np.log(10)*np.arange(-3,3)) >>> for k in K: >>> lambdas = [ (lambda x0=x0,k=k,f=f: f(x0*k)) for f in funcs ] >>> vals = [ f() for f in lambdas ] >>> times = [ timeit.timeit(f,number=1) for f in lambdas ] >>> diffs = scipy.spatial.distance.pdist(np.atleast_2d(vals).T,metric='euclidean') >>> stats.append((k,*vals,*diffs,*times)) >>> >>> func_name_tuples = list(itertools.combinations(func_names,2)) >>> diff_names = [ 'd_{:s}_{:s}'.format(f1,f2) for (f1,f2) in func_name_tuples ] >>> perf_names = [ 't_{:s}'.format(f) for f in func_names ] >>> fields = ['k',*func_names,*diff_names,*perf_names] >>> dtypes = [ (field, '>f4') for field in fields ] >>> labeled_stats = np.array(stats,dtype=dtypes) >>> stats_df = pd.DataFrame(labeled_stats) >>> print(stats_df.to_string(float_format='%8.6g')) k brute numpy scipy d_brute_numpy d_brute_scipy d_numpy_scipy t_brute t_numpy t_scipy 0 0.001 3.1984e+07 3.1984e+07 3.1984e+07 5.58794e-08 6.70552e-08 1.11759e-08 0.212432 0.168858 0.0734278 1 0.01 3.19829e+07 3.19829e+07 3.19829e+07 9.31323e-08 7.82311e-08 1.49012e-08 0.212263 0.16846 0.0791856 2 0.1 3.18763e+07 3.18763e+07 3.18763e+07 7.45058e-09 1.86265e-08 1.11759e-08 0.201706 0.164867 0.0711544 3 1 2.27418e+07 2.27418e+07 2.27418e+07 3.72529e-08 4.84288e-08 1.11759e-08 0.20762 0.166005 0.0724238 4 10 199751 199751 199751 1.16415e-10 2.91038e-11 8.73115e-11 0.202635 0.161932 0.0772684 5 100 252.548 252.548 252.548 3.28555e-11 0 3.28555e-11 0.202512 0.161217 0.0726705
- matscipy.electrochemistry.steric_correction.scipy_distance_based_target_function(x, r=1.0, constraints=None)
Target function. Penalize dense packing for coordinates ||xi-xj||<ri+rj.
- Parameters:
x ((N,dim) ndarray) – particle coordinates
r (float or (N,) ndarray, optional (default=1.0)) – steric radii of particles
- Returns:
float
- Return type:
target function value
- matscipy.electrochemistry.steric_correction.numpy_only_target_function(x, r=1.0, constraints=None)
Target function. Penalize dense packing for coordinates ||xi-xj||<ri+rj.
- Parameters:
x ((N,dim) ndarray) – particle coordinates
r (float or (N,) ndarray, optional (default=1.0)) – steric radii of particles
- Returns:
float
- Return type:
target function value
- matscipy.electrochemistry.steric_correction.neigh_list_based_target_function(x, r=1.0, constraints=None, Dij=None)
Target function. Penalize dense packing for coordinates ||xi-xj||<ri+rj.
- Parameters:
x ((N,dim) ndarray) – particle coordinates
r (float or (N,) ndarray, optional (default=1.0)) – steric radii of particles
constraints (callable, returns (float, (N,dim) ndarray)) – constraint function value and gradient
Dij ((N, N) ndarray, optional (default=None)) – pairwise minimum allowed distance matrix, overrides r
- Returns:
(float, (N,dim) ndarray) –
- f: float, target (or penalty) function, value evaluates to
sum_i sum_{j!=i} max(0,(r_i+r_j)^2-||xi-xj||^2)^2
without double-counting pairs, where r_i and r_j are the steric radii of coordinate points i and j
- df: (N,dim) ndarray of float, the gradient, evaluates to
4*si sj ((r_i’+r_j)^2-||xi-xj||^2)^2)*(xik’-xjk’)*(kdi’j-kdi’i)
for entry (i’,k’), where i subscribes the coordinate point and k the spatial dimension. kd is Kronecker delta, si is sum over i
- Return type:
target function value and gradient
- matscipy.electrochemistry.steric_correction.box_constraint(x, box=array([[0., 0., 0.], [1., 1., 1.]]), r=0.0)
Constraint function. Confine coordinates within box.
- Parameters:
x ((N,dim) ndarray) – coordinates
box ((2,dim) ndarray, optional (default: [[0.,0.,0.],[1.,1.,1.]])) – box corner coordinates
r (float or np.(N,) ndarray, optional (default=0)) – steric radii of particles
- Returns:
float
- Return type:
penalty, positive
- matscipy.electrochemistry.steric_correction.box_constraint_with_gradient(x, box=array([[0., 0., 0.], [1., 1., 1.]]), r=0.0)
Constraint function. Confine coordinates within box.
- Parameters:
x ((N,dim) ndarray) – coordinates
box ((2,dim) ndarray, optional (default: [[0.,0.,0.],[1.,1.,1.]])) – box corner coordinates
r (float or np.(N,) ndarray, optional (default=0)) – steric radii of particles
- Returns:
(float, (N,dim) ndarray of float) –
- g(x), positive:
lik = box[0,k] - xik + ri > 0 for xik out of lower boundaries rik = xik + ri - box[1,k] > 0 for xik out of upper boundaries
where box[0] and box[1] are lower and upper corner of bounding box and k marks spatial dimension.
g(x) = sum_i sum_k ( max(0,lik) + max(0,rik) )^2
- dg(x): gradient, entries accordingly evaluates to
dgdxik(x) = 2*( -max(0,lik) + max(0,rik) )
- Return type:
penalty g(x) and its gradient dg(x)
- matscipy.electrochemistry.steric_correction.apply_steric_correction(x, box=None, r=None, method='L-BFGS-B', options={'disp': True, 'eps': 1e-08, 'gtol': 1e-08, 'maxiter': 100}, target_function=<function neigh_list_based_target_function>, returns_gradient=True, closest_pair_function=<function scipy_distance_based_closest_pair>)
Enforce steric constraints on coordinate distribution within box.
- Parameters:
x ((N,dim) ndarray) – Particle coordinates.
box ((2,dim) ndarray, optional (default: None)) – Box corner coordinates.
r (float or (N,) ndarray, optional (default=None)) – Steric radius of particles. Can be specified particle-wise.
options (dict, optional) – Forwarded to scipy minimzer. https://docs.scipy.org/doc/scipy/reference/optimize.minimize-bfgs.html (default: {‘gtol’:1.e-5,’maxiter’:10,’disp’:True,’eps’:1.e-8})
method (str) – scipy.optimize.minimize method
target_function (func, optional) – One of the target functions within this submodule, or function of same signature. (default: neigh_list_based_target_function)
returns_gradient (bool, optional (default: True)) – If True, then ‘target_function’ is expected to return a tuple (f, df) with f the actual target function value and df its (N,dim) gradient. This flag must be set for ‘neigh_list_based_target_function’.
closest_pair_function (func, optional) – One of the closest pair functions within this submodule, or function of same signature. (default: scipy_distance_based_closest_pair)
- Returns:
float ((N,dim) ndarray) – Modified particle coordinates, meeting steric constraints.
scipy.optimize.OptimizeResult – Forwarded raw object as returned by scipy.optimize.minimize