Import Data, Parse It, and Crunch It

These days I do a lot of data import and massage it one way or another. We are unraveling an interesting technical issue: How does mild electrical stimulation of the spinal cord (aka neuromodulation) stop chronic, debilitating pain signals? No one knows. We are close to presenting the most detailed theory yet elucidated. The project was conceived by Dr. Jeffrey Arle (neurosurgeon) and Dr. Jay Shils (electrical engineer and neurophysiologist), pioneers of neuromodulation. Jeff Arle came up with the key insight into how neuromodulation might work, and it appears in concept he is right (look for our article later this year).

Researchers around the world are exploring neuromodulation for diverse disorders such as chronic pain, movement disorders, depression and drug dependence.

Drs. Arle and Shils lead the way with their new book, Essential Neuromodulation (Academic Press 2011, 504 pp), which compiles critical work by the field's major contributors and provides a springboard for the new applications.

Nick Iftimia did the initial, very detailed research on the spinal cord. I built the initial 2200-connection model and then Nick built the more sophisticated 11,000-connection model. Then I used SolidWorks computer-aided design (CAD) software to build an accurate 3D model of three spinal cord segments and various electrical stimulators, of which the "business end" is an array of platinum-iridium electrodes.

Jay Shils selected COMSOL for our finite-element analysis and I imported the SolidWorks model into it and modeled the physics, that is, the electric potential field generated by the stimulator electrode array. Then I export part of the model, import it into Mathematica, and crunch the data for further processing by Longzhi Mei, our chief programmer, who works in C++ and Java.

All that is a prelude to how I use Mathematica in just this one portion of the project (I use it extensively elsewhere, too). Here is typically how I proceed. First I set a file path.

COMSOLFilePath="C:\\Users\\82V\\Documents\\Neuroscience\\Arle-Shils\\COMSOL\\Dura Thickening-Scar Studies\\Multi-Level-Scar-Studies\\5-Level Scar Model with Revised Geometry\\Calibration to Lee";

And then a file directory path.

SetDirectory@FileNameDrop[COMSOLFilePath,-1]

C:\Users\82V\Documents\Neuroscience\Arle-Shils\COMSOL\Dura Thickening-Scar Studies\Multi-Level-Scar-Studies\5-Level Scar Model with Revised Geometry\Calibration to Lee

I might check the FileNames in the directory.

FileNames@"*.txt"

{New-Control-WhiteMatter-1p5V-300x20x300.txt,New-Control-WhiteMatter-1p65V-300x20x300.txt,New-Control-WhiteMatter-1p6V-300x20x300.txt,New-Control-WhiteMatterV-Cathode-3v_Anodes3v_300x20x300.txt,Old-Control-WhiteMatterV-300x20x300.txt,Old-Control-WhiteMatterV-30x20x30.txt}

Then, for in-Notebook use, I just use FileNameSetter to create a convenient button for navigating and selecting a file, with the familiar Windows dialog box, to process.

FileNameSetter@Dynamic@COMSOLFilePath

FileNameSetter with Dynamic allows using the button repeatedly to select different files.Therefore I might need to just check which filepath is currently assigned to the variable, COMSOLFilePath.

COMSOLFilePath

C:\Users\82V\Documents\Neuroscience\Arle-Shils\COMSOL\Dura Thickening-Scar Studies\Multi-Level-Scar-Studies\5-Level Scar Model with Revised Geometry\Calibration to Lee\VMax@-0p65.txt

Now, obviously all those housekeeping functions are preliminary to doing something interesting with the data in the file. I Import it, and, briefly, I've learned that the "Table" format usually gives me the cleanest Mathematica representation, but note that a great feature of Mathematica is that it has well over a hundred Import and Export formats.

mesh1p12=Import[COMSOLFilePath,"Table"];

I take a quick look at the data file header and first few lines of data, which in this case are {x, y, z, V} where V is electrical potential in volts.

mesh1p12[[1;;15]]

