-
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
51a5762
commit 97cbb84
Showing
333 changed files
with
83,681 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.2.0 | ||
v0.3.0 |
252 changes: 252 additions & 0 deletions
252
v0.3.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": "ed2827e5", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-10-08T19:49:32.363199Z", | ||
"iopub.status.busy": "2024-10-08T19:49:32.363003Z", | ||
"iopub.status.idle": "2024-10-08T19:49:32.577328Z", | ||
"shell.execute_reply": "2024-10-08T19:49:32.576805Z" | ||
} | ||
}, | ||
"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": "ba4ab9d7", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-10-08T19:49:32.579385Z", | ||
"iopub.status.busy": "2024-10-08T19:49:32.579029Z", | ||
"iopub.status.idle": "2024-10-08T19:49:32.583036Z", | ||
"shell.execute_reply": "2024-10-08T19:49:32.582438Z" | ||
} | ||
}, | ||
"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": "1fce33a7", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-10-08T19:49:32.584751Z", | ||
"iopub.status.busy": "2024-10-08T19:49:32.584447Z", | ||
"iopub.status.idle": "2024-10-08T19:49:32.587446Z", | ||
"shell.execute_reply": "2024-10-08T19:49:32.586985Z" | ||
} | ||
}, | ||
"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": "65c544bc", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-10-08T19:49:32.589257Z", | ||
"iopub.status.busy": "2024-10-08T19:49:32.588813Z", | ||
"iopub.status.idle": "2024-10-08T19:49:34.246016Z", | ||
"shell.execute_reply": "2024-10-08T19:49:34.244918Z" | ||
} | ||
}, | ||
"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": "a3e1efe4", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-10-08T19:49:34.248257Z", | ||
"iopub.status.busy": "2024-10-08T19:49:34.248067Z", | ||
"iopub.status.idle": "2024-10-08T19:49:34.251461Z", | ||
"shell.execute_reply": "2024-10-08T19:49:34.250874Z" | ||
} | ||
}, | ||
"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": "cb6da38a", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-10-08T19:49:34.252971Z", | ||
"iopub.status.busy": "2024-10-08T19:49:34.252796Z", | ||
"iopub.status.idle": "2024-10-08T19:49:34.256968Z", | ||
"shell.execute_reply": "2024-10-08T19:49:34.256492Z" | ||
} | ||
}, | ||
"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": "f54e8c06", | ||
"metadata": { | ||
"execution": { | ||
"iopub.execute_input": "2024-10-08T19:49:34.258473Z", | ||
"iopub.status.busy": "2024-10-08T19:49:34.258300Z", | ||
"iopub.status.idle": "2024-10-08T19:49:34.910140Z", | ||
"shell.execute_reply": "2024-10-08T19:49:34.909611Z" | ||
} | ||
}, | ||
"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.6" | ||
} | ||
}, | ||
"nbformat": 4, | ||
"nbformat_minor": 5 | ||
} |
118 changes: 118 additions & 0 deletions
118
v0.3.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.