### Setup 2D Domain and Boundary Conditions

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/EllipticAsymmetric.ipynb

Creates a 2D computational domain using meshgrid and defines boundary conditions. This setup is for solving a 2D PDE problem.

```python
# Create the domain
aX0 = np.linspace(-1,1,101); aX1=aX0
gridScale = aX0[1]-aX0[0]
X0,X1 = np.meshgrid(aX0,aX1,indexing='ij')

# Set the boundary conditions 
domain = (X0**2+X1**2)<0.9**2
bc = 0.3*np.sin(2*np.pi*(X0+X1))
bc[domain] = np.nan
```

--------------------------------

### Setup Domain and Boundary Conditions for Vanishing Viscosity

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/EllipticAsymmetric.ipynb

Creates a 1D domain, calculates grid scale, and sets up boundary conditions for a vanishing viscosity problem. This is a standard setup for 1D PDE simulations.

```python
# Create the domain
X0 = np.linspace(-1,1,101)
gridScale = X0[1]-X0[0]

# Set the boundary conditions 
bc=np.full(X0.shape,np.nan)
bc[0]=0; bc[-1]=1
```

--------------------------------

### Build and Upload Python Wheel Package

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/conda.recipe/steps_to_upload.txt

Build a Python wheel distribution for the project and upload it to a package index using twine. This involves running the setup script and then uploading the contents of the 'dist' directory.

```bash
python conda.recipe/setup.py bdist_wheel
python -m twine upload dist/*
```

--------------------------------

### Integration Setup for Symplectic Schemes (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Algo/Dense.ipynb

Sets up initial conditions and the Hamiltonian for integrating a system using symplectic schemes. This example uses the Celestial Mechanics Hamiltonian and defines initial positions and momenta.

```python
Hamiltonian = H_Celestial
q,p = np.array([1,0]),np.array([0.3,0.7])
```

```python
Hamiltonian = H_Celestial
q,p = np.array([1.,0.]),np.array([0.3,0.7])
```

--------------------------------

### Setup for Comparison of W1 Distances (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/Prox_BeckmanOT.ipynb

Initializes parameters for comparing W1 distances, including the initial condition `x0` and density `ρ`. It also calls a placeholder function `compare_cpu_gpu`.

```python
import numpy as np
# Placeholder for damped_oscillating_test and compare_cpu_gpu for demonstration
def damped_oscillating_test(*args, **kwargs):
    # Dummy implementation
    class MockResult:
        def __init__(self):
            self.y = [np.random.rand(10)]
    class MockInputs:
        def __init__(self):
            self.ξ = np.random.rand(10)
            self.λ = 1.0
    return np.random.rand(), (MockResult(), None, MockInputs())
def compare_cpu_gpu(data):
    pass

x0=5; ρ = np.sqrt(3);
# Assuming Xs and dx are defined elsewhere
# λ,ξ = damped_oscillating_test(x0,ρ,Xs)
# compare_cpu_gpu((λ,ξ,dx))
```

--------------------------------

### Install Package from Anaconda Channel

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/conda.recipe/steps_to_upload.txt

Install the Adaptive Grid Discretizations package from the 'agd-lbr' Anaconda channel for testing purposes. The '--force' flag ensures reinstallation if the package already exists.

```bash
conda install --force agd -c agd-lbr
```

--------------------------------

### Acoustic Wave Equation Parameters Setup (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Repro/WaveAccuracy.ipynb

Sets up initial parameters for acoustic wave equation simulations, including wave speed, density, and domain characteristics. Defines constants for simulation time steps.

```python
ϵ_acoustic = 0.05; ϕ_acoustic = lambda X:ϕ_default(X,ϵ_acoustic)
ρ_flat = 1.
D_flat = Riemann.from_diagonal([3**2, 1]).rotate_by(π/8).m
Nt = 100 # Number of time steps 

```

--------------------------------

### Initialize Nucleation Site Data (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/RandomCells.ipynb

Initializes data structures for nucleation sites, including their number, position, start time, index, and color. This setup is crucial for simulating the initial conditions of a crystallization process.

```python
seeds_tmax = seeds_n**(-1/vdim) # Latest considered nucleation start time.  seeds_t = seeds_tmax*np.random.rand(seeds_n) # Nucleation start time
 seeds_f = np.arange(seeds_n) # Index used for the Voronoi diagram.  seeds_c = np.random.rand(seeds_n,3) # Color of each seed, uniform in [0,1]^3
 seeds_data = {
    'number':seeds_n,
    'position':seeds_p,
    'time':seeds_t,
    'flag':seeds_f, 
    'color':seeds_c,
}
```

--------------------------------

### Setup Simulation Parameters and Cost Function (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_NonDiv/BoatRoutingGeneric_Time.ipynb

Initializes the simulation by creating the domain, defining velocity sets, and calculating the cost function based on intrinsic metric and currents. It prepares all necessary parameters for the time-stepping simulation.

