-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
964c4ff
commit 33aada0
Showing
301 changed files
with
68,974 additions
and
1 deletion.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1 +1 @@ | ||
v0.1.0 | ||
v0.2.0 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,4 @@ | ||
# Sphinx build info version 1 | ||
# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done. | ||
config: 12abd65b70de0084620f26bf724ceadc | ||
tags: 645f666f9bcd5a90fca523b33c5a78b7 |
252 changes: 252 additions & 0 deletions
252
v0.2.0/_downloads/0cd804d4bb1b789ed7d110b2910079bf/jacobian.ipynb
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,252 @@ | ||
{ | ||
"cells": [ | ||
{ | ||
"cell_type": "code", | ||
"execution_count": 1, | ||
"id": "08d80a51", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-04-09T23:05:54.730065Z", | ||
"iopub.status.busy": "2024-04-09T23:05:54.729873Z", | ||
"iopub.status.idle": "2024-04-09T23:05:54.952910Z", | ||
"shell.execute_reply": "2024-04-09T23:05:54.952286Z" | ||
} | ||
}, | ||
"outputs": [], | ||
"source": [ | ||
"import numba\n", | ||
"import numpy as np\n", | ||
"from choclo.prism import gravity_u\n", | ||
"\n", | ||
"@numba.jit(nopython=True, parallel=True)\n", | ||
"def build_jacobian(coordinates, prisms):\n", | ||
" \"\"\"\n", | ||
" Build a sensitivity matrix for gravity_u of a prism\n", | ||
" \"\"\"\n", | ||
" # Unpack coordinates of the observation points\n", | ||
" easting, northing, upward = coordinates[:]\n", | ||
" # Initialize an empty 2d array for the sensitivity matrix\n", | ||
" n_coords = easting.size\n", | ||
" n_prisms = prisms.shape[0]\n", | ||
" jacobian = np.empty((n_coords, n_prisms), dtype=np.float64)\n", | ||
" # Compute the gravity_u field that each prism generate on every observation\n", | ||
" # point, considering that they have a unit density\n", | ||
" for i in numba.prange(len(easting)):\n", | ||
" for j in range(prisms.shape[0]):\n", | ||
" jacobian[i, j] = gravity_u(\n", | ||
" easting[i],\n", | ||
" northing[i],\n", | ||
" upward[i],\n", | ||
" prisms[j, 0],\n", | ||
" prisms[j, 1],\n", | ||
" prisms[j, 2],\n", | ||
" prisms[j, 3],\n", | ||
" prisms[j, 4],\n", | ||
" prisms[j, 5],\n", | ||
" 1.0,\n", | ||
" )\n", | ||
" return jacobian" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": 2, | ||
"id": "c1b2c50d", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-04-09T23:05:54.955476Z", | ||
"iopub.status.busy": "2024-04-09T23:05:54.955085Z", | ||
"iopub.status.idle": "2024-04-09T23:05:54.958839Z", | ||
"shell.execute_reply": "2024-04-09T23:05:54.958383Z" | ||
} | ||
}, | ||
"outputs": [], | ||
"source": [ | ||
"easting = np.linspace(-5.0, 5.0, 21)\n", | ||
"northing = np.linspace(-4.0, 4.0, 21)\n", | ||
"easting, northing = np.meshgrid(easting, northing)\n", | ||
"upward = 10 * np.ones_like(easting)\n", | ||
"\n", | ||
"coordinates = (easting.ravel(), northing.ravel(), upward.ravel())" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": 3, | ||
"id": "e0fc3135", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-04-09T23:05:54.960767Z", | ||
"iopub.status.busy": "2024-04-09T23:05:54.960574Z", | ||
"iopub.status.idle": "2024-04-09T23:05:54.963860Z", | ||
"shell.execute_reply": "2024-04-09T23:05:54.963379Z" | ||
} | ||
}, | ||
"outputs": [], | ||
"source": [ | ||
"prisms = np.array(\n", | ||
" [\n", | ||
" [-10.0, 0.0, -7.0, 0.0, -15.0, -10.0],\n", | ||
" [-10.0, 0.0, 0.0, 7.0, -25.0, -15.0],\n", | ||
" [0.0, 10.0, -7.0, 0.0, -20.0, -13.0],\n", | ||
" [0.0, 10.0, 0.0, 7.0, -12.0, -8.0],\n", | ||
" ]\n", | ||
")" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": 4, | ||
"id": "16350ae8", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-04-09T23:05:54.965756Z", | ||
"iopub.status.busy": "2024-04-09T23:05:54.965566Z", | ||
"iopub.status.idle": "2024-04-09T23:05:56.508938Z", | ||
"shell.execute_reply": "2024-04-09T23:05:56.508360Z" | ||
} | ||
}, | ||
"outputs": [ | ||
{ | ||
"data": { | ||
"text/plain": [ | ||
"array([[-4.49966911e-11, -4.76014375e-11, -3.80054986e-11,\n", | ||
" -2.83231448e-11],\n", | ||
" [-4.49659418e-11, -4.75828152e-11, -3.86898648e-11,\n", | ||
" -2.90433576e-11],\n", | ||
" [-4.48738920e-11, -4.75270202e-11, -3.93578743e-11,\n", | ||
" -2.97540033e-11],\n", | ||
" ...,\n", | ||
" [-3.19387365e-11, -4.58884222e-11, -4.10526880e-11,\n", | ||
" -4.48751701e-11],\n", | ||
" [-3.12924999e-11, -4.52460654e-11, -4.11114162e-11,\n", | ||
" -4.49880072e-11],\n", | ||
" [-3.06338007e-11, -4.45849449e-11, -4.11310225e-11,\n", | ||
" -4.50257165e-11]])" | ||
] | ||
}, | ||
"execution_count": 4, | ||
"metadata": {}, | ||
"output_type": "execute_result" | ||
} | ||
], | ||
"source": [ | ||
"jacobian = build_jacobian(coordinates, prisms)\n", | ||
"jacobian" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": 5, | ||
"id": "386d6302", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-04-09T23:05:56.511716Z", | ||
"iopub.status.busy": "2024-04-09T23:05:56.511507Z", | ||
"iopub.status.idle": "2024-04-09T23:05:56.514816Z", | ||
"shell.execute_reply": "2024-04-09T23:05:56.514379Z" | ||
} | ||
}, | ||
"outputs": [], | ||
"source": [ | ||
"# Define densities for the prisms\n", | ||
"densities = np.array([200.0, 300.0, -100.0, 400.0])\n", | ||
"\n", | ||
"# Compute result\n", | ||
"g_u = jacobian @ densities" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": 6, | ||
"id": "d99cae20", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-04-09T23:05:56.516813Z", | ||
"iopub.status.busy": "2024-04-09T23:05:56.516481Z", | ||
"iopub.status.idle": "2024-04-09T23:05:56.520556Z", | ||
"shell.execute_reply": "2024-04-09T23:05:56.520095Z" | ||
} | ||
}, | ||
"outputs": [], | ||
"source": [ | ||
"@numba.jit(nopython=True, parallel=True)\n", | ||
"def gravity_upward_parallel(coordinates, prisms, densities):\n", | ||
" \"\"\"\n", | ||
" Compute the upward component of the acceleration of a set of prisms\n", | ||
" \"\"\"\n", | ||
" # Unpack coordinates of the observation points\n", | ||
" easting, northing, upward = coordinates[:]\n", | ||
" # Initialize a result array full of zeros\n", | ||
" result = np.zeros_like(easting, dtype=np.float64)\n", | ||
" # Compute the upward component that every prism generate on each\n", | ||
" # observation point\n", | ||
" for i in numba.prange(len(easting)):\n", | ||
" for j in range(prisms.shape[0]):\n", | ||
" result[i] += gravity_u(\n", | ||
" easting[i],\n", | ||
" northing[i],\n", | ||
" upward[i],\n", | ||
" prisms[j, 0],\n", | ||
" prisms[j, 1],\n", | ||
" prisms[j, 2],\n", | ||
" prisms[j, 3],\n", | ||
" prisms[j, 4],\n", | ||
" prisms[j, 5],\n", | ||
" densities[j],\n", | ||
" )\n", | ||
" return result" | ||
] | ||
}, | ||
{ | ||
"cell_type": "code", | ||
"execution_count": 7, | ||
"id": "3fcca485", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-04-09T23:05:56.522616Z", | ||
"iopub.status.busy": "2024-04-09T23:05:56.522194Z", | ||
"iopub.status.idle": "2024-04-09T23:05:57.132985Z", | ||
"shell.execute_reply": "2024-04-09T23:05:57.132357Z" | ||
} | ||
}, | ||
"outputs": [ | ||
{ | ||
"data": { | ||
"text/plain": [ | ||
"True" | ||
] | ||
}, | ||
"execution_count": 7, | ||
"metadata": {}, | ||
"output_type": "execute_result" | ||
} | ||
], | ||
"source": [ | ||
"expected = gravity_upward_parallel(coordinates, prisms, densities)\n", | ||
"np.allclose(g_u, expected)" | ||
] | ||
} | ||
], | ||
"metadata": { | ||
"kernelspec": { | ||
"display_name": "Python 3 (ipykernel)", | ||
"language": "python", | ||
"name": "python3" | ||
}, | ||
"language_info": { | ||
"codemirror_mode": { | ||
"name": "ipython", | ||
"version": 3 | ||
}, | ||
"file_extension": ".py", | ||
"mimetype": "text/x-python", | ||
"name": "python", | ||
"nbconvert_exporter": "python", | ||
"pygments_lexer": "ipython3", | ||
"version": "3.12.2" | ||
} | ||
}, | ||
"nbformat": 4, | ||
"nbformat_minor": 5 | ||
} |
118 changes: 118 additions & 0 deletions
118
v0.2.0/_downloads/21b045eb65658e7f0dbd75b01bedc72a/jacobian.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,118 @@ | ||
#!/usr/bin/env python | ||
# coding: utf-8 | ||
|
||
# In[1]: | ||
|
||
|
||
import numba | ||
import numpy as np | ||
from choclo.prism import gravity_u | ||
|
||
@numba.jit(nopython=True, parallel=True) | ||
def build_jacobian(coordinates, prisms): | ||
""" | ||
Build a sensitivity matrix for gravity_u of a prism | ||
""" | ||
# Unpack coordinates of the observation points | ||
easting, northing, upward = coordinates[:] | ||
# Initialize an empty 2d array for the sensitivity matrix | ||
n_coords = easting.size | ||
n_prisms = prisms.shape[0] | ||
jacobian = np.empty((n_coords, n_prisms), dtype=np.float64) | ||
# Compute the gravity_u field that each prism generate on every observation | ||
# point, considering that they have a unit density | ||
for i in numba.prange(len(easting)): | ||
for j in range(prisms.shape[0]): | ||
jacobian[i, j] = gravity_u( | ||
easting[i], | ||
northing[i], | ||
upward[i], | ||
prisms[j, 0], | ||
prisms[j, 1], | ||
prisms[j, 2], | ||
prisms[j, 3], | ||
prisms[j, 4], | ||
prisms[j, 5], | ||
1.0, | ||
) | ||
return jacobian | ||
|
||
|
||
# In[2]: | ||
|
||
|
||
easting = np.linspace(-5.0, 5.0, 21) | ||
northing = np.linspace(-4.0, 4.0, 21) | ||
easting, northing = np.meshgrid(easting, northing) | ||
upward = 10 * np.ones_like(easting) | ||
|
||
coordinates = (easting.ravel(), northing.ravel(), upward.ravel()) | ||
|
||
|
||
# In[3]: | ||
|
||
|
||
prisms = np.array( | ||
[ | ||
[-10.0, 0.0, -7.0, 0.0, -15.0, -10.0], | ||
[-10.0, 0.0, 0.0, 7.0, -25.0, -15.0], | ||
[0.0, 10.0, -7.0, 0.0, -20.0, -13.0], | ||
[0.0, 10.0, 0.0, 7.0, -12.0, -8.0], | ||
] | ||
) | ||
|
||
|
||
# In[4]: | ||
|
||
|
||
jacobian = build_jacobian(coordinates, prisms) | ||
jacobian | ||
|
||
|
||
# In[5]: | ||
|
||
|
||
# Define densities for the prisms | ||
densities = np.array([200.0, 300.0, -100.0, 400.0]) | ||
|
||
# Compute result | ||
g_u = jacobian @ densities | ||
|
||
|
||
# In[6]: | ||
|
||
|
||
@numba.jit(nopython=True, parallel=True) | ||
def gravity_upward_parallel(coordinates, prisms, densities): | ||
""" | ||
Compute the upward component of the acceleration of a set of prisms | ||
""" | ||
# Unpack coordinates of the observation points | ||
easting, northing, upward = coordinates[:] | ||
# Initialize a result array full of zeros | ||
result = np.zeros_like(easting, dtype=np.float64) | ||
# Compute the upward component that every prism generate on each | ||
# observation point | ||
for i in numba.prange(len(easting)): | ||
for j in range(prisms.shape[0]): | ||
result[i] += gravity_u( | ||
easting[i], | ||
northing[i], | ||
upward[i], | ||
prisms[j, 0], | ||
prisms[j, 1], | ||
prisms[j, 2], | ||
prisms[j, 3], | ||
prisms[j, 4], | ||
prisms[j, 5], | ||
densities[j], | ||
) | ||
return result | ||
|
||
|
||
# In[7]: | ||
|
||
|
||
expected = gravity_upward_parallel(coordinates, prisms, densities) | ||
np.allclose(g_u, expected) | ||
|
Oops, something went wrong.