From d3e99c230be889db3f6e2b39080cd6fcd28b7ce7 Mon Sep 17 00:00:00 2001
From: Brent Huisman <brenthuisman@users.noreply.github.com>
Date: Mon, 5 Oct 2020 10:20:58 +0200
Subject: [PATCH] Update Python examples (#1166)

* Python examples: Updated to use region expressions (quotes)
* Python examples: Uniform plotting
* Python examples: Removed tutorial directory (example roughly matches python/example/single_cell_model.py)
* Python examples: Rewrote flat_cell_builder() to segment_tree()
---
 python/example/{ring.py => network_ring.py}   | 64 ++++++++---------
 python/example/single_cell_model.py           | 41 +++++++++++
 ...builder.py => single_cell_multi_branch.py} | 71 ++++++++++---------
 python/example/single_cell_swc.py             | 19 ++---
 python/tutorial/example1.py                   | 56 ---------------
 5 files changed, 111 insertions(+), 140 deletions(-)
 rename python/example/{ring.py => network_ring.py} (72%)
 create mode 100644 python/example/single_cell_model.py
 rename python/example/{single_cell_builder.py => single_cell_multi_branch.py} (57%)
 delete mode 100644 python/tutorial/example1.py

diff --git a/python/example/ring.py b/python/example/network_ring.py
similarity index 72%
rename from python/example/ring.py
rename to python/example/network_ring.py
index 17d01f10..b7d2a9bc 100644
--- a/python/example/ring.py
+++ b/python/example/network_ring.py
@@ -1,36 +1,45 @@
 import sys
 import arbor
-import matplotlib.pyplot as plt
+import pandas, seaborn
+from math import sqrt
 
 # Construct a cell with the following morphology.
 # The soma (at the root of the tree) is marked 's', and
 # the end of each branch i is marked 'bi'.
 #
-#         b2
+#         b1
 #        /
-# s----b1
+# s----b0
 #        \
-#         b3
+#         b2
 
 def make_cable_cell(gid):
-    b = arbor.flat_cell_builder()
+    # Associate labels to tags
+    labels = arbor.label_dict()
+    labels['soma'] = '(tag 1)'
+    labels['dend'] = '(tag 3)'
+
+    # Build a segment tree
+    tree = arbor.segment_tree()
+
+    # Soma (tag=1) with radius 6 μm, modelled as cylinder of length 2*radius
+    s = tree.append(arbor.mnpos, arbor.mpoint(-12, 0, 0, 6), arbor.mpoint(0, 0, 0, 6), tag=1)
+
+    # Single dendrite (tag=3) of length 100 μm and radius 2 μm attached to soma.
+    b0 = tree.append(s, arbor.mpoint(0, 0, 0, 2), arbor.mpoint(100, 0, 0, 2), tag=3)
 
-    # Soma with radius 6 μm, modelled as cylinder of length 2*radius
-    s  = b.add_cable(length=12, radius=6, name="soma", ncomp=1)
-    # Single dendrite of length 100 μm and radius 2 μm attached to soma.
-    b1 = b.add_cable(parent=s, length=100, radius=2, name="dend", ncomp=1)
-    # Attach two dendrites of length 50 μm to the end of the first dendrite.
+    # Attach two dendrites (tag=3) of length 50 μm to the end of the first dendrite.
     # Radius tapers from 2 to 0.5 μm over the length of the dendrite.
-    b2 = b.add_cable(parent=b1, length=50, radius=(2,0.5), name="dend", ncomp=1)
+    b1 = tree.append(b0, arbor.mpoint(100, 0, 0, 2), arbor.mpoint(100+50/sqrt(2), 50/sqrt(2), 0, 0.5), tag=3)
     # Constant radius of 1 μm over the length of the dendrite.
-    b3 = b.add_cable(parent=b1, length=50, radius=1, name="dend", ncomp=1)
+    b2 = tree.append(b0, arbor.mpoint(100, 0, 0, 1), arbor.mpoint(100+50/sqrt(2), -50/sqrt(2), 0, 1), tag=3)
 
     # Mark location for synapse at the midpoint of branch 1 (the first dendrite).
-    b.add_label('synapse_site', '(location 1 0.5)')
+    labels['synapse_site'] = '(location 1 0.5)'
     # Mark the root of the tree.
-    b.add_label('root', '(root)')
+    labels['root'] = '(root)'
 
-    cell = b.build()
+    cell = arbor.cable_cell(tree, labels)
 
     # Put hh dynamics on soma, and passive properties on the dendrites.
     cell.paint('"soma"', 'hh')
