agd.ExportedCode.Notebooks_Algo.Meissner

View Source

  1# Code automatically exported from notebook Meissner.ipynb in directory Notebooks_Algo
  2# Do not modify
  3from ... import LinearParallel as lp
  4from ... import AutomaticDifferentiation as ad
  5from ... import Plotting
  6norm_infinity = ad.Optimization.norm_infinity
  7
  8import numpy as np
  9from matplotlib import pyplot as plt
 10import scipy
 11import scipy.optimize as sciopt
 12
 13def LinearConstraint_AD(f):
 14    """
 15    Takes a linear constraint f>=0, 
 16    encoded as an ad.Sparse variable, 
 17    and turns it into a scipy compatible constraint.
 18    """
 19    return sciopt.LinearConstraint(f.tangent_operator(), lb = -f.value, ub=np.inf)
 20
 21def QuadraticObjective_AD(f):
 22    """
 23    Takes a quadratic objective function, 
 24    encoded as an ad.Sparse2 variable,
 25    and returns fun, jac and hessian methods.
 26    """
 27    val0 = f.value
 28    grad0 = f.to_first().to_dense().gradient()
 29    n = len(grad0)
 30    hess0 = f.hessian_operator(shape=(n,n))
 31    def fun(x):  return val0+np.dot(grad0,x)+np.dot(x,hess0*x)/2
 32    def grad(x): return grad0+hess0*x
 33    def hess(x): return hess0
 34    return {'fun':fun,'jac':grad,'hess':hess}
 35
 36def NonlinearObjective(f,fargs):
 37    """Returns methods computing the value, gradient and hessian 
 38    of a given objective function f"""
 39    def fun(x):
 40        return f(x,*fargs)
 41    def grad(x): 
 42        return f(ad.Sparse.identity(constant=x),*fargs).to_dense().gradient()
 43    def hess(x): 
 44        return f(ad.Sparse2.identity(constant=x),*fargs).hessian_operator()
 45    return {'fun':fun,'jac':grad,'hess':hess}
 46
 47def NonlinearConstraint(f,fargs):
 48    """
 49    Represents the constraint : np.bincount(indices,weights) >= 0,
 50    where (indices, weights) = f(x,*fargs)
 51    (Indices may be repeated, and the associated values must be summed.)
 52    """
 53    def fun(x):
 54        ind,wei = f(x,*fargs); 
 55        return np.bincount(ind,wei)
 56    def grad(x): 
 57        ind,wei = f(ad.Sparse.identity(constant=x),*fargs)
 58        triplets = (wei.coef.reshape(-1),(ind.repeat(wei.size_ad),wei.index.reshape(-1)))
 59        return scipy.sparse.coo_matrix(triplets).tocsr()
 60    def hess(x,v): # v is a set of weights, provided by the optimizer
 61        ind,wei = f(ad.Sparse2.identity(constant=x),*fargs)
 62        return np.sum(v[ind]*wei).hessian_operator()
 63    return sciopt.NonlinearConstraint(fun,0.,np.inf,jac=grad,hess=hess,keep_feasible=True)
 64
 65def geodesate(points,triangles,n,nosym=False,tol=1e-6):
 66    """
 67    Generate a point set on a sphere, by projecting a regularly refined triangulation.
 68    (This implementation is stupid and inefficient, but will be enough for here.)
 69    Input : 
 70    - points, triangles : a triangulation of the sphere.
 71    - nosym : keep only a single representative among opposite points
 72    Output : the points of a refined triangulation.
 73    """
 74    out = np.zeros((3,1))
 75    def norm2(p): return (p**2).sum(axis=0)
 76    for i,j,k in triangles:
 77        pi,pj,pk = points[i],points[j],points[k]
 78        for u in range(n+1):
 79            for v in range(n+1-u):
 80                p = (u*pi+v*pj+(n-u-v)*pk)/n
 81                p /= np.sqrt(norm2(p))
 82                p = p[:,None]
 83                dist2 = np.min( norm2(out-p) )
 84                if nosym: dist2 = min(dist2, np.min( norm2(out+p) ) )
 85                if dist2<tol: continue 
 86                out = np.concatenate((out,p),axis=1)
 87    if nosym: out *= np.where(out[2]>=0,1,-1)
 88    return out[:,1:]
 89
 90ico_points = np.array([[0., 0., -0.9510565162951536], [0., 0., 0.9510565162951536], [-0.85065080835204, 0., -0.42532540417601994], 
 91  [0.85065080835204, 0., 0.42532540417601994], [0.6881909602355868, -0.5, -0.42532540417601994], 
 92  [0.6881909602355868, 0.5, -0.42532540417601994], [-0.6881909602355868, -0.5, 0.42532540417601994], 
 93  [-0.6881909602355868, 0.5, 0.42532540417601994], [-0.2628655560595668, -0.8090169943749475, -0.42532540417601994], 
 94  [-0.2628655560595668, 0.8090169943749475, -0.42532540417601994], [0.2628655560595668, -0.8090169943749475, 0.42532540417601994], 
 95  [0.2628655560595668, 0.8090169943749475, 0.42532540417601994]])
 96ico_triangles = np.array([[2, 12, 8], [2, 8, 7], [2, 7, 11], [2, 11, 4], [2, 4, 12], [5, 9, 1], [6, 5, 1], [10, 6, 1], [3, 10, 1], [9, 3, 1], 
 97   [12, 10, 8], [8, 3, 7], [7, 9, 11], [11, 5, 4], [4, 6, 12], [5, 11, 9], [6, 4, 5], [10, 12, 6], [3, 8, 10], [9, 7, 3]])-1
 98
 99def sphere_sampling(n,nosym=False):
