## Sparse vs. Dense matrices in Julia
using SparseArrays # This loads the sparse matrix package

## Compare the performance of matrix-vector products for sparse and dense matrices
# We'll do this as a function of how many entries are non-zero in the matrix
using Random, BenchmarkTools, LinearAlgebra 
# n is the size of the sparse matrix
# p is the fraction of entries that are non-zero. 
function trial(n, p)
    Random.seed!(0)
    A = sprand(n, n, p) # Generate a random sparse matrix
    x = rand(n) # Generate a random vector
    y = zeros(n) # Allocate space for the result
    return (@btime(mul!($y, $A, $x))) # this does an in-place multiplication 
end
function dense_trial(n, p)
  Random.seed!(0)
  A = Matrix(sprand(n, n, p)) # Generate a random sparse matrix, but as dense 
  x = rand(n) # Generate a random vector
  y = zeros(n) # Allocate space for the result
  return (@btime(mul!($y, $A, $x))) # this does an in-place multiplication 
end
trial(100, 0.01) # Warm up the JIT compiler
dense_trial(100, 0.01) # Warm up the JIT compiler
trial(1000, 0.1) # Run a trial 
dense_trial(1000, 0.1) # Run a trial 

  129.843 ns (0 allocations: 0 bytes)
  1.350 μs (0 allocations: 0 bytes)
  45.500 μs (0 allocations: 0 bytes)
  89.041 μs (0 allocations: 0 bytes)

1000-element Vector{Float64}:
 24.380572601576063
 29.078132968701016
 26.575361707841537
 21.18796263382051
 29.670362571076172
 36.14012297927813
 26.551888619216744
 31.663151968139935
 23.66509132583557
 24.740730204397682
 20.566221249206468
 24.228092840903123
 25.384055293445066
  ⋮
 22.399432870698334
 26.224067484467533
 22.659537327151032
 29.90517791786851
 27.66856355647173
 22.674037125669766
 25.762337883311627
 21.329546462870077
 29.628815658532766
 21.70315680771119
 27.30104276399619
 23.890768516450965

##
trial(1000, 0.33) # Run a trial 
dense_trial(1000, 0.33) # Run a trial 

  143.000 μs (0 allocations: 0 bytes)
  88.792 μs (0 allocations: 0 bytes)

1000-element Vector{Float64}:
 87.65879414785158
 84.49985208695149
 81.23015351551086
 79.25604914256748
 75.8174366973627
 79.82391725779192
 90.11858231311294
 79.79555198500297
 74.89506475846503
 95.49818148609522
 89.23806108312404
 81.85894282221228
 84.37231956850181
  ⋮
 86.58256983836812
 81.95455475483641
 81.9751324898719
 82.1796160712076
 77.58724245805077
 78.9266037183557
 79.530006257086
 95.41480861614937
 89.98257332957085
 82.86447942309196
 87.72695513617735
 82.06992896076437

## Let's make a little plot of this. 
function get_data(ns, ps)
  Random.seed!(0)
  results_sparse = zeros(length(ns), length(ps))
  results_dense = zeros(length(ns), length(ps))
  for (i, n) in pairs(ns) # pairs is a handy way of getting the index and value... 
    for (j, p) in pairs(ps)
      A = sprand(n, n, p) # Generate a random sparse matrix
      x = rand(n) # Generate a random vector
      y = zeros(n) # Allocate space for the result
      t = @elapsed(mul!(y, A, x)) # this does an in-place multiplication
      A = Matrix(A) # Convert to dense
      t_dense = @elapsed(mul!(y, A, x)) # this does an in-place multiplication
      results_sparse[i,j] = t
      results_dense[i,j] = t_dense
    end
  end
  return results_sparse, results_dense
end 
ns = [100, 200, 400, 800, 1600]
ps = [0.01, 0.05, 0.1, 0.2, 0.33, 0.5]
results_sparse, results_dense = get_data(ns, ps)

([1.375e-6 1.167e-6 … 2.2e-5 2.9875e-5; 1.875e-6 3.833e-6 … 7.375e-6 1.0333e-5; … ; 1.325e-5 1.9625e-5 … 0.000111291 0.000462833; 2.2458e-5 6.7375e-5 … 0.000411 0.000944375], [4.292e-6 2.25e-6 … 6.625e-6 4.2041e-5; 0.0001615 2.75e-5 … 7.834e-6 7.75e-6; … ; 0.00046775 0.000151667 … 0.000385959 0.001171584; 0.000681542 0.001329709 … 0.002080084 0.001563125])

## Now let's plot
using Plots
plot!(size=(800,600))
for (j, ps) in pairs(ps)
  plot!(ns, results_sparse[:,j], label="Sparse, p=$ps", marker=:circle)
  plot!(ns, results_dense[:,j], label="Dense, p=$ps", marker=:circle, linestyle=:dash)
end
xlabel!("Matrix size")
ylabel!("Time (s)")

## Do this with lots of little plots, one for each p
plts = [] 
for (j,p) in pairs(ps)
  plt = plot(ns, results_sparse[:,j], label="Sparse", marker=:circle)
  plot!(plt, ns, results_dense[:,j], label="Dense", marker=:circle, linestyle=:dash)
  plot!(plt, yscale=:log10)
  title!(plt, "p = $p")
  push!(plts, plt) # add pl
end 
plot(plts..., layout=(3,2), size=(800,600)) # plot all the plots

## Scopes in Julia
classvar1 = 5
function f1()
  i = 20 
  classvar1 = 10
  @show classvar1 
  classvar3 = 20 
  @show classvar3 
  for classvar2 in 1:10
    classvar3 = 30
    # println("classvar3 = ", classvar3)
    @show classvar3, classvar2
  end
  @show classvar1 # This will fail because classvar3 is not defined in this scope
  @show classvar3
end
f1()
println(classvar1)

classvar1 = 10
classvar3 = 20
(classvar3, classvar2) = (30, 1)
(classvar3, classvar2) = (30, 2)
(classvar3, classvar2) = (30, 3)
(classvar3, classvar2) = (30, 4)
(classvar3, classvar2) = (30, 5)
(classvar3, classvar2) = (30, 6)
(classvar3, classvar2) = (30, 7)
(classvar3, classvar2) = (30, 8)
(classvar3, classvar2) = (30, 9)
(classvar3, classvar2) = (30, 10)
classvar1 = 10
classvar3 = 30
5

##
function testfunc(xs)
  for i in xs
    @show i 
  end 
end 
testfunc([1,5,10])

i = 1
i = 5
i = 10