~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|  Cornell SSW: Mathematica == Fun(ctions).  |
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Monday March 27, 2017

Author: Varun Gandhi (vg266) - Physics first year grad student

1. What can Mathematica do?

Obligatory feature list:
* Functional programming - first-class functions, lambdas, map, select etc.
* Symbolic manipulation - solving differential equations, nasty integrals etc.
* Rich visualization tools - making pretty graphs is relatively painless (glares at TikZ/PGF)
* Batteries + supercharger included - lots and lots of performant in-built functions
* Good documentation - easy to access with many examples
* Rapid prototyping - consequence of abovementioned features
* Getting help - friendly StackExchange community
and many more (depending on your field). 

(Benchmarks: http://julialang.org/benchmarks/)

Example 1: You can make a fake alien language for a popular science fiction movie 
(which is unfortunately unlike that in the original short story).

Example 2: Lissajous figures (from the Documentation on `Manipulate`).

Manipulate[
ParametricPlot[
{a1 Sin[n1 (x+p1)],a2 Cos[n2 (x+p2)]},{x,0,20Pi},
PlotRange->1,PerformanceGoal->"Quality"],
Style["Vertical",Bold,Medium],
{{n1,1,"Frequency"},1,4},{{a1,1,"Amplitude"},0,1},
{{p1,0,"Phase"},0,2Pi},
Delimiter,Style["Horizontal",Bold,Medium],
{{n2,5/4,"Frequency"},1,4},{{a2,1,"Amplitude"},0,1},
{{p2,0,"Phase"},0,2Pi}, 
ControlPlacement->Left]

Example 3: Mathematica knows Markov processes. MouseOver

2. How do you do X in Mathematica?

Basics: Functions

Mathematica is (mostly) dynamic. You do not need to "declare" variables or bother with types (Mathematica does not have types)1.
Variables are:
* typeless and mutable
* actually just functions with zero arguments
* global by default, so you should avoid mutation unless absolutely necessary
* saved in memory until you quit Mathematica or explicitly cleared

x = 0 (* x should be blue before evaluation as it is unused *)
x = 1; (* ';' can be used to avoid printing *)

x (* x is black because it has been saved. This should print '1' *)
Clear[x];
x (* x was cleared so it will just print x *)

Functions are easy to declare using pattern matching.
Square brackets ‘[ ]’ are reserved for function application.

x = 10; f[x_] := 2^x;
f2[x_,y_]:=(a*x + b*y);
f[2]
f2[1,2]

Two common newbie mistakes are forgetting the ‘:’ or the ‘_’ symbols.

(* Intuitive but incorrect ways *)
x = 5;
g[x] := 2^x; h[x_] = 2^x;
x = 10;
g[x] (* 1024 ? *)
g[5] (* 32 ? *)
g[3] (* 8 ? *)
h[3] (* 8 ? *)

We can declare anonymous functions (== lambdas) using ‘#’ (slot) and ‘&’. 
These are called “pure functions” in Mathematica.
(These can have side-effects, unlike pure functions in other functional languages)

(# + 5/7)&     [2/7] (* CORRECT *)
(# + 5/7)& @ (2/7) (* Prefix notation for unary function application *)
(# + 5/7)&     (2/7) (* INCORRECT. Only square brackets are allowed *)
(# + 5/7)& @   2/7   (* INCORRECT because of operator precedence *)

You can make pure functions by literally calling `Function`.
(The '&' earlier is simply syntactic sugar for `Function`.)

Function[#+5/7][arg]
Function[x,x+5/7][arg]

Pure functions can also be declared with multiple arguments.

(a*#1+b*#2)& [1/4,3/5]
Function[{x,y},a*x+b*y]&[1/4,3/5]

Mathematical equations are relatively easy to write and easy to read. 
You cannot use '_' in function names as it is reserved for pattern matching.

blackBodySpectrum[λ_,T_]:= ((4 π ℏ c^2)/(λ^5))(1/(E^[((2 π ℏ c)/(λ (Subscript[k, B]) T))]-1));
(* Look Ma', no explicit multiplication symbols! *)
ℏ = 6.626 10^-34/(2 π);
c = 3 10^8;
(Subscript[k, B])= 1.38 10^-23;
Plot[blackBodySpectrum[λ,6000],{λ,0,2 10^-6}]

You can use pattern matching to declare functions in a piece-wise manner.

f[2] := 2;
f[x_Integer]:=3;
f[x_String]:=0.5;
f[1] := 1;
{f[1],f[2],f[3],f["123"],f[Cow]}

You can declare multi-expression functions using ‘;’ (`CompoundExpression`). 
The function’s return value is the value of the expression after the last ‘;’.

doSomeSideEffects[x_]:= Print["Hi."];
veryHugeFunction[x_]:=(
doSomeSideEffects[x];
lessHugeFunction[x]
);
veryHugeFunction[z]

Recall that variables and functions are global by default.
You can restrict the scope of temporary variables using `Module`
instead of using `Clear` all the time.

foo1[x_]:=(yy=bar[x]; sideEffects[yy];baz[x]);
foo2[x_]:=Module[{zz},
zz=bar[x];sideEffects[zz];baz[x]];
{yy,zz}
{foo1[arg],foo2[arg]}
{yy,zz}
Clear[yy]

Basics: Lists and associations

Lists are the primary data structure in Mathematica.
Mathematica 10 introduces Associations, which are key-value pairs (dictionaries).
A lot of mathematical functions are ‘Listable’ i.e. they operate on lists automagically.

list1 = {15/2,45/2};
list2 = 15+list1
list3 = (list1 + list2) Degree
Sin[list3]

`Apply` splits a list into elements and feeds them as arguments to a function.

(* Function syntax *)
Apply[myFunc, {a,b,c,d}]
(* Operator syntax *)
myFunc @@ {a,b,c,d}
(* Addition is also a function *)
Plus @@ {a,b,c,d}

Some other useful functions to operate on lists are `Map`, `Thread`, `MapThread`, and `Transpose`.

Clear[x,y,z];
(* Function syntax *)
Map[#^2 Sin[#]&,{x,y,z}]
(* Operator syntax *)
#^2 Sin[#]& /@{x,y,z}

Clear[x1,x2,x3,x4,y1,y2,y3,y4];
Thread[h[{x,y},{a,b}]]
MapThread[h[#1,#2]&,{{x,y},{a,b}}]
Transpose[{{x1,x2,x3,x4}, {y1,y2,y3,y4}}]

(* Format matrix as a table *)
Transpose[{{x1,x2,x3,x4}, {y1,y2,y3,y4}}] // TableForm

// is just postfix notation for unary function application, like @ was prefix notation for the same.

(* Uglier but equivalent *)
TableForm[Transpose[{{x1,x2,x3,x4}, {y1,y2,y3,y4}}]]
(* Also equivalent
TableForm @ Transpose @ {...}
{...} // Transpose // TableForm
(Tranpose @ {...}) // TableForm 
TableForm @ ({...} // Transpose)
and other combinations with [ ]
*)

List indices start at 1. Indexing is done using double brackets '[[ ]]' or 'Part[list, ind]'.

list3 = {3,5,7,9,11,13};
list3[[0]]
list3[[1]]
list3[[1;; 3]]
list3[[-1]]

Since the brackets start increasing quickly with indexing and function application,
it is a good idea to avoid direct indexing for readability.

Associations are key-value pairs. Many functions that operate on lists also work on Associations.

Clear[a,b,c,x,y,z];
assoc = <|a->x,b->y,c->z|>;
assoc[b]
2*# & /@ assoc

Interlude: `Manipulate`, `Animate`, `Plot3D`

`Manipulate` is useful for making interactive plots with sliders and other elements.

Vpoint[q_,qpos_,r_]:=q/Norm[r-qpos];
Manipulate[
ContourPlot[
Vpoint[q1,p[[1]],{x,y}]+Vpoint[q2,p[[2]],{x,y}],
{x,-2,2},{y,-2,2},
Contours->10],
{{q1,-1},-3,3},{{q2,1},-3,3},
{{p,{{-1,0},{1,0}}},{-1,-1},{1,1},Locator},
Deployed->True]

It is also relatively easy to make animations.

Animate[
Plot[{Sin[x-t],Sin[x+t],Sin[x-t]+Sin[x+t]},
{x,0,2π}, PlotRange->{-2,2},PlotLegends->Automatic,
MaxRecursion->0],
{t,0,2π}]

3D plots are easy to make and can be rotated directly without special commands.

Plot3D[{Re[Log[x+I y]],Im[Log[x+I y]]},{x,-2,2},{y,-2,2}, 
ExclusionsStyle->Opacity[0.3],Mesh->None,PlotLegends->Automatic]

Basics: Getting help

Mathematica has a built-in documentation centre accessible using F1 or under the Help Menu.

Say I have code which I copy pasted below which I don’t fully understand.

Show[
Plot[
Sin[x],{x,-3 Pi,3 Pi},
PlotStyle->Red,MaxRecursion->0,
PlotPoints->{None,{-7.3,-5.2,-Pi,Pi,3.8,5}},
PlotRange->{{-3 Pi,3 Pi},{-1,1}}
],
Plot[Sin[x],{x,-5 Pi,5 Pi},PlotStyle->Dashed]
]

There is also a dedicated Mathematica StackExchange which will almost always have an answer
for all your basic questions (use the search function!). People are usually happy to help and fairly polite.

MUST-READ: https://mathematica.stackexchange.com/questions/18393/what-are-the-most-common-pitfalls-awaiting-new-users

Basics: Control flow

Control-flow keywords like ‘if’, ‘do’, ‘for’ etc. are just ordinary functions in Mathematica.
Therefore, unlike many other languages, `If` actually returns a value.
(Recall: All functions return values)

fizzBuzzNum[x_]:=
If[Mod[x,15]==0,"FizzBuzz",
If[Mod[x,3]==0,"Fizz",
If[Mod[x,5]==0,"Buzz",""]]]
Map[fizzBuzzNum,Range[15]]
Array[fizzBuzzNum,15]

Despite that, in most cases you want to avoid writing a loop structure.
If you are writing a `For` loop as a beginner, you are almost surely doing something wrong!

Example 1: Sum all triangular numbers (Tn) below a given M. One definition is Tn = n(n+1)/2.

isTriangular[k_Integer]:=IntegerQ[1/2(-1+√(1+8k))];
(* Naive, functional implementation *)
triangularNumSum1[M_Integer] := Plus@@Array[If[isTriangular[#],#,0] &, M];
triangularNumSum2[M_Integer]:=Plus@@Select[Range[M],isTriangular];

(* Test that the two implementations are equivalent *)
 triangularNumSum1[#]==triangularNumSum2[#]&@ RandomInteger[{1,10^5}]

The implementation is bad because it wastes time doing checks for all numbers. 
If we were writing imperative code, we could change the logic to:
* Start a `For` loop with n=1 (iteration variable), Σ=0
* Compute Tn = n(n+1)/2. 
* If Tn ≤ M, then increment Σ by Tn, else break.

Can this be written in a functional style?

Tn[n_]:=n(n+1)/2;
(* Improved functional implementation *)
triangularNumSum3[M_Integer]:=Plus@@
TakeWhile[Array[Tn,Ceiling[-1+√(2M)]],#≤M&];

a = RandomInteger[{1,10^6}];
TableForm[AbsoluteTiming @#[a]&/@ 
{triangularNumSum1,triangularNumSum2,triangularNumSum3}]
Clear[a]

Example 2: Can we convert the following pseudo-code into a more functional style?

for (x,y) in XYpairs:
    y = foo(x)
    total += bar(x,y)

How can we get rid of the ‘for’ loop? Wouldn’t using `Map` force us to iterate through the list twice?

(*
foo[x_]:=...;
bar[x_,y_]:=...;
*)

(* Decent option. *)
δTotal[XYPairs_]:=Plus@@ (
Function[{x,y},(y=foo[x]; bar[x,y])] @@ #&
/@XYPairs);

XYPairs = {{x1,y1},{x2,y2},{x3,y3}};
newTotal = oldTotal + δTotal[XYPairs]

XYPairs was mutated inside δTotal, which may not be immediately obvious
while refactoring, especially for bigger, intimidating functions.

Instead, try to avoid mutating arguments as far as possible.
Mathematica can memoize function values if the function has no side-effects.
Most library functions, even those named as `Append`, avoid mutating arguments.
(there is an additional function `AppendTo` as an option.)

(* Better option. *)
XYPairs = {{x1,y1},{x2,y2},{x3,y3}};
fooBar[x_]:={#,bar@@#}&@{x,foo[x]};

{newXYPairs, individualδs} = 
Transpose @ (fooBar@#[[1]]&/@XYPairs);

newXYPairs
newTotal = oldTotal+Plus @@ individualδs

Learn to use functions like `Map`, `Array`, `MapThread`, `Transpose`, `Apply`, `Select` (== filter) etc. 
that operate on entire lists instead of writing loops.

Bonus: You can call WolframAlpha from Mathematica

WolframAlpha["Cornell University"]

1	https://stackoverflow.com/questions/4558732/is-mathematica-an-untyped-language