Showing posts with label piecewise. Show all posts
Showing posts with label piecewise. Show all posts

Thursday, February 4, 2016

Wolfram Language Attribute: Flat

Attributes modify the behavior of entire classes of functions. Some Attributes cause functions to have well-known mathematical properties such as commutativity, associativity, or idempotency, notably for valid pattern-matching. Others control the way the Wolfram Language evaluates or does not evaluate a function's arguments, which if you study that you will see is at the very core of Wolfram Language operation.

Listable is a valuable functional language Attribute since it incorporates Thread into a function so it automatically maps over a List and instead of taking a List apart, we apply a function to it as a gestalt.

NumericFunction is an Attribute of all Wolfram Language mathematical functions and is True when all of their arguments are numerical. I suspect it is integral to the interaction of many functions with the underlying knowledge base of Mathematica and therefore to many functions' valid behavior (e.g. Piecewise – always define piecewise mathematical functions with Piecewise).

Attribute: Flat


Flat is the mathematical property of associativity, which dictates that different groupings of numbers to specify the order in which the operation is performed, return the same result. E.g. ((a+b)+c) = (a+(b+c)) and (2x3)x4 = 2x(3x4).

a+(b+c)==(a+b)+c

True

Thomas Bahder uses a neat trick to show how Attributes work and how some affect pattern-matching. Here is the trick revealing the mechanism of Flat. Bahder uses a Condition that will never be True, which causes Wolfram Language to try all permutations of the function's arguments that are valid under Orderless. In the false predicate he inserts a Print statement to reveal how the Wolfram Language evaluator works. HoldForm is necessary to stop the evaluator from putting the arguments into standard order before printing them.

Flat tells the Wolfram Language evaluator, 'If the given pattern x_ + y_ matches any of the following patterns, return True.'

