text:
#=
QR+Cholesky linear system solver
requires precomputed QR factorization of A'

solves linear system in naive.jl by first eliminating s, kap, and tau via the
method in the symindef solver, then reducing the 3x3 symmetric indefinite system
to a series of low-dimensional operations via a procedure similar to that
described by S10.3 of
http://www.seas.ucla.edu/~vandenbe/publications/coneprog.pdf

TODO iterative refinement for 3x3 solution
=#

abstract type QRCholSystemSolver{T <: Real} <: SystemSolver{T} end

function setup_rhs3(
    ::QRCholSystemSolver{T},
    model::Models.Model{T},
    rhs::Point{T},
    sol::Point{T},
    rhs_sub::Point{T},
    ) where {T <: Real}
    @inbounds for (k, cone_k) in enumerate(model.cones)
        rhs_z_k = rhs.z_views[k]
        rhs_s_k = rhs.s_views[k]
        rhs_sub_z_k = rhs_sub.z_views[k]
        if Cones.use_dual_barrier(cone_k)
            z_temp_k = sol.z_views[k]
            @. z_temp_k = -rhs_z_k - rhs_s_k
            Cones.inv_hess_prod!(rhs_sub_z_k, z_temp_k, cone_k)
        else
            Cones.hess_prod!(rhs_sub_z_k, rhs_z_k, cone_k)
            axpby!(-1, rhs_s_k, -1, rhs_sub_z_k)
        end
    end
    return
end