100    """A rather regular sapling of the sphere, obtained by refining an icosahedron mesh."""
101    return geodesate(ico_points,ico_triangles,n,nosym)
102
103def to_mathematica(Z): return str(Z.T.tolist()).replace("[","{").replace("]","}")

def norm_infinity(arr, *args, **kwargs): View Source

46def norm_infinity(arr,*args,**kwargs):
47	"""
48	L-Infinity norm (largest absolute value)
49	"""
50	return norm(arr,np.inf,*args,**kwargs)

L-Infinity norm (largest absolute value)

def LinearConstraint_AD(f): View Source

14def LinearConstraint_AD(f):
15    """
16    Takes a linear constraint f>=0, 
17    encoded as an ad.Sparse variable, 
18    and turns it into a scipy compatible constraint.
19    """
20    return sciopt.LinearConstraint(f.tangent_operator(), lb = -f.value, ub=np.inf)

Takes a linear constraint f>=0, encoded as an ad.Sparse variable, and turns it into a scipy compatible constraint.

def QuadraticObjective_AD(f): View Source

22def QuadraticObjective_AD(f):
23    """
24    Takes a quadratic objective function, 
25    encoded as an ad.Sparse2 variable,
26    and returns fun, jac and hessian methods.
27    """
28    val0 = f.value
29    grad0 = f.to_first().to_dense().gradient()
30    n = len(grad0)
31    hess0 = f.hessian_operator(shape=(n,n))
32    def fun(x):  return val0+np.dot(grad0,x)+np.dot(x,hess0*x)/2
33    def grad(x): return grad0+hess0*x
34    def hess(x): return hess0
35    return {'fun':fun,'jac':grad,'hess':hess}

Takes a quadratic objective function, encoded as an ad.Sparse2 variable, and returns fun, jac and hessian methods.

def NonlinearObjective(f, fargs): View Source

37def NonlinearObjective(f,fargs):
38    """Returns methods computing the value, gradient and hessian 
39    of a given objective function f"""
40    def fun(x):
41        return f(x,*fargs)
42    def grad(x): 
43        return f(ad.Sparse.identity(constant=x),*fargs).to_dense().gradient()
44    def hess(x): 
45        return f(ad.Sparse2.identity(constant=x),*fargs).hessian_operator()
46    return {'fun':fun,'jac':grad,'hess':hess}

Returns methods computing the value, gradient and hessian of a given objective function f

def NonlinearConstraint(f, fargs): View Source