a+b+c/.x_+y_:>u/;Print@FullForm[Plus[x,y]//HoldForm]

HoldForm[Plus[a,Plus[b,c]]]
HoldForm[Plus[b,Plus[a,c]]]
HoldForm[Plus[c,Plus[a,b]]]
HoldForm[Plus[Plus[a,b],c]]
HoldForm[Plus[Plus[a,c],b]]
HoldForm[Plus[Plus[b,c],a]]
HoldForm[Plus[a,b]]
HoldForm[Plus[b,a]]
HoldForm[Plus[a,c]]
HoldForm[Plus[c,a]]
HoldForm[Plus[b,c]]
HoldForm[Plus[c,b]]

a+b+c

Here are related posts on Orderless and OneIdentity, and this post shows how to reveal all built-in functions with Attributes categorize by their Attribute.

This is an excellent out-of-print book but there is plenty to learn from it. I'd compare it to David Wagner's Power Programming in Mathematica in its insight to how the Wolfram Language works.

Bahder, Thomas B. Mathematica for Scientists and Engineers. Reading, MA: Addison-Wesley. 1994.

Saturday, June 20, 2015

Recursive Programming: The General Principle

Recursion, where a function calls itself, is a fundamental programming method. Recursion is a different way of iterating over a List than using a loop. Recurse is from the Latin recurrere, to run back (I picture a dog retrieving a stick over and over).

There are two principles of writing a recursive function, which is piecewise:
  1. Recursive case: Define a function f(n+1) recursively in terms of f(n), or equivalently, f(n) in terms of f(n-1)
  2. Base case: Define one or more base cases, which are stopping points for the recursion, so that it doesn't become infinite. 
An elementary example is exponentiation defined this way:

xn = x•x(n-1)

First we define the base case, so that Mathematica will match it to prevent an infinite recursive loop. Let's just consider positive integers and zero. The base case is any integer raised to the 0th power equals 1.

exponentiation[x_,0]:=1;

The recursive case is:

exponentiation[x_,n_Integer/;x>0]:=x exponentiation[x,n-1]

Start testing a recursive function with the simplest examples:

exponentiation[0,0]
1

exponentiation[1,0]
1

exponentiation[1,1]
1

Let's make sure the Condition requiring x > 0 works; it does.

exponentiation[1,-1]
exponentiation[1,-1]

Trace shows the exponentiation hitting the base case and stopping.

exponentiation[1,1]//Trace
{exponentiation[1,1],1 exponentiation[1,1-1],{{1-1,0},exponentiation[1,0],1},1 1,1}

To compute 22, the function calls itself to computer 22-1 i.e. 21, which in turn calls itself to compute 20, which is the base case, 1. Once it has recursed down to the base case, it has a value to substitute in the preceding computations and the recursion 'unwinds.' All recursive functions work this way.

exponentiation[2,2]//Trace
{exponentiation[2,2],2 exponentiation[2,2-1],{{2-1,1},exponentiation[2,1],2 exponentiation[2,1-1],{{1-1,0},exponentiation[2,0],1},2 1,2},2 2,4}

exponentiation[2,10]
1024

exponentiation[2,100]
1267650600228229401496703205376

exponentiation[2,1000]
10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376

Use Block for Large Recursions and Debugging


You may run into the built-in default recursion limit, $RecursionLimit, of 1024 when performing a large recursion.

In[584]:= exponentiation[2,10000]
During evaluation of In[584]:= $RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>
Out[584]= Hold[exponentiation[2,8978-1]]

The safe way to exceed the limit is using Block. Block temporarily changes the value of a Global or System variable. I've suppressed the output here since it is 3,011 digits long, but you can try this safely at home without "//IntegerDigits//Length".

Block[{$RecursionLimit=10050},exponentiation[2,10000]]//IntegerDigits//Length
3011

Mathematica will warn you and stops an infinite recursion with a default recursion limit specified by $RecursionLimit ($IterationLimit is used similarly to stop infinite loops).

{$RecursionLimit,$IterationLimit}
{1024,4096}

Just as you can use Block to temporarily exceed $RecursionLimit, you can use Block to temporarily 'tighten up' $RecursionLimit while writing a recursive function.

In[592]:= Block[{$RecursionLimit=20},factorial@x_:=x*factorial[x-1];factorial@2]
During evaluation of In[592]:= $RecursionLimit::reclim: Recursion depth of 20 exceeded. >>
Out[592]= Hold[factorial[-14-1]]

Always remember to write the base case rule!

In[593]:= Block[{$RecursionLimit=20},factorial@x_:=x*factorial[x-1];factorial@0=1;factorial@2]
Out[593]= 2

Now we may safely compute the factorial of larger numbers:

factorial@200

788657867364790503552363213932185062295135977687173263294742533244359449963403342920304284011984623904177212138919638830257642790242637105061926624952829931113462857270763317237396988943922445621451664240254033291864131227428294853277524242407573903240321257405579568660226031904170324062351700858796178922222789623703897374720000000000000000000000000000000000000000000000000

Note all the zeroes at the end. They must be from each time a multiple of 10 is incorporated into the factorial.

Tuesday, December 23, 2014

Easy Ways to Create Variations of Function

Easy Ways to Create Variations of Function

You will often need to create more than one version of a function. Here are two different methods, a piecewise function and using an Option. Either are fine, but in general simply for conciseness I use a piecewise function when the function is short (like less than a dozen lines) and an Option when the function is longer.

Here is the dataset for the examples. To see the 6 functions used to create arrays in Mathematica, see Ways to Create Arrays.

dataset=Array[List,{5,3}]

{{{1,1},{1,2},{1,3}},{{2,1},{2,2},{2,3}},{{3,1},{3,2},{3,3}},{{4,1},{4,2},{4,3}},{{5,1},{5,2},{5,3}}}

Here's the first version. It takes a List of data as a first argument and an Integer as its second argument. This simple example uses the index to pick out a Part from the data.

function1[data_List,index_Integer]:=data[[index]]

function1[dataset,4]

{{4,1},{4,2},{4,3}}

The Piecewise Approach

Now we need a second version to take a List of indices in the second argument. We can simply use datatyping by specifying the Head of the second argument. This is one way to do the "piecewise" method. The domain of possible inputs is split into pieces (sub-domains) and each variation of a function is designed to pick out the correct piece (sub-domain) for which it is designed.

function2[data_List,index_List]:=data[[index]]

function2[dataset,{2,4}]

{{{2,1},{2,2},{2,3}},{{4,1},{4,2},{4,3}}}

The Options Approach

You can see that if there were a minor change in a long function, copying the function is not as concise as the following approach using Options. There is more overhead in this example to create the Option and handle it in the function, but in a long function that overhead is less than copying as in the Piecewise approach.

Note that for concise creation of Piecewise mathematical functions, the built-in function Piecewise should be used. For an explanation of how I use Options, see A Template for Optional Arguments [to be published].

ClearAll@function3;
Options@function3="indexOption"->"Integer";

function3[data_List,index_,options___?OptionQ]:=
Module[{indexOption="indexOption"/.{options}/.Options@function3},

If[indexOption=="Integer",dataset[[index]]  ];
If[indexOption=="List",dataset[[{index}]]  ];

(*endModule*)]

Here is the function that prompted me to write this primer. It takes the membrane voltage time series from simulated spiking neuron cells and counts how many spikes occur, with cell range and time series range as arguments. I needed to expand it so I could specify counting spikes in several time series ranges. I didn't want to copy and vary it with the piecewise approach, so I used the options approach. The changes to the original function to implement the new ones are italicized.

ClearAll@batchMeanFiringRateTable;
Options@batchMeanFiringRateTable={"Export"->Off,"MultipleBatch"->Off,"MultipleTimeSeries"->False};

batchMeanFiringRateTable[spikeFileDirectory_String,cellRange_List,timeSeriesRange_List,leftColumnHeading_String:"Neural Group",rightColumnHeading_String:"Average Firing Rate",options___?OptionQ]:=
Module[{exportFileName,firingData,
spikeFiles=getSpikeFiles@spikeFileDirectory,

exportOption="Export"/.{options}/.Options@batchMeanFiringRateTable,
multipleBatchOption="MultipleBatch"/.{options}/.Options@batchMeanFiringRateTable,
multipleTimeSeriesOption="MultipleTimeSeries"/.{options}/.Options@batchMeanFiringRateTable},

(*With "Batch"\[Rule]On, meanFiringRateFromSpikeFile will just return the averageSpikeRate each time it's called*)

(*Need FileBaseName to identify cell group in left column if TableForm; need First@timeSeriesRange to remove the outer List for a single time series range*)

If[multipleTimeSeriesOption==False,firingData=Table[{FileBaseName@batchFile,meanFiringRateFromSpikeFile[batchFile,cellRange,First@timeSeriesRange,"Batch"->On]},{batchFile,spikeFiles}]  (*endIf*)];

If [multipleTimeSeriesOption==True,firingData=Table[{FileBaseName@batchFile,meanFiringRateFromSpikeFile[batchFile,cellRange,timeSeries,"Batch"->On]},{batchFile,spikeFiles},{timeSeries,timeSeriesRange}]/.{{x_,y_},{x_,z_}}->{x,y,z};(*Table will take each spike file and iterate over timeSeriesRange. Don't need First@timeSeriesRange here. *) (*endIf*)];

If[multipleBatchOption==Off,Print@"Multiple batch is off.";Print[firingData//TableForm[#,TableAlignments->Left,TableHeadings->{None,{leftColumnHeading,"Ave AP Rate: "<>ToString@timeSeriesRange}}]&],Return@firingData (*endIf*)];

If[exportOption==On&&multipleBatchOption==Off,exportFileName=FileNameJoin@{DirectoryName@spikeFiles[[1]],ToString@FileNameTake[DirectoryName@spikeFiles[[1]],-1]<>"-BatchFileTable.xlsx"};
Export[exportFileName,firingData];Print@StringForm["File exported to ``.",exportFileName]];
]


batchMeanFiringRateTable["C:\\Users\\Public\\Documents\\UNCuS16_09_2013\\DataFiles\\WDR-Abeta-14-12only\\WDR-Abeta-14-12only-CS8-9-20-RS0p28\\",{1,80},{{1,45},{46,85}},"Neural Group","Method"->"OverSpikingCells","MultipleBatch"->Off,"Export"->Off,"MultipleTimeSeries"->True]