agd.ExportedCode.Notebooks_Algo.SeismicNorm

 1# Code automatically exported from notebook SeismicNorm.ipynb in directory Notebooks_Algo
 2# Do not modify
 3from agd.Metrics.Seismic import Hooke,Thomsen,TTI
 4from agd.Metrics import Riemann
 5from ... import Sphere
 6from ... import LinearParallel as lp
 7from ... import AutomaticDifferentiation as ad
 8from agd.Plotting import SetTitle3D
 9
10import numpy as np
11import copy
12from matplotlib import pyplot as plt
13
14def proj_orthorombic(c):
15    """Project onto the vector space of Hooke tensors corresponding to orthorombic materials in their frame of reference"""
16    to_orthorombic = (c[0,0],c[0,1],c[0,2],c[1,1],c[1,2],c[2,2],c[3,3],c[4,4],c[5,5]) # Seismologists start at 1 ...
17    return Hooke.from_orthorombic(*to_orthorombic).hooke
18
19def frame_score(c,proj,r):
20    """Score for wether c coincides with its projection in the frame defined by r"""
21    c = c.rotate(lp.transpose(r)) # Put in specified frame
22    c.hooke -= proj(c.hooke) # Substract projection 
23    return (c.to_Mandel()**2).sum(axis=(0,1)) # Return frobenius norm squared
24
25def proj_hooke2(c):
26    """Project onto the vector space of Hooke tensors with (2,1) block diagonal structure"""
27    z = np.zeros_like(c[0,0])
28    return ad.array([
29        [c[0,0],c[0,1],     z],
30        [c[1,0],c[1,1],     z],
31        [     z,     z,c[2,2]] ])
32
33def proj_tetragonal(c):
34    """Tetragonal Hooke tensors share the coefficients c11=c22, c13=c23, c44=c55"""
35    c11,c12,c13,c22,c23,c33,c44,c55,c66 = ( # Seismologists start at 1 ...
36        c[0,0],c[0,1],c[0,2],c[1,1],c[1,2],c[2,2],c[3,3],c[4,4],c[5,5])
37    α=(c11+c22)/2
38    γ=(c13+c23)/2
39    δ=(c44+c55)/2
40    return Hooke.from_orthorombic(α,c12,γ,α,γ,c33,δ,δ,c66).hooke
41
42def proj_hexagonal(c):
43    """Hexagonal Hooke tensors are tetragonal, with the additional property that c66=(c11-c12)/2"""
44    c11,c12,c13,c22,c23,c33,c44,c55,c66 = ( # Seismologists start at 1 ...
45        c[0,0],c[0,1],c[0,2],c[1,1],c[1,2],c[2,2],c[3,3],c[4,4],c[5,5])
46    α=(3*(c11+c22)+2*c12+4*c66)/8
47    β=(c11+c22+6*c12-4*c66)/8
48    γ=(c13+c23)/2
49    δ=(c44+c55)/2
50    return Hooke.from_orthorombic(α,β,γ,α,γ,c33,δ,δ,(α-β)/2).hooke
51
52def rotation_from_ball(x):
53    """
54    Produces a rotation matrix from an element of the unit ball
55    B_3 -> S_3 (lower half) -> SO_3, via quaternions
56    B_1 -> SO_2, via rotation of angle pi*x
57    """
58    if   len(x)==3: return Sphere.rotation3_from_ball3(x)[0] 
59    elif len(x)==1: return lp.rotation(np.pi*x[0])
60    else: raise ValueError("Unsupported dimension")
61
62def ball_samples(vdim,dens):
63    """
64    Produce samples of the unit ball of dimension vdim.
65    Approx c(vdim) * dens^vdim elements.
66    """
67    aB,h = np.linspace(-1,1,dens,retstep=True,endpoint=False)
68    B = np.array(np.meshgrid(*(aB,)*vdim,indexing='ij'))
69    B += h*np.random.rand(*B.shape) # Add random perturbations
70    B = B[:,np.linalg.norm(B,axis=0)<=1]
71    return B.reshape(vdim,-1)
72
73def newton_extremize(f,x,niter,max_cond=1e10):
74    """
75    Runs niter steps of Newton's method for extremizing f. 
76    (A.k.a, solve grad(f) = 0)
77    """
78    x_ad = ad.Dense2.identity(constant=x,shape_free=(len(x),))
79    for i in range(niter):
80        f_ad = f(x_ad)
81#        print(np.linalg.eigvalsh(f_ad.coef2[0]))
82        x_ad.value -= lp.solve_AV(f_ad.hessian(),f_ad.gradient())
83    return x_ad.value,f_ad.value
84
85def hooke_frame(c,proj,dens=8,niter=8):
86    """
87    Optimize the frame score, via a Newton method, with a large number of initial seeds.
88    Return the best rotation and associated score.
89    """
90    def f(x): return frame_score(c,proj,rotation_from_ball(x))
91    x = ball_samples({2:1,3:3}[c.vdim],dens) 
92    x,f_x = newton_extremize(f,x,niter)
93    argmin = np.nanargmin(f_x)
94    return lp.transpose(rotation_from_ball(x[:,argmin])),f_x[argmin]