@@ -138,24 +147,9 @@ print('spikes:')
 for sp in spike_recorder.spikes:
     print(' ', sp)
 
-# Plot the voltage trace at the soma of each cell.
-fig, ax = plt.subplots()
+# Plot the recorded voltages over time.
+df = pandas.DataFrame()
 for gid in range(ncells):
-    times = [s.time  for s in samplers[gid].samples(arbor.cell_member(gid,0))]
-    volts = [s.value for s in samplers[gid].samples(arbor.cell_member(gid,0))]
-    ax.plot(times, volts)
-
-legends = ['cell {}'.format(gid) for gid in range(ncells)]
-ax.legend(legends)
-
-ax.set(xlabel='time (ms)', ylabel='voltage (mV)', title='ring demo')
-plt.xlim(0,tfinal)
-plt.ylim(-80,40)
-ax.grid()
-
-plot_to_file=False
-if plot_to_file:
-    fig.savefig("voltages.png", dpi=300)
-    print('voltage samples saved to voltages.png')
-else:
-    plt.show()
+    for s in samplers[gid].samples(arbor.cell_member(gid,0)):
+        df=df.append({'t/ms': s.time, 'U/mV': s.value, 'Cell': f"cell {gid}"}, ignore_index=True)
+seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV",hue="Cell").savefig('network_ring_result.svg')
diff --git a/python/example/single_cell_model.py b/python/example/single_cell_model.py
new file mode 100644
index 00000000..a7ee1ff2
--- /dev/null
+++ b/python/example/single_cell_model.py
@@ -0,0 +1,41 @@
+import arbor
+
+# (1) Create a morphology with a single (cylindrical) segment of length=diameter=6 μm
+tree = arbor.segment_tree()
+tree.append(arbor.mnpos, arbor.mpoint(-3, 0, 0, 3), arbor.mpoint(3, 0, 0, 3), tag=1)
+
+# (2) Define the soma and its center
+labels = arbor.label_dict({'soma':   '(tag 1)',
+                           'center': '(location 0 0.5)'})
+
+# (3) Create cell and set properties
+cell = arbor.cable_cell(tree, labels)
+cell.set_properties(Vm=-40)
+cell.paint('soma', 'hh')
+cell.place('center', arbor.iclamp( 10, 2, 0.8))
+cell.place('center', arbor.spike_detector(-10))
+
+# (4) Make single cell model.
+m = arbor.single_cell_model(cell)
+
+# (5) Attach voltage probe sampling at 10 kHz (every 0.1 ms).
+m.probe('voltage', 'center', frequency=10000)
+
+# (6) Run simulation for 100 ms of simulated activity.
+m.run(tfinal=100)
+
+# (7) Print spike times, if any.
+if len(m.spikes)>0:
+    print('{} spikes:'.format(len(m.spikes)))
+    for s in m.spikes:
+        print('{:3.3f}'.format(s))
+else:
+    print('no spikes')
+
+# (8) Plot the recorded voltages over time.
+import pandas, seaborn # You may have to pip install these.
+df = pandas.DataFrame({'t/ms': m.traces[0].time, 'U/mV': m.traces[0].value})
+seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV").savefig('single_cell_model_result.svg')
+
+# (9) Optionally, you can store your results for later processing.
+df.to_csv('single_cell_model_result.csv', float_format='%g')
diff --git a/python/example/single_cell_builder.py b/python/example/single_cell_multi_branch.py
similarity index 57%
rename from python/example/single_cell_builder.py
rename to python/example/single_cell_multi_branch.py
index 98ce339b..444eb0b3 100644
--- a/python/example/single_cell_builder.py
+++ b/python/example/single_cell_multi_branch.py
@@ -1,10 +1,11 @@
 import arbor
-from arbor import mechanism as mech
-from arbor import location as loc
-import matplotlib.pyplot as plt
+import seaborn
+import pandas
+from math import sqrt
 
 # Make a ball and stick cell model
-b = arbor.flat_cell_builder()
+
+tree = arbor.segment_tree()
 
 # Construct a cell with the following morphology.
 # The soma (at the root of the tree) is marked 's', and
@@ -23,32 +24,42 @@ b = arbor.flat_cell_builder()
 
 # Start with a spherical soma with radius 6 μm,
 # approximated with a cylinder of: length = diameter = 12 μm.
-s  = b.add_cable(length=12, radius=6, name="soma", ncomp=1)
+
+s = tree.append(arbor.mnpos, arbor.mpoint(-12, 0, 0, 6), arbor.mpoint(0, 0, 0, 6), tag=1)
 
 # Add the dendrite cables, labelling those closest to the soma "dendn",
 # and those furthest with "dendx" because we will set different electrical
 # properties for the two regions.
