agd.ExportedCode.Notebooks_FMM.BoatRouting 

 1# Code automatically exported from notebook BoatRouting.ipynb in directory Notebooks_FMM
 2# Do not modify
 3import sys; sys.path.insert(0,"..") # Allow import of agd from parent directory (useless if conda package installed)
 4#from Miscellaneous import TocTools; print(TocTools.displayTOC('BoatRouting','FMM'))
 5
 6from ... import Eikonal
 7from ... import LinearParallel as lp
 8from ... import FiniteDifferences as fd
 9from agd.Metrics import Rander,Riemann
10from ... import AutomaticDifferentiation as ad
11from agd.Plotting import savefig,quiver; #savefig.dirName = 'Images/BoatRouting'
12
13import numpy as np; xp = np
14from copy import copy
15import matplotlib.pyplot as plt
16
17def route_min(z,params):
18    z,μ,ω,M = fd.common_field((z,)+params, depths=(1,0,1,2))
19    z_norm = np.sqrt(lp.dot_VAV(z,M,z))
20    μω_norm = np.sqrt( 2*μ +lp.dot_VAV(ω,M,ω) )
21    cost = z_norm*μω_norm - lp.dot_VAV(z,M,ω)
22    time = z_norm / μω_norm
23    fuel = cost/time
24    rvel = z/time - ω
25    return {
26        'cost':cost, # minimal cost for this travel
27        'time':time, # optimal travel time 
28        'fuel':fuel, # instantaneous fuel consumption
29        'rvel':rvel, # relative velocity, w.r.t current
30    }
31
32def metric(params):
33    μ,ω,M = fd.common_field(params,depths=(0,1,2))
34    return Rander( M*(2*μ + lp.dot_VAV(ω,M,ω)), -lp.dot_AV(M,ω))
35
36def Spherical(θ,ϕ): 
37    """Spherical embedding: θ is longitude, ϕ is latitude from equator toward pole"""
38    return (np.cos(θ)*np.cos(ϕ), np.sin(θ)*np.cos(ϕ), np.sin(ϕ))
39
40def IntrinsicMetric(Embedding,*X):
41    """Riemannian metric for a manifold embedded in Euclidean space"""
42    X_ad = ad.Dense.identity(constant=X,shape_free=(2,)) # First order dense AD variable
43    Embed_ad = ad.asarray(Embedding(*X_ad)) # Differentiate the embedding
44    Embed_grad = Embed_ad.gradient()
45    Embed_M = lp.dot_AA(Embed_grad,lp.transpose(Embed_grad)) # Riemannian metric
46    return Embed_M
47
48def bump(x,y): 
49    """Gaussian-like bump (not normalized)"""
50    return np.exp(-(x**2+y**2)/2)
51
52def Currents(θ,ϕ):
53    """Some arbitrary vector field (water currents)"""
54    bump0 = bump(θ+1,(ϕ+0.3)*2); ω0=(0,1) # intensity and direction of the currents
55    bump1 = 2*bump(2*(θ-0.7),ϕ-0.2); ω1=(1,-1)
56    bump0,ω0,bump1,ω1 = fd.common_field( (bump0,ω0,bump1,ω1), depths=(0,1,0,1))
57    return bump0*ω0+bump1*ω1
58
59def ArrivalTime(hfmIn,params):
60    hfmIn = copy(hfmIn)
61#    if hfmIn.xp is not np: hfmIn['solver']='AGSI' #TODO : why needed ?
62    hfmIn['metric'] = metric(params)
63    hfmIn['exportGeodesicFlow']=1
64    cache = Eikonal.Cache(needsflow=True)
65    hfmOut = hfmIn.Run(cache=cache)
66    
67    flow = hfmOut['flow']
68    no_flow = np.all(flow==0,axis=0)
69    flow[:,no_flow]=np.nan  # No flow at the seed point, avoid zero divide    
70    route = route_min(flow,params)
71    costVariation = route['time']
72    costVariation[no_flow] = 0
73    hfmIn['costVariation'] = np.expand_dims(costVariation,axis=-1)
74    
75    hfmOut2 = hfmIn.Run(cache=cache) # cache avoids some recomputations
76    time = hfmOut2['values'].gradient(0)
77    return time,hfmOut
def route_min(z, params): 
18def route_min(z,params):
19    z,μ,ω,M = fd.common_field((z,)+params, depths=(1,0,1,2))
20    z_norm = np.sqrt(lp.dot_VAV(z,M,z))
21    μω_norm = np.sqrt( 2*μ +lp.dot_VAV(ω,M,ω) )
22    cost = z_norm*μω_norm - lp.dot_VAV(z,M,ω)
23    time = z_norm / μω_norm
24    fuel = cost/time
25    rvel = z/time - ω
26    return {
27        'cost':cost, # minimal cost for this travel
28        'time':time, # optimal travel time 
29        'fuel':fuel, # instantaneous fuel consumption
30        'rvel':rvel, # relative velocity, w.r.t current
31    }
def metric(params): 
33def metric(params):
34    μ,ω,M = fd.common_field(params,depths=(0,1,2))
35    return Rander( M*(2*μ + lp.dot_VAV(ω,M,ω)), -lp.dot_AV(M,ω))
def Spherical(θ, φ): 
37def Spherical(θ,ϕ): 
38    """Spherical embedding: θ is longitude, ϕ is latitude from equator toward pole"""
39    return (np.cos(θ)*np.cos(ϕ), np.sin(θ)*np.cos(ϕ), np.sin(ϕ))