{{%,Model:,5-level-Scar-Model-revised-Geometry-Control-Finer.mph},{%,Version:,COMSOL,4.2.0.150},{%,Date:,May,28,2012,,15:57},{%,Dimension:,3},{%,Nodes:,274151},{%,Expressions:,1},{%,Description:,Electric,potential},{%,x,y,z,V,(V)},{16.3066,4.64461,-1.77884,0.001474},{16.5566,4.47377,-1.76414,0.001851},{16.3066,4.36936,-1.74021,0.00151},{16.3066,4.48662,-1.96791,0.001486},{16.8066,4.36807,-1.73984,0.002314},{16.5556,4.17127,-1.65972,0.001872},{16.6095,4.24711,-1.89929,0.001917}}

I might check the Length of the file; here the Length measures the number of finite element mesh elements plus the header lines (8 of them).

mesh1p12//Length

274159

Now I tell Mathematica to skip the header (start at Part # 9, go to the end) and then find the Minimum of just column 4, the electrical potential. 

Min@mesh1p12[[9;;,4]]

-1.12

Those are just typical Mathematica "one-liners." Now I need to do something more sophisticated, and repeatedly. The resulting function is still concise, thanks to Mathematica incorporating a high-level programming language. Just to give you an idea of what is going on, here is a slice plot of the data at one y-value. This is what the computer in your spinal cord looks like. The gray matter, which are the neurons or processors, is the internal symmetric silhouette with the "horns," and the white matter, which are the wiring or bus lines, is the ellipse surrounding the gray matter. The white matter, that is the wires, receives signals from the peripheral nerves and the brain, feeds them to the gray matter, which is the computer that processes them, and then relays them back to brain or muscle servomotors.

ListPlot[slicePlotXZ,PlotRange->All,AspectRatio->Automatic]

Here's a overview of what the function does, as it's written in my Notebook. "NaN" is COMSOL-speak for "not a number," i.e. no numerical value. Hence I want to first get rid of all those records. Then, since we are only interested in the dorsal (back) side of the spinal cord, I lose all records with x<10, leaving the back side records to Export:

"This will 1) remove all NaN entries and 2) all records in the ventral side of the cord, as determined by the reductionCriterion, which has to be in the form, for example, {x_, ___} /; x<10 . Remember the reduction criterion is taken by DeleteCase, not Cases!"

And here's the function. For my convenience and edification, I ask it to Print various things I want to watch, such as how long it takes to Import these large files, the Length of the initial file, the Length of the reduced file, how long it takes to Export the much smaller crunched file, and the size of the reduced file in bytes.

COMSOL File Reducer

Clear@COMSOLFileReducer2;

COMSOLFileReducer2[COMSOLFilePath_String,reductionCriterion_]:=

Module[{COMSOLOutputFile=COMSOLFilePath,COMSOLwhiteMatterDataFile,COMSOLwhiteMatterDataFileNoNaN,COMSOLwhiteMatterDataFileReducedDomain,exportReducedFileName},

SetDirectory@FileNameDrop[COMSOLFilePath,-1];

Print@Timing[COMSOLwhiteMatterDataFile=Import[COMSOLOutputFile,"Table"];];

Print@Length@COMSOLwhiteMatterDataFile;

COMSOLwhiteMatterDataFileNoNaN=DeleteCases[COMSOLwhiteMatterDataFile,{___,"NaN"}];

COMSOLwhiteMatterDataFileReducedDomain=DeleteCases[COMSOLwhiteMatterDataFileNoNaN, reductionCriterion];

Print@Length@COMSOLwhiteMatterDataFileReducedDomain;

Print@COMSOLwhiteMatterDataFileReducedDomain[[1;;12]];

exportReducedFileName="Reduced"<>FileNameTake@COMSOLOutputFile;

Print@Timing[Export[exportReducedFileName,COMSOLwhiteMatterDataFileReducedDomain,"Table"];];

FileByteCount@exportReducedFileName]

And finally, here is the function in operation, starting with the command line, followed by the Print statement output, which as mentioned is just for my purview. The essential functions reduces the data file to essentials and Exports it in the appropriate format. You can see how much quicker it takes to Export the crunched file of just 1.5MB compared to the Imported file, which typically is 90MB.

In[25]:= COMSOLFileReducer2[COMSOLFilePath,{x_,y_,z_,V_}/;x<10]

During evaluation of In[25]:= {113.881,Null}

During evaluation of In[25]:= 1800008

During evaluation of In[25]:= 36942

During evaluation of In[25]:= {{%,Model:,5-level-Scar-Model-revised-Geometry-Control-Finer.mph},{%,Version:,COMSOL,4.2.0.150},{%,Date:,May,28,2012,,12:29},{%,Dimension:,3},{%,Nodes:,1800000},{%,Expressions:,1},{%,Description:,Electric,potential,in,White,Matter},{%,x,y,z,V,(V)},{10.0067,-55.9804,-6.09667,0.0000541505},{10.1067,-55.9804,-6.09667,0.0000542645},{10.2067,-55.9804,-6.09667,0.0000543788},{10.3067,-55.9804,-6.09667,0.0000544941}}

During evaluation of In[25]:= {3.62,Null}

Out[25]= 1559151

Examples of Postfix Usage and Precedence Analysis

As stated in the Mathematica archives and this blog (two parts), one of the unique aspects of my book is extensive, pervasive usage of Prefix to simplify code and save keystrokes. Like any usage of abbreviated operators, this can require knowledge of operator Precedence.


First, here's a common example of abbreviated operator precedence without Prefix. Do we need parens around the exponential coefficients?

Plot[E^-r t/.{r->.5},{t,-2,10}]

The graph shows something is wrong and Precedence analysis tells us yes, we need parens.

Precedence/@{Power,Minus,Times,ReplaceAll}

{590.,480.,400.,110.}

Plot[E^(-r t)/.{r->.5},{t,-2,10}]

Here's a typical example using Prefix. Will this function work as intended? (It's from S11: What's New in Mathematica 8 by Paul Wellin.)

f1@x_ := 1/x Cos[Log@x/x]; Plot[f1@x, {x, 0, .01}]

What would speed up Precedence analysis would be seeing an operator's Precedence popup on Mouseover, like a Tooltip. And a more minor recommendation to WRI is to make Precedence Listable to avoid needing Map. But that only saves one keystroke and I'm after much bigger game by recommending extensive usage of Prefix.

A tip is to list the operators in the order they appear in your function. Then it's easier to see which are 'stickier' than if they are out of order. Anyway, using Precedence, we see that the function will work as intended.


Unprotect@Precedence; SetAttributes[Precedence,Listable]; Protect@Precedence; Precedence@{Divide,Times,Prefix}

{470.,400.,640.}


f1@x_ := 1/x Cos[Log@x/x]; Plot[f1@x, {x, 0, .01}]

Now that we know the Precedence of Prefix and Divide, we can write the following interesting usage of Range with impunity; we don't need Table or Map. It works because Divide is Listable.

Range@5/6

{1/6,1/3,1/2,2/3,5/6}

Cyclically Selecting Parts Using Mod with Offset = 1

The Doc Center notes this usage of the offset feature in Mod to cyclically select Parts. Setting offset to 1 means Part will cycle from Part #1 through the value of n instead of the expression itself (i.e. Part #0) tbrough n - 1. However I have not seen a useful application of this idiom.

In[1]:= {a, b, c}[[Mod[Range[10], 3, 1]]]


Out[1]= {a, b, c, a, b, c, a, b, c, a}

Mod

The usage of Mod is straightforward, but by "offset" the Doc Center means shifting the cycle of remainders over by the offset along the number line. Normally Mod[m,n] is equivalent to the typical infix mathematical notation m mod n, i.e. it gives the remainder of m/x  (where x is a positive integer) as cycling from 0 through n-1. Offset rotates the cycle on the number line by its value. Here is the usage with default offset of 0:

In[149]:= Mod[Range@7, 3]


Out[149]= {1, 2, 0, 1, 2, 0, 1}

Now using an offset of 1, instead of cycling through {0,1,2}, the remainder cycles through {1,2,3}:

In[151]:= Mod[Range@7, 3, 1]


Out[151]= {1, 2, 3, 1, 2, 3, 1}

Using an offset of 2, instead of cycling between 0 and 2, the remainder cycles through {2,3,4}:

In[152]:= Mod[Range@7, 3, 2]


Out[152]= {4, 2, 3, 4, 2, 3, 4}

So what should happen if we use an offset of -1? The remainder should cycle through {-1,0,1}.

In[153]:= Mod[Range@7, 3, -1]


Out[153]= {1, -1, 0, 1, -1, 0, 1}

Capturing the Remnant of a Partitioned List

When partitioning, you usually don't want to just lose the leftover part at the end of the partitions. Mathematica provides ways to overlap Lists and pad them, but often you just want to capture the remnant leftover elements at the end. Here are two ways to do that. First, we define a function to capture the remnant and then Append it to the partitioned list, for instance in a TableForm usage.

In[77]:= TableForm@Partition[Range@11, 3]
Out[77]//TableForm=

We lost 10 and 11. Let's get them back.

In[78]:= listOfNumbers = Range@11;


remnant[listOfNumbers, 3]


Out[79]= {10, 11}


In[74]:= remnant[aList_List, partitionLength_Integer] :=
 Take[aList, -Mod[Length@aList, partitionLength]]

The most compact way to write a function that requires another short function, such as remnant, is to just protect the variable name and define the function at the same time within Module's local scope list.

In[80]:= unevenListTable[aList_List, partitionLength_Integer] :=
 Module[{remnant = Take[aList, -Mod[Length@aList, partitionLength]]},
  Partition[aList, partitionLength] // Append[#, remnant] & // TableForm]


In[81]:= unevenListTable[Range@11, 3]


Out[81]//TableForm=

There is also a way using a particular Partition syntax. What this says is, "Partition the list into sublists of length 4, take them in blocks of 4 (i.e. do not overlap them), start the first element of the first list at position 1 of the first sublist, and do not pad the last uneven list (i.e. pad it with the empty set)."

In[83]:= Partition[listOfNumbers, 3, 3, 1, {}] // TableForm


Out[83]//TableForm=

To get a glimpse of how this syntax works, by way of contrast, note what happens when we tell Partition to start the first element at position 3 of the first sublist. Since the sublists are only length 3, the first element takes that third slot and then Partition creates sublists of length 3 until it runs into too few elements to do that—the remnant, which is 11—and ignores them.

In[84]:= Partition[Range@11, 3, 3, 3, {}]


Out[84]= {{1}, {2, 3, 4}, {5, 6, 7}, {8, 9, 10}}

RecurrenceTable - Introduction

RecurrenceTable defines a sequence by setting f[n+1] equal to some function of f[n], and an initial condition for f[1]. Recurrence relations are also called difference equations and the method of finite differences uses recurrence relations. These are discrete versions of differential calculus. See the excellent articles in Wolfram Mathworld.

In our computer simulations of neural circuits, while the equations governing the neurons' behavior are phrased as differential equations, in reality we compute them as difference equations. For one thing we add noise to them, and for another, in practice they contain singularities. The "forward Euler method" et cetera won't work. I suspect that many prima facie differential equations are, in practice, unsolvable, and mask an underlying reality actually described by difference equations.

Here is a comparison between Table and RecurrenceTable that simply shows how to use the syntax of RecurrenceTable.

RecurrenceTable

In[21]:= powersOfOneHalf = Table[2^-i, {i, 1, 25}] // N

Out[21]= {0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.0078125, 0.00390625, \
0.00195313, 0.000976563, 0.000488281, 0.000244141, 0.00012207, 0.0000610352, \
0.0000305176, 0.0000152588, 7.62939*10^-6, 3.8147*10^-6, 1.90735*10^-6,
 9.53674*10^-7, 4.76837*10^-7, 2.38419*10^-7, 1.19209*10^-7, 5.96046*10^-8,
 2.98023*10^-8}

In[23]:= RecurrenceTable[{a[n + 1] == .5 a@n, a@1 == 0.5}, a, {n, 25}]

Out[23]= {0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625, 0.0078125, 0.00390625, \
0.00195313, 0.000976563, 0.000488281, 0.000244141, 0.00012207, 0.0000610352, \
0.0000305176, 0.0000152588, 7.62939*10^-6, 3.8147*10^-6, 1.90735*10^-6,
 9.53674*10^-7, 4.76837*10^-7, 2.38419*10^-7, 1.19209*10^-7, 5.96046*10^-8,
 2.98023*10^-8}

We can define f[n] in terms of more than one preceding term. If so, we need to supply the initial conditions for as many terms as specify the recurrence relations. Here is a Fibonacci-type series, a different flavor defined using three preceding terms rather than two. The initial conditions can also be thought of as boundary conditions.

In[4]:= RecurrenceTable[{a[n] == a[n - 1] + a[n - 2] + a[n - 3], a[1] == 1, a[2] == 1,
   a@3 == 1}, a,
   {n, 10}]

Out[4]= {1, 1, 1, 3, 5, 9, 17, 31, 57, 105}

Beginners: Mathematica Front End, Kernel, and Notebooks

Front End and Kernel

There are two components to Mathematica: the Front End and the Kernel. Unless you decide to use a command-line (aka "text-based") interface at a terminal, you are using the graphical user interface (GUI) of Mathematica called the Front End. If you program in Mathematica, the Front End is also your programming environment. You will find the Front End to be a superior GUI and programming environment once you get used to it. However the part of Mathematica that does the computation you request through the Front End is called the Kernel. The Kernel can run on the same computer as the Front End or on remote computers.

Mathematica Notebooks 



Most often you work in a Mathematica Notebook (the filename ends in .nb). Think of a Notebook as an electronic laboratory notebook attached to a very powerful computer (the kernel). You can do almost anything you can think of in a Notebook, such as write text, mathematical formulae, a program, or even a paper, book or webpage. You can communicate with Mathematica's help files, which are in the Documentation Center, request information, perhaps on a function, to be displayed right in the notebook where you are working, or in a separate browser within the Mathematica environment. And you can communicate with the kernel and ask it to do a simple or enormously complicated computation, even if the kernel is located at a remote site.

The Mathematica Notebook consists of a hierarchy of nested cells. You select an entire cell by clicking it's bracket in the right margin. You select contents within a cell in the normal way. Within a cell, the insertion point is a prominent blinking vertical bar, while the cursor is a vertical bar with a capital and base (an "I-beam"). But note that between cells, the insertion point is a prominent horizontal bar running all the way across the notebook, while the cursor is a horizontal I-beam that disappears when not quite between the cells.

You will soon automatically notice whether a cell bracket\[LongDash]the vertical, terminated bar at the right of the cell\[LongDash]is denoted at the top by two small horizontal lines or a triangle. The two small lines indicate an open cell, and the triangle indicates a closed cell. To open the closed cell and view its content, double-click anywhere on the cell bracket. Likewise, in the left margin, open cells are indicated by an upward arrowhead (an angle bracket) and closed cells by a downward arrowhead, and you can toggle open or closed by clicking on the arrowhead.

Clicking Enter on the alphabetic keyboard is equivalent to Return, not Enter. Shift-Enter on the alphabetic keyboard and Enter on the numeric keypad (provided NumLock is on) are equivalent and tell Mathematica to evaluate the current cell if it is an input cell. What "evaluate" means is explicitly explained in Mathematica if you need to understand the details.

The first time you click Enter (in an input cell) after starting a Mathematica session, there is a delay while the Kernel is started and performs its evaluation of your request.

A Mathematica mantra is "Everything is an Expression," which means that the universal, unifying syntax of Mathematica is a normal form expression with a head and argument: head[argument]. Suffice it to say here that the Mathematica Notebooks are entirely rendered via expressions and every aspect of Notebooks are potentially under your control.

<Shift> + <Enter> causes Mathematica to evaluate the contents of the cell containing the input point of that you have selected by clicking its bracket. Remember that the selected cell may not be the last cell in your Notebook, but rather, if you have "skipped around" it is the last cell you have evaluated. Also critical is that the selected cell is processed in light of all previous evaluations, no matter in what order they appear in the notebook. You can always tell what the kernel has seen by the Input and Output line numbers. I.e., that last expression the kernel processed was sent to it in the cell with the highest line number.

You will find it convenient to not always move sequentially down a Notebook but rather to stay in one cell or group of cells, or to skip around a bit. But this will not work unless you pay attention to what the kernel has seen, often using line numbers. Sometimes it is easier to just clear your definitions are start from scratch than to figure out what the kernel has seen.

These mechanics will soon become second nature and you will fall in love with the Mathematica Notebook, your working document in the elegant, unique Mathematica graphical user interface (the "front end").