```python
# Setup the coordinate system and initial condition
X,dx,C_initial = MakeDomain(100 if large_instances else 50)

# Select the set of test velocities
velocity_norm = xp.linspace(1.5,0,12,endpoint=False)
velocity_nangle = 32 if large_instances else 16
velocity_angle = xp.linspace(0,2*np.pi,velocity_nangle,endpoint=False)
velocity = np.moveaxis(ad.array([[ (r*np.cos(θ),r*np.sin(θ)) 
    for r in velocity_norm] for θ in velocity_angle]),-1,0) 

# Quadratic and static model, introduced in 
μ,ω,M = np.ones_like(C_initial),BoatRouting.Currents(*X),BoatRouting.IntrinsicMetric(BoatRouting.Spherical,*X)
μ,ω,M = tuple(e.reshape(e.shape[:-2]+(1,1)+e.shape[-2:]) for e in (μ,ω,M)) # Broadcast ...
v = velocity.reshape(velocity.shape+(1,1)) # Broadcast ...
cost = μ + 0.5* lp.dot_VAV(v-ω,M,v-ω)

# Time steps
t_max = 2.5
dt = 0.5*CFL(velocity,dx)
t_range = xp.arange(0,t_max,dt)
```

--------------------------------

### Declare Seeds and Tips for Front Propagation and Backtracking - Python

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/Curvature.ipynb

Configures the 'seed' and 'tip' parameters within the Eikonal solver's input. It allows for both oriented seeds/tips (with position and orientation) and unoriented seeds/tips (position only), along with their initial values. This setup is crucial for defining the starting and ending points for wave propagation and geodesic calculations.

```python
hfmIn.update({
    # Declare the seeds, from which the front is propagated
    'seed':[-0.5,0.5,0],
    'seed_Unoriented':[0.5,0.5],
    
    # Optionally declare the boundary conditions at the seeds (default is 0)
    'seedValue':0.,
    'seedValue_Unoriented':0.3,
})

# Declare the tips, from which the geodesics are backtracked
hfmIn['tips'] = [[x,y,np.pi/4] 
                for x in Eikonal.CenteredLinspace(-1,1,6) 
                for y in Eikonal.CenteredLinspace(0,1,4)]
hfmIn['tips_Unoriented'] = hfmIn['tips'][:,0:2]
```

--------------------------------

### Setup Parameters for Newton Solver (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_NonDiv/OTBoundary1D.ipynb

Initializes parameters for the Newton solver, including defining the spatial grid 'X', calculating grid scale 'gridScale', setting the forcing function 'f', boundary conditions 'bc', and an initial guess 'guess'. This setup is used for problems with no solution or degenerate solutions.

```python
X = np.linspace(-1,1,101,endpoint=True)
gridScale = X[1]-X[0]
f = 1.+0.9*np.sin(2.*np.pi*X)
bc = np.array((-1.,1.))

guess = np.zeros(X.shape)
```

--------------------------------

### Initialize PDE Parameters and Grid

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/EllipticAsymmetric.ipynb

Sets up initial parameters for a PDE problem, including diffusion coefficients, source terms, and grid-related variables. This is a common setup step before defining the numerical scheme.

```python
# Choose the PDE parameters
omega = fd.as_field(np.array([0.]),X0.shape,depth=1); alpha = 0.; f = 0.
Diff = (1.+0.8*np.sin(2.*np.pi*X0)).reshape((1,1)+X0.shape)
```

--------------------------------

### Add Tips for Backtracking Geodesics

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/DistanceFromBoundary.ipynb

Adds a few 'tips' to the solver configuration, which serve as starting points for backtracking geodesics. This is useful for visualizing the paths computed by the solver.

```python
hfmIn.SetUniformTips([6,6])
```

--------------------------------

### Setup Domain and Boundary Conditions for Critical Viscosity

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/EllipticAsymmetric.ipynb

Sets up a 1D domain and boundary conditions for a case where the viscosity is at its critical stable value. This is a specific scenario to investigate the transition to upwind schemes.

```python
# Create the domain
X0 = np.linspace(-1,1,11)
gridScale = X0[1]-X0[0]

# Set the boundary conditions 
bc=np.full(X0.shape,np.nan)
bc[0]=0; bc[-1]=1
```

--------------------------------

### Setup Grid for Discretization (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/Prox_BeckmanOT.ipynb

Initializes a staggered grid for numerical computations. This function is used to define the spatial domain and resolution for optimal transport calculations.

```python
Xs,Xv,dx = stag_grids(100,[-2,2])
Xs = Xs[0]
```

--------------------------------

### Setup Grid and Ground Cost for Segmentation

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/Prox_MinCut.ipynb

Initializes grid coordinates (`Xϕ`, `Xm`), the differential operator `dx`, and a spatially varying ground cost function `g` using a cosine profile. This setup is typical for image segmentation tasks.

```python
Xϕ,Xm,dx = MinCut.stag_grids((201,),[[-1],[1]]) 
Xm=Xm[0]; Xϕ=Xϕ[0] # Only one coordinate in dimension one
g = np.cos(np.pi*Xϕ)
metric = Metrics.Isotropic(0.1,vdim=1) # Constant cost
```

--------------------------------

### Setup for Two-Dimensional Experiment (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/PorousMinimization.ipynb

Initializes parameters for a two-dimensional experiment, including grid dimensions, spatial coordinates, initial condition, and time parameters. It uses numpy for array manipulation.

```python
nX = 35; nT = 6
#nX = 50; nT = 10 # bit longish
aX,dx = np.linspace(0,1,nX,retstep=True,endpoint=False)
aX+=dx/2
X = np.array(np.meshgrid(aX,aX,indexing='ij'))

Ti = 2e-5
u0 = Barenblatt(Ti,X-0.5)
Tf = 2e-4; T = Tf-Ti
dt = T/nT
```

--------------------------------

### Anaconda Package Upload

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/conda.recipe/steps_to_upload.txt

Log in to Anaconda, build the project package, and upload it to the Anaconda repository. This process requires authentication and specifying the path to the built package.