-b0 = b.add_cable(parent=s,  length=100, radius=2, name="dendn", ncomp=100)
+
+labels = arbor.label_dict()
+labels['soma'] = '(tag 1)'
+labels['dendn'] = '(tag 5)'
+labels['dendx'] = '(tag 6)'
+
+b0 = tree.append(s, arbor.mpoint(0, 0, 0, 2), arbor.mpoint(100, 0, 0, 2), tag=5)
+
 # Radius tapers from 2 to 0.5 over the length of the branch.
-b1 = b.add_cable(parent=b0, length= 50, radius=(2,0.5), name="dendn", ncomp=50)
-b2 = b.add_cable(parent=b0, length= 50, radius=1, name="dendn", ncomp=50)
-b3 = b.add_cable(parent=b1, length= 50, radius=1, name="dendx", ncomp=50)
-b4 = b.add_cable(parent=b1, length= 50, radius=1, name="dendx", ncomp=50)
+
+b1 = tree.append(b0, arbor.mpoint(100, 0, 0, 2), arbor.mpoint(100+50/sqrt(2), 50/sqrt(2), 0, 0.5), tag=5)
+b2 = tree.append(b0, arbor.mpoint(100, 0, 0, 1), arbor.mpoint(100+50/sqrt(2), -50/sqrt(2), 0, 1), tag=5)
+b3 = tree.append(b1, arbor.mpoint(100+50/sqrt(2), 50/sqrt(2), 0, 1), arbor.mpoint(100+50/sqrt(2)+50, 50/sqrt(2), 0, 1), tag=6)
+b4 = tree.append(b1, arbor.mpoint(100+50/sqrt(2), 50/sqrt(2), 0, 1), arbor.mpoint(100+2*50/sqrt(2), 2*50/sqrt(2), 0, 1), tag=6)
 
 # Combine the "dendn" and "dendx" regions into a single "dend" region.
 # The dendrites were labelled as such so that we can set different
 # properties on each sub-region, and then combined so that we can
 # set other properties on the whole dendrites.
-b.add_label('dend', '(join (region "dendn") (region "dendx"))')
+labels['dend'] = '(join (region "dendn") (region "dendx"))'
 # Location of stimuli, in the middle of branch 2.
-b.add_label('stim_site', '(location 1 0.5)')
+labels['stim_site'] = '(location 1 0.5)'
 # The root of the tree (equivalent to '(location 0 0)')
-b.add_label('root', '(root)')
-# The tips of the dendrites (3 locations at b4, b3, b5).
-b.add_label('dtips', '(terminal)')
+labels['root'] = '(root)'
+# The tips of the dendrites (3 locations at b4, b3, b2).
+labels['dtips'] = '(terminal)'
 
 # Extract the cable cell from the builder.
-cell = b.build()
+# cell = b.build()
+cell = arbor.cable_cell(tree, labels)
 
 # Set initial membrane potential everywhere on the cell to -40 mV.
 cell.set_properties(Vm=-40)
@@ -66,12 +77,16 @@ cell.place('"stim_site"', arbor.iclamp( 80, 2, 0.8))
 # Add a spike detector with threshold of -10 mV.
 cell.place('"root"', arbor.spike_detector(-10))
 
+# Discretization: the default discretization in Arbor is 1 compartment per branch.
+# Let's be a bit more precise and make that every 2 micron:
+cell.compartments_length(2)
+
 # Make single cell model.
 m = arbor.single_cell_model(cell)
 
 # Attach voltage probes, sampling at 10 kHz.
-m.probe('voltage', loc(0,0), 10000) # at the soma.
-m.probe('voltage', '"dtips"',   10000) # at the tips of the dendrites.
+m.probe('voltage', '(location 0 0)', 10000) # at the soma.
+m.probe('voltage', '"dtips"',  10000) # at the tips of the dendrites.
 
 # Run simulation for 100 ms of simulated activity.
 tfinal=100
@@ -86,20 +101,8 @@ else:
     print('no spikes')
 
 # Plot the recorded voltages over time.
-fig, ax = plt.subplots()
+df = pandas.DataFrame()
 for t in m.traces:
