Skip to content

Commit

Permalink
fix docs, optimise forward model
Browse files Browse the repository at this point in the history
  • Loading branch information
@FredDeCeuster
FredDeCeuster committed Jul 27, 2023
1 parent a36d7ab commit 1ba337e
Show file tree
Hide file tree
Showing 4 changed files with 156 additions and 56 deletions.
9 changes: 5 additions & 4 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,13 @@ _Probabilistic 3D Reconstruction of Spectral Line Observations._

## About

*p3droslo* is a library for **probabilistic 3D reconstruction** of astronomical **spectral line observations**.
*p3droslo* is a python package that allows you to create **probabilistic 3D reconstructions** of astronomical **spectral line observations**.

Observations of [spectral lines](https://en.wikipedia.org/wiki/Spectral_line) are indespensible in astronomy, since they encode a wealth of information about the physical and chemical conditions of the medium in which they were formed.
Observations of [spectral lines](https://en.wikipedia.org/wiki/Spectral_line) are indespensible in astronomy, since they encode a wealth of information about the physical and chemical conditions of the medium from which they originate.
For instance, their narrow extent in frequency space make them very sensitive to Doppler shifts, such that their shape encodes the motion of the medium along the line of sight.
As a result, given a good model for the line formation process and a good inversion method, these physical and chemical properties can be retrieved from observations.
As a result, given a good model for the line formation process and an inversion method, these physical and chemical properties can be retrieved from observations.
Currently, we mainly focus on retrieving the distributions of the abundance of the chemical species producing the line, the velocity field, and its kinetic temperature.
However, also other parameters can be retrieved.

More information about the model for [spectral line formation](https://p3droslo.readthedocs.io/en/latest/background/spectral_line_formation.html) and the [probabilistic reconstruction methods](https://p3droslo.readthedocs.io/en/latest/background/probabilistic_reconstruction.html) can be found in the [background](https://p3droslo.readthedocs.io/en/latest/background/index.html) pages.

Expand Down Expand Up @@ -57,7 +58,7 @@ We are open to contributions to p3droslo. More information can be found [here](h

## Collaborating

We are always interested in collaborating! If you have a project, an idea, or data that might benefit from probabilistic 3D reconstruction or any other method you can find in this repository, feel free to contact [me](https://freddeceuster.github.io).
We are always interested in collaborating! If you have a project or data that might benefit from probabilistic 3D reconstruction or any other method you can find in this repository, feel free to contact [me](https://freddeceuster.github.io).


## Acknowledgements
Expand Down
98 changes: 98 additions & 0 deletions docs/src/examples/R_Aquilae.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1406,6 +1406,104 @@
"loss.plot()"
]
},
{
"cell_type": "code",
"execution_count": 125,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor([0., 1., 2., 3., 4.], requires_grad=True)"
]
},
"execution_count": 125,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = torch.arange(5.0, requires_grad=True)\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 126,
"metadata": {},
"outputs": [
{
"ename": "RuntimeError",
"evalue": "grad can be implicitly created only for scalar outputs",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mRuntimeError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m/STER/frederikd/p3droslo/docs/src/examples/R_Aquilae.ipynb Cell 50\u001b[0m in \u001b[0;36m<cell line: 12>\u001b[0;34m()\u001b[0m\n\u001b[1;32m <a href='vscode-notebook-cell://ssh-remote%2Benif/STER/frederikd/p3droslo/docs/src/examples/R_Aquilae.ipynb#Y134sdnNjb2RlLXJlbW90ZQ%3D%3D?line=9'>10</a>\u001b[0m optimizer\u001b[39m.\u001b[39mzero_grad()\n\u001b[1;32m <a href='vscode-notebook-cell://ssh-remote%2Benif/STER/frederikd/p3droslo/docs/src/examples/R_Aquilae.ipynb#Y134sdnNjb2RlLXJlbW90ZQ%3D%3D?line=10'>11</a>\u001b[0m \u001b[39m# Backpropagate gradients\u001b[39;00m\n\u001b[0;32m---> <a href='vscode-notebook-cell://ssh-remote%2Benif/STER/frederikd/p3droslo/docs/src/examples/R_Aquilae.ipynb#Y134sdnNjb2RlLXJlbW90ZQ%3D%3D?line=11'>12</a>\u001b[0m loss\u001b[39m.\u001b[39;49mbackward()\n\u001b[1;32m <a href='vscode-notebook-cell://ssh-remote%2Benif/STER/frederikd/p3droslo/docs/src/examples/R_Aquilae.ipynb#Y134sdnNjb2RlLXJlbW90ZQ%3D%3D?line=12'>13</a>\u001b[0m \u001b[39m# Update parameters\u001b[39;00m\n\u001b[1;32m <a href='vscode-notebook-cell://ssh-remote%2Benif/STER/frederikd/p3droslo/docs/src/examples/R_Aquilae.ipynb#Y134sdnNjb2RlLXJlbW90ZQ%3D%3D?line=13'>14</a>\u001b[0m optimizer\u001b[39m.\u001b[39mstep()\n",
"File \u001b[0;32m~/.local/lib/python3.9/site-packages/torch/_tensor.py:487\u001b[0m, in \u001b[0;36mTensor.backward\u001b[0;34m(self, gradient, retain_graph, create_graph, inputs)\u001b[0m\n\u001b[1;32m 477\u001b[0m \u001b[39mif\u001b[39;00m has_torch_function_unary(\u001b[39mself\u001b[39m):\n\u001b[1;32m 478\u001b[0m \u001b[39mreturn\u001b[39;00m handle_torch_function(\n\u001b[1;32m 479\u001b[0m Tensor\u001b[39m.\u001b[39mbackward,\n\u001b[1;32m 480\u001b[0m (\u001b[39mself\u001b[39m,),\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 485\u001b[0m inputs\u001b[39m=\u001b[39minputs,\n\u001b[1;32m 486\u001b[0m )\n\u001b[0;32m--> 487\u001b[0m torch\u001b[39m.\u001b[39;49mautograd\u001b[39m.\u001b[39;49mbackward(\n\u001b[1;32m 488\u001b[0m \u001b[39mself\u001b[39;49m, gradient, retain_graph, create_graph, inputs\u001b[39m=\u001b[39;49minputs\n\u001b[1;32m 489\u001b[0m )\n",
"File \u001b[0;32m~/.local/lib/python3.9/site-packages/torch/autograd/__init__.py:193\u001b[0m, in \u001b[0;36mbackward\u001b[0;34m(tensors, grad_tensors, retain_graph, create_graph, grad_variables, inputs)\u001b[0m\n\u001b[1;32m 189\u001b[0m inputs \u001b[39m=\u001b[39m (inputs,) \u001b[39mif\u001b[39;00m \u001b[39misinstance\u001b[39m(inputs, torch\u001b[39m.\u001b[39mTensor) \u001b[39melse\u001b[39;00m \\\n\u001b[1;32m 190\u001b[0m \u001b[39mtuple\u001b[39m(inputs) \u001b[39mif\u001b[39;00m inputs \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39melse\u001b[39;00m \u001b[39mtuple\u001b[39m()\n\u001b[1;32m 192\u001b[0m grad_tensors_ \u001b[39m=\u001b[39m _tensor_or_tensors_to_tuple(grad_tensors, \u001b[39mlen\u001b[39m(tensors))\n\u001b[0;32m--> 193\u001b[0m grad_tensors_ \u001b[39m=\u001b[39m _make_grads(tensors, grad_tensors_, is_grads_batched\u001b[39m=\u001b[39;49m\u001b[39mFalse\u001b[39;49;00m)\n\u001b[1;32m 194\u001b[0m \u001b[39mif\u001b[39;00m retain_graph \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[1;32m 195\u001b[0m retain_graph \u001b[39m=\u001b[39m create_graph\n",
"File \u001b[0;32m~/.local/lib/python3.9/site-packages/torch/autograd/__init__.py:88\u001b[0m, in \u001b[0;36m_make_grads\u001b[0;34m(outputs, grads, is_grads_batched)\u001b[0m\n\u001b[1;32m 86\u001b[0m \u001b[39mif\u001b[39;00m out\u001b[39m.\u001b[39mrequires_grad:\n\u001b[1;32m 87\u001b[0m \u001b[39mif\u001b[39;00m out\u001b[39m.\u001b[39mnumel() \u001b[39m!=\u001b[39m \u001b[39m1\u001b[39m:\n\u001b[0;32m---> 88\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mRuntimeError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mgrad can be implicitly created only for scalar outputs\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[1;32m 89\u001b[0m new_grads\u001b[39m.\u001b[39mappend(torch\u001b[39m.\u001b[39mones_like(out, memory_format\u001b[39m=\u001b[39mtorch\u001b[39m.\u001b[39mpreserve_format))\n\u001b[1;32m 90\u001b[0m \u001b[39melse\u001b[39;00m:\n",
"\u001b[0;31mRuntimeError\u001b[0m: grad can be implicitly created only for scalar outputs"
]
}
],
"source": [
"optimizer = Adam([m], lr=0.1)\n",
"\n",
"f = 1/m\n",
"mask = (m == 0)\n",
"f[mask] = torch.ones(1)\n",
"\n",
"loss = f**2\n",
"\n",
"# Set gradients to zero\n",
"optimizer.zero_grad()\n",
"# Backpropagate gradients\n",
"loss.backward()\n",
"# Update parameters\n",
"optimizer.step()"
]
},
{
"cell_type": "code",
"execution_count": 121,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 121,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m.requires_grad"
]
},
{
"cell_type": "code",
"execution_count": 127,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tensor([1.0000, 1.0000, 0.5000, 0.3333, 0.2500], grad_fn=<IndexPutBackward0>)"
]
},
"execution_count": 127,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f"
]
},
{
"cell_type": "code",
"execution_count": 112,
Expand Down
9 changes: 5 additions & 4 deletions docs/src/index.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,12 +15,13 @@ Welcome to the *p3droslo* documentation!
About
*****

*p3droslo* is a library for **probabilistic 3D reconstruction** of astronomical **spectral line observations**.
*p3droslo* is a python package that allows you to create **probabilistic 3D reconstructions** of astronomical **spectral line observations**.

Observations of `spectral lines <https://en.wikipedia.org/wiki/Spectral_line>`_ are indespensible in astronomy, since they encode a wealth of information about the physical and chemical conditions of the medium in which they were formed.
Observations of `spectral lines <https://en.wikipedia.org/wiki/Spectral_line>`_ are indespensible in astronomy, since they encode a wealth of information about the physical and chemical conditions of the medium from which they originate.
For instance, their narrow extent in frequency space make them very sensitive to Doppler shifts, such that their shape encodes the motion of the medium along the line of sight.
As a result, given a good model for the line formation process and a good inversion method, these physical and chemical properties can be retrieved from observations.
As a result, given a good model for the line formation process and an inversion method, these physical and chemical properties can be retrieved from observations.
Currently, we mainly focus on retrieving the distributions of the abundance of the chemical species producing the line, the velocity field, and its kinetic temperature.
However, also other parameters can be retrieved.

More information about the model for `spectral line formation <https://p3droslo.readthedocs.io/en/latest/background/spectral_line_formation.html>`_ and the `probabilistic reconstruction methods <https://p3droslo.readthedocs.io/en/latest/background/probabilistic_reconstruction.html>`_ can be found in the `background <https://p3droslo.readthedocs.io/en/latest/background/index.html>`_ pages.

Expand Down Expand Up @@ -67,7 +68,7 @@ We are open to contributions to p3droslo. More information can be found `here <h
Collaborating
*************

We are always interested in collaborating! If you have a project, or idea, or data that might benefit from probabilistic 3D reconstruction or any other method you can find in this repository, feel free to contact `me <https://freddeceuster.github.io>`_.
We are always interested in collaborating! If you have a project or data that might benefit from probabilistic 3D reconstruction or any other method you can find in this repository, feel free to contact `me <https://freddeceuster.github.io>`_.


Acknowledgements
Expand Down
96 changes: 48 additions & 48 deletions src/p3droslo/hydrodynamics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,62 +12,62 @@ def __setup__(self, gamma=1.2, mu=2.381*constants.u.si.value):



gamma = 1.2
mu = 2.381 * constants.u.si.value
# gamma = 1.2
# mu = 2.381 * constants.u.si.value

def steady_state_hydrodynamic_loss(model, f_x=0.0, f_y=0.0, f_z=0.0, heating_m_cooling=0.0):
"""
Loss assuming steady state hydrodynamics, i.e. vanishing time derivatives in Euler's equations.
"""
# def steady_state_hydrodynamic_loss(model, f_x=0.0, f_y=0.0, f_z=0.0, heating_m_cooling=0.0):
# """
# Loss assuming steady state hydrodynamics, i.e. vanishing time derivatives in Euler's equations.
# """

r = model.get_coords(origin=origin_ind)
d = torch.from_numpy(r / np.linalg.norm(r, axis=0)**3)
# r = model.get_coords(origin=origin_ind)
# d = torch.from_numpy(r / np.linalg.norm(r, axis=0)**3)

log_rho = model['log_CO']
log_tmp = model['log_temperature']
log_M = model['log_M']
# log_rho = model['log_CO']
# log_tmp = model['log_temperature']
# log_M = model['log_M']

rho = torch.exp(log_rho)
tmp = torch.exp(log_tmp)
M = torch.exp(log_M)
# rho = torch.exp(log_rho)
# tmp = torch.exp(log_tmp)
# M = torch.exp(log_M)

v_x = v_max * model['velocity_x']
v_y = v_max * model['velocity_y']
v_z = v_max * model['velocity_z']
# v_x = v_max * model['velocity_x']
# v_y = v_max * model['velocity_y']
# v_z = v_max * model['velocity_z']

kBT_o_mu = (constants.k_B.si.value / mu) * tmp
# kBT_o_mu = (constants.k_B.si.value / mu) * tmp

# Energy
eng = 0.5 * (v_x**2 + v_y**2 + v_z**2) + (gamma / (gamma - 1.0)) * kBT_o_mu
# # Energy
# eng = 0.5 * (v_x**2 + v_y**2 + v_z**2) + (gamma / (gamma - 1.0)) * kBT_o_mu

# log rho + log T
log_rho_p_log_tmp = log_rho + log_tmp
# # log rho + log T
# log_rho_p_log_tmp = log_rho + log_tmp

f_x = -constants.G.si.value * M * d[0]
f_y = -constants.G.si.value * M * d[1]
f_x = -constants.G.si.value * M * d[2]

# Continuity equation (steady state): div(ρ v) = 0
loss_cont = model.diff_x(rho * v_x) + model.diff_y(rho * v_y) + model.diff_z(rho * v_z)

# Momentum equation (steady state): v . grad(v) + grad(P) / rho = f
loss_momx = v_x * model.diff_x(v_x) + v_y * model.diff_y(v_x) + v_z * model.diff_z(v_x) + kBT_o_mu * model.diff_x(log_rho_p_log_tmp) - f_x
loss_momy = v_x * model.diff_x(v_y) + v_y * model.diff_y(v_y) + v_z * model.diff_z(v_y) + kBT_o_mu * model.diff_y(log_rho_p_log_tmp) - f_y
loss_momz = v_x * model.diff_x(v_z) + v_y * model.diff_y(v_z) + v_z * model.diff_z(v_z) + kBT_o_mu * model.diff_z(log_rho_p_log_tmp) - f_z

# Energy equation (steady state): div(u v) = 0
loss_engy = rho * (model.diff_x(eng) * v_x + model.diff_y(eng) * v_y + model.diff_z(eng) * v_z) - heating_m_cooling

# Compute the mean squared losses
losses = torch.stack([
((loss_cont/ rho )**2).mean(),
((loss_momx/ v_x )**2).mean(),
((loss_momy/ v_y )**2).mean(),
((loss_momz/ v_z )**2).mean(),
((loss_engy/(rho*eng))**2).mean()
])

# Return losses
return losses
# f_x = -constants.G.si.value * M * d[0]
# f_y = -constants.G.si.value * M * d[1]
# f_x = -constants.G.si.value * M * d[2]

# # Continuity equation (steady state): div(ρ v) = 0
# loss_cont = model.diff_x(rho * v_x) + model.diff_y(rho * v_y) + model.diff_z(rho * v_z)

# # Momentum equation (steady state): v . grad(v) + grad(P) / rho = f
# loss_momx = v_x * model.diff_x(v_x) + v_y * model.diff_y(v_x) + v_z * model.diff_z(v_x) + kBT_o_mu * model.diff_x(log_rho_p_log_tmp) - f_x
# loss_momy = v_x * model.diff_x(v_y) + v_y * model.diff_y(v_y) + v_z * model.diff_z(v_y) + kBT_o_mu * model.diff_y(log_rho_p_log_tmp) - f_y
# loss_momz = v_x * model.diff_x(v_z) + v_y * model.diff_y(v_z) + v_z * model.diff_z(v_z) + kBT_o_mu * model.diff_z(log_rho_p_log_tmp) - f_z

# # Energy equation (steady state): div(u v) = 0
# loss_engy = rho * (model.diff_x(eng) * v_x + model.diff_y(eng) * v_y + model.diff_z(eng) * v_z) - heating_m_cooling

# # Compute the mean squared losses
# losses = torch.stack([
# ((loss_cont/ rho )**2).mean(),
# ((loss_momx/ v_x )**2).mean(),
# ((loss_momy/ v_y )**2).mean(),
# ((loss_momz/ v_z )**2).mean(),
# ((loss_engy/(rho*eng))**2).mean()
# ])

# # Return losses
# return losses

# def loss_cont(self, )

0 comments on commit 1ba337e

Please sign in to comment.