Wednesday, September 23, 2015

How It Works: Union

Here is a recursive version of the set-theoretic function, Union, which returns the common elements of two or more sets (represented as Lists in Mathematica).

There are two base cases: 1) The union of any set with the empty set is the set itself, and 2) The union of any set with a single element set is the single element set appended to the other set. The recursive case is: 3) The union of any two Lists is the union of the first List with the First element of the second List (i.e. the single element base case), which, calling itself, recurses down the Rest of the second List until it hits the empty set base case (an empty List), at which point it unwinds back out to produce the overall union.

Union produces a set of unique elements  — no duplicates — so we need to de-dupe the List. An easy way to do that is Sort it and use a replacement Rule (we could also use the built-in DeleteDuplicates or How It Works: DeleteDuplicates).

I show two ways to control the Precedence of the abbreviate operators for ReplaceRepeated vs Postfix. Since with Precedence of 110 ReplaceRepeated is a little 'stickier' than Postfix at 70 Sort will stick to the ReplaceRepeated function unless we do something. So we can either enclose the compound function in parentheses before applying ReplaceRepeated as in the 2nd base case or we can use Postfix with Slot (#) and pure Function (&) as typical in my extended Postfix syntax as in the recursive case.


myUnion[list1_List,{}]:=list1; (* 1st base case *)

myUnion[list1_List,element_]:=(If[MemberQ[list1,element],list1,list1/.{x___}->{x,element}] //Sort)//.{a___,b_,b_,c___}->{a,b,c}(* 2nd base case *)

] (* recursive case *)







This shows the need to add a delete duplicates function:



A simple replacement Rule does the job, but use ReplaceRepeated or only the first instance of a duplicate will be replaced:




This works without the delete duplicates function.



But this would fail without it. First without delete duplicates.



And now with delete duplicates.



And the base case.



And the recursive case.



Friday, July 10, 2015

How to See the Equivalence of Select and Cases

MatchQ and PatternTest Equate Cases and Select

Select and Cases are essential Mathematica filtering functions. Select filters expressions, usually Lists, with predicates, such as ones that end with Q (evaluate ?*Q to see them) or those that don't, such as Greater, Less, and others that I list here. Cases uses patterns such as testing for a Head (e.g. _Integer) or a more general pattern (e.g. for a 2-element List, {x_,y_}.

Cases and Select are equivalent and a versatile Mathematica programmer knows how to use either one in any circumstance according to the need of the moment. The keys to seeing their equivalence are MatchQ, which translates a pattern into a predicate, and PatternTest (_?test), which translates a predicate into a pattern.





Thus you can take any usage of Cases and translate it into Select.





And you can take any usage of Select and translate it into Cases.





If we built our own predicate for EvenQ using MatchQ, it works the same way.



However note the need for parentheses around the test, which is often the case with PatternTest - the no parens usage fails.





The reason is that Blank and PatternTest both have a higher Precedence ('stickiness') than Function, and so they stick together and PatternTest ends up inside the pure Function instead of using the pure Function as the pattern test. You don't need to worry about that — just as a rule of thumb put parens around any 'homemade' test.



The Q Predicates

This is a handy way to use Partition, especially using {} so a final short sublist won't be dropped from a 'ragged' List. The syntax says

Partition[list, sublist length, offset length, start Position, padding]

and here means divide the List into Partitions of length 4, moving the window that same length 4 with each Partition so there is no overlap, start the List at Position 1,  and the empty set {} for padding tells Partition to return the final sublist even if it's shorter than what is specified by the 2nd argument, and don't pad that final short list. Evaluate this code in your Notebook:


Here's the same result in a List.



Two excellent introductions to Mathematica:

Wednesday, July 8, 2015

New in Mathematica 10: Associations as General-Purpose Associative Arrays

Associations are the most important data structure introduced in Mathematica in years. As the Guide says:

Along with lists, associations are fundamental constructs in the Wolfram Language. They associate keys with values, allowing highly efficient lookup and updating even with millions of elements. Associations provide generalizations of symbolically indexed lists, associative arrays, dictionaries, hashmaps, structs, and a variety of other powerful data structures. 

They of course will underlie databases and artificial intelligence data structures. They are intended to replace less structured expressions such as pairs in Lists or Rules. The code written for Association functions has been highly optimized for speed. I suspect that Associations are the core data structure of the Wolfram Alpha natural language processor, which is the front end for the curated data collections in Alpha, just as the Notebook is the front end for the kernel.

So when I needed a lookup table today I decided to try Associations. Here are some simple functions that define an axon's threshold in terms of the stimulating pulse seen at the axon, according to the Lapicque equation Vth = Vrh(1+τch/pw).

(You don't need to understand these!)

absoluteRheobase@fiberDiameter_:=0.000589+0.01518 Exp[-fiberDiameter/2.3477]

chronaxie@fiberDiameter_:=102.764 +162.46744 Exp[-fiberDiameter/5.53435];


Now since fiber threshold is driven by their diameter and the pulse width of the stimulating pulse, we make sample tables of those. Esc + m + esc gives us the mu in microseconds and micrometers; no need to open a Palette.

pulseWidthTable50to500=Table[pw,{pw,100,500,100}] (* in μs *);

fiberDiameters1to15=Table[diameters,{diameters,5,10}] (* in μm *);

Now here is the function that creates the Associations. It took a bit of fiddling, but by making a Table of Rules in one step as should be done with good functional programming, it is then easy to Apply Association to the entire Table, simply replacing the Head, List, with Association, which converts the List of Rules into an Association. In the Table pulseWidth is the "i" iterator and fiberDiameter is the "j" iterator, and notice that we don't need to iterate numerically from 1 to Length@pulseWidthTable50to500, for example, we just directly iterate over the List itself, which is simpler.



Here are several versions of lookup functions. The first one just uses the same syntax as we'd use to access elements of an associative array made with function@value, such as we'd construct with a memo function.

thresholdLookup[lookupTable_Association, pulseWidth_Integer, diameter_] := lookupTable@{pulseWidth, diameter}

absoluteThresholdTable@{100, 7}


Likely it's safer and more elegant to use the new Lookup function instead.

thresholdLookup2[lookupTable_Association, pulseWidth_Integer, diameter_] := Lookup[lookupTable, {{pulseWidth, diameter}}]

thresholdLookup2[absoluteThresholdTable, 100, 7]


We see that using a List as the Key results in returning a List as the Value, like in solutions for equations returned by Solve, etc. So to access the Value we can use First@key to remove the List brackets, as we often do with solutions to equations.

thresholdLookup3[lookupTable_Association, pulseWidth_Integer, diameter_] := Lookup[lookupTable, {{pulseWidth, diameter}}] // First

thresholdLookup3[absoluteThresholdTable, 100, 7]


Years ago I talked to Stephen Wolfram at an artificial life conference about artificial intelligence and his views were kind of naive. Things have changed. I expect that Mathematica and its Wolfram Language are headed inexorably toward large-scale AI applications,  approaching, and exceeding in some respects, human intelligence over the next two decades.

Sunday, June 28, 2015

How to Find Memory Used in Computations

At some point you will need to — or just want to — push the limits of Mathematica and it may be helpful to measure memory consumed during a computation. Here are methods inspired by Wagon, Mathematica in Action (2nd ed., pp 431 et seq.). (See also Memory Management Tools.)

The first method is simple; store the initial MemoryInUse in a variable, do your calc, then subtract it from the current MemoryInUse.



Second, suppose you want to monitor memory in several places in your code. Initialize the variable currentMemory, then sprinkle


through your code like Print statements, and the enclosing Reap will cause each one to output the net memory consumed in bytes between the previous use and the current one. I say 'net' memory consumed because Mathematica is designed to optimize memory use and so clears memory in use every chance it gets. Following the pattern of Timing, I reverse the output to give the memory results first, then the function's output.

(*initialize s*); s=0;


Did you notice, as Wagon points out, that Do uses minimal memory during a calculation? It cleans up after every iteration.

Like Timing, this third method of checkMemory requires you to wrap your computation in checkMemory, access the result of your computation with First and if desired suppress the output with a semi-colon. Mimicking the structure of a built-in function is an example of a good programming technique. The designers of Mathematica are master programmers and they have created and debugged templates or paradigms of function structures. Often you can follow the structure of built-in functions to create sound ones of your own.

preComputationMemory=MemoryInUse[];Print["Memory in use at start is ",preComputationMemory];
Print["Memory in use at end is ",MemoryInUse[]];{memoryUsedInComputation,computation1}

Once it's done, this computation uses 32 bytes of memory and the result is 14.3927.


Memory in use at start is 31286024
Memory in use at end is 31286056

We find the memory use with the function output for a small Table.


Memory in use at start is 31289472
Memory in use at end is 31289752

Then suppress output, with a semi-colon as is done in Timing, for a large Table.


Memory in use at start is 31340688
Memory in use at end is 111344696

I left off the summation of the Table to show MemoryInUse during the calc. If we add Total, we see the max MemoryInUse reclaimed (-80,000,024 bytes).


Memory in use at start is 31359736
Memory in use at end is 31359864

Compare memory used in a Do loop. Pretty interesting!


Memory in use at start is 31364400
Memory in use at end is 31364432



So there you have it, Do takes 30 times as long as the functional command Table[...]//Total, but uses a minute fraction of the latter's memory.

Friday, June 26, 2015

How to Query 11 Different Search Engines

You can query a search engine with URLFetch. For the URL you need a String that includes a standard prefix and a suffix that begins with something like "/?q=querykeywords". I have removed some unnecessary junk from most of the query Strings to make them as short and simple as possible.

There are several ways to functionally compose a query for a search engine. This one uses StringJoin. I use Which to handle the cases of 1) a single keyword, 2) multiple keywords, or 3) any other data type, which is an error. Using True as the third test in conditional statements is a standard way to handle 'all else' cases.