```bash
anaconda login --username $ANACONDA_LOGIN --password $ANACONDA_PASSWORD
conda-build AdaptiveGridDiscretizations/
anaconda upload --force $PATH_TO_AGD.tar.bz2 --user agd-lbr
```

--------------------------------

### Problem Setup: Isotropic Metric on a Square

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_NonDiv/EikonalEulerian.ipynb

Defines the problem parameters for testing the numerical scheme, specifically an isotropic metric on a square domain with a single seed point. This setup is used to initialize the grid, metric, boundary conditions, and decompose the tensors.

```python
#Define the square [-1,1]^2, sampled on a cartesian grid
aX0 = np.linspace(-1,1,51); aX1 = aX0
gridScale=aX0[1]-aX0[0]
X0,X1 = np.meshgrid(aX0,aX1,indexing='ij')

# Define the domain, and the problem parameters
metric = fd.as_field(np.eye(2),X0.shape)
bc=np.full(X0.shape,np.nan)
bc[X0.shape[0]//2,X0.shape[1]//2] = 0

# Decompose the tensors dual to the Riemannian metric
coefs, offsets = Selling.Decomposition(lp.inverse(metric))
cofs[:,np.logical_not(np.isnan(bc))] = np.nan
```

--------------------------------

### Setup 1D Experiment Parameters

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/PorousMinimization.ipynb

Initializes the parameters for a one-dimensional experiment, including spatial and temporal discretization, initial time, initial condition (u0) derived from the Barenblatt solution, and time step size.

```python
nX = 200; nT = 20
#nX = 5; nT=1
aX,dx = np.linspace(0,1,nX,retstep=True,endpoint=False)
aX+=dx/2
X = aX[None]

Ti = 1e-4
u0 = Barenblatt(Ti,X-0.5)
Tf = 1e-3; T = Tf-Ti
dt = T/nT
```

--------------------------------

### Execute Dispersion Comparison and Plotting Setup

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Repro/WaveAccuracy.ipynb

Executes the dispersion comparison for 2D acoustic wave propagation with specified parameters and sets up the plotting function for visualizing results.

```python
%%time
def mk_needle(uθ,κ): return Riemann.needle(uθ,1,κ).m
result_2d = dispersions(ppw=[4,6,10],mk_aniso=mk_needle,name='Two dimensional')
```

--------------------------------

### Artistic Effect Setup (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/AnisotropicDiffusion.ipynb

Loads a synthetic image for demonstrating artistic effects using diffusion. The image is read and normalized, serving as the input for generating artistic patterns.

```python
image0 = imread("Notebooks_Div/TestData/CoherenceEnhancingDiffusion2D_TestImage.png")
image0 = xp.asarray(image0/255)
```

--------------------------------

### Example Usage of Rotation Parameters

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Repro/WaveAccuracy.ipynb

Demonstrates the usage of the `rotation_params` function with 2D and 3D unit vectors, verifying the computed rotation using linear parallelism functions.

```python
u = np.array([-1.,2])
θ = rotation_params(u)
assert np.allclose(lp.rotation(θ) @ np.array([0,1.]), u/np.linalg.norm(u))
u = np.array([1.,2,-3])
θ,ax = rotation_params(u)
assert np.allclose(lp.rotation(θ,ax) @ np.array([0,0,1.]),u/np.linalg.norm(u))
```

--------------------------------

### Initialize Grid and Initial Function

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/AnisotropicDiffusion.ipynb

This code initializes a grid using `xp.linspace` and creates a 2D Gaussian function `u0` as the initial state. This setup is used to demonstrate the effects of non-linear diffusion over time.

```python
aX,dx = xp.linspace(-1,1,retstep=True)
X = ad.array(np.meshgrid(aX,aX,indexing='ij')) 
u0 = 2*np.exp( - 4*(X**2).sum(axis=0))
```

--------------------------------

### Create and Activate Anaconda Environment

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/conda.recipe/steps_to_upload.txt

Set up a dedicated Conda environment for the project using a provided YAML file and activate it for subsequent operations. This ensures all dependencies are isolated.

```bash
conda env create --file AdaptiveGridDiscretizations/Miscellaneous/agd-hfm_dev.yaml
conda activate agd-hfm_dev
```

--------------------------------

### Simulate Scheme Evolution Over Time (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_NonDiv/BoatRouting_Time.ipynb

Numerically simulates the scheme's evolution over a range of time steps, starting from an initial condition. It uses the precomputed linear structure for efficiency and accumulates the results.

```python
%%time
C_solution = ad.array(list(accumulate(
    t_range[1:],
    initial=C_initial, 
    func=lambda C_t,t: Scheme(C_t,dt,params,scheme0)
)))

```

--------------------------------

### Import HFM Library and Utilities (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/Illusion.ipynb

Imports necessary modules from the HFM library, including the Eikonal solver and plotting utilities. It also configures the import path to allow importing from a parent directory, which is useful when the conda package is not installed. This setup is crucial for running the subsequent examples.

```Python
import sys; sys.path.insert(0,"..") # Allow import of agd from parent directory (useless if conda package installed)
#from Miscellaneous import TocTools; print(TocTools.displayTOC('Illusion','FMM'))
```

```Python
from agd import Eikonal
from agd.Plotting import savefig; #savefig.dirName = 'Figures/Illusion'
```

```Python
import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy
```

--------------------------------