-    ax.plot(t.time, t.value)
-
-legend_labels = ['{}: {}'.format(s.variable, s.location) for s in m.traces]
-ax.legend(legend_labels)
-ax.set(xlabel='time (ms)', ylabel='voltage (mV)', title='cell builder demo')
-plt.xlim(0,tfinal)
-plt.ylim(-80,50)
-ax.grid()
-
-# Set to True to save the image to file instead of opening a plot window.
-plot_to_file=False
-if plot_to_file:
-    fig.savefig("voltages.png", dpi=300)
-else:
-    plt.show()
+    df=df.append( pandas.DataFrame({'t/ms': t.time, 'U/mV': t.value, 'Location': t.location, "Variable": t.variable}) )
+
+seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV",hue="Location",col="Variable").savefig('single_cell_multi_branch_result.svg')
diff --git a/python/example/single_cell_swc.py b/python/example/single_cell_swc.py
index f50768bb..bb9429d3 100644
--- a/python/example/single_cell_swc.py
+++ b/python/example/single_cell_swc.py
@@ -11,7 +11,7 @@
 import arbor
 from arbor import mechanism as mech
 from arbor import location as loc
-import matplotlib.pyplot as plt
+import pandas, seaborn
 
 # Load a cell morphology from an swc file.
 tree = arbor.load_swc('../../test/unit/swc/pyramidal.swc')
@@ -72,19 +72,8 @@ else:
     print('no spikes')
 
 # Plot the recorded voltages over time.
-fig, ax = plt.subplots()
+df = pandas.DataFrame()
 for t in m.traces:
-    ax.plot(t.time, t.value)
+    df=df.append( pandas.DataFrame({'t/ms': t.time, 'U/mV': t.value, 'Location': t.location, "Variable": t.variable}) )
 
-legend_labels = ['{}: {}'.format(s.variable, s.location) for s in m.traces]
-ax.legend(legend_labels)
-ax.set(xlabel='time (ms)', ylabel='voltage (mV)', title='swc morphology demo')
-plt.xlim(0,tfinal)
-plt.ylim(-80,80)
-ax.grid()
-
-plot_to_file=False
-if plot_to_file:
-    fig.savefig("voltages.png", dpi=300)
-else:
-    plt.show()
+seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV",hue="Location",col="Variable").savefig('single_cell_swc.svg')
diff --git a/python/tutorial/example1.py b/python/tutorial/example1.py
deleted file mode 100644
index e698363d..00000000
--- a/python/tutorial/example1.py
+++ /dev/null
@@ -1,56 +0,0 @@
-import arbor
-
-# Create a morphology with a single segment of length=diameter=5 μm
-tree = arbor.segment_tree()
-tree.append(arbor.mnpos, arbor.mpoint(-3, 0, 0, 3), arbor.mpoint(3, 0, 0, 3), tag=1)
-
-labels = arbor.label_dict({'soma': '(tag 1)', 'center': '(location 0 0.5)'})
-
-cell = arbor.cable_cell(tree, labels)
-
-# Set initial membrane potential everywhere on the cell to -40 mV.
-cell.set_properties(Vm=-40)
-# Put hh dynamics on soma, and passive properties on the dendrites.
-cell.paint('"soma"', 'hh')
-# Attach stimuli with duration of 100 ns and current of 0.8 nA.
-cell.place('"center"', arbor.iclamp( 10, duration=0.1, current=0.8))
-# Add a spike detector with threshold of -10 mV.
-cell.place('"center"', arbor.spike_detector(-10))
-
-# Make single cell model.
-m = arbor.single_cell_model(cell)
-
-# Attach voltage probes, sampling at 1 MHz.
-m.probe('voltage', '"center"',  1000000)
-
-# Run simulation for tfinal ms with time steps of 1 μs.
-tfinal=30
-m.run(tfinal, dt=0.001)
-
-# Print spike times.
-if len(m.spikes)>0:
-    print('{} spikes:'.format(len(m.spikes)))
-    for s in m.spikes:
-        print('  {:7.4f}'.format(s))
-else:
-    print('no spikes')
-
-# Plot the recorded voltages over time.
-import matplotlib.pyplot as plt
-fig, ax = plt.subplots()
-for t in m.traces:
-    ax.plot(t.time, t.value)
-
-legend_labels = ['{}: {}'.format(s.variable, s.location) for s in m.traces]
-ax.legend(legend_labels)
-ax.set(xlabel='time (ms)', ylabel='voltage (mV)', title='cell builder demo')
-plt.xlim(0,tfinal)
-plt.ylim(-80,50)
-ax.grid()
-
-# Set to True to save the image to file instead of opening a plot window.
-plot_to_file=False
-if plot_to_file:
-    fig.savefig("voltages.png", dpi=300)
-else:
-    plt.show()
-- 
GitLab