searchEngineQuery[searchEnginePrefix_String, keywords_] := 
   Head@keywords === String, 
   url = StringJoin[searchEnginePrefix, keywords],
   Head@keywords === List, 
   url = StringJoin[searchEnginePrefix, 
     Table[i <> "+", {i, Most@keywords}], Last@keywords],
   True, Print@"Error: wrong data type."];

We test it on a single keyword (I omit the full output, which you can try at home).


We try it on a List of keywords.


We try it on a wrong data type.

searchEngineQuery[queryTemplateDuckDuckGo, 5]

Error: wrong data type.

And here are sample fetches for the search engines. Again, I omit the full output but you can try them yourself. I precede each sample with a rating the quality of 'cleanliness' of results by StringLength, the idea being that the shorter length results Strings have less junk in them and more of the actual results. By that measure DuckDuckGo, Blekko and Alhea give the 'cleanest' output for further processing. But don't misconstrue that measure as a quality rating of the results themselves.



11 Search Engines Queried

Google: 42,192


DuckDuckGo: 11,496


Bing: 80,722

URLFetch[""] 94,714


AOL: 164,072


Wow: 117,642


InfoSpace: 19,161


Info: 42,104


Blekko: 14,660


DogPile: 104,627


Alhea: 21,239


Thursday, June 25, 2015