### Problem Parameter Setup and Scheme Initialization (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Repro/WaveAccuracy.ipynb

Sets up problem parameters such as domain resolution (Nxs), points per wavelength (ppws), and time parameters. It then builds the necessary coefficients (M, C, S) and initializes various wave propagation Hamiltonians for different numerical schemes (e.g., Selling, Lebedev) and orders of accuracy.

```python
Nxs = [150] # Considered domain resolutions
ppws = [7,8,10,12,15,20] # ppw in flat domain. 
Nt=100

ϵ_elastic = 0.02
ϕ_elastic = lambda X:ϕ_default(X,ϵ_elastic)
M_flat = np.eye(2)
C_flat = Hooke.olivine[0].extract_xz().rotate_by(π/12).hooke; C_flat/=np.max(C_flat)
cfl_mult = 0.5

results_Nx = []
for Nx in Nxs:
    # Build the domain, and coefficients
    X,dx = how.make_domain(Nx,2)
    M,C,S,make_sol = make_test_e(M_flat,C_flat,X,dx,ϕ_elastic)
    # Now, the Lebedev scheme. Since point density is doubled, we take this into account. 
    # In three dimensions, point density would be quadrupled (->replace sqrt(2) with 2**(2/3) in next line)
    Nx_Leb = int(np.round(Nx/np.sqrt(2)))
    X_Leb,dx_Leb = how.make_domain(Nx_Leb,2)
    _,X_10,X_01,_ = ec.shifted_grids(X_Leb,dx_Leb)
    M_Leb,C_Leb,S_Leb,_ = make_test_e(M_flat,C_flat,X_Leb,dx_Leb,ϕ=ϕ_elastic)
    _,_,_,make_sol_10   = make_test_e(M_flat,C_flat,X_10, dx_Leb,ϕ=ϕ_elastic)
    _,_,_,make_sol_01   = make_test_e(M_flat,C_flat,X_01, dx_Leb,ϕ=ϕ_elastic)
    

    # Build the hamiltonians
    stag2 = ec.staggered(dx_Leb,2,'Periodic')
    stag4 = ec.staggered(dx_Leb,4,'Periodic')
    λ,e = Eikonal.VoronoiDecomposition(C); sel_C = λ,Hooke.moffset(e) 
    λ,e = ts.conv_decomp(*ts.deconv(C,σ=5,ϵ=5e-4,depth=2)); conv_C = λ,Hooke.moffset(e)
    λ,e = Eikonal.VoronoiDecomposition(C,smooth=True); smooth_C = λ,Hooke.moffset(e)
#    lin_D = ts.linear_decomp(D) # Non-positive for this anisotropy
    # Set a timestep
    WaveH = aw.ElasticHamiltonian_Sparse(M,C,dx,2,S) # TODO : sel_C
    dt = cfl_mult*WaveH.dt_max()
    T = Nt*dt
    WaveHs_e = [ # All Hamiltonians
        ("Selling2",WaveH),
        ("CorrSelling2",ec.SellingCorrelated2H_ext(M,sel_C, X,dx,2,S)),
        ("ConvSelling2",ec.SellingCorrelated2H_ext(M,conv_C,X,dx,2,S)),
        ("SmoothSelling2",ec.SellingCorrelated2H_ext(M,smooth_C,X,dx,2,S)),
        ("Lebedev2",ec.LebedevH2_ext(M_Leb,C_Leb,stag2,X_Leb,S_Leb)),
        ("Selling4",aw.ElasticHamiltonian_Sparse(M,C,dx,4,S)), # TODO : sel_C
        ("CorrSelling4",ec.SellingCorrelated2H_ext(M,sel_C, X,dx,4,S)),
        ("ConvSelling4",ec.SellingCorrelated2H_ext(M,conv_C,X,dx,4,S)),
        ("SmoothSelling4",ec.SellingCorrelated2H_ext(M,smooth_C,X,dx,4,S)),
        ("Lebedev4",ec.LebedevH2_ext(M_Leb,C_Leb,stag4,X_Leb,S_Leb)),
#		("Sellinv6",aw.ElasticHamiltonian_Sparse(M,C,dx,6,S)),
#		("CorrSelling6",ec.SellingCorrelated2H_ext(M,C,X,dx,6,S)),		
        ]
    results_ppw = []
    for ppw in ppws:
        print(f"{Nx=}, {ppw=}")
        planewaves = [planewave_with_spatial_frequency(Nx/ppw,planewave) for planewave in planewaves_default_e]
        # Exact colocated solution
        (q_exact,p_exact) = make_sol(planewaves)        
        sol_Col = (q_exact(  dt/2),p_exact(0),  
                   q_exact(T+dt/2),p_exact(T))
        # Exact solution on the Lebedev grids
        q_10,p_10 = make_sol_10(planewaves)
        q_01,p_01 = make_sol_01(planewaves)
        sol_Leb = np.array(((q_10(  dt/2),q_01(  dt/2)), (p_10(0),p_01(0)),
                            (q_10(T+dt/2),q_01(T+dt/2)), (p_10(T),p_01(T))))
		
        # Additional parameters
        results_scheme = []
        for name,WaveH in WaveHs_e:
            q0,p0,qf,pf = sol_Leb if name.startswith("Lebedev") else sol_Col # 
            q1,p1 = WaveH.Euler_p(q0,p0,dt,niter=Nt)
            results_ord = []
            for norm_p in (1,2,np.inf):
                err = norm(q1-qf,norm_p)/norm(q1,norm_p) # relative error
                results_ord.append(err)
            results_scheme.append(results_ord)
        results_ppw.append(results_scheme)
    results_Nx.append(results_ppw)

results_e = np.array(results_Nx) # Nx,ppw,scheme,ord
```

