{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3319aa2f",
   "metadata": {},
   "source": [
    "# Lab 3: Introduction to Julia\n",
    "\n",
    "There are traditionally two categories of computer languages: Compiled (C/C++, Fortran) and Interpreted (Python, R).\n",
    "\n",
    "Julia language has inherited the advantage from both categories.\n",
    "\n",
    "Some helpful links:\n",
    " * Julia documentation (https://julialang.org/)\n",
    " * Julia analysis demo at PyHEP 2023 (https://github.com/Moelf/PyHEP_2023_JuliaDemo/blob/main/Julia_analysis_demo.ipynb)\n",
    " * Talk (https://jiling.web.cern.ch/jiling/dump/2021_Harvard_JuliaHEP.html)\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c3cca61",
   "metadata": {},
   "source": [
    "## Motivation\n",
    "\n",
    "![](julia_comparison.png)\n",
    "Ref: https://arxiv.org/abs/2109.09973\n",
    "\n",
    "Julia is a solution to the \"two-language problem,\" e.g. using writing performant C++ code and a high-level Python wrapper.\n",
    "\n",
    "Consider this function taken from Numba's introduction. Numba is great because:\n",
    "- you don't write C++\n",
    "- you can write fast loop!\n",
    "\n",
    "#### Python (Numba)\n",
    "```python\n",
    "@jit(nopython=True)\n",
    "def go_fast(a):\n",
    "    trace = 0.0\n",
    "    for i in range(a.shape[0]):\n",
    "        trace += np.tanh(a[i, i])\n",
    "    return a + trace\n",
    "\n",
    "In [2]: x = np.arange(100).reshape(10, 10)\n",
    "In [3]: %timeit go_fast(x)\n",
    "814 ns ± 5.2 ns (mean ± σ)    # Numba\n",
    "16500 ns ± 70.5 ns (mean ± σ) # CPython\n",
    "```\n",
    "\n",
    "but it does have draw back (all JIT x Python solutions do): it's not compatible with all Python.\n",
    "\n",
    "the Julia \"sale\" is that:\n",
    "- you still avoid C++\n",
    "- you can write fast anything!*\n",
    "- there's no compatibility issue and \"unhandled\" construct, it's just the language itself\n",
    "\n",
    "#### Julia\n",
    "```julia\n",
    "function go_faster(a)\n",
    "    trace = 0.0\n",
    "    for i in axes(a, 1)\n",
    "        trace += tanh(a[i, i])\n",
    "    end\n",
    "    return a .+ trace\n",
    "end\n",
    "julia> x = reshape(0:99, 10, 10)\n",
    "julia> @benchmark go_faster($x)\n",
    "158.947 ns ± 100.451 ns (mean ± σ) # σ due to GC\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d74b149f",
   "metadata": {},
   "source": [
    "## Array basics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54bf6db0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A function makerandom with two methods is declared.\n",
    "# This looks just like two functions with the same name but with different arguments.\n",
    "# An array r with n Float64 elements is declared; element i is referred to as r[i].\n",
    "# The Base function round(T,m) rounds the number m to the closest integer and converts to type T (here Int)\n",
    "# The Base function rand() generates (when given no arguments) a random Float64 in [0,1).\n",
    "# The functional difference between the methods is Float64 (first) vs Float32 (second) numbers returned.\n",
    "# In the second method r[i]=rand() implies a loss of precision from 64 to 32 bit Float.\n",
    "\n",
    " function makerandom(n::Int)\n",
    "    r=Array{Float64}(undef,n)\n",
    "    for i=1:n\n",
    "       r[i]=rand()\n",
    "    end\n",
    "    return r\n",
    " end\n",
    "\n",
    " function makerandom(m::Float64)\n",
    "    n=round(Int,m)\n",
    "    r=Array{Float32}(undef,n)\n",
    "    for i=1:n\n",
    "       r[i]=rand()\n",
    "    end\n",
    "    return r\n",
    " end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59df592d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The function is called twice, with different arguments (example of multiple dispatch).\n",
    "# In each case 5 times; n is assigned an integer-constant value, type will be Int (Int64).\n",
    "# The Base function convert(Float64,n) converts the integer n to a Float64.\n",
    "\n",
    " n=5\n",
    " m=convert(Float64,n)\n",
    "\n",
    "# The function is called with Int argument (first method dispatched)\n",
    "# a will be the 5-element Float64 array returned by the function\n",
    "\n",
    " a=makerandom(n)\n",
    " for i=1:n\n",
    "    println(i,\"  \",a[i])\n",
    " end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1372f69f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The function is called with Float64 argument (second method dispatched)\n",
    "# b will be the 5-element Float32 array returned by the function\n",
    "\n",
    " a=makerandom(m)\n",
    " for i=1:n\n",
    "    println(i,\"  \",a[i])\n",
    " end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6212a221",
   "metadata": {},
   "source": [
    "## Matrix multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87f0cb54",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function creates an n*n random matrix (similar to randomarray.jl but 2-dim array)\n",
    "\n",
    " function randmatrix(n::Int)\n",
    "    mat=Array{Float64}(undef,n,n)\n",
    "    for j=1:n\n",
    "       for i=1:n\n",
    "          mat[i,j]=rand()\n",
    "       end\n",
    "    end\n",
    "    return mat\n",
    " end\n",
    "\n",
    "# Ask for matrix size and read it (stdin = standard input) using readline\n",
    "# Text input has to be parsed to a type (here Int)\n",
    "\n",
    " println(\"Give matrix size\")\n",
    " n=3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5a4a217",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create two random matrices and their product\n",
    "# Note that \"*\" is actual matrix multiplication\n",
    "\n",
    " a=randmatrix(n)\n",
    " b=randmatrix(n)\n",
    " c=a*b\n",
    "\n",
    "# print the results\n",
    "\n",
    " println()\n",
    " println(\"A, B, A*B\")\n",
    " for i=1:n\n",
    "    println(a[i,:],\"  \",b[i,:],\"  \",c[i,:])\n",
    " end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fff06d8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Element-by-element operations can be done with \".\" added to operator\n",
    "# Here element-by-element multiplication\n",
    "\n",
    " c=a.*b\n",
    " println()\n",
    " println(\"A, B, A elements multiplied by B elements\")\n",
    " for i=1:n\n",
    "    println(a[i,:],\"  \",b[i,:],\"  \",c[i,:])\n",
    " end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6be65ee3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Matrix inverse using Base function inv()\n",
    "# Multiply inverse of A by A to check if OK\n",
    "\n",
    " b=inv(a)\n",
    " c=a*b\n",
    " println()\n",
    " println(\"A, 1/A, A*(1/A)\")\n",
    " for i=1:n\n",
    "    println(a[i,:],\"  \",b[i,:],\"  \",c[i,:])\n",
    " end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14e12e5d",
   "metadata": {},
   "source": [
    "## Intro to Monte Carlo Simulation\n",
    "In physics, the term ”Monte Carlo” refers to the use of random numbers. Monte\n",
    "Carlo integration is the simplest of a wide range of ”Monte Carlo methods”, where averages are\n",
    "calculated using uniform random sampling.\n",
    "\n",
    "Monte Carlo simulation methods are related to the elementary Monte Carlo integration methods that we are discussing here, but are based on more efficient non-uniform sampling schemes. We will learn more about the detail in the later part of this course."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a183920",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculates pi using a circle with radius r=1/2 (slightly faster than r=1) by MC sampling.\n",
    "\n",
    "function pionebin(n) # samples and computes the average for one bin with n points\n",
    "   sum=0\n",
    "   for i=1:n\n",
    "      x2=(rand()-0.5)^2\n",
    "      y2=(rand()-0.5)^2\n",
    "      if x2+y2 <= 0.25 \n",
    "         sum=sum+1\n",
    "      end\n",
    "   end\n",
    "   return 4.0*sum/n\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0b099c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "function samplepi(nbin,nsamp)  # Doing MC sampling of nbin bins, each sampling nsamp points\n",
    "   av=0.\n",
    "   er=0.    \n",
    "   for j=1:nbin\n",
    "      p=pionebin(nsamp)\n",
    "      av=av+p\n",
    "      er=er+p^2\n",
    "      println(j,\"   \",p)\n",
    "   end\n",
    "   av=av/nbin\n",
    "   er=er/nbin\n",
    "   er=((er-av^2)/(nbin-1))^0.5\n",
    "   return av,er\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12baec0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "nbin=10\n",
    "nsamp=10000\n",
    "\n",
    "ap,ep=samplepi(nbin,nsamp)\n",
    "\n",
    "println()\n",
    "println(\"Final result and error:   \",ap,\"   \",ep)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9225cd0b",
   "metadata": {},
   "source": [
    "## Under the hood\n",
    "\n",
    "Julia is able to show us what is happening \"under the hood\"\n",
    "\n",
    "![](julia_compiler.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa131b23",
   "metadata": {},
   "outputs": [],
   "source": [
    "function f(x)\n",
    "    cos(x) + 2*x\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de2395b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "@code_typed f(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3aafe040",
   "metadata": {},
   "outputs": [],
   "source": [
    "@code_llvm f(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42f52cdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "@code_native f(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cd0a138",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.1",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}