using InstantiateFromURL
# optionally add arguments to force installation: instantiate = true, precompile = true
github_project("QuantEcon/quantecon-notebooks-julia", version = "0.7.0")

using LinearAlgebra, Statistics

a = [10, 20, 30]

a = [1.0, 2.0, 3.0]

typeof(randn(100))

ndims(a)

size(a)

Array{Int64, 1} == Vector{Int64}
Array{Int64, 2} == Matrix{Int64}

[1, 2, 3] == [1; 2; 3]  # both column vectors

[1 2 3]  # a row vector is 2-dimensional

zeros(3)

zeros(2, 2)

fill(5.0, 2, 2)

x = Array{Float64}(undef, 2, 2)

fill(0, 2, 2)  # fills with 0, not 0.0

fill(false, 2, 2)  # produces a boolean matrix

x = [1, 2, 3]
y = x
y[1] = 2
x

x = [1, 2, 3]
y = copy(x)
y[1] = 2
x

x = [1, 2, 3]
y = similar(x)
y

x = [1, 2, 3]
y = similar(x, 4)  # make a vector of length 4

x = [1, 2, 3]
y = similar(x, 2, 2)  # make a 2x2 matrix

a = [10, 20, 30, 40]

a = [10 20 30 40]  # two dimensional, shape is 1 x n

ndims(a)

a = [10 20; 30 40]  # 2 x 2

a = [10; 20; 30; 40]

ndims(a)

a = [10 20 30 40]'

ndims(a)

a = [10 20 30 40]
a[end-1]

a[1:3]

a = randn(2, 2)
a[1, 1]

a[1, :]  # first row

a[:, 1]  # first column

a = randn(2, 2)

b = [true false; false true]

a[b]

a = zeros(4)

a[2:end] .= 42

a

a = [1 2; 3 4]
b = a[:, 2]
@show b
a[:, 2] = [4, 5] # modify a
@show a
@show b;

a = [1 2; 3 4]
@views b = a[:, 2]
@show b
a[:, 2] = [4, 5]
@show a
@show b;

@views b = a[:, 2]
view(a, :, 2) == b

a = [1 2; 3 4]
b_slice = a[:, 2]
@show typeof(b_slice)
@show typeof(a)
@views b = a[:, 2]
@show typeof(b);

a = [1 2; 3 4]
b = a'   # transpose
typeof(b)

a = [1 2; 3 4]
b = a' # transpose
c = Matrix(b)  # convert to matrix
d = collect(b) # also `collect` works on any iterable
c == d

d = [1.0, 2.0]
a = Diagonal(d)

@show 2a
b = rand(2,2)
@show b * a;

b = [1.0 2.0; 3.0 4.0]
b - Diagonal([1.0, 1.0])  # poor style, inefficient code

b = [1.0 2.0; 3.0 4.0]
b - I  # good style, and note the lack of dimensions of I

typeof(I)

x = [1 2 3]
y = x  # name `y` binds to whatever value `x` bound to

x = [1 2 3]
y = x        # name `y` binds to whatever `x` bound to
z = [2 3 4]
y = z        # only changes name binding, not value!
@show (x, y, z);

x = [1 2 3]
y = x       # name `y` binds to whatever `x` bound to
z = [2 3 4]
y .= z      # now dispatches the assignment of each element
@show (x, y, z);

function f(x)
    return [1 2; 3 4] * x  # matrix * column vector
end

val = [1, 2]
f(val)

function f(x)
    return [1 2; 3 4] * x # matrix * column vector
end

val = [1, 2]
y = similar(val)

function f!(out, x)
    out .= [1 2; 3 4] * x
end

f!(y, val)
y

function f(x)
    return [1 2; 3 4] * x  # matrix * column vector
end

val = [1, 2]
y = similar(val)

function f!(out, x)
    out = [1 2; 3 4] * x   # MISTAKE! Should be .= or [:]
end
f!(y, val)
y

y = [1 2]
y .-= 2    # y .= y .- 2, no problem

x = 5
# x .-= 2  # Fails!
x = x - 2  # subtle difference - creates a new value and rebinds the variable

x = 2

function f(x)
    x = 3     # MISTAKE! does not modify x, creates a new value!
end

f(x)          # cannot modify immutables in place
@show x;

using StaticArrays
xdynamic = [1, 2]
xstatic = @SVector [1, 2]  # turns it into a highly optimized static vector

f(x) = 2x
@show f(xdynamic)
@show f(xstatic)

# inplace version
function g(x)
    x .= 2x
    return "Success!"
end
@show xdynamic
@show g(xdynamic)
@show xdynamic;

# g(xstatic) # fails, static vectors are immutable

a = [-1, 0, 1]

@show length(a)
@show sum(a)
@show mean(a)
@show std(a)      # standard deviation
@show var(a)      # variance
@show maximum(a)
@show minimum(a)
@show extrema(a)  # (mimimum(a), maximum(a))

b = sort(a, rev = true)  # returns new array, original not modified

b = sort!(a, rev = true)  # returns *modified original* array

b == a  # tests if have the same values

b === a  # tests if arrays are identical (i.e share same memory)

a = ones(1, 2)

b = ones(2, 2)

a * b

b * a'

A = [1 2; 2 3]

B = ones(2, 2)

A \ B

inv(A) * B

ones(2) * ones(2)  # does not conform, expect error

ones(2)' * ones(2)

@show dot(ones(2), ones(2))
@show ones(2) ⋅ ones(2);

b = ones(2, 2)
b * ones(2)

ones(2)' * b

ones(2, 2) * ones(2, 2)   # matrix multiplication

ones(2, 2) .* ones(2, 2)   # element by element multiplication

A = -ones(2, 2)

A.^2  # square every element