Spherical embedding: θ is longitude, ϕ is latitude from equator toward pole

def IntrinsicMetric(Embedding, *X): 
41def IntrinsicMetric(Embedding,*X):
42    """Riemannian metric for a manifold embedded in Euclidean space"""
43    X_ad = ad.Dense.identity(constant=X,shape_free=(2,)) # First order dense AD variable
44    Embed_ad = ad.asarray(Embedding(*X_ad)) # Differentiate the embedding
45    Embed_grad = Embed_ad.gradient()
46    Embed_M = lp.dot_AA(Embed_grad,lp.transpose(Embed_grad)) # Riemannian metric
47    return Embed_M

Riemannian metric for a manifold embedded in Euclidean space

def bump(x, y): 
49def bump(x,y): 
50    """Gaussian-like bump (not normalized)"""
51    return np.exp(-(x**2+y**2)/2)

Gaussian-like bump (not normalized)

def Currents(θ, φ): 
53def Currents(θ,ϕ):
54    """Some arbitrary vector field (water currents)"""
55    bump0 = bump(θ+1,(ϕ+0.3)*2); ω0=(0,1) # intensity and direction of the currents
56    bump1 = 2*bump(2*(θ-0.7),ϕ-0.2); ω1=(1,-1)
57    bump0,ω0,bump1,ω1 = fd.common_field( (bump0,ω0,bump1,ω1), depths=(0,1,0,1))
58    return bump0*ω0+bump1*ω1

Some arbitrary vector field (water currents)

def ArrivalTime(hfmIn, params): 
60def ArrivalTime(hfmIn,params):
61    hfmIn = copy(hfmIn)
62#    if hfmIn.xp is not np: hfmIn['solver']='AGSI' #TODO : why needed ?
63    hfmIn['metric'] = metric(params)
64    hfmIn['exportGeodesicFlow']=1
65    cache = Eikonal.Cache(needsflow=True)
66    hfmOut = hfmIn.Run(cache=cache)
67    
68    flow = hfmOut['flow']
69    no_flow = np.all(flow==0,axis=0)
70    flow[:,no_flow]=np.nan  # No flow at the seed point, avoid zero divide    
71    route = route_min(flow,params)
72    costVariation = route['time']
73    costVariation[no_flow] = 0
74    hfmIn['costVariation'] = np.expand_dims(costVariation,axis=-1)
75    
76    hfmOut2 = hfmIn.Run(cache=cache) # cache avoids some recomputations
77    time = hfmOut2['values'].gradient(0)
78    return time,hfmOut