Recursive Programming on Lists (Arrays): Map

Here is a simple recursive mapping function. This is a recursion on Lists (arrays) with arbitrary argument types, not numbers.

Recall that Map is a classic functional-style, 'structural' command since it leads us to think in terms of applying a function to an entire structure, such as an array (a List in Mathematica). The Map function saves us the trouble of 'deconstructing' the array with a Do or For loop as in procedural programming.

We define a base case and a recursive case:

1. Base case: If the source List is empty, just return it

2. Recursive case: Call itself on the Rest of each successively smaller sub-List


We test the base case:



We test some recursive cases:





How does it work? You can go through the Trace, but essentially by calling itself repeatedly on the Rest of list (via the first, recursive 'clause' of recursiveMap, recursiveMap[function,Rest@list]), it constructs successively smaller sub-Lists and eventually hits the empty List, which is the base case. 

Then it begins 'unwinding' and Prepend and the second clause of recursiveMap is repeatedly invoked (function@First@list). Starting with the last List of a single element (c in this Trace), it applies the function f to the First element of each sub-List and Prepends it to a growing List of results until the First element of the entire list is Prepended.



Source: David Wagner

Wednesday, June 24, 2015

Using Recursion to Define a Periodic Function

I'm working my way through Heikki Ruskeepaa's book, Mathematica Navigator and found this elegant one-liner that shows how to define a periodic function with recursion. Note that the base case required to stop the recursion is not a single value but the condition 0 <= t <= 2.

Clear@sawtooth;sawtooth[time_,scale_:1]:=If[0<=time<=2,scale time,sawtooth[scale (time-2)]]

How does the recursion work? If t => 2, sawtooth keeps subtracting 2 from t recursively until t is back in the range 0 <= t <= 2 and uses the computes that value for y. We can see sample values of a function plotted with DiscretePlot.


Of course a faster way of defining the domain is to use Mod, but there is a lesson here. We don't need the speed of Mod, so either implementation is fine - using Mod for the most compact and fast function, or using an elegant new recursive technique to learn it. This function omits the values at multiples of 2, but is effectively the same function:

Clear@sawtooth2;sawtooth2[time_,scale_:1]:=scale Mod[time,2]

Mathematica has built-in functions to produce various canonical waveforms  such as SawtoothWave. We are working with sawtooth and triangle waves in the lab since, interestingly, there is a lot of nonlinearity at the electrode-tissue interface and in the kHz+ frequency range, a rectangular wave is transformed by what is called partial half-wave rectification into something like a triangle offset in the negative (cathodic) direction from V = 0. We believe this is a key to understanding how biphasic neuromodulation works in the kHz+ frequency range.

Plot[SquareWave[{-50,50},t],{t,0,10},ExclusionsStyle->Dotted,AxesLabel->{"time (mS)","mV"},PlotLabel->"Electric potential at the electrode"]

Plot[TriangleWave[{-40,10},t],{t,0,10},AxesLabel->{"time (mS)","mV"},PlotLabel->"Electric potential seen by the axon"]