48def NonlinearConstraint(f,fargs):
49    """
50    Represents the constraint : np.bincount(indices,weights) >= 0,
51    where (indices, weights) = f(x,*fargs)
52    (Indices may be repeated, and the associated values must be summed.)
53    """
54    def fun(x):
55        ind,wei = f(x,*fargs); 
56        return np.bincount(ind,wei)
57    def grad(x): 
58        ind,wei = f(ad.Sparse.identity(constant=x),*fargs)
59        triplets = (wei.coef.reshape(-1),(ind.repeat(wei.size_ad),wei.index.reshape(-1)))
60        return scipy.sparse.coo_matrix(triplets).tocsr()
61    def hess(x,v): # v is a set of weights, provided by the optimizer
62        ind,wei = f(ad.Sparse2.identity(constant=x),*fargs)
63        return np.sum(v[ind]*wei).hessian_operator()
64    return sciopt.NonlinearConstraint(fun,0.,np.inf,jac=grad,hess=hess,keep_feasible=True)

Represents the constraint : np.bincount(indices,weights) >= 0, where (indices, weights) = f(x,*fargs) (Indices may be repeated, and the associated values must be summed.)

def geodesate(points, triangles, n, nosym=False, tol=1e-06): View Source

66def geodesate(points,triangles,n,nosym=False,tol=1e-6):
67    """
68    Generate a point set on a sphere, by projecting a regularly refined triangulation.
69    (This implementation is stupid and inefficient, but will be enough for here.)
70    Input : 
71    - points, triangles : a triangulation of the sphere.
72    - nosym : keep only a single representative among opposite points
73    Output : the points of a refined triangulation.
74    """
75    out = np.zeros((3,1))
76    def norm2(p): return (p**2).sum(axis=0)
77    for i,j,k in triangles:
78        pi,pj,pk = points[i],points[j],points[k]
79        for u in range(n+1):
80            for v in range(n+1-u):
81                p = (u*pi+v*pj+(n-u-v)*pk)/n
82                p /= np.sqrt(norm2(p))
83                p = p[:,None]
84                dist2 = np.min( norm2(out-p) )
85                if nosym: dist2 = min(dist2, np.min( norm2(out+p) ) )
86                if dist2<tol: continue 
87                out = np.concatenate((out,p),axis=1)
88    if nosym: out *= np.where(out[2]>=0,1,-1)
89    return out[:,1:]

Generate a point set on a sphere, by projecting a regularly refined triangulation. (This implementation is stupid and inefficient, but will be enough for here.) Input :

points, triangles : a triangulation of the sphere.
nosym : keep only a single representative among opposite points Output : the points of a refined triangulation.

ico_points = array([[ 0. , 0. , -0.95105652], [ 0. , 0. , 0.95105652], [-0.85065081, 0. , -0.4253254 ], [ 0.85065081, 0. , 0.4253254 ], [ 0.68819096, -0.5 , -0.4253254 ], [ 0.68819096, 0.5 , -0.4253254 ], [-0.68819096, -0.5 , 0.4253254 ], [-0.68819096, 0.5 , 0.4253254 ], [-0.26286556, -0.80901699, -0.4253254 ], [-0.26286556, 0.80901699, -0.4253254 ], [ 0.26286556, -0.80901699, 0.4253254 ], [ 0.26286556, 0.80901699, 0.4253254 ]])

ico_triangles = array([[ 1, 11, 7], [ 1, 7, 6], [ 1, 6, 10], [ 1, 10, 3], [ 1, 3, 11], [ 4, 8, 0], [ 5, 4, 0], [ 9, 5, 0], [ 2, 9, 0], [ 8, 2, 0], [11, 9, 7], [ 7, 2, 6], [ 6, 8, 10], [10, 4, 3], [ 3, 5, 11], [ 4, 10, 8], [ 5, 3, 4], [ 9, 11, 5], [ 2, 7, 9], [ 8, 6, 2]])

def sphere_sampling(n, nosym=False): View Source

100def sphere_sampling(n,nosym=False):
101    """A rather regular sapling of the sphere, obtained by refining an icosahedron mesh."""
102    return geodesate(ico_points,ico_triangles,n,nosym)

A rather regular sapling of the sphere, obtained by refining an icosahedron mesh.

def to_mathematica(Z): View Source

104def to_mathematica(Z): return str(Z.T.tolist()).replace("[","{").replace("]","}")