ones(2, 2) + ones(2, 2)  # same as ones(2, 2) .+ ones(2, 2)

A = ones(2, 2)

2 * A  # same as 2 .* A

x = [1, 2]
x .+ 1     # not x + 1
x .- 1     # not x - 1

a = [10, 20, 30]

b = [-100, 0, 100]

b .> a

a .== b

b

b .> 1

a = randn(4)

a .< 0

a[a .< 0]

a = [10, 20, 30, 40]

b = reshape(a, 2, 2)

b

b[1, 1] = 100  # continuing the previous example

b

a

a = [1 2 3 4]  # two dimensional

dropdims(a, dims = 1)

log(1.0)

log.(1:4)

[ log(x) for x in 1:4 ]

A = [1 2; 3 4]

det(A)

tr(A)

eigvals(A)

rank(A)

a = 10:12        # a range, equivalent to 10:1:12
@show Vector(a)  # can convert, but shouldn't

b = Diagonal([1.0, 2.0, 3.0])
b * a .- [1.0; 2.0; 3.0]

a = 0.0:0.1:1.0  # 0.0, 0.1, 0.2, ... 1.0

maxval = 1.0
minval = 0.0
stepsize = 0.15
a = minval:stepsize:maxval # 0.0, 0.15, 0.3, ...
maximum(a) == maxval

maxval = 1.0
minval = 0.0
numpoints = 10
a = range(minval, maxval, length=numpoints)
# or range(minval, stop=maxval, length=numpoints)

maximum(a) == maxval

t = (1.0, "test")
t[1]            # access by index
a, b = t        # unpack
# t[1] = 3.0    # would fail as tuples are immutable
println("a = $a and b = $b")

t = (val1 = 1.0, val2 = "test")
t.val1      # access by index
# a, b = t  # bad style, better to unpack by name with @unpack
println("val1 = $(t.val1) and val1 = $(t.val1)") # access by name

t2 = (val3 = 4, val4 = "test!!")
t3 = merge(t, t2)  # new tuple

function f(parameters)
    α, β = parameters.α, parameters.β  # poor style, error prone if adding parameters
    return α + β
end

parameters = (α = 0.1, β = 0.2)
f(parameters)

using Parameters

function f(parameters)
    @unpack α, β = parameters  # good style, less sensitive to errors
    return α + β
end

parameters = (α = 0.1, β = 0.2)
f(parameters)

using Parameters
paramgen = @with_kw (α = 0.1, β = 0.2)  # create named tuples with defaults

# creates named tuples, replacing defaults
@show paramgen()  # calling without arguments gives all defaults
@show paramgen(α = 0.2)
@show paramgen(α = 0.2, β = 0.5);

typeof(nothing)

function f(y)
    x = nothing
    if y > 0.0
        # calculations to set `x`
        x = y
    end

    # later, can check `x`
    if isnothing(x)
        println("x was not set")
    else
        println("x = $x")
    end
    x
end

@show f(1.0)
@show f(-1.0);

function f(x)
    if x > 0.0
        return sqrt(x)
    else
        return nothing
    end
end
x1 = 1.0
x2 = -1.0
y1 = f(x1)
y2 = f(x2)

# check results with isnothing
if isnothing(y1)
    println("f($x2) successful")
else
    println("f($x2) failed");
end

function f(x)
    x > 0.0 ? sqrt(x) : nothing  # the "a ? b : c" pattern is the ternary
end

f(1.0)

x = [1.0, nothing]

x = [0.1, -1.0, 2.0, -2.0]
y = f.(x)

# presumably check `y`

function f(x)
    @assert x > 0.0
    sqrt(x)
end

f(1.0)

function f(x; z = nothing)

    if isnothing(z)
        println("No z given with $x")
    else
        println("z = $z given with $x")
    end
end

f(1.0)
f(1.0, z=3.0)

function f(x)
    if x > 0.0
        return x
    else
        return NaN
    end
end

f(0.1)
f(-1.0)

@show typeof(f(-1.0))
@show f(-1.0) == NaN  # note, this fails!
@show isnan(f(-1.0))  # check with this

# throws exception, turned off to prevent breaking notebook
# sqrt(-1.0)

# to see the error
try sqrt(-1.0); catch err; err end  # catches the exception and prints it

# throws exception, turned off to prevent breaking notebook
# convert(Int64, 3.12)

# to see the error
try convert(Int64, 3.12); catch err; err end  # catches the exception and prints it.

function f(x)
    try
        sqrt(x)
    catch err                # enters if exception thrown
        sqrt(complex(x, 0))  # convert to complex number
    end
end

f(0.0)
f(-1.0)

x = [3.0, missing, 5.0, missing, missing]

f(x) = x^2

@show missing + 1.0
@show missing * 2
@show missing * "test"
@show f(missing);      # even user-defined functions
@show mean(x);

x = missing

@show x == missing
@show x === missing  # an exception
@show ismissing(x);

x = [1.0, missing, 2.0, missing, missing, 5.0]

@show mean(x)
@show mean(skipmissing(x))
@show coalesce.(x, 0.0);  # replace missing with 0.0;

function compute_asymptotic_var(A, Σ;
                                S0 = Σ * Σ',
                                tolerance = 1e-6,
                                maxiter = 500)
    V = Σ * Σ'
    S = S0
    err = tolerance + 1
    i = 1
    while err > tolerance && i ≤ maxiter
        next_S = A * S * A' + V
        err = norm(S - next_S)
        S = next_S
        i += 1
    end
    return S
end

A = [0.8  -0.2;
     -0.1  0.7]

Σ = [0.5 0.4;
     0.4 0.6]

maximum(abs, eigvals(A))

our_solution = compute_asymptotic_var(A, Σ)

using QuantEcon

norm(our_solution - solve_discrete_lyapunov(A, Σ * Σ'))