using RCall

R"""
library(Matrix)
Rgibbs <- function(N, thin) {
  mat <- matrix(0, nrow=N, ncol=2)
  x <- y <- 0
  for (i in 1:N) {
    for (j in 1:thin) {
      x <- rgamma(1, 3, y * y + 4) # 3rd arg is rate
      y <- rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))
    }
    mat[i,] <- c(x, y)
  }
  mat
}
"""

R"""
system.time(Rgibbs(10000, 500))
"""

using Distributions

function jgibbs(N, thin)
    mat = zeros(N, 2)
    x = y = 0.0
    for i in 1:N
        for j in 1:thin
            x = rand(Gamma(3, 1 / (y * y + 4)))
            y = rand(Normal(1 / (x + 1), 1 / sqrt(2(x + 1))))
        end
        mat[i, 1] = x
        mat[i, 2] = y
    end
    mat
end

jgibbs(100, 5); # warm-up
@elapsed jgibbs(10000, 500)

;ls ~/.julia/packages

using Distributions

pathof(Distributions)

using RCall

x = randn(1000)
R"""
hist($x, main="I'm plotting a Julia vector")
"""

R"""
library(ggplot2)
qplot($x)
"""

x = R"""
rnorm(10)
"""

# collect R variable into Julia workspace
y = collect(x)

# an integer, same as int in R
y = 1
typeof(y) 

# a Float64 number, same as double in R
y = 1.0
typeof(y) 

# Greek letters:  `\pi<tab>`
π

typeof(π)

# Greek letters:  `\theta<tab>`
θ = y + π

# emoji! `\:kissing_cat:<tab>`
😽 = 5.0

# `\alpha<tab>\hat<tab>`
α̂ = π

# vector of Float64 0s
x = zeros(5)

# vector Int64 0s
x = zeros(Int, 5)

# matrix of Float64 0s
x = zeros(5, 3)

# matrix of Float64 1s
x = ones(5, 3)

# define array without initialization
x = Matrix{Float64}(undef, 5, 3)

# fill a matrix by 0s
fill!(x, 0)

# initialize an array to be constant 2.5
fill(2.5, (5, 3))

# rational number
a = 3//5

typeof(a)

b = 3//7

a + b

# uniform [0, 1) random numbers
x = rand(5, 3)

# uniform random numbers (in single precision)
x = rand(Float16, 5, 3)

# random numbers from {1,...,5}
x = rand(1:5, 5, 3)

# standard normal random numbers
x = randn(5, 3)

# range
1:10

typeof(1:10)

1:2:10

typeof(1:2:10)

# integers 1-10
x = collect(1:10)

# or equivalently
[1:10...]

# Float64 numbers 1-10
x = collect(1.0:10)

# convert to a specific type
convert(Vector{Float64}, 1:10)

using Random # standard library
Random.seed!(123) # seed
x = rand(1_000_000) # 1 million random numbers in [0, 1)

@time sum(x) # first run includes compilation time

@time sum(x) # no compilation time after first run

# just the runtime
@elapsed sum(x)

# just the allocation
@allocated sum(x)

using BenchmarkTools

bm = @benchmark sum($x)

using Statistics # standard library
benchmark_result = Dict() # a dictionary to store median runtime (in milliseconds)
benchmark_result["Julia builtin"] = median(bm.times) / 1e6

using Libdl

C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()   # make a temporary file

# compile to a shared library by piping C_code to gcc
# (works only if you have gcc installed):

