# using package
using LinearAlgebra
using Plots

# assign variables
x = 55699
x = 1 + 1 * im    # complex number
x = [1, 2, 3]   # mutable list
x = [1 2; 3 4]  # matrix

2.0^2   # not **
sqrt(9)     # give Float64
sqrt(-1)    # give error
sqrt(-1.0 + 0.0 * im)     # need complex number input to get complex result

# ? to get help
# ] to enter pkg mode
# add QuadGK
# ; to enter shell

# include("./myfile.jl") to load a file

# divisor and remainder

div(8,5)   # should be 1
8 ÷ 5   # \div<tab>

rem(8,5)
8 % 5

# function:

f = x -> 2 * x  # in-line
# or
f(x) = 2 * x    # lousy

# general way to construct a function
# std. form
function f(x)
    ret = 2 * x
    return ret
end

# with optional arguments 
# required, optional + sequential; optional keywords (arbitrary order)
function f(x, y = 2, z = 3; a = 4.0, b = 5.0)
    ret = x + y + z^a + b
    println([x, y, z, a, b])
    return ret
end
f(1)
f(1, 10, b = 173, a = 77)
# have to specify y before we can specify z
f(1, 2, 4, b = 173)

# use of NamedTuples
options = (a=5.0, b=3.0)
f(1, 2, 4; options...)

# evaluate a function
f = x->x^2
f(2)
# pipe it
2 |> f


rand()  # uniform, 0 to 1
rand(10)    # array of 10 random number
rand(2, 2)   # 2 x 2 random matrix (real)

zeros(x)  # suitable zero for the variable type of x
ones(x) 

# code internal: multiple dispatch and type stability
methods(*, [Int, Int])
@less 2 + 2.0
@code_llvm 2 + 4.0
@code_native 1 + 1
# check type
typeof(x)
eltype(x)

# arrays and generators
N = 50
xmin = 0.0
xmax = 1.0
xarr = [xmin + (i1 - 1) * (xmax - xmin) / (N - 1) for i1 in 1:N]

# create an array with similar dimension and fill it with values
x1arr = similar(xarr)
fill!(x1arr, pi)

# similar construction but with a generator (iterator)
xarr = range(0.0, 1.0, length = 50)
xarr = range(0.0, 1.0, step = 1/49)

gen = (x^2 for x in xarr)

# convert a generator to an array
xarr = collect(gen)

# "." operation: on each element
xarr .^ 2
sin.(xarr)
# more convenient
@. sin(xarr)^4

[sin(xx)^4 for xx in range(0.0, 1.0, length = 50)]

# summation
# by standard loop
N = 10
ret = 0.0
for i1 in 1:N
    ret += 1 / i1^2
end

# with list
ret = sum([1 / i1^2 for i1 in 1:N])

# really condensed notation
ret = sum(@. 1 / (1:N)^2)

# cool way to write recursive function
# 3-way if statement
fac(N) = N == 1 ? 1 : N * fac(N - 1)

function fac1(N::Int)
    if N == 1
        ret = 1
    else
        ret = N * fac1(N - 1)
    end
    return ret
end

# what if N is not an integer?
using SpecialFunctions
gamma(N + 1)  # :D
methods(gamma)

# list construction

data = []
for i1 in 1:10
    push!(data, i1)
end
data = [x for x in 1:10]

# matrix operation
x = [1 2; 3 4]
y = rand(2, 2)

x * y
inv(x)
sqrt(x)
exp(x)
log(x)
tr(x)
diag(x)

eigen(x)
eigvals(x)
eigvecs(x)

# A * x = b
A = rand(2,2)
b = rand(2)
A \ b

# asking user for input
char = readline()
ff = parse(Float64, readline())
i1 = parse(Int, readline())

# simple I/O (write to & read from file)

# prepare your data in whatever way
# e.g. 10 rows, 3 columns
df = DataFrame(rand(10, 3), :auto)

# DataFrames & CSV (recommended)
using CSV

# to write
CSV.write("result.dat", df)

# to read
df = CSV.File("result.dat")

# nice looking dataset
df
push!(df, [1, 2, 3])
df.x1
df.x2

# in REPL, ; to shell for further manipulation of the datafile.

# useful options: 
# comment with #, spaces as delimiters, no header, etc.
dataf = CSV.File(
	"result.dat",
	comment = "#",
	delim = " ",
  header = false
)

# DelimitedFiles
# see delimited files https://docs.julialang.org/en/v1/stdlib/DelimitedFiles/

using DelimitedFiles
# to write
f(x, y) = x * y
data = []
for i1 in 1:10
    x = i1
    y = rand()
    push!(data, [x, y, f(x, y)]') # entry as column then transposed
    # push!(data, [x y f(x, y)])  # entry as row
end
data = vcat(data...) # turn data into a matrix

# to get the 2nd column
data[:, 2]

filename = "result.dat"
open(filename, "w") do f
    Base.write(f, "# some captions \n")
    writedlm(f, data)
end

# to read, comments like #... are ignored
file = filename
data2 = readdlm(file, comments = true)

data2 == data
# to fetch, e.g. xvec
xvec = data[:, 1]

# quick plot
rval = 2.5
p1 = plot(xarr, x -> x^2, label="r=$rval")
title("title of plot")
legend(["line 1", "line 2"], "northwest")
savefig(p1, "fig.png")  # or pdf

# ways to get the same matrix
mat1 = [1 2 3; 4 5 6; 7 8 9]
mat1[1, 3]

# column based thinking
mat2 = [[1, 4, 7] [2, 5, 8] [3, 6, 9]]

# reshape
mat3 = [i1 for i1 in 1:9]
mat3 = reshape(mat3, (3, 3))
mat3'