Skip to content

Commit

Permalink
add geometries for geant4 (#177)
Browse files Browse the repository at this point in the history
Co-authored-by: natl <[email protected]>
  • Loading branch information
natl and natl authored Jun 19, 2022
1 parent 1e62811 commit 057f9cc
Show file tree
Hide file tree
Showing 9 changed files with 494 additions and 267 deletions.
1 change: 1 addition & 0 deletions docs/source/examples.rst
Original file line number Diff line number Diff line change
Expand Up @@ -11,4 +11,5 @@ Examples
examples/multiple-solenoids
examples/modifying-solenoids
examples/generating-geometries-for-geant
examples/generating-structures-for-geant

30 changes: 20 additions & 10 deletions docs/source/examples/generating-geometries-for-geant.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -64,13 +64,17 @@
" \"\".join(np.random.choice([\"G\", \"A\", \"T\", \"C\"], num_basepairs_turned))\n",
")\n",
"\n",
"chain_straight.to_frame().to_csv(\"results/50nm_straight.csv\", sep=' ', index=False)\n",
"chain_turned.to_frame().to_csv(\"results/50nm_turn.csv\", sep=' ', index=False)\n",
"chain_turned_twisted.to_frame().to_csv(\"results/50nm_turn_twist.csv\", sep=' ', index=False)\n",
"chain_straight.to_frame().to_csv(\"results/50nm_straight.csv\", sep=\" \", index=False)\n",
"chain_turned.to_frame().to_csv(\"results/50nm_turn.csv\", sep=\" \", index=False)\n",
"chain_turned_twisted.to_frame().to_csv(\n",
" \"results/50nm_turn_twist.csv\", sep=\" \", index=False\n",
")\n",
"\n",
"chain_straight.to_plot().savefig(\"results/50nm_straight.png\", sep=' ', index=False)\n",
"chain_turned.to_plot().savefig(\"results/50nm_turn.png\", sep=' ', index=False)\n",
"chain_turned_twisted.to_plot().savefig(\"results/50nm_turn_twist.png\", sep=' ', index=False)"
"chain_straight.to_plot().savefig(\"results/50nm_straight.png\", sep=\" \", index=False)\n",
"chain_turned.to_plot().savefig(\"results/50nm_turn.png\", sep=\" \", index=False)\n",
"chain_turned_twisted.to_plot().savefig(\n",
" \"results/50nm_turn_twist.png\", sep=\" \", index=False\n",
")"
]
},
{
Expand Down Expand Up @@ -119,7 +123,9 @@
"\n",
"chain4_straight.to_frame().to_csv(\"results/50nm_4_straight.csv\", sep=\" \", index=False)\n",
"chain4_turned.to_frame().to_csv(\"results/50nm_4_turn.csv\", sep=\" \", index=False)\n",
"chain4_turned_twisted.to_frame().to_csv(\"results/50nm_4_turn_twist.csv\", sep=\" \", index=False)\n",
"chain4_turned_twisted.to_frame().to_csv(\n",
" \"results/50nm_4_turn_twist.csv\", sep=\" \", index=False\n",
")\n",
"\n",
"chain4_straight.to_plot().savefig(\"results/50nm_4_straight.png\")\n",
"chain4_turned.to_plot().savefig(\"results/50nm_4_turn.png\")\n",
Expand Down Expand Up @@ -168,7 +174,9 @@
"solenoid_turned.translate([0, 0, -side_length / 2.0])\n",
"solenoid_turned_twisted.translate([0, 0, -side_length / 2.0])\n",
"\n",
"solenoid_straight.to_frame().to_csv(\"results/solenoid_straight.csv\", sep=\" \", index=False)\n",
"solenoid_straight.to_frame().to_csv(\n",
" \"results/solenoid_straight.csv\", sep=\" \", index=False\n",
")\n",
"solenoid_turned.to_frame().to_csv(\"results/solenoid_turned.csv\", sep=\" \", index=False)\n",
"solenoid_turned_twisted.to_frame().to_csv(\n",
" \"results/solenoid_turned_twisted.csv\", sep=\" \", index=False\n",
Expand Down Expand Up @@ -213,7 +221,7 @@
"source": [
"side_length = 1000 # angstrom\n",
"radius_solenoid = 100 # angstrom\n",
"nhistones = 51 # histones\n",
"nhistones = 51 # histones\n",
"separation = 250 # angstroms\n",
"\n",
"solenoid4_straight = dnachain.MultiSolenoidVolume(\n",
Expand Down Expand Up @@ -251,7 +259,9 @@
"solenoid4_turned.translate([0, 0, -side_length / 2.0])\n",
"solenoid4_turned_twisted.translate([0, 0, -side_length / 2.0])\n",
"\n",
"solenoid4_straight.to_frame().to_csv(\"results/solenoid4_straight.csv\", sep=\" \", index=False)\n",
"solenoid4_straight.to_frame().to_csv(\n",
" \"results/solenoid4_straight.csv\", sep=\" \", index=False\n",
")\n",
"solenoid4_turned.to_frame().to_csv(\"results/solenoid4_turned.csv\", sep=\" \", index=False)\n",
"solenoid4_turned_twisted.to_frame().to_csv(\n",
" \"results/solenoid4_turned_twisted.csv\", sep=\" \", index=False\n",
Expand Down
222 changes: 222 additions & 0 deletions docs/source/examples/generating-structures-for-geant.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,222 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Making large scale structures in Geant4\n",
"\n",
"Geant4 simulations use a large scale structure "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import sys\n",
"from pathlib import Path\n",
"\n",
"try:\n",
" # The voxelisation library produces the cubic voxelisation that\n",
" # can be used to build DNA\n",
" from fractaldna.structure_models import voxelisation as v\n",
"\n",
" # The hilbert module produces and handles L-Strings\n",
" from fractaldna.structure_models import hilbert as h\n",
"except (ImportError, ModuleNotFoundError):\n",
" sys.path.append(str(Path.cwd().parent.parent.parent))\n",
" from fractaldna.structure_models import voxelisation as v\n",
" from fractaldna.structure_models import hilbert as h\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"from mpl_toolkits.mplot3d import Axes3D"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fractal Cell Shapes\n",
"\n",
"### Generating a square shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Start with the initial L-String X for a Hilbert Curve\n",
"initial_string = \"X\"\n",
"# Iterate it as required (here, nn=3)\n",
"# for nn = 8, this takes about two hours on my 16GB RAM Mac\n",
"nn = 3\n",
"iterated_lstring = h.iterate_lstring(initial_string)\n",
"for _ in range(nn - 1):\n",
" iterated_lstring = h.iterate_lstring(iterated_lstring)\n",
"\n",
"vf = v.VoxelisedFractal.fromLString(iterated_lstring, pbar=True)\n",
"vf.center_fractal()\n",
"# fig = vf.to_plot()\n",
"# fig.savefig('results/fractal-X-3-centred.svg')\n",
"\n",
"# If you are saving a BIG fractal, try using the to_text() method instead\n",
"# as large dataframes are very slow to generate beyond 6 iterations.\n",
"# (Very slow as in multiple hours slow!)\n",
"# with open(f'results/fractal-X-{nn}-centred.csv', 'w') as ff:\n",
"# ff.write(vf.to_text())\n",
"vf.to_frame().to_csv(f\"results/fractal-X-{nn}-centred.csv\", index=False, sep=\" \")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<img src='results/fractal-X-3-centred.svg'>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Generating a Rectangular Shape\n",
"The seed `XFXFX` will generate a rectangular shape with aspect ratio 1:1:3"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Start with the initial L-String XFXFX for a Hilbert Curve\n",
"initial_string = \"XFXFX\"\n",
"# Iterate it as required (here, nn=4)\n",
"nn = 4\n",
"iterated_lstring = h.iterate_lstring(initial_string)\n",
"for _ in range(nn - 1):\n",
" iterated_lstring = h.iterate_lstring(iterated_lstring)\n",
"\n",
"vf = v.VoxelisedFractal.fromLString(iterated_lstring, pbar=True)\n",
"vf.center_fractal()\n",
"# fig = vf.to_plot()\n",
"# fig.savefig(f'results/fractal-XFXFX-{nn}-centred.svg')\n",
"\n",
"vf.to_frame().to_csv(f\"results/fractal-XFXFX-{nn}-centred.csv\", index=False, sep=\" \")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Note that the x-axis is compressed relative to the others in the below image*\n",
"\n",
"<img src='results/fractal-XFXFX-2-centred.svg'>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Importing a shape from a path\n",
"\n",
"The `voxelisation` model can convert the path of this curve to a voxelised representation, of straight\n",
"and curved boxes.\n",
"\n",
"In this example we perform this on a text file with X/Y/Z columns:\n",
"```\n",
"X\tY\tZ\n",
"-22\t-106\t216\n",
"-22\t-107\t216\n",
"-22\t-107\t215\n",
"-22\t-108\t215\n",
"-22\t-108\t214\n",
"-23\t-108\t214\n",
"-23\t-109\t214\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%%capture\n",
"df = pd.read_csv(\"results/example-path.csv\", sep=\"\\t\")\n",
"fig = plt.figure()\n",
"ax = fig.add_subplot(111, projection=\"3d\")\n",
"ax.plot(df.X, df.Y, df.Z)\n",
"fig.savefig(\"example-path.svg\")\n",
"\n",
"vf = v.VoxelisedFractal.from_path(df.values)\n",
"fig_fractal = vf.to_plot()\n",
"fig_fractal.savefig(\"example-path-voxels.svg\")\n",
"\n",
"vf.to_frame().to_csv(\"results/example-path-voxels.csv\", sep=\" \")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Left: Source Plot, Right:Voxelised Plot\n",
"\n",
"<img src='results/example-path.svg' width='45%'> <img src='results/example-path-voxels.svg' width='45%'>\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generating Random Volumes\n",
"\n",
"It can be useful to generate randomised volumes for testing a simulation. \n",
"This was the subject of (this article)[https://doi.org/10.1016/j.ejmp.2018.02.011].\n",
"\n",
"To generate a randomised volume, the `fractaldna.structure_models.random_placements`\n",
"is available.\n",
"\n",
"In that paper, 200,000 non overlapping prisms were simulated in a r=3000nm \n",
"ball.\n",
"The prisms had dimensions 30x30x100nm"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from fractaldna.structure_models import random_placements as rp\n",
"\n",
"# Generating 200,000 prisms can take around 4-5 hours\n",
"prisms = rp.generate_non_overlapping_prisms(\n",
" n_prisms=200_000,\n",
" size=[30, 30, 100], # nanometres\n",
" rad=3000, # nanometres\n",
" early_exit=-1,\n",
" verbose=True,\n",
")\n",
"\n",
"df = prisms.to_frame()\n",
"df.to_csv(\"results/prisms_200k_r3000.csv\", sep=\" \", index=False)"
]
}
],
"metadata": {
"language_info": {
"name": "python",
"pygments_lexer": "ipython3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
13 changes: 11 additions & 2 deletions fractaldna/structure_models/hilbert.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D # NOQA
from tqdm import tqdm

# interpret F as DrawForward(1);
# interpret + as Yaw(90);
Expand Down Expand Up @@ -120,7 +121,11 @@ def iterate_lstring(inString: str):


def generate_path(
lstring: str, n: int = 2, distance: float = 10.0, rounding: int = 0
lstring: str,
n: int = 2,
distance: float = 10.0,
rounding: int = 0,
pbar: bool = False,
) -> List[np.array]:
"""
Generate a path from an l-string
Expand All @@ -134,6 +139,7 @@ def generate_path(
:param n: steps on path between forward movements
:param distance: distance between points forward movements
:param rounding: rounding to apply to each position
:param pbar: display a progress bar
:return: list of XYZ points
"""
if distance < 10 ** (-1 * rounding):
Expand Down Expand Up @@ -161,7 +167,10 @@ def forward(axis, n=n, distance=distance):
r"|": lambda axis: np.dot(rotu(axis, np.pi), axis),
}

for char in lstring:
pbar_lstring = tqdm(
lstring, total=len(lstring), disable=(not pbar), desc="Generating Path"
)
for char in pbar_lstring:
if char == r"F":
lastpos = pos[len(pos) - 1]
for point in forward(axis):
Expand Down
7 changes: 2 additions & 5 deletions fractaldna/structure_models/random_placements.py
Original file line number Diff line number Diff line change
Expand Up @@ -477,9 +477,6 @@ def new_prism(rad=rad):
if verbose is True:
pbar = tqdm(iterable=None, total=n_prisms)
while n_prisms_placed < n_prisms:
if (n_prisms // 10000) * 10000 == n_prisms_placed:
sys.stderr.write(f"Generating prism {n_prisms_placed} of {n_prisms}")

if attempts >= early_exit and early_exit >= 0:
sys.stderr.write(
f"Too many failed prism placements, stopped after {n_prisms_placed} placements."
Expand All @@ -494,8 +491,8 @@ def new_prism(rad=rad):
n_prisms_placed += 1
else:
attempts += 1

pbar.close()
if verbose is True:
pbar.close()
return prisms


Expand Down
Loading

0 comments on commit 057f9cc

Please sign in to comment.