open(`gcc -std=c99 -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code) 
end

# define a Julia function that calls the C function:
c_sum(X::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)

# make sure it gives same answer
c_sum(x)

bm = @benchmark c_sum($x)

# store median runtime (in milliseconds)
benchmark_result["C"] = median(bm.times) / 1e6

using RCall

R"""
library(microbenchmark)
y <- $x
rbm <- microbenchmark(sum(y))
"""

# store median runtime (in milliseconds)
@rget rbm # dataframe
benchmark_result["R builtin"] = median(rbm[:time]) / 1e6

using RCall

R"""
sum_r <- function(x) {
  s <- 0
  for (xi in x) {
    s <- s + xi
  }
  s
}
library(microbenchmark)
y <- $x
rbm <- microbenchmark(sum_r(y))
"""

# store median runtime (in milliseconds)
@rget rbm # dataframe
benchmark_result["R loop"] = median(rbm[:time]) / 1e6

using PyCall
PyCall.pyversion

# get the Python built-in "sum" function:
pysum = pybuiltin("sum")
bm = @benchmark $pysum($x)

# store median runtime (in miliseconds)
benchmark_result["Python builtin"] = median(bm.times) / 1e6

using PyCall

py"""
def py_sum(A):
    s = 0.0
    for a in A:
        s += a
    return s
"""

sum_py = py"py_sum"

bm = @benchmark $sum_py($x)

# store median runtime (in miliseconds)
benchmark_result["Python loop"] = median(bm.times) / 1e6

# bring in sum function from Numpy 
numpy_sum = pyimport("numpy")."sum"

bm = @benchmark $numpy_sum($x)

# store median runtime (in miliseconds)
benchmark_result["Python numpy"] = median(bm.times) / 1e6

benchmark_result

x = randn(5, 3)

size(x)

size(x, 1) # nrow() in R

size(x, 2)

# total number of elements
length(x)

# 5 × 5 matrix of random Normal(0, 1)
x = randn(5, 5)

# first column
x[:, 1]

# first row
x[1, :]

# sub-array
x[1:2, 2:3]

# getting a subset of a matrix creates a copy, but you can also create "views"
z = view(x, 1:2, 2:3)

# same as
@views z = x[1:2, 2:3]

# change in z (view) changes x as well
z[2, 2] = 0.0
x

# y points to same data as x
y = x

# x and y point to same data
pointer(x), pointer(y)

# changing y also changes x
y[:, 1] .= 0
x

# create a new copy of data
z = copy(x)

pointer(x), pointer(z)

# 1-by-3 array
[1 2 3]

# 3-by-1 vector
[1, 2, 3]

# multiple assignment by tuple
x, y, z = randn(5, 3), randn(5, 2), randn(3, 5)

[x y] # 5-by-5 matrix

[x y; z] # 8-by-5 matrix

x = randn(5, 3)

y = ones(5, 3)

x .* y # same x * y in R

x .^ (-2)

sin.(x)

x = randn(5)

using LinearAlgebra
# vector L2 norm
norm(x)

# same as
sqrt(sum(abs2, x))

y = randn(5) # another vector
# dot product
dot(x, y) # x' * y

# same as
x'y

x, y = randn(5, 3), randn(3, 2)
# matrix multiplication, same as %*% in R
x * y

x = randn(3, 3)

# conjugate transpose
x'

b = rand(3)
x'b # same as x' * b

# trace
tr(x)

det(x)

rank(x)

using SparseArrays

# 10-by-10 sparse matrix with sparsity 0.1
X = sprandn(10, 10, .1)

# convert to dense matrix; be cautious when dealing with big data
Xfull = convert(Matrix{Float64}, X)

# convert a dense matrix to sparse matrix
sparse(Xfull)

# syntax for sparse linear algebra is same as dense linear algebra
β = ones(10)
X * β

# many functions apply to sparse matrices as well
sum(X)

function myfunc(x)
    return sin(x^2)
end

x = randn(5, 3)
myfunc.(x)

map(x -> sin(x^2), x)

map(x) do elem
    elem = elem^2
    return sin(elem)
end

# Mapreduce
mapreduce(x -> sin(x^2), +, x)

# same as
sum(x -> sin(x^2), x)

[sin(2i + j) for i in 1:5, j in 1:3] # similar to Python

typeof(1.0), typeof(1)

supertype(Float64)

subtypes(AbstractFloat)

# Is Float64 a subtype of AbstractFloat?
Float64 <: AbstractFloat

# On 64bit machine, Int == Int64
Int == Int64

# convert to Float64
convert(Float64, 1)

# same as
Float64(1)

# Float32 vector
x = randn(Float32, 5)

# convert to Float64
convert(Vector{Float64}, x)

# same as
Float64.(x)

# convert Float64 to Int64
convert(Int, 1.0)

convert(Int, 1.5) # should use round(1.5)

round(Int, 1.5)

g(x) = x + x

g(1.5)

g("hello world")

g(x::Float64) = x + x

g(x::Number) = x + x

methods(g)

# an Int64 input
@which g(1)

# a Vector{Float64} input
@which g(randn(5))

g(2), g(2.0)

@code_lowered g(2)

@code_warntype g(2)

@code_warntype g(2.0)

@code_llvm g(2)

@code_llvm g(2.0)

@code_native g(2)

@code_native g(2.0)

# a function defined earlier
function tally(x::Array)
    s = zero(eltype(x))
    for v in x
        s += v
    end
    s
end

using Random
Random.seed!(123)
a = rand(20_000_000)
@time tally(a) # first run: include compile time

@time tally(a)

using BenchmarkTools

@benchmark tally($a)

using Profile

Profile.clear()
@profile tally(a)
Profile.print(format=:flat)

;cat bar.jl

;julia --track-allocation=user bar.jl

;cat bar.jl.21740.mem