Actual tree.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append(\"/home/roozbeh/Morphology\")\n",
    "import numpy as np\n",
    "from numpy import genfromtxt\n",
    "import pandas as pd\n",
    "import McNeuron\n",
    "import matplotlib.pyplot as plt\n",
    "from copy import deepcopy\n",
    "from numpy.linalg import inv\n",
    "from collections import Counter\n",
    "import pickle\n",
    "from sympy.solvers import solve\n",
    "from sympy import nsolve\n",
    "from sympy import Symbol\n",
    "from numpy import sqrt as sqrt\n",
    "from McNeuron import asymmetric as asy\n",
    "import scipy\n",
    "import scipy.io as spio\n",
    "%matplotlib inline\n",
    "import Atree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Making swc files from cynalder model of the trees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# tree_data = Atree.data_util.pandas_non_missing()\n",
    "# pickle.dump(tree_data, open('./trees4roozbeh/tree_data.csv', 'wb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "UnicodeDecodeError",
     "evalue": "'ascii' codec can't decode byte 0xf0 in position 6: ordinal not in range(128)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mUnicodeDecodeError\u001b[0m                        Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-93a5510e4038>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mtree_data\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpickle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'./trees4roozbeh/tree_data.csv'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'rb'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mUnicodeDecodeError\u001b[0m: 'ascii' codec can't decode byte 0xf0 in position 6: ordinal not in range(128)"
     ]
    }
   ],
   "source": [
    "tree_data = pickle.load(open('./trees4roozbeh/tree_data.csv', 'rb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "result = tree_data.sort_values(['location'])\n",
    "result = result.reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mat = spio.loadmat(\"./trees4roozbeh/1.mat\", squeeze_me=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "mat['Rad'].shape, mat['Len'].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = 30\n",
    "tree = McNeuron.Neuron(result['swc'][index])\n",
    "McNeuron.visualize.plot_2D(tree,projection=np.array([[0,1,0], [0,0,1], [1,0,0]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ten = np.load('/data/roozbeh_arch/All Data/ten_sample.npy')\n",
    "\n",
    "meta_data_path = '/home/roozbeh/Neuron Morphology/All Data/meta.csv'\n",
    "regular_swc_path = '/home/roozbeh/Neuron Morphology/All Data/regular_swc.npy'\n",
    "meshed_5_swc_path = '/home/roozbeh/Neuron Morphology/All Data/meshed_5_swc.npy'\n",
    "raw_swc_path = '/home/roozbeh/Neuron Morphology/All Data/raw_swc.npy'\n",
    "morph_data  = McNeuron.data_util.read_data(meta_data_path=meta_data_path,\n",
    "                                               swc_path=raw_swc_path)\n",
    "morph_data_5  = McNeuron.data_util.read_data(meta_data_path=meta_data_path,\n",
    "                                               swc_path=meshed_5_swc_path)\n",
    "morph_data_r  = McNeuron.data_util.read_data(meta_data_path=meta_data_path,\n",
    "                                               swc_path=regular_swc_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def rall_power(self, power=1.5, branching_degree='all'):\n",
    "    \"\"\"\n",
    "    Computing the Rall power. It looks at the branching and compute the ratio of\n",
    "    Sum(diameter of children to the power) over the diameter of the mother. Notice\n",
    "    that it if the number of children can be more than two.\n",
    "    \n",
    "    Parameters:\n",
    "    -----------\n",
    "    power: float\n",
    "        the power to compute Rall power\n",
    "        \n",
    "    Returns:\n",
    "    --------\n",
    "    rall_power: numpy vector\n",
    "        A vector with the size of the number of branching nodes\n",
    "        \n",
    "    \"\"\"\n",
    "    self.set_branch_order()\n",
    "    if branching_degree is 'all':\n",
    "        branch_index = np.where(self.features['branch order'][1:]>1)[0]\n",
    "    else:\n",
    "        branch_index = np.where(self.features['branch order'][1:]==branching_degree)[0]\n",
    "    branch_index += 1\n",
    "    n_branch = len(branch_index)\n",
    "    rall_power = np.zeros(n_branch)\n",
    "    for i in range(n_branch):\n",
    "        index_branch = branch_index[i]\n",
    "        children = np.where(self.parent_index==index_branch)[0]\n",
    "        children_diam = self.diameter[children]\n",
    "        mother_diam = self.diameter[index_branch]\n",
    "        # print children_diam, mother_diam\n",
    "        if power is 'max':\n",
    "            rall_power[i] = mother_diam/np.max(children_diam)\n",
    "        else:\n",
    "            rall_power[i] = np.power(mother_diam, power)/np.power(children_diam, power).sum()\n",
    "    rall_power[np.isnan(rall_power)] = -1\n",
    "    rall_power[np.isinf(rall_power)] = -1\n",
    "    return rall_power\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "McNeuron.visualize.plot_2D(tree)\n",
    "plt.plot(tree.diameter)\n",
    "rall_power(tree, power=2, branching_degree='all')[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = 92\n",
    "power = 1./3.\n",
    "power = 2\n",
    "hist = np.zeros(19)\n",
    "for index in range(1000,2000):\n",
    "    tree = McNeuron.Neuron(morph_data_r['swc'][index])\n",
    "    #McNeuron.visualize.plot_2D(tree)\n",
    "    h = np.histogram(rall_power(tree, power=power, branching_degree='all'), bins=np.arange(0,2,.1));\n",
    "    hist = h[0] + hist\n",
    "plt.plot(hist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = 50\n",
    "tree = McNeuron.Neuron(result['swc'][index])\n",
    "tree = McNeuron.Neuron(ten[0])\n",
    "tree = McNeuron.Neuron(morph_data_r['swc'][index])\n",
    "plt.hist(rall_power(tree, power='max', branching_degree='all')[:], bins=np.arange(0,2,.1))\n",
    "plt.show()\n",
    "plt.hist(rall_power(tree, power=1, branching_degree='all')[:], bins=np.arange(0,2,.1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = 2\n",
    "tree = McNeuron.Neuron(result['swc'][index])\n",
    "plt.hist(rall_power(tree,power=.5), bins=50)\n",
    "plt.show()\n",
    "plt.hist(rall_power(tree,power=1), bins=50)\n",
    "plt.show()\n",
    "plt.hist(rall_power(tree,power=1.5), bins=50)\n",
    "plt.show()\n",
    "plt.hist(rall_power(tree,power=2.5), bins=50)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#McNeuron.visualize.plot_dendrogram(mat['CylData'][:,8].astype(int))\n",
    "#plt.plot(mat['CylData'][:,8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for i in range(0,80,10):\n",
    "    index = i\n",
    "    tree = McNeuron.Neuron(result['swc'][index])\n",
    "    print result['location'][index]\n",
    "    McNeuron.visualize.plot_2D(tree,projection=np.array([[0,1,0], [0,0,1], [1,0,0]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tree.basic_features()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(tree.discrepancy(np.arange(.1,4,.1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "L = np.array([])\n",
    "for i in range(100):\n",
    "    try:\n",
    "        mat = spio.loadmat(\"./trees4roozbeh/\"+str(i)+\".mat\", squeeze_me=True)\n",
    "        parent = mat['CPar']\n",
    "    except:\n",
    "        L = np.append(L,i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "reg = McNeuron.subsample.regular_subsample(A)\n",
    "n2 = McNeuron.Neuron(reg)\n",
    "n2.motif_features()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n2.n_node, bb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dd = n2.features['die die depth']\n",
    "bd = n2.features['branch die depth']\n",
    "bb = n2.features['branch branch depth']\n",
    "plt.plot(bd**2/(4*dd*bb))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n.features['branch die'], n.features['branch branch'], n.features['die die']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n.features['branch order']>2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "McNeuron.visualize.plot_dendrogram(n2.parent_index[:400])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(n.branch_order(n.parent_index))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n.features.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n.features['die die']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "branch_order = n.branch_order(n.parent_index)\n",
    "branch_angle, side_angle = n.branch_angle(branch_order)\n",
    "branch_angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(n.features['branch depth'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bins = np.array([18, 18, 18, 28, 30, 100])\n",
    "colors = ['red', 'blue', 'green', '#3CC2D9', '#FFA500', '#800080']\n",
    "features = ['branch angle', 'local angle', 'global angle', 'distance from root',\n",
    "          'segmental neural length', 'neuronal/euclidean for segments']\n",
    "i = 5\n",
    "a = np.histogram(n.features[features[i]],\n",
    "                 bins =  np.arange(1,1.23,.01))[0].astype(float)\n",
    "b = a/sum(a)\n",
    "plt.bar(np.arange(b.shape[0]), \n",
    "        b, \n",
    "        width=.8, \n",
    "        color=colors[i])\n",
    "#plt.axis('off')\n",
    "plt.savefig('actual-tree-neuronal-euclidean for segments.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n.features['initial segments']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.hist(n.features[features[i]], 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.hist(n.features['branch angle'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n.features['discrepancy space']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.hist(n.features['branch order'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "neuron = McNeuron.Neuron(data['X'][0])\n",
    "vector = neuron.location - neuron.location[:, neuron.parent_index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "BBdepth = n.features['branch branch depth']\n",
    "BDdepth = n.features['branch die depth']\n",
    "DDdepth = n.features['die die depth']\n",
    "asy_index = ((BDdepth+1)**2).astype(float)/(4*(BBdepth+1)*(DDdepth+1)).astype(float)\n",
    "plt.plot(asy_index[:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n.features['asymmetric ratio']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n.features['branch die']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n.features['die die']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n.features['branch branch']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = pickle.load(open(\"Primary Cell Class.p\", \"rb\" ))\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py34]",
   "language": "python",
   "name": "conda-env-py34-py"
  },
  "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}