--------------------------------

### Build and Serve Project Documentation with pdoc

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Miscellaneous/Notes.md

Builds the project documentation using pdoc and places the output in the `../AdaptiveGridDiscretizations_help/docs` directory. It also provides a command to view the documentation interactively.

```console
pip install pdoc
pdoc -t Miscellaneous/ -o ../AdaptiveGridDiscretizations_help/docs agd
```

```console
pdoc -t Miscellaneous/ agd
```

--------------------------------

### Python: Initializing Input Arrays

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Algo/Reverse.ipynb

Initializes NumPy arrays `a0` and `b0` for use in the AD loop examples. `a0` contains powers of 2, and `b0` is its reciprocal. These arrays serve as the starting point for the arithmetico-geometric algorithm.

```python
a0=np.array([1.,2.,4.,8.])
b0=1/a0
```

--------------------------------

### Acoustic Wave Equation Test Case Setup (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Repro/WaveAccuracy.ipynb

Generates parameters and the exact solution for the acoustic wave equation in a deformed homogeneous domain. It takes flat parameters, domain information, and a deformation function as input, returning the deformed density, dissipation, and a function to generate the exact solution.

```python
def make_test_a(ρ_flat,D_flat,X,dx, ϕ=ϕ_default):
    """Parameters and exact solution of the acoustic wave equation in a deformed homogeneous domain"""
    from agd.ExportedCode.Notebooks_Div.HighOrderWaves import ExactSol_a,tq_a,tp_a,tρ_a,tD_a
    ϕ,dϕ,inv_dϕ,Jϕ,d2ϕ = how.differentiate(lambda x:ϕ(x),X)
    ρ = tρ_a(lambda x:ρ_flat,ϕ,Jϕ)
    D = tD_a(lambda x:D_flat,ϕ,inv_dϕ,Jϕ)
    def make_sol(planewaves=planewaves_default_a):
        """Exact homogeneous planewave solution, deformed by ϕ."""
        q_flat,p_flat = ExactSol_a(ρ_flat,D_flat,*zip(*planewaves))
        def q(t): return tq_a(lambda x:q_flat(t,x),ϕ)
        def p(t): return tp_a(lambda x:p_flat(t,x),ϕ,Jϕ)
        return q,p
    return ρ,D,make_sol
```

--------------------------------

### Import Libraries for Adaptive Grid Discretizations

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Algo/VoronoiVectors.ipynb

Imports necessary libraries for the adaptive grid discretization project. It includes setting up the path for local imports and then importing specific modules for metrics, selling, linear parallelism, and finite differences. This setup is crucial for running the subsequent examples.

```python
import sys; sys.path.insert(0,"..") # Allow import of agd from parent directory (useless if conda package installed)
#from Miscellaneous import TocTools; print(TocTools.displayTOC('VoronoiVectors','Algo'))
```

```python
from agd.Metrics import Riemann
from agd import Selling
from agd import LinearParallel as lp
from agd import FiniteDifferences as fd
```

```python
import numpy as np
from matplotlib import pyplot as plt
import itertools
```

--------------------------------

### Initialize Solver for Tubular Structure Segmentation (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/AsymmetricQuadratic.ipynb

Initializes the Eikonal solver with an 'AsymmetricQuadratic2' model for segmenting a tubular structure. It specifies seed and tip points, which are crucial for defining the start and end points of the path to be segmented. This setup is tailored for applications like vessel segmentation.

```python
hfmIn = Eikonal.dictIn({
    'model':'AsymmetricQuadratic2',
    'exportValues':1,
    'seed':[0.75,-0.5],
    'tip':[-0.75,-0.5]
})
```

--------------------------------

### Install AGD Library using pip

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/README.md

Installs the AGD library using pip, the standard package installer for Python. This is the recommended method for installation.

```console
pip install agd
```

--------------------------------

### Run Preliminary Tests

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/conda.recipe/steps_to_upload.txt

Execute test scripts to verify the functionality of versioning, table of contents generation, and markdown processing within the project. These scripts are located in the Miscellaneous directory.

```python
python Miscellaneous/TestVersion.py # --new_version=???
python Miscellaneous/TestTocs.py # --update 
python Miscellaneous/TestMarkdown.py
```

--------------------------------

### Setup for 2D Convex Body Example (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Algo/Meissner.ipynb

Initializes parameters for demonstrating 2D convex body calculations. It sets a random seed for reproducibility, defines a number of angular points (nθ), generates sorted, somewhat random angles, and creates unit vectors X representing points on a circle. It also sets a height function 'u' close to 1 to ensure non-empty facets.

```python
np.random.seed(42)
nθ=10

# Use strictly increasing but somewhat random angles
θ = np.sort(np.pi*np.random.rand(nθ))+np.linspace(0,np.pi,nθ,endpoint=False)
X = np.array([np.cos(θ),np.sin(θ)]) # Unit vectors

# Choose u close to 1, so that all facets are non-empty
u = np.ones_like(θ)+0.2*np.random.rand(nθ)
```

--------------------------------

### Setting Up Numerical Grid and Scheme Parameters (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_NonDiv/BoatRoutingGeneric_Time.ipynb

