Tag Archives: #MonthOfJulia

#MonthOfJulia Day 20: Calculus

Julia-Logo-Calculus

Mathematica is the de facto standard for symbolic differentiation and integration. But many other languages also have great facilities for Calculus. For example, R has the deriv() function in the base stats package as well as the numDeriv, Deriv and Ryacas packages. Python has NumPy and SymPy.

Let’s check out what Julia has to offer.

Numerical Differentiation

First load the Calculus package.

julia> using Calculus

The derivative() function will evaluate the numerical derivative at a specific point.

julia> derivative(x -> sin(x), pi)
-0.9999999999441258
julia> derivative(sin, pi, :central)       # Options: :forward, :central or :complex
-0.9999999999441258

There’s also a prime notation which will do the same thing (but neatly handle higher order derivatives).

julia> f(x) = sin(x);
julia> f'(0.0)                             # cos(x)
0.9999999999938886
julia> f''(0.0)                            # -sin(x)
0.0
julia> f'''(0.0)                           # -cos(x)
-0.9999977482682358

There are functions for second derivatives, gradients (for multivariate functions) and Hessian matrices too. Related packages for derivatives are ForwardDiff and ReverseDiffSource.

Symbolic Differentiation

Symbolic differentiation works for univariate and multivariate functions expressed as strings.

julia> differentiate("sin(x)", :x)
:(cos(x))
julia> differentiate("sin(x) + exp(-y)", [:x, :y])
2-element Array{Any,1}:
 :(cos(x))    
 :(-(exp(-y)))

It also works for expressions.

julia> differentiate(:(x^2 * y * exp(-x)), :x)
:((2x) * y * exp(-x) + x^2 * y * -(exp(-x)))
julia> differentiate(:(sin(x) / x), :x)
:((cos(x) * x - sin(x)) / x^2)

Have a look at the JuliaDiff project which is aggregating resources for differentiation in Julia.

Numerical Integration

Numerical integration is a snap.

julia> integrate(x -> 1 / (1 - x), -1 , 0)
0.6931471805602638

Compare that with the analytical result. Nice.

julia> diff(map(x -> - log(1 - x), [-1, 0]))
1-element Array{Float64,1}:
 0.693147

By default the integral is evaluated using Simpson’s Rule. However, we can also use Monte Carlo integration.

julia> integrate(x -> 1 / (1 - x), -1 , 0, :monte_carlo)
0.6930203819567551

Sympy-logo

Symbolic Integration

There is also an interface to the Sympy Python library for symbolic computation. Documentation can be found here. You might want to restart your Julia session before loading the SymPy package.

julia> using Sympy

Revisiting the same definite integral from above we find that we now have an analytical expression as the result.

julia> integrate(1 / (1 - x), (x, -1, 0))
log(2)
julia> convert(Float64, ans)
0.6931471805599453

To perform symbolic integration we need to first define a symbolic object using Sym().

julia> x = Sym("x");                       # Creating a "symbolic object"
julia> typeof(x)
Sym (constructor with 6 methods)
julia> sin(x) |> typeof                    # f(symbolic object) is also a symbolic object
Sym (constructor with 6 methods)

There’s more to be said about symbolic objects (they are the basis of pretty much everything in SymPy), but we are just going to jump ahead to constructing a function and integrating it.

julia> f(x) = cos(x) - sin(x) * cos(x);
julia> integrate(f(x), x)
     2            
  sin (x)         
- ─────── + sin(x)
     2     

What about an integral with constant parameters? No problem.

julia> k = Sym("k");
julia> integrate(1 / (x + k), x)
log(k + x)

We have really only grazed the surface of SymPy. The capabilities of this package are deep and broad. Seriously worthwhile checking out the documentation if you are interested in symbolic computation.

newton_and_leibniz

I’m not ready to throw away my dated version of Mathematica just yet, but I’ll definitely be using this functionality often. Come back tomorrow when I’ll take a look at solving differential equations with Julia.

The post #MonthOfJulia Day 20: Calculus appeared first on Exegetic Analytics.

#MonthOfJulia Day 19: Units

Julia-Logo-SIUnits

The packages we’ll be looking at today should bring joy to the hearts of all Physical Scientists. Actually they should make any flavour of Scientist happy.

It is natural for man to relate the units of distance by which he travels to the dimensions of the globe that he inhabits. Thus, in moving about the earth, he may know by the simple denomination of distance its proportion to the whole circuit of the earth. This has the further advantage of making nautical and celestial measurements correspond. The navigator often needs to determine, one from the other, the distance he has traversed from the celestial arc lying between the zeniths at his point of departure and at his destination. It is important, therefore, that one of these magnitudes should be the expression of the other, with no difference except in the units. But to that end, the fundamental linear unit must be an aliquot part of the terrestrial meridian. … Thus, the choice of the metre was reduced to that of the unity of angles.
Pierre-Simon Laplace

SIUnits

The SIUnits package provides unit-checked operations for quantities expressed in SI units.

julia> using SIUnits
julia> using SIUnits.ShortUnits

It supports both long and short forms of units and all the expected arithmetic operations.

julia> 1KiloGram + 2kg
3 kg 
julia> 4Meter - 2m
2 m
julia> 4m / 2s
2.0 m s⁻¹

Note that it only recognises the American spelling of “meter” and not the (IMHO correct) “metre”! But this is a small matter. And I don’t want to engage in any religious wars.

Speaking of small matters, it’s possible to define new units of measure. Below we’ll define the micron and Angstrom along with their conversion functions.

julia> import Base.convert
julia> Micron = SIUnits.NonSIUnit{typeof(Meter),:µm}()
µm
julia> convert(::Type{SIUnits.SIQuantity},::typeof(Micron)) = Micro*Meter
convert (generic function with 461 methods)
julia> Angstrom = SIUnits.NonSIUnit{typeof(Meter),:Å}()
Å
julia> convert(::Type{SIUnits.SIQuantity},::typeof(Angstrom)) = Nano/10*Meter
convert (generic function with 462 methods)

And now we can freely use these new units in computations.

julia> 5Micron
5 µm
julia> 1Micron + 1m
1000001//1000000 m
julia> 5200Angstrom                        # Green light
5200 Å

Physical

The Physical package is documented here. Apparently it’s not as performant as SIUnits but it does appear to have a wider scope of functionality. We’ll use it to address an issue raised on Day 17: converting between Imperial and Metric units.

Let’s kick off by loading the package.

using Physical

There’s a lot of functionality available, but we are going to focus on just one thing: converting pounds and inches into kilograms and metres. First we define a pair of derived units. To do this, of course, we need to know the appropriate conversion factors!

julia> Inch = DerivedUnit("in", 0.0254*Meter)
1 in 
julia> Pound = DerivedUnit("lb", 0.45359237*Kilogram)
1 lb 

We can then freely change the average heights and weights that we saw earlier from Imperial to Metric units.

julia> asbase(66Inch)
1.6764 m 
julia> asbase(139Pound)
63.04933943 kg 

On a related note I’ve just put together a package of physical constants for Julia.

julia> using PhysicalConstants
julia> PhysicalConstants.MKS.SpeedOfLight
2.99792458e8
julia> PhysicalConstants.MKS.Teaspoon
4.92892159375e-6

Did you know that a teaspoon was 4.92892 millilitres? There I was, wallowing in my ignorance, thinking that it was 5 millilitres. Pfffft. Silly me. There are

teaspoon-volume

Units can be a contentious issue. Watch the video below to see what Richard Feynman had to say about the profusion of units used by Physicists to measure energy. Also check out the full code for today along with the index to the entire series of #MonthOfJulia posts on github.

For those who want some proof that physicists are human, the proof is in the idiocy of all the different units which they use for measuring energy.
Richard P. Feynman

The post #MonthOfJulia Day 19: Units appeared first on Exegetic Analytics.

#MonthOfJulia Day 18: Plotting

Julia-Logo-Plotting

There’s a variety of options for plotting in Julia. We’ll focus on those provided by Gadfly, Bokeh and Plotly and.

Gadfly

Gadfly is the flavour of the month for plotting in Julia. It’s based on the Grammar of Graphics, so users of ggplot2 should find it familiar.

gadfly-logo

To start using Gadfly we’ll first need to load the package. To enable generation of PNG, PS, and PDF output we’ll also want the Cairo package.

julia> using Gadfly
julia> using Cairo

You can easily generate plots from data vectors or functions.

julia> plot(x = 1:100, y = cumsum(rand(100) - 0.5), Geom.point, Geom.smooth)
julia> plot(x -> x^3 - 9x, -5, 5)

Gadfly plots are by default rendered onto a new tab in your browser. These plots are mildly interactive: you can zoom and pan across the plot area. You can also save plots directly to files of various formats.

julia> dampedsin = plot([x -> sin(x) / x], 0, 50)
julia> draw(PNG("damped-sin.png", 800px, 400px), dampedsin)

damped-sin

Let’s load up some data from the nlschools dataset in R’s MASS package and look at the relationship between language score test and IQ for pupils broken down according to whether or not they are in a mixed-grade class.

julia> using RDatasets
julia> plot(dataset("MASS", "nlschools"), x="IQ", y="Lang", color="COMB",
            Geom.point, Geom.smooth(method=:lm), Guide.colorkey("Multi-Grade"))

nlschools

Those two examples just scratched the surface. Gadfly can produce histograms, boxplots, ribbon plots, contours and violin plots. There’s detailed documentation with numerous examples on the homepage.

Watch the video below (Daniel Jones at JuliaCon 2014) then read on about Bokeh and Plotly.

Bokeh

Bokeh is a visualisation library for Python. Bokeh, like D3, renders plots as Javascript, which is viewable in a web browser. In addition to the examples on the library homepage, more can be found on the homepage for Julia’s Bokeh package.

The first thing you’ll need to do is install the Bokeh library. If you already have a working Python installation then this is easily done from the command line:

$ pip install bokeh

Next load up the package and generate a simple plot.

julia> using Bokeh
julia> autoopen(true);
julia> x = linspace(0, pi);
julia> y = cos(2 * x);
julia> plot(x, y, title = "Cosine")
Plot("Cosine" with 1 datacolumns)

The plot will be written to a file bokeh_plot.html in the working directory, which will in turn be opened by the browser. Use plotfile() to change the name of the file. The plot is interactive, with functionality to pan and zoom as well as resize the plot window.

bokeh-plot

Plotly

The Plotly package provides a complete interface to plot.ly, an online plotting service with interfaces for Python, R, MATLAB and now Julia. To get an idea of what’s possible with plot.ly, check out their feed. The first step towards making your own awesomeness with be loading the package.

using Plotly

Next you should set up your plot.ly credentials using Plotly.set_credentials_file(). You only need to do this once because the values will be cached.

Data series are stored in Julia dictionaries.

julia> p1 = ["x" => 1:10, "y" => rand(0:20, 10), "type" => "scatter", "mode" => "markers"];
julia> p2 = ["x" => 1:10, "y" => rand(0:20, 10), "type" => "scatter", "mode" => "lines"];
julia> p3 = ["x" => 1:10, "y" => rand(0:20, 10), "type" => "scatter", "mode" => "lines+markers"];
julia> Plotly.plot([p1, p2, p3], ["filename" => "basic-line", "fileopt" => "overwrite"])
Dict{String,Any} with 5 entries:
  "error"    => ""
  "message"  => ""
  "warning"  => ""
  "filename" => "basic-line"
  "url"      => "https://plot.ly/~collierab/17"

You can either open the URL provided in the result dictionary or do it programatically:

julia> Plotly.openurl(ans["url"])

plotly-scatter

By making small jumps through similar hoops it’s possible to create some rather intricate visualisations like the 3D scatter plot below. For details of how that was done, check out my code on github.
plotly-3d-scatter

Obviously plotting and visualisation in Julia are hot topics. Other plotting packages worth checking out are PyPlot, Winston and Gaston. Come back tomorrow when we’ll take a look at using physical units in Julia.




The post #MonthOfJulia Day 18: Plotting appeared first on Exegetic Analytics.