function solve_subsystem3(
    syssolver::QRCholSystemSolver{T},
    solver::Solver{T},
    sol::Point{T},
    rhs::Point{T},
    ) where {T <: Real}
    model = solver.model
    copyto!(sol.vec, rhs.vec)
    x = sol.x
    y = sol.y
    z = sol.z

    copyto!(syssolver.QpbxGHbz, x)
    mul!(syssolver.QpbxGHbz, model.G', z, true, true)
    lmul!(solver.Ap_Q', syssolver.QpbxGHbz)

    if !iszero(model.p)
        ldiv!(solver.Ap_R', y)
        @views copyto!(sol.vec[1:model.p], y)

        if !isempty(syssolver.Q2div)
            mul!(syssolver.GQ1x, syssolver.GQ1, y)
            block_hess_prod!.(syssolver.HGQ1x_k, syssolver.GQ1x_k, model.cones)
            mul!(syssolver.Q2div, syssolver.GQ2', syssolver.HGQ1x, -1, true)
        end
    end

    if !isempty(syssolver.Q2div)
        @views x_sub2 = sol.vec[(model.p + 1):model.n]
        ldiv!(x_sub2, syssolver.fact, syssolver.Q2div)
    end

    lmul!(solver.Ap_Q, x)

    mul!(syssolver.Gx, model.G, x)
    block_hess_prod!.(syssolver.HGx_k, syssolver.Gx_k, model.cones)

    axpby!(true, syssolver.HGx, -1, z)

    if !iszero(model.p)
        copyto!(y, syssolver.Q1pbxGHbz)
        mul!(y, syssolver.GQ1', syssolver.HGx, -1, true)
        ldiv!(solver.Ap_R, y)
    end

    return sol
end

function block_hess_prod!(
    prod_k::AbstractVecOrMat{T},
    arr_k::AbstractVecOrMat{T},
    cone_k::Cones.Cone{T},
    ) where {T <: Real}
    if Cones.use_dual_barrier(cone_k)
        Cones.inv_hess_prod!(prod_k, arr_k, cone_k)
    else
        Cones.hess_prod!(prod_k, arr_k, cone_k)
    end
    return
end

#=
direct dense
=#

mutable struct QRCholDenseSystemSolver{T <: Real} <: QRCholSystemSolver{T}
    lhs_sub::Symmetric{T, Matrix{T}}
    lhs_sub_fact::Symmetric{T, Matrix{T}}
    fact::Factorization{T}

    rhs_sub::Point{T}
    sol_sub::Point{T}
    sol_const::Point{T}
    rhs_const::Point{T}

    use_sqrt_hess_cones::Vector{Bool}
    GQ1::AbstractMatrix{T}
    GQ2::AbstractMatrix{T}
    QpbxGHbz::Vector{T}
    Q1pbxGHbz::SubArray
    Q2div::SubArray
    GQ1x::Vector{T}
    HGQ1x::Vector{T}
    HGQ2::Matrix{T}
    Gx::Vector{T}
    HGx::Vector{T}
    HGQ1x_k::Vector
    GQ1x_k::Vector
    HGQ2_k::Vector
    GQ2_k::Vector
    HGx_k::Vector
    Gx_k::Vector

    function QRCholDenseSystemSolver{T}() where {T <: Real}
        syssolver = new{T}()
        return syssolver
    end
end

function load(
    syssolver::QRCholDenseSystemSolver{T},
    solver::Solver{T},
    ) where {T <: Real}
    model = solver.model
    (n, p, q) = (model.n, model.p, model.q)
    nmp = n - p
    cone_idxs = model.cone_idxs

    syssolver.lhs_sub = Symmetric(zeros(T, nmp, nmp), :U)
    syssolver.lhs_sub_fact = zero(syssolver.lhs_sub)

    syssolver.use_sqrt_hess_cones = falses(length(cone_idxs))

    # very inefficient method used for sparse G * QRSparseQ
    # see https://github.com/JuliaLang/julia/issues/31124#issuecomment-501540818
    GQ = model.G * solver.Ap_Q

    syssolver.GQ2 = GQ[:, (p + 1):end]
    syssolver.HGQ2 = zeros(T, q, nmp)
    syssolver.QpbxGHbz = zeros(T, n)
    syssolver.Q2div = view(syssolver.QpbxGHbz, (p + 1):n)
    syssolver.Gx = zeros(T, q)
    syssolver.HGx = zeros(T, q)
    syssolver.HGQ2_k = [view(syssolver.HGQ2, idxs, :) for idxs in cone_idxs]
    syssolver.GQ2_k = [view(syssolver.GQ2, idxs, :) for idxs in cone_idxs]
    syssolver.HGx_k = [view(syssolver.HGx, idxs, :) for idxs in cone_idxs]
    syssolver.Gx_k = [view(syssolver.Gx, idxs, :) for idxs in cone_idxs]

    if !iszero(p)
        syssolver.GQ1 = GQ[:, 1:p]
        syssolver.Q1pbxGHbz = view(syssolver.QpbxGHbz, 1:p)
        syssolver.GQ1x = zeros(T, q)
        syssolver.HGQ1x = zeros(T, q)
        syssolver.HGQ1x_k = [view(syssolver.HGQ1x, idxs, :) for idxs in cone_idxs]
        syssolver.GQ1x_k = [view(syssolver.GQ1x, idxs, :) for idxs in cone_idxs]
    end

    setup_point_sub(syssolver, model)

    return syssolver
end

function update_lhs(
    syssolver::QRCholDenseSystemSolver{T},
    solver::Solver{T},
    ) where {T <: Real}
    model = solver.model

    # update LHS and factorization
    isempty(syssolver.Q2div) || update_lhs_fact(syssolver, solver)

    # update solution for fixed c,b,h part
    rhs_const = syssolver.rhs_const
    @inbounds for (k, cone_k) in enumerate(model.cones)
        @views h_k = model.h[model.cone_idxs[k]]
        block_hess_prod!(rhs_const.z_views[k], h_k, cone_k)
    end
    solve_subsystem3(syssolver, solver, syssolver.sol_const, rhs_const)

    return syssolver
end

function update_lhs_fact(
    syssolver::QRCholDenseSystemSolver{T},
    solver::Solver{T},
    ) where {T <: Real}
    model = solver.model
    lhs = syssolver.lhs_sub.data
    cones = model.cones
    use_sqrt_hess_cones = syssolver.use_sqrt_hess_cones
    HGQ2 = syssolver.HGQ2
    GQ2_k = syssolver.GQ2_k
    HGQ2_k = syssolver.HGQ2_k

    nmp = size(lhs, 1)
    @inbounds for k in eachindex(cones)
        use_sqrt_hess_cones[k] = Cones.use_sqrt_hess_oracles(nmp, cones[k])
    end

    # sqrt cones
    if any(use_sqrt_hess_cones)
        idx = 1
        @inbounds for k in eachindex(cones)
            use_sqrt_hess_cones[k] || continue
            cone_k = cones[k]
            arr_k = GQ2_k[k]
            q_k = size(arr_k, 1)
            @views prod_k = HGQ2[idx:(idx + q_k - 1), :]
            if Cones.use_dual_barrier(cone_k)
                Cones.inv_sqrt_hess_prod!(prod_k, arr_k, cone_k)
            else
                Cones.sqrt_hess_prod!(prod_k, arr_k, cone_k)
            end
            idx += q_k
        end
        @views outer_prod!(HGQ2[1:(idx - 1), :], lhs, true, false)
    else
        fill!(lhs, 0)
    end

    # not sqrt cones
    @inbounds for k in eachindex(cones)
        use_sqrt_hess_cones[k] && continue
        arr_k = GQ2_k[k]
        prod_k = HGQ2_k[k]
        block_hess_prod!(prod_k, arr_k, cones[k])
        mul!(lhs, arr_k', prod_k, true, true)
    end

    # TODO try equilibration, iterative refinement etc like posvx/sysvx
    solver.time_upfact += @elapsed syssolver.fact =
        posdef_fact_copy!(syssolver.lhs_sub_fact, syssolver.lhs_sub)

    if !issuccess(syssolver.fact)
        println("positive definite linear system factorization failed")
    end

    return
end


text:
"""
    TransferOperatorTriangulationExact

A transfer operator estimated over a triangulation of a set of points. Transition
probabilities are computed from exact simplex intersections.

The return type for the `TransferOperatorEstimatorTriangulationExact` estimator.

## Fields
- **`transfermatrix`**. The estimated transfer matrix.

"""
struct TransferOperatorTriangulationExact <: AbstractTransferOperator
    transfermatrix::AbstractArray{Float64, 2}
end

export TransferOperatorTriangulationExact

ext: jl


text:
using ColorSchemesInvert
using Base.Test
using Colors

@testset "invert tests" begin
cs = linspace(RGB(0,0,0), RGB(1,1,1),5);
@test invert(cs, cs[3]) == 0.5

# Note that invert() takes the first closest match.
cs = [RGB(0,0,0), RGB(1,1,1),
      RGB(0,0,0),
      RGB(0,0,0), RGB(1,1,1)];
@test invert(cs, cs[3]) == 0

# Note that invert() takes the left index when two are identical.
cs = [RGB(1,1,1), RGB(0,0,0), RGB(0,0,0), RGB(0,1,1), RGB(1,1,0)];
@test invert(cs, cs[2]) == 0.25
cs = [RGB(0,0,0), RGB(0,0,0)];
@test invert(cs, cs[1]) == 0

cs = [RGB(0,0,0)];
@test_throws InexactError invert(cs, cs[1])
# (The above line throws for the same reason the below line does.
#  If this behavior ever changes, so should `invert`.)
@test_throws InexactError get(cs, 1.0, (1,1))
end

@testset "convertToScheme tests" begin
# Add color to a grayscale image.
red_cs = linspace(RGB(0,0,0), RGB(1,0,0))
gray_img = linspace(RGB(0,0,0), RGB(1,1,1))
new_img = convertToScheme(red_cs, gray_img)
@test all(.≈(new_img, red_cs, atol=0.5))  # This is broken.. It should be way more specific. See next test.

# Should be able to uniquely match each increasing color with the next
# increasing color in the new scale.
red_cs = [linspace(RGB(0,0,0), RGB(1,0,0))...]
blue_scale_img = [linspace(RGB(0,0,0), RGB(0,0,1))...]
new_img = convertToScheme(red_cs, blue_scale_img)
@test_broken unique(new_img) == new_img
end


text:
#=
Find the unique positive integer whose square has the form 1_2_3_4_5_6_7_8_9_0,
where each “_” is a single digit.
=#

#TODO: Try iterating by calculating the squares as permutations then grabbing sqrts.
function calc()
  # This is kinda icky.
  for i in 1238902700:3238902700
    s = i^2
    if length("$s") != 19
      continue
    end
    if occursin(r"1\d2\d3\d4\d5\d6\d7\d8\d9\d0", "$s")
      return round(Int64, sqrt(s))
    end
  end
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
for n = 1:5, m = 1:5
    println("$n * $m = $(n * m)")
end

arr = rand(1:10, 10)
println(arr)
searched = 4
for (ix, curr) in enumerate(arr)
    if curr == searched
        println("The searched elemente $searched occurs on indexs $ix ")
        break
    end
end

for n in 1:10
    if 3 <= n <= 6
        continue #skip current iteration
    end
    println(n)
end

#Syntax error usually result in LoadError
if code in codes
    println( "This is an accdeptable code")
else
    throw(DomainError())
end

rethrow() is useful to hadn the current exception toa higher calling code level
warn("Something is not right here")

info("Did you know htis?")

type CustomException <: Exception end

#Dangerous code comes in the try bl9ock
#try-catch-finally
#You do htis so that your programs can contintue to run even if the do something bokners

a = []
try
    pop!(a)
catch ex
    println(typeof(ex))
    showerror(STDOUT,ex)
end

ret = try
    a = 4 * 2
catch ex
end

f = open("fi1.txt") #returns an IOStream(<file file1.txt)
try
    #operate on file f
catch ex
finally close(f)
end
x =
9
function funscope(n)
    x = 0 # x ins in the local scope of the function
    for i = 1:n
        x # x is local to the for loop
        x = i + 1
        if (x == 7)
            println("This is the local x in  for : $x ") # = > 7
        end
    end
x
println("this is the local x in fucnscope: $x ") # = > 0
 x = 15
end

funscope(10)
println("This is th eglboal x: $x") #

#if you need local bindings for a variable, use the let block
anon = cell(2)
for i = 1:2
    anon[i] () -> println(i)
    i += i
end
begin
    local x = 1
    let local x = 2
        println(x) # 2
    end
    x
    println(x) # 1
end
begin
    local	x	=	1
    let
        local	x	=	2
        println(x)
    end
    x
    println(x)	#>	1 end
end

# for loops update a variable
i = 0
for i = 1:10
end
println(i)

#comprehensions creat a temporary one
i = 0
[i for i = 1:10]
println(i)

function fib_producer(n)
    a, b = (0,1)
    for i = 1:n
        produce(b)
        a, b = (b, a + b)
    end
end
fib_producer(5)
tsk1 = Task( () -> fib_producer(10))
consume(tsk1)

tsk1 = @task fib_producer(10)

# Chapter 5
Collection types!

c = [1,2,3]
d = [1 2 3]

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

matrix = [1 2; 3 4]

[1 2 ] * [3 ; 4]
Array{Float64, 2}

ma1 = rand(3,5)
ndims(ma1)

size(ma1)

size(ma1,1)

size(ma1,2)

length(ma1)

idm = eye(3)

idm[1:end, 2]

idm[:, 2]

idm[2, :]

idmc = idm[2:end, 2:end]

idm[2, :]

idmc = idm[2:end, 2:end]

idm[2,:] = 0

idm[2:end, 2:end] = [ 5 7 ; 9 11 ]

?sub

ma1 = rand(3,3)

ma2 = rand(3,3)

X = inv(ma1) * ma2

X

X = ma1 \ ma2

X = ma2 / ma1

Tuples are a fixed sized gorup of values separated by commans and optionally surrounded by parentheses
tuples are heterogenous containers
arrays are homogenous containers
('z', 3.13) :: (Char,Float64)

t3 = (5, 6, 7, 8)

for i in t3
    println(i)
end

a, b = t3

a

b

# Dictionariessssssssssss
Keys are unique
Basically a collection of tw-element tuples of the form (key,value)
d1 = [1 => 4.2, 2=> 5.3]
d2 = ["a" => 1, (2,3) => true]

d3 = [:A => 100, :B => 200]

d3[2]
dmus = [:first_name => "Luis", :surname => "Armstrong", :occupation => "musician", :date_of_birth => "4/8/1901"]

in((:Z, 999), d3)

(:Z, 999) in d3

in((:A, 100), d3)

(:A, 100) in d3

d3[:A]
d4 = Dict()

d5 = Dict{Float64, Int64}()

d5["c"]

delete!(d3, :B)

ki = keys(d3)
for k in keys(d3)
    println(k)
end

:A in keys(d3)

haskey(d3,:A)

collect(keys(d3))

d5 = ["JS Bach " => 1685, "Woody Allen" => 1935, "Barack Obama" => 1961]

for (k,v) in d5
    println("$k was born in $v")
end

dpairs = ["a", 1, "b", 2, "c", 3]

d6 = [dpairs[i] => dpairs[i+1] for i in 1:2:length(dpais)]

d6 = Dict{String => Int64}[dpairs[i] => dpairs[i+1] for i in 1:2:length(dpairs)]

function showfactor(n)
    d = factors(n)
    println("factors for $n")
    for (k,v) in d
        print("$k^$v\t")
    end
end

@time showfactor(3752)

@apropos factor

arrkey = [key for (key, value) in dict]
#this is the same as collect(keys(dict))
arrval = [value for (key,value) in dict]
#htis is the same as collect(values(dict))
arrn = [dict[i] for i = 1:10]

arrn = map((i) -> dict[i], [1:n])

d1 = Dict( 1 => 4.2, 2=> 5.3)

d2 = Dict{Any,Any}("a" => 1, (2,3) => true)

capitals = Dict{String, String}("France" => "Paris", "China" => "Beijing")

d5 = Dict(zip(keys1, values1))

# Dictionaries have keys that are unique, but values are not, and the keys are not ordered
# Sets are collections where the order does not matter
# Checking set ∈ in a set is independent of size
x = Set([1:100])
@time 2 in x
x2 = Set([1:1000000])
@time 2 in x2

using JLD
words1 = "To be or not to be, that is the question!"
@save "words1.txt" words1
str = readall("words1.txt")
nonalpha = r"(\W\s?)"
test = replace(words1, nonalpha, ' ')

digits = r"(\d+)"
test2 = replace(words1, digits, ' ')
word_list = split(words1, ' ')
word_freq = Dict{String, Int64}()
for word in word_list
    word = strip(word)
    if isempty(word) continue end
    haskey(word_freq, word) ?
        word_freq[word] += 1 :
        word_freq[word] = 1
    end
println("Word : freq \n")
words = sort!(collect(keys(word_freq)))
for word in words
    println("$word : $(word_freq[word])")
end
###############################
#Chapter 6. More on Types, Methods, and
#Modules
###############################

#Type annotations and conversions
(31+42) :: Float64

# Type conversions and promotions
convert(Int64, 7.0)

# In general, convert(Type, x) tries to shove x into an instane of Type.
# generally, type(x) also works.
Int64(7.0)
# But it doesn't work with 1. lost precision, 2. incompatible source value
# We can also define our own conversions and promotions

# Because of fallback promotion, Julia is super clever and instead of
# defining the + operator for shittons of combos, you can just
(x:: Number, y :: Number) = +(promote(x,y)...)
promote(1.5, im)
promote(true, 'c', 1.0)
#  promote_rule(::Type{Int8}, ::Type{Uint16}) = Int64.
# check out shit at base/promotion.jl

# TYYYYYYYYYPE Hierarchy
1.07^10
1.10^7

TYPE    Instantiate     Subtypes
Concrete    Y               N
Abstract    N               Y

Abstract such as Number or Real. Groups multiple subtypes together.
Abstract types cannot have any fields.
Any is the supertype of all types, and all the objects are also instances of Any
All types are supertypes of None and no object is an instance of None.
Now, None is Union() and Nothing becomes Void

# USER DEFINED AND COMPOSITE types
type Point
    x :: Float64
    y :: Float64
    z :: Float64
end

# Point is a concrete type.

p1 = Point(1,2,3)
p1
typeof(p1)

Composite type is most similar to a struct in C or a class without methods in Java

Instead of Python or Java's object.func(args), Julia uses function(object,args)

Julia has no classes.

Data and functions are separate.

Methods are not typed to a single type, because multiple dispatch connects them with different types.

When adding a new method for a type, you don't have to change the code of the type itself, as with other oolangs.

names(p1)

fields(p1)
methods(point)
methods(Point)
p1.y

p1.z = 3.13
p1.z = "A"
Objects as arguments to functions are passed by REFERENCE.
if you don't want your objects to be changeable, use immutable.

immutable Vector3D
    x::Float64
    y::Float64
    z::Float64
end
p = Vector3D(1,2,3)
Immutables are thread safe!
Hugely optimized for Performance.
p
p.x = 3.4

is(p, Vector3D)
p == Vector3D
abstract type Employee end
type Developer <: Employee
    name :: String
    iq
    favorite_lang :: String
end

Employee() # not going to work because yo ucan't instantiate an abstract type
The type Developer has two implicit constructors, but we can define another *outer constructor*
Developer(name, iq) = Developer(name,iq, "Java")
Outer constructors help default object creation.
devel1 = Developer("Bob", 110)

devel2 = Developer("William", 145, "Julia")
type Manager
    name :: String
    iq
    department :: String
end

man1 = Manager("Julia", 120, "ICT")
function cleverer(d:: Developer, e:: Employee)
    println("The developer $(d.name) is cleverer I think!")
end

cleverer(devel1, devel2)
cleverer(devel1, man)

function cleverer(e:: Employee, d:: Developer)
    if e.iq <= d.iq
        println("The dev $(d.name) is cleverer!")
    else
        println("The employee $(e.name) is cleverer!")
    end
end

function cleverer(d1:: Developer, d2:: Developer)
    if d1.iq <= d2.iq
        println("The developer $(d2.name) is cleverer!")
    else
        println("The deeveloper $(d1.name) is clverer!")
    end
end

cleverer(devel1, devel2)

cleverer(devel2, devel1)

TYPES and Collections - inner constructors
type Person
    firstname :: String
    lastname :: String
    sex :: Char
    age :: Float64
    children :: Array{String, 1}
end
p1 = Person("Alan", "Bates", 'M', 45.5, ["Jeff", "Stephan"])
people = Person[]
push!(people, p1)

push!(people, Person("Julia", "Smith", 'F', 27, ["Viral"]))
import Base.fullname
fullname(p::Person) = "$(p.firstname) $(p.lastname)"

show(people)

fullname(p::Person) = "$(p.firstname) $(p.lastname)"

Or the much better performant

fullname(p::Person) = string(p.firstname, " ", p.lastname)
print(fullname(p1))
type Family
    name :: String
    members :: Array{String, 1}
    big :: Bool
    Family(name :: String) = new(name, String[], false)
    Family(name :: String, members) = new(name, members, length(members)> 4)
end

fam = Family("Bates-Smith", ["Alan", "Julia", "Jeff", "Stephan", "Viral"])
new is exclusive to inner constructors

When you use inne constructors, there are no default constructors anymore. Outer Constructors calling a limited set of inner constructors is often
the best practice.

TYPE UNIONS
================
type VPoint
    x:: Float64
    y:: Float64
end
type Vector2D
    x ::Float64
    y ::Float64
end

p = VPoint(2.0,5.0)
v = Vector2D(1.0,3.0)

+(p, v)
import Base.+
+(p :: VPoint, q :: VPoint) = VPoint( p.x + q.x, p.y + q.y )
+(u :: Vector2D, v :: Vector2D) = VPoint(u.x + v.x, u.y + v.y)
+(u :: Vector2D,  p :: VPoint) = VPoint(u.x + p.x, u.y + p.y)
+(p:: VPoint, v:: Vector2D) = VPoint(p.x + v.x, p.y + v.y)

+(p,v)

VecOrPoint = Union{Vector2D, VPoint}

isa(p, VecOrPoint)

isa(v, VecOrPoint)

+(u :: VecOrPoint, v :: VecOrPoint) = VecOrPoint( u.x + v.x, u.y + v.y)
This is amazing because we just skipped defining 4 methods and did just 1.

PARAMETRIC TYPES AND METHODS
============================
You can have...
Array{Int64,1},
Array{Int16,1},
Array{Int32,1},
Array{ASCIISTRING,1},
Array{Float64,1}...

Array is thus a Parametric type
Array{T, 1}

In general, types can take type parameters, so that declarations can take a whole bunch of values

type TPoint{T}
    y :: T
    x :: T
end

Now, all of these are subtypes of TPoint.
issubtype(TPoint{String}, TPoint)

HOWEVERRRRRRRRRRRR
issubtype(TPoint{Float64}, TPoint{Real})
This is called type contravariance.

You can indicate the T in the constructor...

test = TPoint{Int64}(2,5)
or shorten to...
test = TPoint(2,5)

now for example,

teststring = TPoint("London", "Great-Britan!")

Now, we want to restrict the parameter type to T and only subtypes of Real...

type RPoint{T <: Real}
    x :: T
    y :: T
end

now if we try this next line it should error.

test2 = RPoint("London", "Great-Britain!")

Just as well, methods can have type parameters (in essence, restrictions on the mathematial domain)
add{T}(x:: T, y :: T) = y + x

add(2,3)

add(2,3.0)

add{T <: Number}(x :: T, y :: T) = y+x

function vecfloat{T <: AbstractFloat}(x :: Vector{T}) end

PAckages/Libraries are usually contained in a module, whose name starts with uppercase letter.

module Package1
println("FMLMLIA")
end

import Winston, Gadfly
Winston.plot(rand(4))
Gadfly.plot(x=[1:10],y=rand(10))

The global scope is in the Main module
The Core Module contains all built-in identifiers and is always available.
The standard library is also available (all the contents of /base is cointained in
    Base, Pkg, Collections, Graphics, Test, Profile)
A Module has a type: Module
typeof(Base)

#This tells you all the vars you have defined!
names(Main)

methods(names)

@doc names
#whos() tells you all the vars used and the corresponding types and allocation in memory
whos()

whos(Colors)

You can make internal definitions such as constants, variables, types, functions and so on visible to other modules (make them public)
by declaring them with export.

e.g.
export Type1, perc

If module LibA is defined in modules_ext.jl, then require("modules_ext.jl") will load this in the current code.
using LibA ---- Makes its exported names available in the current namespace.

Only the exported functions remain visible. The others remain private.

eg
# The below module save in file temperature_converter.jl in same directory
# module TemperatureConverter
#     export as_celsius
#
#     function as_celsius(temperature, unit)
#         if unit == :Celsius
#             return temperature
#         elseif unit == :Kelvin
#             return kelvin_to_celsius(temperature)
#         end
#     end
#
# function kelvin_to_celsius(temperature)
#     # private function
#     return temperature + 273
# end
#
# end
println("$(as_celsius(100,:Celsius))")
println("$(as_celsius(100, :Kelvin))")
println("#(kelvin_to_celsius(0))")

IN GENERAL, there are different ways of importin definitions from another module LibA.

1. using LibA --- LibA gets searched for exported definitions if needed. You can use funcitons in LibA without qualifying it.
2. using LibA.varB or, using LibC: varC, funcC
3. import LibD.funcD only imports one name na dcan also be useed if funcD was not exported.
    now must use funcD as LibD.funcD.
4. importall LibE --- import all the exported names in LibE

Imported vars are read-only, current module cannot create vars with the same names as the imported ones.
A source file can cointain many modules, or one module can be defined in several source files.
if a module contains a function __init__(), this will be executed when the module is first loaded.

include("file.jl") is, to the compiler, no different to copying file1.jl and pasting it directly in the current file or the REPL.

The variable LOAD_PATH contains a list of dirs where Julia looks for module files when running using, import, or include.
    you can exted it with push!(LOAD_PATH, "new/path/to/search")
 modules are compiled on load, which slows shit down if not for precompilation.

CHAPTER 7 - METAPROGRAMMING IN JULIAAAAAAAAAAAAA
================================================

Everything is an expression that returns a value when executed.
Everypice of the code is internally represented as an ordinary Julia data structure, also called an expression.
By working on expressions, a Julia Program can transform and generate new code.

It inherits homoiconicity from lisp: "code is data and data is code"
1. Expressions and symbols
2. eval and interpolation
3. defining macros
4. built-in macros
5. reflection capabilities

Abstract syntax trees and grammars...
an expression is simly an object that represents Julia code.
eg
2 + 3 is

    +
  2    3

To make Julia see this as an expression and block its evaluation, we have to quote it.
Or rather, to preced it by a colon(:) like
:(2 + 3). When you eval :(2 + 3), it just returns :( 2 + 3)
typeof( :(2 + 3))

quote
    a = 42
    b = a^2
    a - b
end

We can give expressions a name
d1 = :(2+3)

d1.head

d1.args

d1.typ

e2 = :(2 + a * b - c)

e2.head

e2.args

e2.typ

expressions work recursively.
We can simplify every subexpression in the same way until we end up with elementary symbols and literals.
In the content of an expression, symbols are used to singnify access to variables.
They represent the variable in the structure of the tree code.
The prevent evaluation chatacter is :.
:x returns :x

dump(d1)

dump() helps to show the args, type, and head of an expression. (works way better in REPL.
Eval and interpolation.

Because we know th etype of expressions, we can build new ones like this:

e3 = Expr(:call, :*, 3, 4 )

and this should be allowed.
We can evaluate this expression with the function eval()
eval(e4)

but we should be careful that these expressions are dfined when evaluated in the global scopeeg
f1 = :(3 - zzz)
eval(f1)

which should error because zzz was not defined. Similarly, we can have evals that change the current state of the global scope.
f2 = eval(:(b = 1))

To make writing expressions easier, we can interpolate with the $ sign.
It evaluates immediately what the expression means.
a = 4
b = 1

e4 = :(a+b)

e5 = :($a+b)

There are two kinds of evaluations here:
1. expression interpolation with the $ --- evaluates when the expression is consturcted (at parse time)
2. Quotation ( with : or quote) evaluates only when the expression is passed to eval at runtime.

This is huge for reducing boilerplates in external libraries, or to generate lots of similar routines in other situations.
    In robotics, being able to generate a program and then run it is very useful.

Functions on expressions

function math_expr(op, op1, op2)
    expr = Expr(:call,op, op1,op2)
    return expr
end

ex = math_expr(:+, 1, Expr(:call, :*, 4,5))

eval(ex)

And now an example that doubles any numeric argument but leaves the expressions alone.

function make_expr2(op, opr1, opr2)
    oprf1, oprf2 = map(x -> isa(x, Number) ? 2*x : x, (opr1, opr2))
    retexpr = Expr(:call, op, oprf1, oprf2)
    return retexpr
end

make_expr2(:+, 1, 2)

ex = make_expr2(:+, 1, Expr(:call, :*, 5,8))
eval(ex)

We move on to macros!

MACROOOOOOOOOOOOS
=================
like functions, but take expressions instead of values.
The expansion occurs at parse time when the syntax tree is being built, not when the code is actually executed.
A macro maps a tuple of arguments to a return expression, and the resulting expression is compiled directly rather than
requiring a runtime eval().

macro sayhello()
    return :( println("Hello world!") )
end

@sayhello

FUNCTION: it takes the input values and returns the computed values at runtime.
MACRO: It takes the input expressions and returns the modified expressions at parse time.

ie, a macro is a custom program transformation.
macros are defined as such
macro mname()
    #code
    println("YO!")
end
: is used for
1. creating a symbol via :x, :(:), :foo == Symbol("Foo"), or Symbol(:var, '_', "sym")
2. creating expression with :(2+3), x = Expr(:call, :+, 2,3)
ex = quote
    foo
    bar
    baz
end

@mname

@mname sum

macro macint(ex)
    quote
        println("start")
        $ex
        println("after")
    end
end

@macint println("Where am I?")
macro assert1(ex)
    :($ex ? nothing : error("Assertion failed: ", $(string(ex))))
end
@assert1 1 == 1

@assert 1 == 42

REMEMBER macroexpand!
macroexpand(:(@assert 1 == 42))
macro unless(test, branch)
    quote
        if !$test
            $branch
        end
    end
end

arr = [3.14, 42, 'b']

@unless 42 in arr println("arr does not contain 42")

macroexpand(:(@unless 41 in arr println("arr does not contain 41")))

quote
    if !(41 in arr)
        println("arr does not cointain 41")
    end
end
arr = ["a", "b", "c"]

macro convarr(arr, T)
    :(reshape($T[$arr...], size($arr)...))
end

@convarr arr Symbol

Because you are injecting code into the module they are called, gotta watch nameclashes.
Hygienic macro respect this.
    1. Declare vars in the macro as local, so as to not conflict with outer variables.
    2. Use the escape function to make sure that an interpolated expression is not expanded, but instead is used literally.
    3. Don't call eval inside a macro (likiely the vars you are evaluating don't even exist at that point)

macro timeit(ex)
    quote
        local t0 = time()
        local val = $(esc(ex))
        local t1 = time()
        print("elapsed time in seconds: ")
        @printf "%.9f" t1 - t0
        val
    end
end

@timeit factorial(10)
@timeit a^3

Gotta differentiate between calling context and macro context.

Allow you to create DSLs.

Built-in macros
===============

TESTING!
@assert
using Base.Test
@test 1==3

@test 1 ≈ 1.1

DEBUGGING!
@which sort(arr)

@show

BENCHMARKING!
REMEMBER @timed
@timed [x^2 for x in 1:1000]
@elapsed [x^2 for x in 1:1000]

@allocated [x^2 for x in 1:1000]

TASKS!
a = @async 1+2
consume(a)

REFLECTION CAPABILITIES!
Code is represented by expressions that are datastructures of type Expr.
The structure of a program and its types can be explored programmatically, just like any other data.
    This means a running program can dynamically discover its own properties, which is called REFLECTION.

Some functions that do this are...
typeof(), subtypes(),
methods(),
names and types.

eg

type person
    name :: String
    height :: Float64
end

fields(person)
code_lowered(+, (Int, Int))
code_typed(+, (Int, Int))

CHAPTER 8 - I/O, Networking, ParallelComputing
==============================================
1. Basico IO
2. working with CSVs
3. using DataFrames
4. working with TCP sockets and servers
5. Interacting with databases
6. Parallel operations nad ocmputing

Basico IO
read(STDIN, Char)

write(STDOUT, "Julia")
STDIN and STDOUT are simply streams
readbytes(STDIN,3)

REMEMBER eachline
If you need to read all the lines in a n input stream, remember to use the readline.

stream = STDIN
for line in eachline(stream)
    print("Found $line")
    print("processed")
end

we use eof(stream) with a while loop...

while !eof(stream
    x = readline(line, Char)
    println("Found $x")
    #process
end

WORKING WITH FILES
To work with files, we use the IOStream type.
it has the supertype IO and the following characteristics...
    1. fields given by

fields(IOStream)
Array{Symbol,1}: :handle :ios :name :mark

2. the types are given by IOStream.types

IOStream.types

The file handler is a pointer of the type Ptr, which is a reference to the file object.

touch("empty.txt")

fname = "empty.txt"
f1 = open(fname)
typeof(f1)

data = readlines(f1)

for line in data
    print(line)
end

close(f1)

REMEMBER do format
open(fname) do file
    process(file)
end

Writing a file requires the "w" flag, then writing to it with print, println or write, then closing.

type Wine
    fixed_acidity :: Array{Float64}
    volatile_acidity :: Array{Float64}
    citric_acid :: Array{Float64}
    quality :: Array{Float64}
end

#wine1 = Wine(data[1,:]...)
The elements of the row are splatted with the ... operator into the Wine constructor.

using DataFrames

dv = @data([7, 3, NA, 5, 42])
df = DataFrame()
df[:Col1] = 1:4
df[:Col2] = [ e, pi, sqrt(2), 42]
df[:Col3] = [ true, false, true, false]
show(df)

Julia can work with other human-readable file formats through specialized packages:
For JSON, use the JSON package. The parse method converts the JSON strings into
Dictionaries, and the json method turns any Julia object into a JSON string.
For XML, use the LightXML package
For YAML, use the YAML package
For HDF5 (a common format for scientific data), use the HDF5 package
For working with Windows INI files, use the IniFile package

We want to be able to handle multiple clients!

server = listen(8081)
while true
    conn = accept(server)
    @async begin
        try
            while true
                line = readline(conn)
                println(line)
                write(conn, line)
            end
        catch ex
            print("connection ended with error $ex")
        end
    end
end

The listen function has some variants, for example, listen(IPv6(0),2001) creates a TCP
server that listens on port 2001 on all IPv6 interfaces. Similarly, instead of readline, there are
also simpler read methods:
read(conn, Uint8): This method blocks until there is a byte to read from conn, and then
returns it. Use convert(Char, n) to convert a Uint8 value into Char. This will let you see
the ASCII letter for Uint8 you read in.
read(conn, Char): This method blocks until there is a byte to read from conn, and then
returns it.

DATABASESSSSSSSSSS
==================
Open Database Connectivity

PARALLEL OPERATIONS AND COMPUTINGGGGGGGGG
=========================================
Julia implements the actor model like Erlang, Elixir, or Dart.

Creating Processes
julia -p n to start a REPL with n workers
These are diferent processes, not threads, so they do not share memory.

workers()

for pid in workers()
    print("sup $pid here!")
end
addprocs(1)

workers()
myid()

nprocs()

All these workers communicate via TCP ports and run on the same machine, which is why it is
called a local cluster. To activate workers on a cluster of computers, start Julia as follows:
julia --machinefile machines driver.jl
Here, machines is a file that contains the names of the computers you want to engage, like this:
node01
node01
node02
node02
node03
Here node01, node02, and node03 are the three names of computers in the cluster, and we want
to start two workers each on node01 and node02, and one worker on node03.
The driver.jl file is the script that runs the calculations and has the process identifier 1. This
command uses a password-less SSH login to start the worker processes on the specified
machines. The following screenshot shows all the eight processors on an eight core machine
when engaged in a parallel operation:

REVIEW THE CLUSTERMANAGER

USING LOW-LEVEL COMMS

f(x) = x^2

r1 = remotecall(x -> x^2, 2, 1000)

fetch(r1)

r2 = @spawnat 2 sqrt(2) #lets worker 2 calculate sqrt(2)
fetch(r2)
r3 = @spawn sqrt(5)

typeof(r3)

fetch(r3)

#TO EXECUTE ON ALlLL THE WORKERS
r = [@spawnat w sqrt(5) for w in workers()]

fetch(r[1])

To execute the same statement on all hte workers, we can use the @everywhere
@everywhere println(myid())

x = 5

@everywhere println(x)
@everywhere w = 8

@everywhere using Colors

@everywhere function fib(n)
    if (n<2) then
        return n
    else return fib(n-1) + fib(n-2)
    end
end

TIP!
Separate code into files: functions.jl as the functions and parameters tha need to be run in parallel,
and the driver.jl that manages the processing and collecting the results.

Other alternative is to separate the load files on the command line

julia -p n -L file2.jl -L file2.jl driver.jl

reduce data movement to improve performance.

if every worker needs a var, use a global.

for pid in workers()
    remotecall(pid, x -> (global d; d = x; nothing), d)
end

PARALLEL LOOPS AND MAPS
=======================
function buffon(n)
hit = 0
for i = 1:n
    mp = rand()
    ϕ = (rand() * π) - π/2
    xright = mp + cos(ϕ)/2
    xleft = mp - cos(ϕ)/2
    p = (xright >= 1 || xleft <= 0) ? 1 : 0
    hit += p
end
miss = n - hit
πapprox = n /hit * 2
end

@time buffon(10_000_000)

function buffon_par(n)
    hit = @parallel (+) for i = 1:n
    mp = rand()
    ϕ = (rand() * π) - π/2
    xright = mp + cos(ϕ)/2
    xleft = mp - cos(ϕ)/2
        (xright >= 1 || xleft <= 0) ? 1 : 0
    end
    miss = n - hit
    πapprox = n/ hit * 2
end

@time buffon_par(10_000_000)
addprocs(1)

FORMAT:
@parallel ("reducing operator") "the for loop"
Calculations must be independent.

So that means this will fail...

arr = zeros(10_000_000)
@parallel for i = 1: 10_000_000
    arr[i] = i
end

because each processor has a copy of it.

If the computational task is to apply a function to all elements in some collection, you can use a
parallel map operation through the pmap function. The pmap function takes the following form:

pmap(f, coll),

applies a function f on each element of the collection coll in parallel, but
preserves the order of the collection in the result. Suppose we have to calculate the rank of a
number of large matrices. We can do this sequentially as follows:

function rank_marray()

    marr = [rand(1000,1000) for i=1:10]

    for arr in marr
        println(rank(arr))
    end
end

function prank_marray()
    marr = [rand(1000,1000) for i = 1:10]
    println(pmap(rank,marr))
end

@time prank_marray()

These are typical map-reduce problems.

DISTRIBUTED ARRAYS

@everywhere using DistributedArrays

arr = drand((100,100), workers()[1:1], [1,1])

whos(arr)

arr.pmap

arr.indexes

arr.cuts

arr.chunks

da = @parallel [2i for i = 1:10]

DArray((10,100)) do I
    println(I)
    return zeros(length(I[1]), length(I[2]))
end

Julia's model for building a large parallel application works by means of a global distributed
address space. This means that you can hold a reference to an object that lives on another machine
participating in a computation. These references are easily manipulated and passed around
between machines, making it simple to keep track of what's being computed where. Also,
machines can be added in mid computation when needed.

Chapter 9 - Running External Programs
=====================================

Running shell stuff
C and Fortran
Calling Python
Performance tips

RUNNING SHELL COMMANDS
----------------------

pwd()
cd("d:\\test\\week1")
; ls
; mkdir folder
; cd folder

Wanna run a shell command? Use the run() function!

cmd = `echo "Julia is smart"`
typeof(cmd)
@doc Cmd
run(cmd)
run("date")
@doc run

1.02^300

Julia forks commands as child processes from the Julia process. Instead of immediately running
the command in the shell, backticks create a Cmd object to represent the command, which can then
be run, connected to other commands via pipes, and read or write to it.

INTERPOLATION
file = "file1.txt"
cmd = `cat $file`
run(cmd) # prints the ocntents of file1.txt

PIPELINING
Use the |> symbol!

run(`cat $file` |> "test.txt")

and you can chain it thus

run(`echo $("\n hi \n Julia")` |> `cat` |> `grep -n J`)

Another example is to search for the word "is" in all the text files in the current folder; use the
following command:
run(`grep is $(readdir())`)

To capture the result of a command in Julia, use readall or readline:
a = readall(`cat "tosort.txt"` |> `sort`)

Now a has the value "A\r\nB\r\nC\n".

Multiple commands can be run in parallel with the & operator:
run(`cat "file1.txt"` & `cat "tosort.txt"`)
This will print the lines of the two files intermingled, because the printing happens concurrently

CALLING C AND FORTRAN
---------------------

For calling C, we need to work with pointer types.
A native Ptr{T} is the memory address for a variable of type T.
The term bitstype is also used. This is for concrete type - data consists of bits:
    Int8, Uint8, Int32, Float64, Bool, Char
To pass a string to C, it is converted to a contiguous type byte array represenatiton with the function
    bytestring(); given Ptr to a C string, it returns a Julia string.

typeof(Ptr)
@doc Ptr

lang = ccall( (:getenv, "libc"), Ptr{UInt8}, (Ptr{UInt8}, (Ptr{UInt8},), "LANGUAGE"))

@doc ccall

Base.Libdl.find_library(["libc"])

When calling FORTRAN, all inputs must be passed by reference.
Arguments to C are in general automatically converted.
The returned values in C types are also converted to Julia types.
Arrays of Booleans are handled differently in C and J, so gotta handle manually.
Some system dependent types applies too.

ccal() auto ensures args are preserved from GC until call returns.

C types are mapped to JUlia types.
short => Int16
double => Float64

Embedding is possible

you can use C with Callbacks

If you have existing C code, you must compile it as a shared lib to call it from Julia.
With GCC, use =shared -fPIC
Support for Cpp is through Cpp and Clang

Calling Python
==============

using PyCall
pyeval("10*10")

@pyimport math

math.sin(math.pi / 2)

PERFORMANCE TIPS!
=================

1. Don't use globals, or make them constant.
2. Structure code around functions that work on local vars via the funciton args - this returns the results rather than mutating global objects
3. Type stability. Don't change types over time.
4. If you need large arrays, use a sizehint from the start.
5. if arr is very large and you don't need it, kill it.
6. sometimes disabling gc via gc_disable() can be useful.
7. Anonymous funcitons yayyyy
8. use small functions.
9. code different versions of a funciton according to the types, so that multiple dispatch applies.
10. use types for keyword arguments, avoid the splat ... for dynamic lists of keyword args.
11. using mutating APIS like push!() is useful to avoid copying.
12. prefer array operations to comprehensions.
13. don't use a try/catch in a tight loop.
14. Use immutables
15. avoid Any
16. Avoid abstract types in a collection
17. don't copy unless necessary, avoid using a large number of variables or temp arrays or collections
18. devectorize your code
19. use @parallel if appropriate
20. reduce data movement between workers in parallel execution
21. fix deprecation warnings
22. use the @inbounds macro
23. don't use eval at runtime

Tools to use
============

check speed and allocations with @time

linter tool is sweet for performance profiling

consider using TypeCheck

using ProfileView

CHAPTER 10 - THE STDLIB and Packages
====================================

Stdlib, package manager, publish a package, graphics in Julia, GAdfly on data.asdf

THE STDLIB
----------
filter(f, coll) applies f to all the elements of coll

fold, min, max, sum, any, all, map, filter

filter(x -> iseven(x), 1:10)

mapreduce(f, op, coll) applies f to all elements in col, then reduces with op()
so
mapreduce(x -> sqrt(x), +, 1:10) == sum(map(x -> sqrt(x), 1:10))

The pipeline operator lets you write very functionally styled code
x |> f applies f to the function argument x.

1:10 |> (x -> sqrt(x)) |> sum

REMEMBER clipboard()
a = 42

clipboard(a)

GRAPHICS!
=========
using Gadfly

draw(SVG("gadfly.svg", 6inch, 3inch), plot([x -> x^2], 0, 25))
p1 = plot([x -> cos(x)/x], 5, 25)
draw(PNG("gadfly.png",300,100), pl)

APPENDIX AAAAAAAAAAAAAAAAAAAAAA
===============================
Useful.

===============================================================
PART II. Module 2
===============================================================

How fast can Julia be?

function pisum()
    sum = 0.0
    for j = 1:500
        sum = 0.0
        for k = 1:10_000
            sum += 1.0/(k*k)
        end
    end
    sum
end

Analyzing Julia Performance

@timev sqrt.(rand(1000))

The Julia Profiler
inspects call stack every 1ms or 10ms, counts how often each line of ocde is encountered,
    helps prioritize optimization efforts.

using Base.Profile

function testfunc()
    x = rand(10000,10000)
    y = std(x,1)
    return y
end

testfunc()

@profile testfunc()

Profile.print()

Profile.init()

There are a few profiler options that are sometimes useful, although the defaults are a good
choice for most use cases. Primary among them is the sampling interval. This can be provided as
keyword arguments to the Profile.init() method. The default delay is 1 millisecond on Linux,
and should be increased for very long-running programs through the following line of code:
julia> Profile.init(delay=.01)
The delay may be reduced as well, but the overhead of profiling can increase significantly if it is
lowered too much.
Finally, you may have realized that the profiler stores its samples in memory to be viewed later.
In order to profile a different program during an existing Julia session, it may be necessary to
clear the stored profile from memory. The Profile.clear() function does this, and must
therefore be run between any two invocations of @profile within the same Julia process

using ProfileView

ProfileView.view()

test()
@profile test()
ProvileView.view()

Naaaais

Using the Memory Allocation tracker
track -alloation=user

Profile.clear_malloc_data()

To get Julia to track memory allocation, start the julia process with the –trackallocation=user
option as follows:
julia> track -allocation=user
This will start a normal Julia session in which you can run your code as usual. However, in the
background, Julia will track all the memory used, which will be written to .mem files when Julia
exits. There will be a new .mem file for each .jl file that is loaded and executed. These files
will contain the Julia code from their corresponding source files, with each line annotated with
the total amount of memory that was allocated as a result of executing this line.
As we discussed before, when running Julia code, the compiler will compile user code at
runtime. Once again, we do not want to measure the memory allocation due to the compiler. To
achieve this, first run the code under measurement once, after starting the Julia process. Then run
the Profile.clear_malloc_data() function to restart the allocation measurement counters.
Finally, run the code under measurement once again, and then exit the process. This way, we will
get the most accurate memory measurements.

CHAPTER 3 - Types in Julia
types are tag-on values that restrict the range of potential values

foo(x::Integer) = "an integer"
foo(x::String) = "a string"
function bar(a, b)
    x :: Int64 = 0
    y = a + b
    return y
end

foo(1)

foo("1")

foo(1.5)

Dispatch is the process of selecting a funciton to be executed at runtime
Dispatch is a runtime process, method overloading is a compile-time concept.
In traditional OOP, dispatch at runtime occurs only on the receiver of the method ie., the object before the dot -- and hence
why it is called "single dispatch".

Julia progs, ∴, cointain many small function definitions for differnt types of arguments.
God pratice is to constrain arg types to the widest level possible.
Use only tight methods when failure is assured with wider types.

consider:

sumsqr(x,y) = x^2 + y^2

sumsqr(1,2)

sumsqr(1.5,2.5)

sumsqr(1+2im, 2+3im)

ABSTRACT TYPES
Abstract types cannot have any instantiated values == Cabnn only be nodes, not leaves of the type hierarchy.
Represent sets of related types.
Int32 and Int64 both inherit from the Integer abstract type

abstract Number

abstract Real <: Number

abstract FloatingPoint <: Real

abstract Integer <: Real

abstract Unsigned <: Integer

Number has no explicit super type - thus a subtype of Any.

Concrete types can be instantiated to values.
Concrete types cannot have any subtypes

COMPOSTIE and immutable types
=============================

type Pixel
    x :: Int64
    y :: Int64
    color :: Int64
end

p = Pixel(4,4,100)

p.x = 10

p.x

immutable IPixel
    x :: Int64
    y :: Int64
    color :: Int64
end

p = IPixel(4,4,100)

p.x = 10

type TPixel{T}
    x :: Int64
    y :: Int64
    color :: T
end

Parametros de tipos are the powerful shiznit.
This now accepts methods for TPixel{Int64}, or TPixel{String},
    and TPixel is now an abstract type itself.

Array{T,N} ==== T is a type and N is a numberrr

TYPE INFERENCE
==============
The type inference algo is a forward dataflow analysis - NOT the Hindley-Milner algo in
    ML langs.

check this
[x for x=1:5]

No type annotations, but we get an Int64.

TYPE STABILITY
==============
Be it. Use it.
Return value of a funciton is dependent only on the types of its arguments and not the values.

Thus type stability can continue across chains of function invocations wihtout actually running the code,
    even though the language is fully dynamic.

function mytrunc(x)
if x < 0
    return 0
else
    return x
end
end

mytrunc(-1)

mytrunc(-2.5)

mytrunc(2.5)

Works for Ints and Floats.

typeof(mytrunc(2.5))

typeof(mytrunc(-2.5))

@code_warntype mytrunc(-2.5)

function trunc_fixed(x)
    if x < 0
        return zero(x)
    else
        return x
    end
end

trunc_fixed(-2.4)

trunc_fixed(-2)

typeof(trunc_fixed(-2))

typeof(trunc_fixed(-2.4))

REMEMBER move the type variable part of the code into a function.

@code_llvm trunc_fixed(2)

@code_llvm trunc_fixed(0)
@code_llvm trunc_fixed(-2.5)

@code_llvm mytrunc(2)
@code_llvm mytrunc(0)
@code_llvm mytrunc(2.5)

using BenchmarkTools

@benchmark mytrunc(2.5)
@benchmark trunc_fixed(2.5)

Identifiying Type-Stability
===========================
look for Union returns and %jl_value_t% in the assembly code

LOOP VARS
===========================

Don't change types of vars from one iteration of the loop to the other.

function sumsqrtn(n)
    r = 0
    for i = 1:n
        r = r + sqrt(i)
    end
    return r
end

function fixsumsqrtn(n)
    r = 0.0
    for i = 1:n
        r = r + sqrt(i)
    end
    return r
end

@benchmark sumsqrtn(1_000_000)

@benchmark fixsumsqrtn(1_000_000)

KERNEL METHODS
==============

function string_zeros(s :: AbstractString)
    x = Array(s == "Int64" ? Int64 : Float64, 1_000_000)
    for i in 1:length(x)
        x[i] = 0
    end
    return x
end

@benchmark string_zeros("Int64")

But here, we are inferring in one part and filling a loop in another -
this calls for 2 functions.

function string_zeros_stable(s :: AbstractString)
    x = Array(s == "Int64" ? Int64 : Float64, 1000000)
    return fill_zeros(x)
end

function fill_zeros(x)
    for i in 1:length(x)
        x[i] = 0
    end
    return x
end

@benchmark string_zeros_stable("Int64")

TYPES IN STORAGE LOCATIONS
==========================

When defining a storage locations for data, it is importnat to specify a concrete type.
    So fo things that hold data, like arrays, dictionaries, or fields in composite types, it is
    best to explicitly define what type it will hold.

Arrays
======

experiment with 2 arrays

a = Int64[1,2,3,4,5,6,7,8,9,10]
bbb = UInt128[1,2,3,4,5,6,7,8,9,10]

function arr_sumsqr{ T <: Number}(x :: Array{T,1})
    r = zero(T)
    for i = 1:length(x)
        r = r + x[i] ^ 2
    end
    return r
end

@benchmark arr_sumsqr(a)

COMPOSITE TYPES
===============

immutable BPoint
    x
    y
end

Compiler does not know actual size, so there might be pointer misdirection.
waaaay better to do this.

immutable ConcretePoint
    x :: Float64
    y :: Float64
end

PARAMETRIC COMPOSITE TYPES ARE THE SHIZNIT
==========================================
What if we wanna work with Float32 and Float 16?

We could try....

immutable PointWithAbstract
    x :: AbstractFloat
    y :: AbstractFloat
end

But this would still mean some pointer misdirections, as compiler does not know
exact layout of the data in memory. Way better to try this.

immutable ParametricPoint{ T <: AbstractFloat}
    x :: T
    y :: T
end

When we write the type in this manner, our code remains generic. We can write our methods with
the confidence that the ParametricPoint type can hold values for any type of a floating point
number. Yet, at runtime, when an instance of this type is created, it is instantiated with a particular
type of float. In other words, once an instance is created, T becomes known. At this point, all the
benefits of specifying the concrete type discussed before are applicable. Both storage and type
inferences are efficient now.

CHAPTER 4 - FUNCTIONS AND MACROSSSSSSSSS
========================================

Write shit tons of small functions based on type.
Learn about the problem with globals.

The trouble with goals.
-----------------------
To keep track of all writes to global variables is akin to solving the halting problem.
Therefore, the data-flow algo sucks.

p = 2

function pow_array(x :: Vector{Float64})
    s = 0.0
    for y in x
        s = s + y^p
    end
    return s
end

t = rand(100000)

@benchmark pow_array(t)

@code_warntype pow_array(t)

const p2 = 2

function pow_array2(x :: Vector{Float64})
    s = 0.0
    for y in x
        s = s + y^p2
    end
    return s
end

@benchmark pow_array2(t)

in JULIA, a const var can change value, but not its type.

@code_warntype pow_array2(2)

INLINING
===========
As we've mentioned before, Julia code consists of many small functions. Unlike most other
language implementations, some of the core primitives in the base library are also implemented in
Julia. This means that the function call overhead has the potential to be a bottleneck in a Julia
program. This is mitigated using some aggressive inlining performed by the Julia compiler.
Inlining is an optimization performed by a compiler, where the contents of a function or method is
inserted directly into the body of the caller of that function. Thus, instead of making a function
call, execution continues directly by executing the operations of the callee within the caller's
body.
In addition, many compiler optimization techniques work within the body of a single function.
Inlining, therefore, allows many more optimization to be effective within the program.

Inlining is done by the Julia compiler, the rest is LLVM magic.

DEFAULT-INLINING
===================

Inline worthiness is done by Julia compiler.
Essentially, small functions with deterministic types.
However, it bloats code size and ∴ balance is needed.

testtrunc(x) = x < 0 ? zero(x) : x

function sqrt_sin(x)
    y = testtrunc(x)
    return sin(sqrt(y)+1)
end

Now we look at the processed AST after the compiler has run its type inference and inlining passes.

@code_typed sqrt_sin(1)

CONTROLLING INLINING
====================

REMEMBER @inline
needs to be in front of a function definition.

No call-site annotation to force inlining.

We define f(x) that performs some numeric ops on args, then g(x) calls f
after transforming its argument.

function f(x)
    a = x*5
    b = a+3
    c = a-4
    d = b/c
end

This is too long to be inlined by default.

g(x) = f(x)

@code_typed g(3)

@inline function f_inline(x)
    a = x*5
    b = a + 3
    c = a-4
    d = b/c
end

g_inline(x) = f_inline(2*x)

@code_typed g_inline(3)

@code_llvm g_inline(3)

Disabling inlining is useful and necessary for debugging.

REMEMBER -inline=no

CLOSURES AND ANONYMOUS FUNCTIONS
================================

. Anonymous functions, as the name suggests, are
functions without a name, and they are usually created at the point where they are passed in to
another function as an argument. In Julia, they are created with the -> operator separating the
arguments from the function body. These, and named functions created within the scope of another
function, and referring to variables from this outer scope, are called closures. This name arises
from the idea of these functions "closing over" the outer scope.

sqr(x) = x ^ 2

@benchmark map(sqr, rand(100000))

@benchmark map(x -> x ^ 2, rand(100000))

USING MACROS FOR PERFORMANCE
============================
macros can help do stuff in compile time and not in run time.

Raw Source
    Parse
        Macro Expansion - macroexpand()
            Lowering        - @code_lowered
                Type inference
                    Inlining
                        Generated Function expansion - @code_typed
                            Code Gen
                                Native Compilation      - @code_native
                                    RUN!

USING Macros

We gonna try evaluating a polynomial

function poly_naive(x, a...)
    p = zero(x)
    for i = 1:length(a)
        p = p + a[i] * x ^ (i-1)
    end
    return p
end

f_naive(x) = poly_naive(x, 1,2,3,4,5)

f_naive(3.5)

1 + 2*3.5 + 3*3.5^2 + 4*3.5^3 + 5*3.5^4

@benchmark f_naive(3.5)

REMEMBER peakflops()

peakflops()

function poly_horner(x, a...)
    b = zero(x)
    for i = length(a) : -1 : 1
        b = a[i] + b*x
    end
    return b
end

f_horner(x) = poly_horner(x,1,2,3,4,5)

@benchmark f_horner(3.5)

REMEMBER
The Horner macro
Improving the speed of this computation starts with realizing that the coefficients of the
polynomial are constants. They do not change and are known when writing the program. In other
words, they are known at compile time. So, maybe we can expand and write out the expression
for the Horner's rule for our polynomial. This will take the following form, for the polynomial
that we used previously:

muladd(x,muladd(x,muladd(x,muladd(x,5,4),3),2),1)

This is likely to be the fastest way to compute our polynomial. However, writing this out for
every polynomial that we might want to use will be extremely annoying. We loose the benefit of
having a general library function that can compute any polynomial.

This is exactly what macros are for!

macro horner(x,p...)
    ex = esc(p[end])
    for i = length(p) -1 : -1 : 1
        ex = :(muladd(t, $ex, $(esc(p[1]))))
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

f_horner_macro(x) = @horner(x,1,2,3,4,5)

@benchmark f_horner_macro(3.5)

REMEMBER Strategy of generating customized code for particualr problems.

GENERATED FUNCTIONS also called staged functions
===================

Macros for after type inference.

Generated functions come in two parts.
They are invoked once for each unique type of its arguments.
    Here, the arguments themselves take the values of their types.
    The return value of this execution MUST be an expression that is used as the body of the function when called with valeus of these types.

This cycle is executed each time the function is called with new types.
    The function is called with types as values once, and then the returned expression is used for all invocations with argument values of htis type.

Using generated functions for performance
    ==============================

REMEMBER use GenFuncs for performance by moving constant parts of the computation earlier into the compilation stage

Toy problem: number of cells in a multidim array.

function prod_dim{T, N}(x :: Array{T,N})
    s = 1
    for i = 1:N
        s = s*size(x,i)
    end
    return s
end

prod_dim(rand(10,5,5))

Optimizing this computation with a generated function starts with the observation that the number
of iterations of the loop is equal to the number of dimensions of the array, which is encoded as a
type parameter for arrays. In other words, for a particular type of input (and array of a particular
dimension), the loop size is fixed. So, what we can try to do in a generated function is move the
loop to the compile time:

@generated function prod_dim_gen_impl{T, N}(x :: Array{T,N})
    ex = :(1)
    for i = 1:N
        ex  = :(size(x, $i) * $ex)
    end
    return ex
end

prod_dim_gen_impl(rand(10,5,5))

Double checking
function prod_dim_gen_impl{T, N}(x::Array{T, N})
    ex = :(1)
        for i = 1:N
            ex = :(size(x, $i) * $ex)
        end
    return ex
end

x = rand(2,2,2)

prod_dim_gen_impl(x)

x = rand(2,2,2,2);

prod_dim_gen_impl(x)

REMEMBER this is a form of loop unrolling

Also, note that this function looks much simpler without the loop. The number of tokens in this
function is significantly reduced. This might make the function inline-worthy and cause the
compiler to inline this function, making this code even faster.

USING NAMED PARAMETERS
======================

REMEMBER Useful for creating APIs
However, compiler cannot infer the types of named parameters effectively

named_param(x; y = 1, z = 1) = x^y + x^z

pos_param(x,y,z) = x^y + x^z

@time for i in 1:100000; named_param(4; y = 2, z = 3); end

@time for i in 1:100000; pos_param(4,2,3); end

CHAPTER 5 - FAST NUMBERS
========================

Default numeric types are as close to metal as possible - so that we get C like speed.
Default sizes depend on the OS.

WORDS_SIZE

REMEMBER bits(3)

bits(3)

isbits(Int64)

REMEMBER isbits

isbits(String)

As a Julia value, basic numeric types can be boxed.

bits(-3)

myadd(x, y) = x + y
Looking at the output of of the following compiled code, (even if, like me, you are not an expert at
reading assembly), it should be apparent that, other than the function overhead to set the stack and
return the result, the generated code simply consists of the CPU instruction to add two machine
integers, addq. There is no boxing/unboxing operation remaining in the native code when the
function is called. Take a look at the following:
julia> @code_native myadd(1,2)

.section __TEXT,__text,regular,pure_instructions
Filename: none
Source line: 1
pushq %rbp
movq %rsp, %rbp
Source line: 1
addq %rsi, %rdi
movq %rdi, %rax
popq %rbp
ret

There is an even bigger advantage to storing numbers using the machine representation. Arrays of
these numbers can be stored using contiguous storage. A type tag is stored once at the start.
Beyond this, data in numerical arrays is stored in a packed form. This not only means that these
arrays can be passed to C libraries as-is (minus the type tag) but also that the compiler can
optimize computations on these arrays easily. There is no need for pointer dereferencing when
operating on numerical arrays of bit types.

Integer Overflow
================

There are no overflow checks present within any base mathematical op in Julia

typemax(Int64)

bits(typemax(Int64))

typemin(Int64)

bits(typemin(Int64))

typemax(Int64) + 1

Pg. 307

bits(TPixel)

40*18

.001 *45000

46440/17.77

.001*3300

1.12^30

bits(.3)

bits(.5)

bits(.7)

bits(.7*.5)

bits(.7*5)

bits(.7*.5)

bits(round(.7*.5,3))

round(.7*.5,3)

fibbi(n) = n <= 1 ? 1 : n*fibbi(n-1)

for i in 1:10
    println("fibbi($i) = ", fibbi(i))
end

max_stars_repo_path: Dojo/C++ Programming Practices and Principles - Stroustrup/HFTAldridge/juliahpc.jl

max_stars_repo_name: miguelraz/PathToPerformance

max_issues_repo_name: miguelraz/PathToPerformance

max_forks_repo_name: miguelraz/PathToPerformance


text:
using NQCBase: au_to_ang, eV_to_au, eV_per_ang_to_au
using .JuLIP: JuLIP

export JuLIPModel

"""
    struct JuLIPModel{T} <: AdiabaticModel

Model for interfacing with JuLIP potentials.
"""
struct JuLIPModel{T} <: AdiabaticModel
    atoms::JuLIP.Atoms{T}
end

NQCModels.ndofs(::JuLIPModel) = 3

function NQCModels.potential(model::JuLIPModel, R::AbstractMatrix)
    JuLIP.set_positions!(model.atoms, au_to_ang.(R))
    return eV_to_au(JuLIP.energy(model.atoms))
end

function NQCModels.derivative!(model::JuLIPModel, D::AbstractMatrix, R::AbstractMatrix)
    JuLIP.set_positions!(model.atoms, au_to_ang.(R))
    force = JuLIP.forces(model.atoms)
    D .= -eV_per_ang_to_au.(JuLIP.mat(force))
    return D
end


text:
function DubinsCar(scenario=:three_obstacles; N=101)
    if scenario == :three_obstacles
        opts = SolverOptions(
            cost_tolerance_intermediate=1e-2,
            penalty_scaling=10.,
            penalty_initial=10.0,
        )

        #  Car w/ obstacles
        model = RobotZoo.DubinsCar()
        n,m = size(model)

        N = 101 # number of knot points
        tf = 5.0
        dt = tf/(N-1)

        x0 = @SVector [0., 0., 0.]
        xf = @SVector [3., 3., 0.]

        Q = Diagonal(@SVector [1., 1., 1.])
        R = Diagonal(@SVector [0.5, 0.5])
        Qf = 10.0*Diagonal(@SVector ones(n))
        obj = LQRObjective(Q,R,Qf,xf,N)

        # create obstacle constraints
        r_circle_3obs = 0.25

        circle_x = 3*@SVector [0.25, 0.5, 0.75]
        circle_y = 3*@SVector [0.25, 0.5, 0.75]
        circle_r = @SVector fill(r_circle_3obs+model.radius, 3)

        obs = CircleConstraint(n, circle_x, circle_y, circle_r)
        bnd = BoundConstraint(n,m, u_min=[0,-3],u_max=[3,3])
        goal = GoalConstraint(xf)

        conSet = ConstraintList(n,m,N)
        add_constraint!(conSet, obs, 2:N-1)
        add_constraint!(conSet, bnd, 1:N-1)
        add_constraint!(conSet, goal, N:N)

        # Create problem
        U = [@SVector fill(0.01,m) for k = 1:N-1]
        car_3obs_static = Problem(model, obj, xf, tf, constraints=conSet, x0=x0)
        initial_controls!(car_3obs_static, U)
        rollout!(car_3obs_static)
        return car_3obs_static, opts

    elseif scenario==:turn90
        opts = SolverOptions(
            cost_tolerance_intermediate=1e-3,
            active_set_tolerance=1e-4
        )

        # model
        model = RobotZoo.DubinsCar()
        n,m = size(model)
        tf = 3.

        # cost
        d = 1.5
        x0 = @SVector [0., 0., 0.]
        xf = @SVector [d, d,  deg2rad(90)]
        Qf = 100.0*Diagonal(@SVector ones(n))
        Q = (1e-2)*Diagonal(@SVector ones(n))
        R = (1e-2)*Diagonal(@SVector ones(m))

        # problem
        U = [@SVector fill(0.1,m) for k = 1:N-1]
        obj = LQRObjective(Q,R,Qf,xf,N)

        # constraints
        cons = ConstraintList(n,m,N)
        add_constraint!(cons, GoalConstraint(xf), N)

        prob = Problem(model, obj, xf, tf, x0=x0, U0=U)
        rollout!(prob)

        return prob, opts

    elseif scenario==:parallel_park
        opts = SolverOptions(
            cost_tolerance_intermediate=1e-3,
            active_set_tolerance_pn=1e-4
        )

        # model
        model = RobotZoo.DubinsCar()
        n,m = size(model)
        tf = 3.

        # cost
        d = 1.5
        x0 = @SVector [0., 0., 0.]
        xf = @SVector [0., d,  0.]
        Qf = 100.0*Diagonal(@SVector ones(n))
        Q = (1e-2)*Diagonal(@SVector ones(n))
        R = (1e-2)*Diagonal(@SVector ones(m))

        # constraints
        u_bnd = 2.
        x_min = @SVector [-0.25, -0.001, -Inf]
        x_max = @SVector [0.25, d + 0.001, Inf]
        bnd_x = BoundConstraint(n,m, x_min=x_min, x_max=x_max)
        bnd_u = BoundConstraint(n,m, u_min=-u_bnd, u_max=u_bnd)
        # bnd = BoundConstraint(n,m,x_min=x_min,x_max=x_max,u_min=-u_bnd,u_max=u_bnd)
        goal = GoalConstraint(xf)

        conSet = ConstraintList(n,m,N)
        # add_constraint!(conSet, bnd, 1:N-1)
        add_constraint!(conSet, bnd_u, 1:N-1)
        add_constraint!(conSet, bnd_x, 2:N-1)
        add_constraint!(conSet, goal, N:N)

        # problem
        U = [@SVector fill(0.1,m) for k = 1:N-1]
        obj = LQRObjective(Q,R,Qf,xf,N)

        prob = Problem(model, obj, xf, tf, constraints=conSet, x0=x0, U0=U)
        rollout!(prob)

        return prob, opts

    elseif scenario==:escape
        opts = SolverOptions(
            constraint_tolerance=1e-5,
            cost_tolerance_intermediate=1e-3,
            penalty_scaling=100.,
            penalty_initial=10.,
            active_set_tolerance_pn=1e-4,
            active_set_tolerance_al=1e-4
        )

        # Static Car Escape
        T = Float64;

        # model
        model = RobotZoo.DubinsCar()
        n,m = size(model)
        x0 = @SVector [2.5,2.5,0.]
        xf = @SVector [7.5,2.5,0.]
        N = 101
        tf = 3.0

        # cost
        Q = (1e-3)*Diagonal(@SVector ones(n))
        R = (1e-2)*Diagonal(@SVector ones(m))
        Qf = 100.0*Diagonal(@SVector ones(n))
        obj = LQRObjective(Q,R,Qf,xf,N)

        # constraints
        r = 0.5
        s1 = 30; s2 = 50; s3 = 15

        circles_escape = NTuple{3,Float64}[]

        for i in range(0,stop=5,length=s1)
            push!(circles_escape,(0.,i,r))
        end
        for i in range(0,stop=5,length=s1)
            push!(circles_escape,(5.,i,r))
        end
        for i in range(0,stop=5,length=s1)
            push!(circles_escape,(10.,i,r))
        end
        for i in range(0,stop=10,length=s2)
            push!(circles_escape,(i,0.,r))
        end
        for i in range(0,stop=3,length=s3)
            push!(circles_escape,(i,5.,r))
        end
        for i in range(5,stop=8,length=s3)
            push!(circles_escape,(i,5.,r))
        end

        n_circles_escape = 3*s1 + s2 + 2*s3

        circles_escape
        x,y,r = collect(zip(circles_escape...))
        x = SVector{n_circles_escape}(x)
        y = SVector{n_circles_escape}(y)
        r = SVector{n_circles_escape}(r)

        obs = CircleConstraint(n,x,y,r)
        bnd = BoundConstraint(n,m,u_min=-5.,u_max=5.)
        goal = GoalConstraint(xf)

        conSet = ConstraintList(n,m,N)
        add_constraint!(conSet, obs, 2:N-1)
        add_constraint!(conSet, bnd, 1:N-1)
        add_constraint!(conSet, goal, N:N)

        # Build problem
        U0 = [@SVector ones(m) for k = 1:N-1]

        car_escape_static = Problem(model, obj, xf, tf;
            constraints=conSet, x0=x0)
        initial_controls!(car_escape_static, U0);

        X_guess = [2.5 2.5 0.;
                   4. 5. .785;
                   5. 6.25 0.;
                   7.5 6.25 -.261;
                   9 5. -1.57;
                   7.5 2.5 0.]
        X0_escape = Altro.interp_rows(N,tf,Array(X_guess'))
        initial_states!(car_escape_static, X0_escape)

        return car_escape_static, opts
    end
end


text:
const allpermorder = [:NATURAL, :ENTROPIC, :REV_ENTROPIC, :RANDOM]
"""
    ardca(Z::Array{Ti,2},W::Vector{Float64}; kwds...)
Auto-regressive analysis on the L×M alignment `Z` (numerically encoded in 1,…,21), and the `M`-dimensional normalized
weight vector `W`.

Return two `struct`: `::ArNet` (containing the inferred hyperparameters) and `::ArVar`

Optional arguments:
* `lambdaJ::Real=0.01` coupling L₂ regularization parameter (lagrange multiplier)
* `lambdaH::Real=0.01` field L₂ regularization parameter (lagrange multiplier)
* `epsconv::Real=1.0e-5` convergence value in minimzation
* `maxit::Int=1000` maximum number of iteration in minimization
* `verbose::Bool=true` set to `false` to stop printing convergence info on `stdout`
* `method::Symbol=:LD_LBFGS` optimization strategy see [`NLopt.jl`](https://github.com/JuliaOpt/NLopt.jl) for other options
* `permorder::Union{Symbol,Vector{Ti}}=:ENTROPIC` permutation order. Possible values are `:NATURAL,:ENTROPIC,:REV_ENTROPIC,:RANDOM` or a custom permutation vector

# Examples
```
julia> arnet, arvar= ardca(Z,W,lambdaJ=0,lambdaH=0,permorder=:REV_ENTROPIC,epsconv=1e-12);
```
"""
function ardca(Z::Array{Ti,2},W::Vector{Float64};
                lambdaJ::Real=0.01,
                lambdaH::Real=0.01,
                epsconv::Real=1.0e-5,
                maxit::Int=1000,
                verbose::Bool=true,
                method::Symbol=:LD_LBFGS,
                permorder::Union{Symbol,Vector{Int}}=:ENTROPIC
                ) where Ti <: Integer

    checkpermorder(permorder)
    all(x -> x > 0, W) || throw(DomainError("vector W should normalized and with all positive elements"))
    isapprox(sum(W), 1) || throw(DomainError("sum(W) ≠ 1. Consider normalizing the vector W"))
    N, M = size(Z)
    M = length(W)
    q = Int(maximum(Z))
    aralg = ArAlg(method, verbose, epsconv, maxit)
    arvar = ArVar(N, M, q, lambdaJ, lambdaH, Z, W, permorder)
    θ,psval = minimize_arnet(aralg, arvar)
    Base.GC.gc() # something wrong with SharedArrays on Mac
    ArNet(θ,arvar),arvar
end
"""
    ardca(filename::String; kwds...)
Run [`ardca`](@ref) on the fasta alignment in `filename`

Return two `struct`: `::ArNet` (containing the inferred hyperparameters) and `::ArVar`

Optional arguments:
* `max_gap_fraction::Real=0.9` maximum fraction of insert in the sequence
* `remove_dups::Bool=true` if `true` remove duplicated sequences
* `theta=:auto` if `:auto` compute reweighint automatically. Otherwise set a `Float64` value `0 ≤ theta ≤ 1`
* `lambdaJ::Real=0.01` coupling L₂ regularization parameter (lagrange multiplier)
* `lambdaH::Real=0.01` field L₂ regularization parameter (lagrange multiplier)
* `epsconv::Real=1.0e-5` convergence value in minimzation
* `maxit::Int=1000` maximum number of iteration in minimization
* `verbose::Bool=true` set to `false` to stop printing convergence info on `stdout`
* `method::Symbol=:LD_LBFGS` optimization strategy see [`NLopt.jl`](https://github.com/JuliaOpt/NLopt.jl) for other options
* `permorder::Union{Symbol,Vector{Ti}}=:ENTROPIC` permutation order. Possible
  values are `:NATURAL,:ENTROPIC,:REV_ENTROPIC,:RANDOM` or a custom permutation
  vector

# Examples
```
julia> arnet, arvar =  ardca("pf14.fasta", permorder=:ENTROPIC)
```
"""
function ardca(filename::String;
                theta::Union{Symbol,Real}=:auto,
                max_gap_fraction::Real=0.9,
                remove_dups::Bool=true,
                kwds...)
    W, Z, N, M, q = read_fasta(filename, max_gap_fraction, theta, remove_dups)
    W ./= sum(W)
    ardca(Z, W; kwds...)
end

function checkpermorder(po::Symbol)
    po ∈ allpermorder || error("permorder :$po not iplemented: only $allpermorder are defined");
end

(checkpermorder(po::Vector{Ti}) where Ti <: Integer) = isperm(po) || error("permorder is not a permutation")

function minimize_arnet(alg::ArAlg, var::ArVar{Ti}) where Ti
    @extract var : N q q2
    @extract alg : epsconv maxit method
    vecps = Vector{Float64}(undef,N - 1)
    θ = Vector{Float64}(undef, ((N*(N-1))>>1)*q2 + (N-1)*q)
    Threads.@threads for site in 1:N-1
        x0 = zeros(Float64, site * q2 + q)
        opt = Opt(method, length(x0))
        ftol_abs!(opt, epsconv)
        xtol_rel!(opt, epsconv)
        xtol_abs!(opt, epsconv)
        ftol_rel!(opt, epsconv)
        maxeval!( opt, maxit)
        min_objective!(opt, (x, g) -> optimfunwrapper(x, g, site, var))
        elapstime = @elapsed  (minf, minx, ret) = optimize(opt, x0)
        alg.verbose && @printf("site = %d\tpl = %.4f\ttime = %.4f\t", site, minf, elapstime)
        alg.verbose && println("status = $ret")
        vecps[site] = minf
        offset = div(site*(site-1),2)*q2 + (site-1)*q + 1
        θ[offset:offset+site * q2 + q - 1] .= minx
    end
    return θ, vecps
end

function optimfunwrapper(x::Vector, g::Vector, site, var)
    g === nothing && (g = zeros(Float64, length(x)))
    return pslikeandgrad!(x, g, site,  var)
end

function pslikeandgrad!(x::Vector{Float64}, grad::Vector{Float64}, site::Int, arvar::ArVar)
    @extract arvar : N M q q2 lambdaJ lambdaH Z W IdxZ
    LL = length(x)
    for i = 1:LL - q
        grad[i] = 2.0 * lambdaJ  * x[i]
    end
    for i = (LL - q + 1):LL
        grad[i] = 2.0 * lambdaH * x[i]
    end
    pseudolike = 0.0
    vecene = zeros(Float64, q)
    expvecenesumnorm = zeros(Float64, q)
    @inbounds for m in 1:M
        izm = view(IdxZ, :, m)
        zsm = Z[site+1,m] # the i index of P(x_i|x_1,...,x_i-1) corresponds here to i+1
        fillvecene!(vecene, x, site, izm, q, N)
        lnorm = logsumexp(vecene)
        expvecenesumnorm .= @. exp(vecene - lnorm)
        pseudolike -= W[m] * (vecene[ zsm ] - lnorm)
        sq2 = site * q2
        @avx for i in 1:site
            for s in 1:q
                grad[ izm[i] + s ] += W[m] * expvecenesumnorm[s]
            end
            grad[ izm[i] + zsm ] -= W[m]
        end
        @avx for s = 1:q
            grad[ sq2 + s ] += W[m] * expvecenesumnorm[s]
        end
        grad[ sq2 + zsm ] -= W[m]
    end
    pseudolike += l2norm_asym(x, arvar)
end

function fillvecene!(vecene::Vector{Float64}, x::Vector{Float64}, site::Int, IdxSeq::AbstractArray{Int,1}, q::Int, N::Int)
    q2 = q^2
    sq2 = site * q2
    @inbounds for l in 1:q
        scra = 0.0
        @avx for i in 1:site
            scra += x[IdxSeq[i] + l]
        end
        scra += x[sq2 + l] # sum H
        vecene[l] = scra
    end
end

function logsumexp(X::Vector)
    u = maximum(X)
    isfinite(u) || return float(u)
    return u + log(sum(x -> exp(x - u), X))
end

function l2norm_asym(vec::Array{Float64,1}, arvar::ArVar)
    @extract arvar : q N lambdaJ lambdaH
    LL = length(vec)
    mysum1 = 0.0
    @inbounds @avx for i = 1:(LL - q)
        mysum1 += vec[i] * vec[i]
    end
    mysum1 *= lambdaJ
    mysum2 = 0.0
    @inbounds @avx for i = (LL - q + 1):LL
        mysum2 += vec[i] * vec[i]
    end
    mysum2 *= lambdaH
    return mysum1 + mysum2
end

max_issues_repo_name: pagnani/ArDCA

max_forks_repo_name: pagnani/ArDCA


text:
using Plots
using DifferentialEquations

function seirmodel(du, u, p, t)
    S = u[1]
    E = u[2]
    I = u[3]
    R = u[4]
    dS = -p[1]*S*I
    dE = p[1]*S*I-p[2]*E
    dI = p[2]*E-p[3]*I
    dR = p[3]*I
    du[1] = dS
    du[2] = dE
    du[3] = dI
    du[4] = dR
end

function solveSEIR(u0, β, γ, δ, time_span)

    prob = ODEProblem(seirmodel, u0, time_span, [β, δ, γ])

    sol = solve(prob)

    return sol

end

R0=2.0
N=10000
γ=0.1
β=γ*R0/N
δ=0.25

sol = solveSEIR(Array{Float64}([N-1,0,1,0]), β, γ, δ, (0.0,360.0))

plot_graph = plot(sol, vars=(0,1), title="SEIR model/Runge-Kutta R0=$R0 γ=$γ δ=$δ", label="Suscepible",
                  color=:blue, seriestype=:line, legend=:left, size=(1200,900))
plot!(plot_graph, sol, vars=(0,2), label="Exposed", color=:darkmagenta)
plot!(plot_graph, sol, vars=(0,3), label="Infectious", color=:red)
plot!(plot_graph, sol, vars=(0,4), label="Removed", color=:green)

png(plot_graph, "seir_ode.png")

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
using QuNet
using DataFrames
using CSV
using Statistics

"""
Convert heatmap data to a Pandas Dataframe
"""
function coord_list_to_dataframe(coord_list)
    # List comprehension to extract data
    e = [coord_list[i][1] for i in 1:length(coord_list)]
    f = [coord_list[i][2] for i in 1:length(coord_list)]

    # # Convert cood_list to data_frame and write to csv
    df = DataFrame(Efficiency = e, Fidelity = f)
    return df
end

# Graph + sampling parameters
grid_size = 100
num_trials = 500

edge_costs = Dict("loss"=>0.05, "Z"=>0.05)
net = GridNetwork(grid_size, grid_size, edge_costs=edge_costs)

#usage:
# (network::Union{QNetwork, QuNet.TemporalGraph},
#    num_trials::Int64, num_pairs::Int64; max_paths=3, src_layer::Int64=-1,
#    dst_layer::Int64=-1, edge_perc_rate=0.0)

pair_sizes = [1, 50, 500, 2500]
max_paths = [1, 2, 3, 4]
for num_pairs in pair_sizes
    println("Collecting for $num_pairs userpairs")
    for num_paths in max_paths
        println("\tCollecting for $num_paths paths")
        coord_list = QuNet.heat_data(net, num_trials, num_pairs, max_paths=num_paths)
        df = coord_list_to_dataframe(coord_list)
        CSV.write("data/heatmap_data/$(num_pairs)pair$(num_paths)path.csv", df)
    end
end

# NOTE: Also do one data collection for much larger sample size:
grid_size = 100
num_trials = 5000
num_pairs = 50
max_paths = 4
println("Collecting for 50 userpair")
coord_list = QuNet.heat_data(net, num_trials, num_pairs, max_paths=max_paths)
df = coord_list_to_dataframe(coord_list)
CSV.write("data/HiRes_50pair4path.csv", df)

"""
This function gets the end-to-end failure rate (P₀) of the graph in the heatmap
plot. By my calculations, this is about 0.20 for the specified parameters
"""
function get_P0_for_graph()
    num_trials = 100
    num_pairs = 50
    max_paths = 4

    failure_counts = []
    for i in 1:num_trials
        println("num_trial $i")
        edge_costs = Dict("loss"=>0.05, "Z"=>0.05)
        net = GridNetwork(grid_size, grid_size, edge_costs=edge_costs)
        user_pairs = QuNet.make_user_pairs(net, num_pairs)
        pathset, pur_paths, pathuse_count = QuNet.greedy_multi_path!(net, QuNet.purify, user_pairs, 4)
        push!(failure_counts, pathuse_count[1])
    end
    return(mean(failure_counts./num_pairs))
end

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
function movements_to_side(r, side :: HorizonSide)
    count = 0
    while isborder(r,side) == false
        move!(r,side)
        count = count + 1
    end
    return count
end

function ugol(r)

    count_s_n = 0
    count_w_o = 0
    home = Int[]

    while (isborder(r, Sud) == false) || (isborder(r, West) == false)
        s = movements_to_side(r, Sud)
        w = movements_to_side(r, West)
        count_s_n = count_s_n + s
        count_w_o = count_w_o + w
        push!(home, s)
        push!(home, w)
    end
    return home, count_s_n, count_w_o
end

function count_move!(r, side :: HorizonSide, count)

    for i = 1:count
        move!(r, side)
    end
end

function go_home(r, home)
    reverse!(home)

    for i in 1: length(home)
        if i % 2 == 0
            count_move!(r, Nord, home[i])
        else
            count_move!(r, Ost, home[i])
        end
    end
end

function Crest(r)

    home = Int[]
    home, count_s_n , count_w_o = ugol(r)

    for side in (Nord, Ost, Sud, West)
        if (side == Nord) || (side == Sud)
            count_move!(r, side, count_s_n)
            putmarker!(r)
            count_s_n = movements_to_side(r,side)
        else
            count_move!(r, side, count_w_o)
            putmarker!(r)
            count_w_o = movements_to_side(r,side)
        end
    end
    go_home(r,home)

end

Crest(r)
show(r)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
function get_possible_decisions(state::DynamicIndividualState,
     static_state::StaticIndividualState;
     allow_nothing::Bool = false, allow_overfarming::Bool = false)
     # 00000001 - shield
     # 00000010 - nothing
     # 00000100 - fast move
     # 00001000 - charged move
     # 00010000 - charged move 1
     # 00100000 - charged move 2

     fast_moves_pending = get_fast_moves_pending(state)

     if get_hp(state[0x01]) == 0x0000 || get_hp(state[0x02]) == 0x0000
          return 0x00, 0x00
     end

     cmp = get_cmp(state)
     if isodd(cmp) # if team 1 is using a charged move and has cmp
          return (get_energy(state[0x01]) >=
               get_energy(static_state[0x01].charged_move_2) ?
               0x30 : 0x10, has_shield(state, 0x02) ? 0x03 : 0x02)
     elseif !iszero(cmp) # if team 2 is using a charged move and has cmp
          return (has_shield(state, 0x01) ? 0x03 : 0x02,
               get_energy(state[0x02]) >=
               get_energy(static_state[0x02].charged_move_2) ?
               0x30 : 0x10)
     end

     d = 0x00, 0x00
     for i = 0x01:0x02
          if fast_moves_pending[i] <= 0x0001
               if allow_overfarming || get_energy(state[i]) != 0x0064
                    d = i == 0x01 ? (d[1] + 0x04, d[2]) : (d[1], d[2] + 0x04)
               end
               if get_energy(state[i]) >= get_energy(static_state[i].charged_move_1)
                    d = i == 0x01 ? (d[1] + 0x08, d[2]) : (d[1], d[2] + 0x08)
               end
               if allow_nothing
                    d = i == 0x01 ? (d[1] + 0x02, d[2]) : (d[1], d[2] + 0x02)
               end
          else
               d = i == 0x01 ? (d[1] + 0x02, d[2]) : (d[1], d[2] + 0x02)
          end
     end

     return d
end


text:
# This file was generated, do not modify it. # hide
X = permutedims(reshape([1, 2, 3, 4], 2, 2))


text:
using Plots, Printf, DelimitedFiles, Random, Statistics, ProgressMeter
using MDToolbox

using BSON: @save, @load
@load "afm_generate_test_data.bson" afms models models_template qs myseed sigma_noise nframe

configs = [AfmizeConfig(10.0 * (pi / 180),
    r,
    MDToolbox.Point2D(-250, -200),
    MDToolbox.Point2D(250, 200),
    MDToolbox.Point2D(6.25, 6.25),
    MDToolbox.defaultParameters())
    for r in [10.0, 15.0, 20.0, 25.0, 30.0]]

# single frame
r = []
@time @showprogress 1 "Computing..." for i = 1:20
    r2 = getafmposterior(afms[i], models_template, qs, configs)
    push!(r, r2)
end

# all frames with parallel map
#r = pmap(x -> getafmposterior(x, models_template, qs, configs), afms[1:20])

@save "afm_analyze_test_data_sequential.bson" r configs

max_stars_repo_name: matsunagalab/notebooks

max_issues_repo_name: matsunagalab/notebooks

max_forks_repo_name: matsunagalab/notebooks


text:
using Test;

# Learning:
# - my first "multiple dispatch"... I added a fallback here so that
#   minimum will still work even though I'm giving one Float (±Inf) and one Int
# - Also, apparently I can type `\pm` and VSCode will convert it to ± for Julia magic. Cool!
function minimum(a, b)
    return a < b ? a : b
end

function maximum(a, b)
    return a > b ? a : b
end

function find_pair(prev_lines, val)
    for i in 1:length(prev_lines)
        for j in (i + 1):length(prev_lines)
            if prev_lines[i] + prev_lines[j] == val
                return true
            end
        end
    end
    return false
end

function parse_line(s)
    return parse(Int64, s)
end

function run(fname, preamble_length)
    # read file
    f = open(fname, "r")
    lines = map(parse_line, readlines(f))
    close(f)

    result = -1
    for (idx, val) in enumerate(lines)
        if idx > preamble_length
            if ! find_pair(lines[idx - preamble_length:idx - 1], val)
                result = val
                break
            end
        end
    end

    println("Ran ", fname, " got ", result)
    return result
end

function run2(fname, target)
    # read file
    f = open(fname, "r")
    lines = map(parse_line, readlines(f))
    close(f)

    result = -1
    len = length(lines)
    for idx in 1:len
        total = 0
        low = +Inf
        high = -Inf
        inner_idx = idx
        while total < target && inner_idx <= len
            val = lines[inner_idx]
            low = minimum(low, val)
            high = maximum(high, val)
            total += val
            inner_idx += 1
        end

        if total == target
            println("Total = ", total, " Low = ", low, " High = ", high)
            result = low + high
            break
        end
    end

    println("[part 2] Ran ", fname, " got ", result)
    return result
end

@test run("9ex.txt", 5) == 127
run("9.txt", 25)

@test run2("9ex.txt", 127) == 62
run2("9.txt", 3199139634)


text:
using Test
using Dates

using ScottishTaxBenefitModel

using .ModelHousehold:
    Household,
    Person,
    People_Dict,
    default_bu_allocation,
    get_benefit_units,
    get_head,
    get_spouse,
    is_single,
    pers_is_carer,
    pers_is_disabled,
    search

using .Definitions

using .BenefitCap:
    apply_benefit_cap!

using .STBIncomes

using .Intermediate:
    MTIntermediate,
    make_intermediate

using .STBParameters:
    UniversalCreditSys,
    BenefitCapSys,
    LegacyMeansTestedBenefitSystem

using .Results:
    BenefitUnitResult,
    HouseholdResult,
    init_household_result,
    to_string

using .Utils:
    eq_nearest_p,
    to_md_table

using .ExampleHelpers

sys = get_system( scotland=true )

@testset "Benefit Cap Shakedown" begin
    examples = get_all_examples()
    hbs = collect(100:100:1000)
    for (hht,hh) in examples
        bus = get_benefit_units( hh )
        intermed = make_intermediate(
            hh,
            sys.hours_limits,
            sys.age_limits,
            sys.child_limits )
        head = get_head( hh )
        spouse = get_spouse( hh )
        spid :: BigInt = spouse === nothing ? head.pid : spouse.pid
        for hb in 200.0:200:1_000
            for cb in [20,50,100]
                res = init_household_result( hh )
                res.bus[1].pers[head.pid].income[HOUSING_BENEFIT] = hb
                res.bus[1].pers[spid].income[HOUSING_BENEFIT] = cb
                for buno in eachindex(bus)
                    apply_benefit_cap!(
                        res.bus[buno],
                        hh.region,
                        bus[buno],
                        intermed.buint[buno],
                        sys.bencap,
                        legacy_bens
                    )
                    println( "on family $hht bu $buno hb=$hb cb = $cb")
                    println( res.bus[buno].bencap )
                    if buno == 1
                        println( "res.bus[$buno].pers[\$(head.pid)].income[HOUSING_BENEFIT] = $(res.bus[buno].pers[head.pid].income[HOUSING_BENEFIT])\n\n" )
                    end
                end # bus
            end # cbs
        end # hbs
    end # hhs
end # testset

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
# MATLAB related functionality
using MATLAB

##############################################################################
##
## Matlab setup
##
##############################################################################

# Give access to BSR Benchmark code
@matlab addpath("/ais/gobi3/u/jsnell/projects/eccv2014/BSR/bench/benchmarks")

# Give access to BSR lib code for gPb and UCM
@matlab addpath("/u/jsnell/gobi3/projects/eccv2014/BSR/grouping/lib")

function mstatementquiet(ex::Expr)
    ss = IOBuffer()
    MATLAB.write_mstatement(ss, ex)
    print(ss, ";")
    bytestring(ss)
end

macro mquiet(ex)
    :( MATLAB.eval_string($(mstatementquiet(ex))) )
end

##############################################################################
##
## Algorithms for computing superpixels
##
##############################################################################

function slic{T<:Real}(img::Array{T}, regionsize::Int=20, regularizer::Float64=0.1)
    @mput img
    eval_string("img = im2single(img);")
    eval_string("seg = vl_slic(img, $regionsize, $regularizer);")
    convert(Matrix{Int}, @mget seg) + 1     # minimum segment should be 1
end

##############################################################################
##
## Conversion functions
##
##############################################################################

function thinboundary(bmap::BitMatrix)
    @mput bmap
    @mquiet bmapthin = bwmorph(bmap, "thin", inf)
    convert(BitArray, @mget bmapthin)
end

function correspondpixels(bmap1::BitMatrix, bmap2::BitMatrix, maxdist::Float64=0.0075, outliercost::Float64=100.0)
    @mput bmap1  bmap2 maxdist outliercost
    @mquiet (match1, match2, cost, oc) = correspondPixels(double(bmap1), double(bmap2), maxdist, outliercost)
    @mget match1 match2 cost oc
    ([int(match1)], [int(match2)], cost, oc)
end

##############################################################################
##
## Compute gPb
##
##############################################################################

function gPb{T<:Real}(img::Array{T})
    @mput img
    @mquiet gPb_orient = globalPb_im(img)
    @mget gPb_orient
end

##############################################################################
##
## UCM
##
##############################################################################

function ucm_regions{T<:Real}(img::Array{T}, scale::T=0.0)
    gpb_res = gPb(img)
    @mput gpb_res
    eval_string("ucm2 = contours2ucm(gpb_res, 'doubleSize');")
    eval_string("labels2 = bwlabel(ucm2 <= $scale);")
    eval_string("labels = labels2(2:2:end, 2:2:end);")
    convert(Matrix{Int}, @mget labels)
end

function ucm_regions{T<:Real}(img::Array{T}, scales::Vector{T})
    gpb_res = gPb(img)
    @mput gpb_res
    eval_string("ucm2 = contours2ucm(gpb_res, 'doubleSize');")
    alllabels = map(scales) do scale
        eval_string("labels2 = bwlabel(ucm2 <= $scale);")
        eval_string("labels = labels2(2:2:end, 2:2:end);")
        convert(Matrix{Int}, @mget labels)
    end
    Matrix{Int}[l for l in alllabels]
end

function ucm_strength{T<:Real}(img::Array{T}, scale::T=0.0)
    # returns the result of a 0.0 strength ucm,
    # along with the corresponding horizontal and vertical strengths
    gpb_res = gPb(img)
    @mput gpb_res
    eval_string("ucm2 = contours2ucm(gpb_res, 'doubleSize');")
    eval_string("labels2 = bwlabel(ucm2 <= $scale);")
    eval_string("labels = labels2(2:2:end, 2:2:end);")
    labels = convert(Matrix{Int}, @mget labels)

    ucm2 = @mget ucm2

    hstrength = ucm2[2:2:end, 3:2:end]
    vstrength = ucm2[3:2:end, 2:2:end]

    labels, hstrength, vstrength
end


text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

# Options
@enum(LogLevels::Int,
      TRACE  = 1,
      DEBUG  = 2,
      INFO   = 3,
      NOTICE = 4,
      WARN   = 5,
      ERROR  = 6)

@enum(FixedVariableTreatments::Int,
      RELAX_BOUND = 1,
      MAKE_PARAMETER = 2)

@enum(InertiaCorrectionMethod::Int,
      INERTIA_AUTO = 1,
      INERTIA_BASED = 2,
      INERTIA_FREE = 3)

@enum(KKTLinearSystem::Int,
      SPARSE_KKT_SYSTEM = 1,
      SPARSE_UNREDUCED_KKT_SYSTEM = 2,
      DENSE_KKT_SYSTEM = 3)

@enum(Status::Int,
      SOLVE_SUCCEEDED = 1,
      SOLVED_TO_ACCEPTABLE_LEVEL = 2,
      SEARCH_DIRECTION_BECOMES_TOO_SMALL = 3,
      DIVERGING_ITERATES = 4,
      INFEASIBLE_PROBLEM_DETECTED = 5,
      MAXIMUM_ITERATIONS_EXCEEDED = 6,
      MAXIMUM_WALLTIME_EXCEEDED = 7,
      INITIAL = 11,
      REGULAR = 12,
      RESTORE = 13,
      ROBUST  = 14,
      RESTORATION_FAILED = -1,
      INVALID_NUMBER_DETECTED = -2,
      ERROR_IN_STEP_COMPUTATION = -3,
      NOT_ENOUGH_DEGREES_OF_FREEDOM = -4,
      USER_REQUESTED_STOP = -5,
      INTERNAL_ERROR = -6)

const STATUS_OUTPUT_DICT = Dict(
    SOLVE_SUCCEEDED => "Optimal Solution Found.",
    SOLVED_TO_ACCEPTABLE_LEVEL => "Solved To Acceptable Level.",
    SEARCH_DIRECTION_BECOMES_TOO_SMALL => "Search Direction is becoming Too Small.",
    DIVERGING_ITERATES => "Iterates divering; problem might be unbounded.",
    MAXIMUM_ITERATIONS_EXCEEDED => "Maximum Number of Iterations Exceeded.",
    MAXIMUM_WALLTIME_EXCEEDED => "Maximum wall-clock Time Exceeded.",
    RESTORATION_FAILED => "Restoration Failed",
    INFEASIBLE_PROBLEM_DETECTED => "Converged to a point of local infeasibility. Problem may be infeasible.",
    INVALID_NUMBER_DETECTED => "Invalid number in NLP function or derivative detected.",
    ERROR_IN_STEP_COMPUTATION => "Error in step computation.",
    NOT_ENOUGH_DEGREES_OF_FREEDOM => "Problem has too few degrees of freedom.",
    USER_REQUESTED_STOP => "Stopping optimization at current point as requested by user.",
    INTERNAL_ERROR => "Internal Error.")


text:
"""
Every node is in its own module.
"""
function partition_single(A::Array{Int64, 2})
  # Every species is its own partition
  S = eye(Int64, sum(size(A)))
  M = Modular(A, S)
  return M
end


text: include("init.jl")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "P4est"
version = v"2.2.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://p4est.github.io/release/p4est-2.2.tar.gz", "1549cbeba29bee2c35e7cc50a90a04961da5f23b6eada9c8047f511b90a8e438"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd p4est-2.2/
if [[ "${target}" == *-freebsd* ]]; then
  export LIBS="-lm"
fi
CFLAGS="-I${prefix}/include" LDFLAGS="-L${prefix}/lib" ./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --without-blas
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64)
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libp4est", :libp4est),
    LibraryProduct("libsc", :libsc)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of row numbers.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Show row number" begin
    expected = """

      Row
      Col. 1
      Col. 2
      Col. 3
      Col. 4

      1
      1
      false
      1.0
      1

      2
      2
      true
      2.0
      2

      3
      3
      false
      3.0
      3

      4
      4
      true
      4.0
      4

      5
      5
      false
      5.0
      5

      6
      6
      true
      6.0
      6

"""
    result = pretty_table(
        String,
        data;
        alignment = [:l, :r, :c, :r],
        backend = Val(:html),
        standalone = false,
        show_row_number = true
    )

    @test result == expected

    expected = """

      Row number
      Col. 1
      Col. 2
      Col. 3
      Col. 4

      1
      1
      false
      1.0
      1

      2
      2
      true
      2.0
      2

      3
      3
      false
      3.0
      3

      4
      4
      true
      4.0
      4

      5
      5
      false
      5.0
      5

      6
      6
      true
      6.0
      6

"""
    result = pretty_table(
        String,
        data;
        alignment = [:l, :r, :c, :r],
        backend = Val(:html),
        standalone = false,
        row_number_column_title = "Row number",
        show_row_number = true
    )

    @test result == expected
end


text:
using Test, Random
using Julog, PDDL
using ParallelPlanners

include("heuristics.jl")
include("planners.jl")