Initializes the computational grid, defines time parameters, and sets up the initial cost function `C_initial` and the numerical scheme `scheme0` for solving the HJB PDE. It accounts for grid size, time step calculation based on CFL condition, and initial cost distribution.

```python
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

# Coordinate system
nX = 100 if large_instances else 50
aX,dx = xp.linspace(-1,1,nX,retstep=True)
X = ad.array(np.meshgrid(aX,aX,indexing='ij'))

four_point = False

# Wait until all domain is reached even at smallest speed
t_max = 1./np.min(velocity_norm)
dt = 0.5*CFL(velocity,dx,four_point) 
dt = float(dt) # cupy compatibility
t_range = np.arange(0,t_max,dt)

C_initial = 300 * norm(X,axis=0) 

scheme0 = MakeScheme0(C_initial,dx,velocity,four_point)
```

--------------------------------

### Setup and Run Eikonal Solver (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/HighAccuracy.ipynb

Configures and runs the Eikonal solver using the 'agd' library for an isotropic model in the Poincare hyperbolic plane. It sets up the input parameters, defines the grid, calculates the cost function, and executes the solver.

```python
dimx=100
hfmIn = Eikonal.dictIn({
    'model':'Isotropic2',
    'seed': [0.,1.5],
    'exportValues':1.,
})
hfmIn.SetRect(sides=[[-0.5,0.5],[1,2]],dimx=dimx,sampleBoundary=True)
hfmIn.SetUniformTips((6,6))
    
X,Y = hfmIn.Grid()
hfmIn['cost'] = PoincareCost((X,Y))
```

```python
hfmOut = hfmIn.Run()
```

--------------------------------

### Example Usage of Automatic Differentiation Checks (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/WaveExamples.ipynb

Provides example calls to the `check_ad` function to demonstrate its usage with different parameters for Acoustic and Elastic wave types. These examples cover various configurations of dimensions, orders, boundary conditions, and iteration counts. Dependencies include the `check_ad` function itself.

```python
def some_check_ad():
    check_ad(1,10,'Acoustic')
    check_ad(2,7,'Acoustic',order_t=4,bc='Neumann')
    check_ad(2,6,'Acoustic',order_t=4,bc='Dirichlet',niter=10)
    check_ad(3,5,'Acoustic',order_x=4)
```

--------------------------------

### Setup Domain and Initial Wave Conditions (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/WaveExamples.ipynb

Initializes the simulation domain, kinetic energy matrix, elastic constants, and generates initial conditions for pressure and shear waves. It also sets up boundary conditions based on the computation backend (CPU or GPU).

```python
dom,X,dx = make_domain(3,vdim=2)
M = xp.ones((1,1)) # Kinetic energy. Could use anisotropic matrix field, in case of coordinate changes.
C = mica.extract_xz().rotate_by(π/6).hooke
q0p = explosion(X) # pressure wave trigger
q0s = torsion(X) # Shear wave trigger
p0 = np.zeros_like(q0p)
bc = 'Neumann' if xp is np else 'Dirichlet' # Neumann is only implemented on CPU
```

--------------------------------

### One-Dimensional Grid and Cost Function Setup (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/Prox_MinCut.ipynb

Initializes a one-dimensional grid, defines a ground cost function 'g', and sets up a variable metric for a segmentation experiment. This prepares the input data for demonstrating mincut segmentation in one dimension.

```python
Xϕ,Xm,dx = MinCut.stag_grids((101,),[[-1],[1]]) 
Xm=Xm[0]; Xϕ=Xϕ[0] # Only one coordinate in dimension one
g = (np.abs(Xϕ)<=0.1) - (np.abs(Xϕ)>=0.9).astype(float)
metric = Metrics.Isotropic( 0.1*(1.+0.5*np.cos(2*np.pi*Xm)) )
```

--------------------------------

### Install AGD Library using Conda (Deprecated)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/README.md

Installs the AGD library using Conda. This version is deprecated and does not include GPU codes.

```console
conda install agd -c agd-lbr
```

--------------------------------

### Initialize Seismic Fast Marching Solver (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/HighAccuracy.ipynb

Sets up the input parameters for a Fast Marching solver, specifically configured for a 'Seismic2' model. It defines the initial seed points and specifies that computed values should be exported.

```python
hfmIn = Eikonal.dictIn({
    'model':'Seismic2',
    'seeds':np.array([[2.,0.]]),
    'exportValues':1.,
})

# Set the coordinate system
hfmIn.SetRect([[1,3],[-1,1]],dimx=200)
X = hfmIn.Grid()
```

--------------------------------

### Initialize Constant Metric Grid and Run Fast Marching Solver (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/Seismic.ipynb

Sets up a constant metric grid for fast marching, defines grid dimensions, and runs the solver. The 'factoringRadius' set to -1 indicates factorization over the entire domain. The output includes solver completion time and default field values.

```python
n=101
hfmIn_Constant.SetRect(sides=[[-1,1],[-1,1]],dimx=n)
hfmIn_Constant.update({
    'factoringRadius':-1, # factorization is used over all the domain 
    'verbosity':1,
})
X = hfmIn_Constant.Grid()
```

```python
exactSolution = hfmIn_Constant['metric'].norm(X)
```

```python
hfmOut = hfmIn_Constant.Run()
```

--------------------------------

### Install HFM Library using Conda

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/Summary.ipynb

Installs the HamiltonFastMarching (HFM) library and its dependencies using the conda package manager. This is a straightforward method for setting up the library for use in Python environments.

```console
conda install -c agd-lbr hfm
```

