Lecture 38: Introduction to Julia
=================================

Julia is a programming language which aims to combine the high productivity
of scripting languages with the high performance of compiled languages.
The paper [BEKS17]_ describes the language.

.. _figjulialang.png:

.. figure:: ./figjulialang.png
    :align: center

    The logo of Julia, from the Software Engineering Daily web site.

Getting Started
---------------

In the cloud, at 
:index:`CoCalc`: https://cocalc.com
you can open a Jupyter notebook with a Julia kernel.

To install Julia on your own computer, follow the instructions
at <https://julialang.org>.  To execute Julia programs in a Jupyter notebook,
do ``using IJulia`` followed by ``notebook()`` in the REPL of Julia,
that is: when running in a Terminal or PowerShell window.
With ``using`` one can install ``Plots``, ``SymPy`` and ``Symbolics``.
The ``SymPy.jl`` is a wrapper to the Python package ``SymPy``,
whereas ``Symbolics`` [GMCGSER22]_ is a computer algebra system 
written in Julia.

To demonstrate symbolic computation with Julia,
consider the problem of computing the slope of the tangent line
at a point on a circle.

Implicit Differentiation with SymPy.jl
--------------------------------------

In the equation of the circle \ :math:`x^2 + y^2 - 1 = 0`,
the variable \ :math:`x` is the independent variable and
the variable \ :math:`y` depends on \ :math:`x`.
We will obtain the slope via implicit differentiation,
as \ :math:`\displaystyle \frac{d}{d~\!x} y(x)`.


::

   using SymPy

If ``SymPy`` is not yet installed, then you will be prompted to install.

The declarations of the variables result in ``dy`` to be assigned the
unknown slope \ :math:`\displaystyle \frac{d}{d~\!x} y(x)`.

::

   x = Sym("x")
   y = SymFunction("y")
   dy = diff(y(x), x)

Observe that ``y`` is declared as a symbolic function,
in *any* symbol.  If we want to use ``y`` when we define
the equation, then we *must* use ``y(x)`` to indicate that
the function ``y`` depends on ``x``.

::

   c = x^2 + y(x)^2 - 1
   dc = diff(c, x)
   dydx = solve(dc, dy)

The ``solve`` returns an ``Array`` which starts in Julia with index one.
The slope is then obtained as

::

   symbolicslope = dydx[1]

and equals ``-x/y(x)``.

Let us now take a random point on the circle, generated by a random angle,
using the polar representation.

::

   angle = 2*pi*rand()
   (a, b) = (cos(angle), sin(angle))
   sd = Dict(x => a, y(x) => b)

The coordinates of the point ``(a, b)`` are stored in a dictionary,
which is a convenient data structure for substitution.
Observe that, as Julia is geared towards numerical computing,
the ``pi`` is a 64-bit floating-point approximation for \ :math:`\pi`.
The numerical value for the slope is then obtained via substitution
as follows:

::

   slope = symbolicslope.subs(sd)

The numerical value in ``slope`` serves as a reference to check
the alternative computation of the slope via the Taylor series.

The ``series`` of ``SymPy`` works for one variable only,
but it works well on purely symbolic expressions, so we can compute
the expression for the tangent line by running ``series`` twice.
Because we continue the session, we declare new symbols ``X`` and ``Y``
for the variables in the circle equation and ``A`` and ``B`` for the
coordinates of the point.

::

   X, Y, A, B = Sym("X, Y, A, B")
   nc = X^2 + Y^2 - 1
   cx = SymPy.series(nc, X, A, 2)

In ``cx`` is the Taylor series expansion of the expression
stored in ``nc`` about ``X = A``, of order ``2``.
To convert the ``cx`` into a symbolic expression on which
we can apply ``series`` again, we select the coefficient
to make the expression ``px`` (a polynomial as truncated series).

::

   px = cx.coeff(X, 0) + cx.coeff(X - A, 1)*(X - A)
   cy = series(px, Y, B, 2)
   cY = cy.subs(Dict(A => a, B => b)

Then the slope is obtained as

::

   - cY.coeff(X, 1)/cY.coeff(Y, 1)

which agrees with the value earlier assigned to ``slope``.

Implicit Differentiation with Symbolics.jl
------------------------------------------

We will now compute the value of the slope a third time.

::

   using Symbolics.jl
   @variables x, y(x)

Observe that with ``Symbolics`` (opposite to what we did with ``SymPy``),
we declare ``y(x)`` so ``y`` is explicitly defined as depending on ``y``.
This implies that typing ``y`` suffices and, unlike with ``SymPy``,
we do not need to type ``y(x)`` when we define the equation of the circle.

In order not to confuse with the ``Differential`` of ``SymPy``,
which we used earlier in this notebook,
we declare that we use the ``Differential`` of ``Symbolics``.

::

   D = Symbolics.Differential(x)
   dy = D(y)

As before, ``dy``
is \ :math:`\displaystyle \frac{d}{d~\!x} y(x)`.

::

   c = x^2 + y^2 - 1
   dc = D(c)

The differential operator ``D`` does not expand automatically,
we have to ask for it.

::

   eq = expand_derivatives(dc)
   SBslope = Symbolics.solve_for(eq,dx)

In ``SBslope`` is the expression ``-x/y(x)``
so the symbolic part of this problem is thus indeed solved.
What is left is the substitution, again using a dictionary,
mapping the coordinates ``(a, b)`` of the random point
to the symbols ``x`` and ``y``.

::

   nsd = Dict(x => a, y => b)
   substitute(SBslope, nsd)

and we see again the value that was earlier assigned
to ``slope``.

.. rubric:: References


.. [BEKS17] J. Bezanson, A. Edelman, S. Karpinski, V.B. Shah.
   A fresh approach to numerical computing.
   *SIAM Review* 59(1): 65-98, 2017.

.. [GMCGSER22] S. Gowda, Y. Ma, A. Cheli, M. Gwozdz, V.B. Shah,
   A. Edelman, and C. Rackauckas.
   *High-performance symbolic-numerics via multiple dispatch*,
   arxiv:2105.03949v3 [cs.CL] 5 Feb 2022.