Multi-Simulation-Setup.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup for Multiple Parameter Simulations\n",
    "\n",
    "This notebook shows how to set up multiple simulation runs for systematic parameter studies. \n",
    "\n",
    "## 1) Generate directories for parameter scan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start by importing some modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "import numpy as np\n",
    "import tifffile as tif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define which input parameter should be varied. Possible keys are *memory*, *aperture*, *branching_angle*, \n",
    "*bifurcation_angle*,*branchingprob*,*bifurcationprob*,*deathprob* and *reactivationprob*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_key = \"memory\";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The range of values that it should take are set in this numpy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_range = [1,2,3,4,5];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now set and create the directory for this simulation if it does not yet exist. This notebook needs to run from the base directory of the framework:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "home = os.getcwd();\n",
    "basedir = \"Examples/Osteocytes\"\n",
    "simdir = os.path.join(home,basedir,var_key);\n",
    "\n",
    "if not os.path.exists(simdir):\n",
    "        os.makedirs(simdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, create one subdirectoy for each value of the parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for val in var_range:\n",
    "    new_path = basedir + \"/\" + str(var_key) + \"/\" + str(val);\n",
    "    if not os.path.exists(new_path):\n",
    "        os.makedirs(new_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This cell creates subdirectories for each parameter value corresponding to the number of identical runs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "runs = 10\n",
    "\n",
    "for v in var_range:\n",
    "    for r in range(runs):\n",
    "        save_dir = os.path.join(home,simdir,str(v),str(r),);\n",
    "        if not os.path.exists(save_dir):\n",
    "            os.mkdir(save_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## 2) Copy the parameter files to the new simulation directories\n",
    "\n",
    "Make copies of all configuration files with the standard parameters for each parameter configuration.\n",
    "\n",
    "Then, in the last line, `sed` is used to actually set the varying parameter to the predefined values in each directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for val in var_range:\n",
    "    new_path = basedir + \"/\" + str(var_key) + \"/\" + str(val);\n",
    "    if not os.path.exists(new_path):\n",
    "        os.makedirs(new_path)\n",
    "    os.system(\"cp \"+basedir+\"/growth_parameters.csv \"+new_path)\n",
    "    os.system(\"cp \"+basedir+\"/starting_positions_ocy.csv \"+new_path)\n",
    "    os.system(\"cp \"+basedir+\"/structured_image_dir.csv \"+new_path)\n",
    "    os.system(\"cp \"+basedir+\"/internal_parameters.csv \"+new_path)\n",
    "    os.system(\"sed -i \\'s/\"+str(var_key)+\".*/\"+str(var_key)+\",\"+str(val)+\"/g\\' \"+str(new_path)+\"/internal_parameters.csv\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "The following cells only need to be run if new starting positions and environments should be used, or if slurm job files should be created. \n",
    "\n",
    "***\n",
    "\n",
    "## 3) Define new starting positions (optional)\n",
    "Run this cell if you want to create a new set of starting positions and directions. Otherwise, the existing file will be used and copied to the new simulation directories.\n",
    "\n",
    "Here, we create the cell array for the Osteocyte Network scenario:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cells = np.array([[128,64,128],[128,128,128],[128,192,128],\n",
    "                  [64,96,128], [64,160,128],\n",
    "                  [192,96,128],[192,160,128],\n",
    "                  [128,96,96], [128,160,96],\n",
    "                  [64,64,96],  [64,128,96],  [64,192,96],\n",
    "                  [192,64,96], [192,128,96], [192,192,96],\n",
    "                  [128,96,160],[128,160,160],\n",
    "                  [64,64,160], [64,128,160], [64,192,160],\n",
    "                  [192,64,160],[192,128,160],[192,192,160]]);\n",
    "\n",
    "N_cells = cells.shape[0]\n",
    "\n",
    "# offset of cell surface from cell center\n",
    "w = 7\n",
    "l = 12\n",
    "h = 4\n",
    "\n",
    "# number of outgrowing processes per cell\n",
    "N = 20\n",
    "\n",
    "# equally distributed angles\n",
    "idx = np.arange(0, N, dtype=float) + 0.5;\n",
    "phi = np.arccos(1 - 2*idx/N);\n",
    "theta = np.pi * idx * (1 + 5**0.5);\n",
    "\n",
    "# azimuth and polar angles of initial direction\n",
    "az = np.mod(np.rad2deg(phi),180);\n",
    "pol = np.mod(np.rad2deg(theta),360);\n",
    "\n",
    "# positions on cell surface relative to center\n",
    "x = w * np.cos(theta) * np.sin(phi);\n",
    "y = l * np.sin(theta) * np.sin(phi);\n",
    "z = h * np.cos(phi)\n",
    "\n",
    "# generate column vectors for positions of all cells\n",
    "xc = np.squeeze(np.repeat(cells[:,0],N) + np.tile(x,(1,N_cells)));\n",
    "yc = np.squeeze(np.repeat(cells[:,1],N) + np.tile(y,(1,N_cells)));\n",
    "zc = np.squeeze(np.repeat(cells[:,2],N) + np.tile(z,(1,N_cells)));\n",
    "\n",
    "azc = np.squeeze(np.tile(az,(1,N_cells)));\n",
    "polc = np.squeeze(np.tile(pol,(1,N_cells)));\n",
    "\n",
    "# index column to identify processes from the same cell\n",
    "idc = np.repeat(np.arange(1, N_cells+1, dtype=float),N);\n",
    "\n",
    "# save starting positions to text file\n",
    "np.savetxt(os.path.join(basedir,'starting_positions_ocy.csv'),np.column_stack((xc,yc,zc,azc,polc,idc)), fmt='%.2f', delimiter=',', header='Starting Positions\\nX,Y,Z,Az,Pol,CN')\n",
    "\n",
    "# now copy again into all the subdirectories to get the new file:\n",
    "for val in var_range:\n",
    "    new_path = basedir + \"/\" + str(var_key) + \"/\" + str(val);\n",
    "    os.system(\"cp \"+basedir+\"/starting_positions_ocy.csv \"+new_path)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## 4) Create the structured grid for the environment (optional)\n",
    "\n",
    "Run this cell to re-create the structured images and feature maps, common to all simulation runs. Otherwise, the already existing images will be used.\n",
    "\n",
    "In this example, we create an anisotropic bone matrix for the Osteocyte Network scenario:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Initialize_Network_Growing as ing\n",
    "import Start_Conditions.StructField_Maker as StructField_Maker# Path to this notebook\n",
    "\n",
    "path_structured_images = os.path.join(basedir,'structured_image_dir.csv');\n",
    "\n",
    "# generate and save empty volume \n",
    "fe = np.zeros((256,256,256));\n",
    "tif.imsave(os.path.join(basedir,'empty.tif'),np.float32(fe))\n",
    "\n",
    "# create 256³ structured image procedurally\n",
    "f = np.zeros((4,4));\n",
    "\n",
    "# some collagen fibers\n",
    "f[0:2,0:2] = 300;\n",
    "\n",
    "# tile to fill volume\n",
    "f2 = np.tile(f, [256,64,64]);\n",
    "\n",
    "# introduce some lamellae with perpendicular fibers\n",
    "f2[64:256:64,:,::16] = 1;\n",
    "f2[64:256:64,::16,:] = 1;\n",
    "\n",
    "# save at the appropriate position\n",
    "tif.imsave(os.path.join(basedir,'fiberz.tif'),np.float32(f2))\n",
    "\n",
    "# Generate feature maps from image data of growth environment\n",
    "features = ing.StructField_Maker.structured_field(path_structured_images, home, sigma_divd=2, sigma_divt1=2, sigma_divt2=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## 5) Create slurm job file for each parameter value (optional)\n",
    "\n",
    "If you run this cell, a slurm .job file will be created that defines a job array for each parameter value, corresponding to the number of runs defined above. \n",
    "\n",
    "This example is for the Julia HPC cluster at JMU. Settings for partition, number of cores, starting the environment etc. may have to be adapted to your cluster.\n",
    "\n",
    "From the slurm head node, the job can then be submitted by `sbatch jobfile.job`.\n",
    "\n",
    "The simulations can also be run individually by calling:\n",
    "\n",
    "`python run_growth.py /Examples/Osteocytes/$var_key/$var_value $run`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path and conda environment on the compute nodes (may differ from where this notebook runs)\n",
    "slurmhome = \"/home/phk57mf/pythrahyper_net/\";\n",
    "env = \"pythra-mayavi\"\n",
    "\n",
    "for v in var_range:\n",
    "\n",
    "    job_file = os.path.join(home,simdir,str(v),var_key+\"_\"+str(v)+\".job\")\n",
    "\n",
    "    with open(job_file,'w') as f:\n",
    "\n",
    "            f.writelines(\"#!/bin/bash\\n\")\n",
    "            f.writelines(\"#SBATCH --job-name=OCY_\"+str(var_key)+\"_\"+str(v)+\"_%a_\"+\".job\\n\")\n",
    "            f.writelines(\"#SBATCH --output=.out\"+\"/OCY_\"+str(var_key)+\"_\"+str(v)+\"_%a.out\\n\")\n",
    "            f.writelines(\"#SBATCH --error=.out\"+\"/OCY_\"+str(var_key)+\"_\"+str(v)+\"_%a.err\\n\")\n",
    "            f.writelines(\"#SBATCH --time=2-00:00\\n\")\n",
    "            f.writelines(\"#SBATCH --mem=16000\\n\")\n",
    "            f.writelines(\"#SBATCH --partition=standard\\n\")\n",
    "            f.writelines(\"#SBATCH -c 8\\n\")\n",
    "            f.writelines(\"#SBATCH -n 1\\n\")\n",
    "            f.writelines(\"#SBATCH --array=0-\"+str(runs-1)+\"\\n\")\n",
    "            f.writelines(\"source activate \"+env+\"\\n\")\n",
    "            f.writelines(\"cd \"+slurmhome+\"\\n\")\n",
    "            f.writelines(\"python run_growth.py /\"+os.path.join(basedir,var_key,str(v))+\"/ $SLURM_ARRAY_TASK_ID\\n\")\n",
    "            \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}