--------------------------------

### Elastic Wave Equation Test Case Setup (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Repro/WaveAccuracy.ipynb

Generates parameters and the exact solution for the elastic wave equation in a deformed homogeneous domain. It takes flat material properties, domain information, and a deformation function, returning the deformed stiffness tensors and a function to generate the exact solution.

--------------------------------

### Install CuPy for CUDA 10.2

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Miscellaneous/Notes.md

Installs the CuPy library, a NumPy-compatible array library accelerated by CUDA, for a specific CUDA toolkit version (10.2 in this case). This command should be run after the `agd-hfm_gpu` environment is activated.

```console
pip install cupy-cuda102
```

--------------------------------

### Setup for Two-Dimensional Test Case (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/Prox_BeckmanOT.ipynb

Initializes parameters for a two-dimensional test case, including grid dimensions, spatial step `dx`, and initial probability distributions `μ` and `ν`. It defines the difference `ξ` and regularization parameter `λ`.

```python
import numpy as np
# Assuming stag_grids is defined
# Placeholder for stag_grids for demonstration
def stag_grids(shape, bounds):
    n = shape[0]
    dx = np.array([(bounds[1][i] - bounds[0][i]) / shape[j] for i, j in enumerate([0, 1])])
    X = [np.linspace(bounds[0][i], bounds[1][i], shape[i]) for i in [0, 1]]
    Xv = np.meshgrid(*X)
    return (X, Xv, dx)

n = 100
Xs, Xv, dx = stag_grids((n, n + 1), [[-1, -1], [1, 1]])
μ = np.zeros(Xs[0].shape); ν = μ.copy()
μ[n // 6, n // 3] = 1 / np.prod(dx); ν[3 * n // 4, 5 * n // 6] = 0.8 / np.prod(dx)

ξ = ν - μ; λ = 1
```

--------------------------------

### Get Reeds-Shepp Model Shape

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/Curvature.ipynb

Retrieves and displays the shape of the input configuration for the Reeds-Shepp model. This is a simple call to the 'shape' attribute.

```python
hfmIn.shape
```

--------------------------------

### Setup Homogeneous Domain and Initial Conditions (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/WaveExamples.ipynb

Initializes a homogeneous domain, sets physical parameters like density and metric, and defines initial wave states (q0, p0). It also calculates the time step based on the maximum wave speed and grid spacing.

```python
dom,X,dx = make_domain(3,vdim=2)
ρ = 1.
#D = np.eye(2)
D = Riemann.from_diagonal((1,0.5)).rotate_by(π/6).m
q0 = churp(X)
p0 = np.zeros_like(q0)

dt_max = dx # Speed = 1
T = 1
niter = (int)(3*T/dt_max); dt = T/niter
```

--------------------------------

### Initialize Fast Marching Input (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/HighAccuracy.ipynb

Sets up the input dictionary for the Fast Marching (FM) solver, specifying parameters like the model type, seed point, and grid dimensions. It configures the spatial domain and initializes uniform tips for the solver.

```Python
dimx=200
hfmIn = Eikonal.dictIn({
    'model':'Riemann2',
    'seed': [0.,0.],
    'exportValues':1.,
    'geodesicSolver':'Discrete'
})
hfmIn.SetRect(sides=[[-1,1],[-0.5,0.5]],dimx=dimx)
hfmIn.SetUniformTips((6,6))
```

--------------------------------

### Example Euler Angle Conversion (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/Curvature3.ipynb

Demonstrates the usage of the EulerAngles function by converting a scaled unit vector back to its Euler angles.

```python
EulerAngles(0.1*UnitVector(1.9,4.1))
```

--------------------------------

### Initialize and Run Acoustic Dispersion Simulation (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/ElasticComparisons.ipynb

This script initializes parameters and runs the acoustic dispersion simulation for a specific numerical order. It sets up the figure, defines parameters like ppw, kappa, theta, order_x, and kmax, calculates dx and ks, and then calls the frame_acoustic function to generate the plot. Finally, it saves the figure.

```python
fig = plt.figure(figsize=[6,6]); plt.axis('equal')
ppw = 4; κ=3; θ=π/6; order_x=2; kmax = 2.1 # parameters
dx = 2*π/(ppw*np.sqrt(κ)); ks = mk_ks(kmax); title = f"Acoustic dispersions, {κ=}, {order_x=}, {ppw=}"
frame_acoustic(κ,θ,dx,order_x,ks,title)
savefig(fig,to_filename(title))
```

--------------------------------

### Get Anisotropy Bound in Python

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/HighAccuracy.ipynb

Retrieves the upper bound of the anisotropy for the defined metric. This is a single numerical value representing a property of the metric.

```python
np.max(hfmIn['metric'].anisotropy_bound())
```

--------------------------------

### Setting Up Eikonal Solver Parameters (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/ClosedPaths.ipynb

Initializes the Eikonal solver with specific parameters for the 'Dubins2' model, including curvature radius, cost, seed, and tip points. It also sets the rectangular domain and the number of orientation angles for the computation.

```python
hfmIn = Eikonal.dictIn({
    'model':'Dubins2',
    'xi':0.3, # Minimal curvature radius
    'cost':1,
    'seed_Unoriented':[0.2,0.5],
    'tip_Unoriented':[1.75,0.5],
    'exportValues':True,
})
hfmIn.SetRect([[0,2],[0,1]],dimx=200)
hfmIn.nTheta = 96
```

--------------------------------

