# I downloaded Julia from here:
# 
# https://juliacomputing.com/products/juliapro.html
#
# And I found these Julia-for-Python-Programmers-style cheatsheets useful:
# 
# https://juliadocs.github.io/Julia-Cheat-Sheet/
# https://cheatsheets.quantecon.org/
# http://math.mit.edu/~stevenj/Julia-cheatsheet.pdf
#
# I needed to call this all of this code once (and only once) to get PyPlot
# to work. The Julia interpreter gave me an error message about
# mpl_toolkits.mplot3d not being installed, but afterwards "using PyPlot"
# seemed to work anyway.
#
# ENV["PYTHON"] = ""
# Pkg.build("PyCall")
# Pkg.add("PyCall")
# Pkg.add("PyPlot")
#

using PyPlot

x_min       = -5.0
x_max       = 5.0
num_samples = 101
dx          = (x_max - x_min) / (num_samples - 1)

# Create array with evenly spaced samples between min and max values.
x = linspace(x_min,x_max,num_samples)

mu            = 0.0 # mean
sigma_sqr     = 1.0 # variance
sigma_sqr_inv = 1.0/sigma_sqr

# Non-normalized 1D Gaussian function; can easily construct using
# Matlab-style expression syntax, i.e., using an array as a scalar.
G = exp.(-(1.0/2.0)*sigma_sqr_inv*(x - mu).^2.0)

plot(x,G,label=L"$\mu = 0.0, \sigma^2 = 1.0$");
legend();

# Simplest possible numerical integration.
println(sum(G*dx))

# meshgrid function adapted from https://github.com/JuliaLang/julia/blob/master/examples/ndgrid.jl

function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    return (repmat(vx, m, 1), repmat(vy, 1, n))
end

x1_min         = -5.0
x2_min         = -5.0
x1_max         = 5.0
x2_max         = 5.0
num_samples_x1 = 11
num_samples_x2 = 11
dx1            = (x1_max - x1_min) / (num_samples_x1 - 1)
dx2            = (x2_max - x2_min) / (num_samples_x2 - 1)
half_dx1       = 0.5*dx1
half_dx2       = 0.5*dx2

# Matlab-style meshgrid function.
X2,X1 = meshgrid(linspace(x2_min,x2_max,num_samples_x2), linspace(x1_min,x1_max,num_samples_x1))

pcolormesh(linspace(x2_min-half_dx2,x2_max+half_dx2,num_samples_x2+1),
           linspace(x1_min-half_dx1,x1_max+half_dx1,num_samples_x1+1),
           X1);
plt[:gca]()[:set_aspect]("equal");
colorbar();

pcolormesh(linspace(x2_min-half_dx2,x2_max+half_dx2,num_samples_x2+1),
           linspace(x1_min-half_dx1,x1_max+half_dx1,num_samples_x1+1),
           X2);
plt[:gca]()[:set_aspect]("equal");
colorbar();

# Can easily use the output arrays from meshgrid in mathematical expressions.
F = X1.^2.0 + X2.^2.0

pcolormesh(linspace(x2_min-half_dx2,x2_max+half_dx2,num_samples_x2+1),
           linspace(x1_min-half_dx1,x1_max+half_dx1,num_samples_x1+1),
           F);
plt[:gca]()[:set_aspect]("equal");
colorbar();

x1_min         = -5.0
x2_min         = -5.0
x1_max         = 5.0
x2_max         = 5.0
num_samples_x1 = 101
num_samples_x2 = 101
dx1            = (x1_max - x1_min) / (num_samples_x1 - 1)
dx2            = (x2_max - x2_min) / (num_samples_x2 - 1)
half_dx1       = 0.5*dx1
half_dx2       = 0.5*dx2

X2,X1 = meshgrid(linspace(x2_min,x2_max,num_samples_x2), linspace(x1_min,x1_max,num_samples_x1))

F = X1.^2.0 + X2.^2.0

plot_surface(X2,X1,F,cmap="viridis");

# Non-normalized 2D Gaussian function.
function G_func_2d(mu, Sigma_sqr_inv, x)
    G = exp.(-(1.0/2.0)*(x - mu)'*Sigma_sqr_inv*(x - mu))
    return G[1,1]
end

mu            = [0 0]'
Sigma_sqr_inv = [1 2; 2 10]
# Sigma_sqr_inv = [1 0; 0 1]
G             = zeros(num_samples_x1,num_samples_x2)

# Can also easily construct functions with explicit for-loops
# noticeably faster than Python or Matlab.
for i1 = 1:num_samples_x1
    for i2 = 1:num_samples_x2
        x = [ X1[i1,i2] X2[i1,i2] ]'
        G[i1,i2] = G_func_2d(mu, Sigma_sqr_inv, x)
    end
end

plot_surface(X2,X1,G,cmap="viridis");

pcolormesh(linspace(x2_min-half_dx2,x2_max+half_dx2,num_samples_x2+1),
           linspace(x1_min-half_dx1,x1_max+half_dx1,num_samples_x1+1),
           G);
plt[:gca]()[:set_aspect]("equal");
colorbar();

# Can slice arrays as in Matlab and Python.
plot(X1[:,61],G[:,61]);