### Plot Solution (First Example)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/EllipticAsymmetric.ipynb

Sets the title for a plot and displays the solution obtained from solving the weak form of the PDE. This is used for visualizing the simulation results.

```python
plt.title("Solution to -(du')'=0")
plt.plot(X0,variational.solve_weakform())
```

--------------------------------

### Setup for Optimization (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/Monopolist.ipynb

Initializes grid parameters, performs Delaunay triangulation, identifies boundary triangles, defines the convex region Y, and sets up an initial guess function for the solution.

```python
nX=11
X,DX = square_X(nX-1)
T = Delaunay(X)
Tbd = boundary_triangulation(T,DX)
Y = np.array([(1,0),(1,1),(0,1),(0,0)]).T
def guess(X): return 0.1+0.25*((X[0]**2+X[1]**2)+X[0]+X[1])
```

--------------------------------

### Initialize and Run Numerical Scheme (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_NonDiv/BoatRoutingGeneric_Time.ipynb

This snippet demonstrates the initialization and execution of a numerical scheme for a given cost function and initial condition. It includes setting up parameters like time step, grid, and initial state, then iteratively applying the scheme over a time range. The CFL condition is used to determine a stable time step.

```python
four_point=False

# Quadratic parameters : unit weight for time, unit horizontal drift, euclidean geometry
model_params = (1.,xp.array((1.,0.)),xp.eye(2)) # μ, ω, M

# Velocities sample the unit ball
velocity_norm = xp.linspace(1,0,8,endpoint=False)
velocity_angle = xp.linspace(0,2*np.pi,16,endpoint=False)
velocity = np.moveaxis(ad.array([[ (r*np.cos(θ),r*np.sin(θ)) 
    for r in velocity_norm] for θ in velocity_angle]),-1,0) 
cost = quadratic_cost(model_params,velocity)

# Coordinate system
aX,dx = xp.linspace(-1,1,retstep=True)
X = ad.array(np.meshgrid(aX,aX,indexing='ij'))

# Initial condition
K_Lip = 5.
C_initial = K_Lip * norm(X,axis=0)
dt = CFL(velocity,dx)

# Time discretization
t_max = 1.5
dt = 0.5*CFL(velocity,dx,four_point)
t_range = xp.arange(0,t_max,dt)
```

--------------------------------

### Initialize Scheme Parameters and Initial Condition (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_NonDiv/BoatRouting_Time.ipynb

Sets up parameters for a numerical scheme, including model parameters, coordinate system, and initial conditions. It defines unit weight for time, horizontal drift, geometry, grid, and an initial concentration based on a Lipshitz constant.

```python
# Parameters : unit weight for time, unit horizontal drift, euclidean geometry
model_params = (1.,xp.array((1.,0.)),xp.eye(2))
params = scheme_params(model_params) 

# Coordinate system
aX,dx = xp.linspace(-1,1,retstep=True)
X = np.meshgrid(aX,aX,indexing='ij')

# Initial condition
K_Lip = 3.
C_initial = K_Lip * ad.Optimization.norm(X,axis=0)
dt = CFL(K_Lip,dx,params)
```

--------------------------------

### Visualize Geodesics (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_FMM/Seismic.ipynb

Plots the calculated geodesics on top of the domain. Since the metric is constant in this example, the geodesics appear as straight lines, illustrating the shortest paths.

```python
# The geodesics are straight lines, since the metric is constant
fig = plt.figure(figsize=[4,4]); plt.title('Minimal geodesics, Seismic metric (constant)'); plt.axis('equal'); 
for geo in hfmOut['geodesics']:  plt.plot(*geo)
```

--------------------------------

### Build and Serve Jupyter Notebook Slides

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Miscellaneous/Notes.md

Converts a Jupyter notebook to HTML slides and serves them locally. This command is useful for presenting code or project findings directly from a notebook.

```console
jupyter nbconvert InteractiveStencils.ipynb --to slides --post serve
```

--------------------------------

### Import Plotting Libraries (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_Div/AnisotropicDiffusion.ipynb

Imports Matplotlib for plotting and SciPy's datasets for sample data, enabling the visualization of results and the use of standard datasets for examples.

```python
import matplotlib
from matplotlib import pyplot as plt
import scipy.datasets
```

--------------------------------

### Initialize Random Data and Boundary Conditions (Python)

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_NonDiv/NonlinearMonotoneSecond2D.ipynb

Initializes random data for testing and sets up a mock Dirichlet boundary condition. This is a prerequisite for running the scheme comparisons.

```python
np.random.seed(42)
u_random = np.random.uniform(-1,1,guess.shape)

bc_unit = Domain.MockDirichlet(guess.shape,1,padding=0.)
```

--------------------------------

### Backtrack Geodesics from Specific Time

Source: https://github.com/mirebeau/adaptivegriddiscretizations/blob/main/Notebooks_NonDiv/BoatRoutingGeneric_Time.ipynb

Computes and visualizes geodesics starting from a specific time that minimizes travel cost. This accounts for potential time-dependent costs and delays in the simulation.

```python
best_step = np.argmin(interp(X,C_solution)(y0),axis=0) # Step of the simulation minimizing cost
t_delay = len(t_range) - best_step-1 # Account for reversed time
geodesics = odeint_array(C_velocity[::-1],y0,t_range[::-1],X, t_delay=t_delay)
for geo,delay in zip(geodesics.T,t_delay): plt.plot(*geo[:,int(delay):]) # Drop the first points according to delay
```