text:
# ------------------------------------------------------------------------------------------
# # Tools - Using arrays to store data
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Introduction to arrays
#
# **Arrays are collections of boxes that we can use to store data.** In the last notebook,
# we saw an image that can help us to picture a 1-dimensional array:
#
#
#
# *Why do we want an object like this to store our data?*
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# An alternative to using an array in some contexts would be to name every individual piece
# of data, as follows:
#
# ```julia
# a = 1.1
# b = 2.2
# c = 3.3
# ```
#
# We can visualize how this data is stored:
#
#
#
#
#
# This is like having a separate box for every piece of data, rather than a series of
# connected boxes for all our data.
#
# The more data we have, the more annoying it becomes to keep track of all these boxes and
# their names. Furthermore, if we want to do the same thing with many pieces of data, it's
# much easier to put all of these pieces of data in one place to work with them at once.
#
# For example, we may want to multiply `a`, `b`, and `c` by `2`. We could multiply three
# times:
#
# ```julia
# a * 2
# b * 2
# c * 2
# ```
#
# Or, instead, we could create one array (let's call it `numbers`) and multiply that array
# by `2`:
# ```julia
# numbers * 2
# ```
#
# The syntax for creating this array, `numbers`, is
#
# ```julia
# numbers = [a, b, c]
# ```
#
# Or, we could have just written
#
# ```julia
# numbers = [1.1, 2.2, 3.3]
# ```
#
# It's worth noting that in Julia, 1-dimensional arrays are also called "vectors".
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Creating arrays
#
# In the last section, we saw that we could create the array `numbers` by typing our
# elements, `a`, `b`, and `c` (or `1.1`, `2.2`, and `3.3`), inside square brackets.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
#
#
# Create an array called `first_array` that stores the numbers 10 through 20.
#
# Reminder: ESC+b to open a box below.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Array comprehensions
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Alternatively we can create arrays via *"array comprehensions"*. This is a nice way to
# automate creating an array, if you don't want to have to type long lists of numbers inside
# square brackets.
#
# In an array comprehension, you write code inside square brackets that will generate the
# array as that code is executed. Let's see an example.
#
# ```julia
# counting = [2 * i for i in 1:10]
# ```
#
# The above line of code is saying that we want to create an array called `counting` that
# stores two times each of the integers between 1 and 10. This array comprehension has a few
# different parts, so let's dissect them:
#
#
#
# If we wanted to create another array, `roots`, that stored the square roots of all
# integers between 1 and 10, we could modify the code above to
#
# ```julia
# roots = [sqrt(i) for i in 1:10]
# ```
#
# #### Exercise
#
# Create an array, `squares`, that stores the squares of all integers between 1 and 100.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Looking inside an array
#
# We can "index" into an array to grab contents inside the array at a particular position.
# To index into our `counting` array, we place square brackets after the name of the array
# and put the position number of the element/data we want inside those brackets.
#
# For example, we might say
#
# ```julia
# counting[3]
# ```
#
# to grab the third element in the array called `counting`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 2
#
# Execute the code in the next cell. It will generate an array called `myprimes` that stores
# the first 168 prime numbers. Index into `myprimes` to grab the 89th smallest prime. What
# is this prime?
# ------------------------------------------------------------------------------------------

using Primes

myprimes = primes(1000); # The semicolon suppresses the output, try removing it

# ------------------------------------------------------------------------------------------
# ### Slicing
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Instead of grabbing a single number in an array, we can also take a **slice** of an array,
# i.e. a subset of the array, which can include multiple values. To take a slice of
# `counting` that includes the 3rd, 4th, and 5th entries, we use the following syntax with a
# colon (`:`):
#
# ```julia
# counting[3:5]
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 3
#
# Index into `myprimes` to grab the 89th through the 99th smallest primes (inclusive).
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Modifying an array
#
# In this section, we'll see how to edit an array we've already created. We'll see how to
#
# 1) Update an item at an existing position in an array;
#
# 2) Add an item to the end of an array;
#
# 3) Remove an item from the end of an array.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Update an item at an existing position in an array
#
# First off, we can update an item at an existing position in an array by indexing into the
# array and changing the value at the desired index. For example, let's say we have an array
# called `myfriends`:
#
# ```julia
# myfriends = ["Ted", "Robin", "Barney", "Lily", "Marshall"]
# ```
#
# We can grab Barney's name from `myfriends` via
#
# ```julia
# myfriends[3]
# ```
#
# and we can change "Barney" to something else by reassigning `myfriends[3]`:
#
# ```julia
# myfriends[3] = "Baby Bop"
# ```
#
# Note that a single `=` assigns a new variable to the value on the left of the `=` sign.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 4
#
# Use an array comprehension to create an array, `mysequence`, that stores the numbers 4
# through 10. Index into `mysequence` and update it to replace the last element, `10`, with
# `11`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Add an item to the end of an array
#
# We can add an item to the end of an array with the `push!` function. Don't worry about the
# exclamation mark at the end of `push!`; we'll talk about why that's there later, when we
# discuss functions.
#
# For now, note that when you call `push!` on an array and a value, you add that value to
# the end of the array. For example,
#
# ```julia
# push!(counting, 1000)
# ```
#
# will turn `counting`, which was formerly
#
# ```julia
# [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
# ```
# into
# ```julia
# [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1000]
# ```
#
# Prove to yourself this works!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 5
#
# Copy the following code to declare the array `fibonacci`:
#
# ```julia
# fibonacci = [1, 1, 2, 3, 5, 8, 13]
# ```
#
# Use `push!` to add `21` to `fibonacci` after the number `13`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Remove an item from the end of an array
#
# To remove an item from the end of an array, use the `pop!` function. When using `pop!` on
# arrays, `pop!` only needs one input argument, namely the array you want to change. For
# example, if
#
# ```julia
# counting == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1000]
# ```
#
# then
# ```julia
# pop(counting)
# ```
# will remove `1000` from the end of `counting`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 6
#
# Use `pop!` to remove `21` from `fibonacci`. What does this function call return?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 7
#
# What is the last element of `fibonacci`? Try
#
# ```julia
# fibonacci[end]
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Arrays of arrays and multidimensional arrays
#
# So far we've only seen 1D arrays of scalars, but arrays can also store other arrays and
# can have an arbitrary number of dimensions.
#

# For example, the following are arrays of arrays:
# ------------------------------------------------------------------------------------------

favorites = [ ["koobideh", "chocolate", "eggs"], ["penguins", "cats", "sugargliders"] ]

numbers = [ [1, 2, 3], [4, 5], [6, 7, 8, 9] ]

# ------------------------------------------------------------------------------------------
# One way to create a multidimensional array is to use the `rand` function. `rand` takes an
# arbitrary number of arguments where the number of arguments determines the number of
# dimensions the array created by `rand` will have!
#
# For example, if we pass `rand` two integer inputs, it will generate a 2D array, while
# three integer inputs will generate a 3D array:
# ------------------------------------------------------------------------------------------

rand(4, 3)

rand(4, 3, 2)

# ------------------------------------------------------------------------------------------
# If we want to grab a value from a 2D array, we index into the array, specifying the row
# and column of the element of interest! For example, if we have the array `A` given by
#
# ```julia
# A = [ 1  2  3
#       4  5  6
#       7  8  9 ]
#  ```
#  we could grab the number 6 by saying
#
#  ```
#  A[2, 3]
#  ```
#
#  since 6 is in the 2nd row and 3rd column of `A`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 8
#
#  Copy and execute the following code to get an array, `myprimematrix`, which stores the
# first 100 primes
#
#  ```
#  myprimematrix = reshape(primes(541), (10, 10))
#  ```
#
#  Grab the prime in the 8th row and 5th column via indexing. What is the value of this
# prime?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Copying arrays
#
# Be careful when you want to copy arrays!
#
# Let's say you want to copy an existing array, like `fibonacci`, to an array called
# `somenumbers`. Remember how we defined `fibonacci`:
#
# ```
# fibonacci = [1, 1, 2, 3, 5, 8, 13]
# ```julia
#
# What if we try to copy fibonacci by saying
#
# ```
# somenumbers = fibonacci
# ```julia
# ?
# Execute this code to try to copy fibonacci. Look at `somenumbers` to see if it stores what
# you want.
# ------------------------------------------------------------------------------------------

fibonacci

somenumbers = fibonacci

somenumbers

# ------------------------------------------------------------------------------------------
# Now let's update `somenumbers` by changing its first element.
# ------------------------------------------------------------------------------------------

somenumbers[1] = 404

# ------------------------------------------------------------------------------------------
# Now let's look inside fibonacci.
# ------------------------------------------------------------------------------------------

fibonacci

# ------------------------------------------------------------------------------------------
# #### Exercise 9
#
# What is the first element in `fibonacci`?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Copying or not?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Did copying `fibonacci` like this work?
#
# No, unfortunately not. When we tried to copy, all we did was give `fibonacci` a new name,
# `somenumbers`. Now when we update `somenumbers`, we're also updating `fibonacci` because
# they are the same object!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# If we actually want to make a *copy* of the array bound to `fibonacci`, we can use the
# `copy` function:
# ------------------------------------------------------------------------------------------

# First, restore fibonacci

fibonacci[1] = 1
fibonacci

somemorenumbers = copy(fibonacci)

somemorenumbers[1] = 404

fibonacci

# ------------------------------------------------------------------------------------------
# In this last example, fibonacci was not updated. Therefore we see that the arrays bound to
# `somemorenumbers` and `fibonacci` are distinct.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 10
#
# Copy `myprimematrix` to `mynewprimematrix`. Update `mynewprimematrix[3,3]` to `1234`.
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/02. Tools - Using arrays to store data.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# ------------------------------------------------------------------------------------------
# Often used in machine learning, a "one hot" vector is a vector of all zeros, except for a
# single `1` entry.
# Representing it as a standard vector is memory-inefficient, so it cries out for a special
# implementation.
# ------------------------------------------------------------------------------------------

struct OneHotVector <: AbstractVector{Int}
    idx::Int
    len::Int
end

Base.size(v::OneHotVector) = (v.len,)

Base.getindex(v::OneHotVector, i::Integer) = Int(i == v.idx)

OneHotVector(3, 10)

A = rand(5,5)

A * OneHotVector(3, 5)

Vector(OneHotVector(3,5))

# ------------------------------------------------------------------------------------------
# ## Exercise
#
# Generalize it to any element type.
# ------------------------------------------------------------------------------------------

ext: jl


text:
export LBFGS!, gradients_GPU
using Random

function gradients_GPU(loss::PyObject, vars = missing;use_locking::Bool = false, kwargs...)
    opt = tf.train.Optimizer(use_locking, "default_"*randstring(8))
    vars = coalesce(vars, get_collection())
    grads = opt.compute_gradients(loss, vars, colocate_gradients_with_ops=true)
    [x[1] for x in grads]
end

@doc raw"""
    LBFGS!(sess::PyObject, loss::PyObject, max_iter::Int64=15000;
    vars::Array{PyObject}=PyObject[], callback::Union{Function, Nothing}=nothing, kwargs...)
`BFGS!` is a simplified interface for BFGS optimizer. See also [`ScipyOptimizerInterface`](@ref).
`callback` is a callback function with signature
```julia
callback(vs::Array, iter::Int64, loss::Float64)
```
`vars` is an array consisting of tensors and its values will be the input to `vs`.
# Example 1
```julia
a = Variable(1.0)
loss = (a - 10.0)^2
sess = Session(); init(sess)
BFGS!(sess, loss)
```
# Example 2
```julia
θ1 = Variable(1.0)
θ2 = Variable(1.0)
loss = (θ1-1)^2 + (θ2-2)^2
cb = (vs, iter, loss)->begin
    printstyled("[#iter $iter] θ1=$(vs[1]), θ2=$(vs[2]), loss=$loss\n", color=:green)
end
sess = Session(); init(sess)
cb(run(sess, [θ1, θ2]), 0, run(sess, loss))
BFGS!(sess, loss, 100; vars=[θ1, θ2], callback=cb)
```
# Example 3
Use `bounds` to specify upper and lower bound of a variable.
```julia
x = Variable(2.0)
loss = x^2
sess = Session(); init(sess)
BFGS!(sess, loss, bounds=Dict(x=>[1.0,3.0]))
```
"""->
function LBFGS!(sess::PyObject, loss::PyObject, max_iter::Int64=15000;
    vars::Array{PyObject}=PyObject[], callback::Union{Function, Nothing}=nothing, kwargs...)
    __cnt = 0
    __loss = 0
    __var = nothing
    out = []
    function print_loss(l, vs...)
        if !isnothing(callback); __var = vs; end
        if mod(__cnt,1)==0
            println("iter $__cnt, current loss=",l)
        end
        __loss = l
        __cnt += 1
    end
    __iter = 0
    function step_callback(rk)
        if mod(__iter,1)==0
            println("================ STEP $__iter ===============")
        end
        if !isnothing(callback)
            callback(__var, __iter, __loss)
        end
        push!(out, __loss)
        __iter += 1
    end
    kwargs = Dict(kwargs)
    if haskey(kwargs, :bounds)
        kwargs[:var_to_bounds] = kwargs[:bounds]
        delete!(kwargs, :bounds)
    end
    if haskey(kwargs, :var_to_bounds)
        desc = "`bounds` or `var_to_bounds` keywords of `BFGS!` only accepts dictionaries whose keys are Variables"
        for (k,v) in kwargs[:var_to_bounds]
            if !(hasproperty(k, "trainable"))
                error("The tensor $k does not have property `trainable`, indicating it is not a `Variable`. $desc\n")
            end
            if !k.trainable
                @warn("The variable $k is not trainable, ignoring the bounds")
            end
        end
    end
    opt = ScipyOptimizerInterface(loss, method="L-BFGS-B",options=Dict("maxiter"=> max_iter, "ftol"=>1e-12, "gtol"=>1e-12); kwargs...)
    @info "Optimization starts..."
    ScipyOptimizerMinimize(sess, opt, loss_callback=print_loss, step_callback=step_callback, fetches=[loss, vars...])
    out
end

"""
    LBFGS!(value_and_gradients_function::Function, initial_position::Union{PyObject, Array{Float64}}, max_iter::Int64=50, args...;kwargs...)
Applies the BFGS optimizer to `value_and_gradients_function`
"""
function LBFGS!(value_and_gradients_function::Function,
    initial_position::Union{PyObject, Array{Float64}}, max_iter::Int64=50, args...;kwargs...)
    tfp.optimizer.bfgs_minimize(value_and_gradients_function,
        initial_position=initial_position, args...;max_iterations=max_iter, kwargs...)[5]
end

@doc raw"""
    LBFGS!(sess::PyObject, loss::PyObject, grads::Union{Array{T},Nothing,PyObject},
    vars::Union{Array{PyObject},PyObject}; kwargs...) where T<:Union{Nothing, PyObject}
Running BFGS algorithm
``\min_{\texttt{vars}} \texttt{loss}(\texttt{vars})``
The gradients `grads` must be provided. Typically, `grads[i] = gradients(loss, vars[i])`.
`grads[i]` can exist on different devices (GPU or CPU).
# Example 1
```julia
import Optim # required
a = Variable(0.0)
loss = (a-1)^2
g = gradients(loss, a)
sess = Session(); init(sess)
BFGS!(sess, loss, g, a)
```
# Example 2
```julia
import Optim # required
a = Variable(0.0)
loss = (a^2+a-1)^2
g = gradients(loss, a)
sess = Session(); init(sess)
cb = (vs, iter, loss)->begin
    printstyled("[#iter $iter] a = $vs, loss=$loss\n", color=:green)
end
BFGS!(sess, loss, g, a; callback = cb)
```
"""
function LBFGS!(sess::PyObject, loss::PyObject, grads::Union{Array{T},Nothing,PyObject},
        vars::Union{Array{PyObject},PyObject}; callback::Union{Function, Missing}=missing, kwargs...) where T<:Union{Nothing, PyObject}

    if isa(grads, PyObject); grads = [grads]; end
    if isa(vars, PyObject); vars = [vars]; end
    if length(grads)!=length(vars); error("Length of grads and vars do not match"); end

    idx = ones(Bool, length(grads))
    for i = 1:length(grads)
        if isnothing(grads[i])
            idx[i] = false
        end
    end
    grads = grads[idx]
    vars = vars[idx]

    sizes = []
    for v in vars
        push!(sizes, size(v))
    end
    grds = vcat([tf.reshape(g, (-1,)) for g in grads]...)
    vs = vcat([tf.reshape(v, (-1,)) for v in vars]...); x0 = run(sess, vs)
    pl = placeholder(x0)
    n = 0
    assign_ops = PyObject[]
    for (k,v) in enumerate(vars)
        push!(assign_ops, assign(v, tf.reshape(pl[n+1:n+prod(sizes[k])], sizes[k])))
        n += prod(sizes[k])
    end

    __loss = 0.0
    __losses = Float64[]
    __iter = 0
    __value = nothing
    function f(x)
        run(sess, assign_ops, pl=>x)
        __loss = run(sess, loss)
        return __loss
    end

    function g!(G, x)
        run(sess, assign_ops, pl=>x)
        __value = x
        G[:] = run(sess, grds)
    end

    function callback1(x)
        __iter = x.iteration
        __loss = x.value
        push!(__losses, __loss)
        if !ismissing(callback)
            callback(__value, __iter, __loss)
        end
        false
    end

    Optim.optimize(f, g!, x0, Optim.LBFGS(), Optim.Options(show_trace=true, callback=callback1, kwargs...))
    return __losses
end


text:
#--------------------------------------------------------------------#
#     Input data used for heat and moisture transfer calculation     #
#--------------------------------------------------------------------#

#DODATI RUBNE UVJETE

#########
#Timestep if necessary#
##########

#dt=1        #Timestep in h
#T=8760      #Time domain

###########
#Mesh size#
###########

dl=0.0005 #defined in m

########################
#Geometry and materials#
########################

#Elements -- matrix of vectors E_1=[material_key,x_0,y_0,x_1,y_1] to E_n

ne=12 #number of elements

E=zeros(Float64,ne,9) #Geometry and material matrix + boundary conditions

#Element 1
E[1,1]=12
E[1,2]=0
E[1,3]=0
E[1,4]=0.5
E[1,5]=0.0015
E[1,6]=1
E[1,7]=230
E[1,8]=1
E[1,9]=1

#Element 2
E[2,1]=12
E[2,2]=0
E[2,3]=0.0015
E[2,4]=0.0015
E[2,5]=0.0365
E[2,6]=2
E[2,7]=230
E[2,8]=1
E[2,9]=1

#Element 3
E[3,1]=12
E[3,2]=0.0015
E[3,3]=0.0015
E[3,4]=0.5
E[3,5]=0.035
E[3,6]=3
E[3,7]=0.029
E[3,8]=1
E[3,9]=1

#Element 4
E[4,1]=12
E[4,2]=0
E[4,3]=0.035
E[4,4]=0.015
E[4,5]=0.0365
E[4,6]=4
E[4,7]=230
E[4,8]=1
E[4,9]=1

#Element 5
E[5,1]=12
E[5,2]=0.015
E[5,3]=0.035
E[5,4]=0.5
E[5,5]=0.0365
E[5,6]=5
E[5,7]=0.029
E[5,8]=1
E[5,9]=1

#Element 6
E[6,1]=12
E[6,2]=0
E[6,3]=0.0365
E[6,4]=0.015
E[6,5]=0.0415
E[6,6]=6
E[6,7]=0.12
E[6,8]=1
E[6,9]=1

#Element 7
E[7,1]=12
E[7,2]=0.015
E[7,3]=0.0365
E[7,4]=0.5
E[7,5]=0.0415
E[7,6]=7
E[7,7]=0.029
E[7,8]=1
E[7,9]=1

#Element 8
E[8,1]=12
E[8,2]=0
E[8,3]=0.0415
E[8,4]=0.5
E[8,5]=0.0475
E[8,6]=8
E[8,7]=1.15
E[8,8]=1
E[8,9]=1

##############################
#Boundary condition materials#
##############################

#Element 17 (dolje)
E[9,1]=100
E[9,2]=0
E[9,3]=-dl
E[9,4]=0.5
E[9,5]=0
E[9,6]=17
E[9,7]=(9.09)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[9,8]=20
E[9,9]=1

#Element 18 (lijevo)
E[10,1]=101
E[10,2]=-dl
E[10,3]=0
E[10,4]=0
E[10,5]=0.0475
E[10,6]=18
E[10,7]=(1)*2/dl #Assign as (h)*dl/2 if h is convection coef.
E[10,8]=0
E[10,9]=-1

#Element 19 (gore)
E[11,1]=102
E[11,2]=0
E[11,3]=0.0475
E[11,4]=0.5
E[11,5]=0.0475+dl
E[11,6]=19
E[11,7]=(16.67)*dl/2 #Assign as (h)*dl/2 if h is convection coef.
E[11,8]=0
E[11,9]=1

#Element 20 (desno)
E[12,1]=103
E[12,2]=0.5
E[12,3]=0
E[12,4]=0.5+dl
E[12,5]=0.0475
E[12,6]=20
E[12,7]=(1)*2/dl #Assign as (h)*dl/2 if h is convection coef.
E[12,8]=20
E[12,9]=-1


text:
using FranklinCLI
using Test

@testset "FranklinCLI.jl" begin
    # Write your tests here.
end


text:
using Documenter
using Distances

makedocs(sitename="Distances.jl")


text:
@testset "inplace keyword" begin
    rosenbrock = MultivariateProblems.UnconstrainedProblems.examples["Rosenbrock"]
    f = MVP.objective(rosenbrock)
    g! = MVP.gradient(rosenbrock)
    h! = MVP.hessian(rosenbrock)
    function g(x)
        G = similar(x)
        g!(G, x)
        G
    end
    function h(x)
        n = length(x)
        H = similar(x, n, n)
        h!(H, x)
        H
    end
    initial_x = rosenbrock.initial_x

    inp_res = optimize(f, g, h, initial_x; inplace = false)
    op_res  = optimize(f, g!, h!, initial_x; inplace = true)

    for op in (Optim.minimizer, Optim.minimum, Optim.f_calls,
                Optim.g_calls, Optim.h_calls, Optim.iterations, Optim.converged)
        @test all(op(inp_res) .=== op(op_res))
    end
end


text:
### Guru Interfaces
finufft_plan{T} = Ptr{T} where T <: Union{Float32,Float64}

function finufft_makeplan(type::Integer,
                          dim::Integer,
                          n_modes::T,
                          iflag::Integer,
                          ntrans::Integer,
                          eps::Float64,
                          opts::nufft_opts=finufft_default_opts()) where T <: Union{Array{BIGINT}, NTuple{3,BIGINT}}

    # see https://stackoverflow.com/questions/40140699/the-proper-way-to-declare-c-void-pointers-in-julia for how to declare c-void pointers in julia
    plan_p = Ref{finufft_plan{Float64}}()

    ret = ccall( (:finufft_makeplan, libfinufft),
            Cint,
            (Cint,
            Cint,
            Ref{BIGINT},
            Cint,
            Cint,
            Cdouble,
            Ptr{finufft_plan{Cdouble}},
            Ref{nufft_opts}),
            type,dim,n_modes,iflag,ntrans,eps,plan_p,opts
            )
    check_ret(ret)

    plan = plan_p[]
    return plan
end
function finufft_makeplan(type::Integer,
    dim::Integer,
    n_modes::T,
    iflag::Integer,
    ntrans::Integer,
    eps::Float32,
    opts::nufft_opts=finufft_default_opts()) where T <: Union{Array{BIGINT}, NTuple{3,BIGINT}}

    # see https://stackoverflow.com/questions/40140699/the-proper-way-to-declare-c-void-pointers-in-julia for how to declare c-void pointers in julia
    plan_p = Ref{finufft_plan{Float32}}()

    ret = ccall( (:finufftf_makeplan, libfinufft),
            Cint,
            (Cint,
            Cint,
            Ref{BIGINT},
            Cint,
            Cint,
            Cfloat,
            Ptr{finufft_plan{Cfloat}},
            Ref{nufft_opts}),
            type,dim,n_modes,iflag,ntrans,eps,plan_p,opts
            )
    check_ret(ret)

    plan = plan_p[]
    return plan
end

function finufft_setpts(plan::finufft_plan{Float64},
                        xj::Array{Float64},
                        yj::Array{Float64}=Float64[],
                        zj::Array{Float64}=Float64[],
                        s::Array{Float64}=Float64[],
                        t::Array{Float64}=Float64[],
                        u::Array{Float64}=Float64[])

    nj = length(xj)
    nk = length(s)

    ret = ccall( (:finufft_setpts, libfinufft),
                Cint,
                (finufft_plan{Cdouble},
                BIGINT,
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{Cdouble},
                BIGINT,
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{Cdouble}),
                plan,nj,xj,yj,zj,nk,s,t,u
                )

    check_ret(ret)
end
function finufft_setpts(plan::finufft_plan{Float32},
    xj::Array{Float32},
    yj::Array{Float32}=Float32[],
    zj::Array{Float32}=Float32[],
    s::Array{Float32}=Float32[],
    t::Array{Float32}=Float32[],
    u::Array{Float32}=Float32[])

    nj = length(xj)
    nk = length(s)

    ret = ccall( (:finufftf_setpts, libfinufft),
                Cint,
                (finufft_plan{Cfloat},
                BIGINT,
                Ref{Cfloat},
                Ref{Cfloat},
                Ref{Cfloat},
                BIGINT,
                Ref{Cfloat},
                Ref{Cfloat},
                Ref{Cfloat}),
                plan,nj,xj,yj,zj,nk,s,t,u
                )

    check_ret(ret)
end

function finufft_exec(plan::finufft_plan{Float64}, cj::Array{ComplexF64}, fk::Array{ComplexF64})

    ret = ccall( (:finufft_execute, libfinufft),
                Cint,
                (finufft_plan{Cdouble},
                Ref{ComplexF64},
                Ref{ComplexF64}),
                plan,cj,fk
                )

    check_ret(ret)
end
function finufft_exec(plan::finufft_plan{Float32}, cj::Array{ComplexF32}, fk::Array{ComplexF32})

    ret = ccall( (:finufftf_execute, libfinufft),
                Cint,
                (finufft_plan{Cfloat},
                Ref{ComplexF32},
                Ref{ComplexF32}),
                plan,cj,fk
                )

    check_ret(ret)
end

function finufft_destroy(plan::finufft_plan{Float64})
    ret = ccall( (:finufft_destroy, libfinufft),
                Cint,
                (finufft_plan{Cdouble},),
                plan
                )

    check_ret(ret)
end
function finufft_destroy(plan::finufft_plan{Float32})
    ret = ccall( (:finufftf_destroy, libfinufft),
            Cint,
            (finufft_plan{Cfloat},),
            plan
            )

    check_ret(ret)
end


text:
@reexport module CapacityExpansion

# Standard library
using Statistics
using Printf
using Dates

using RecipesBase
using Distributions
using Parameters
using StochasticPrograms
using Flux
using HydroModels
using HydroModels: AbstractHydroModel, StochasticHydroModel, River, Plant, Area, Scenario, Q̄, μ, %

import HydroModels: modelindices

using Plots: font, text, Shape

export
    CapacityExpansionData,
    CapacityExpansionScenario,
    RecurrentCapacityExpansionSampler,
    CapacityExpansionModel

include("data.jl")
include("scenarios.jl")
include("indices.jl")
include("model.jl")

end


text:
# --- list frequency of 2-node and 3-node motifs in a network
# ---  don't care about multiple representations of the same motif

# --- counts 2-node motif frequency in nw
function motifcount2(nw)
    n = length(nw[:,1])
    ret = zeros(Float64,9)
    for i=1:n
        for j=1:n
            ret[1+3*(1+nw[i,j])+(1+nw[j,i])] += 1
        end
    end
    ret /= n^2
end

# --- returns motif id from the edge list
function motifmap(ij,ji,ik,ki,jk,kj)
    nij = 3*(1+ij)+(1+ji)
    nik = 3*(1+ik)+(1+ki)
    njk = 3*(1+jk)+(1+kj)
    return 1+81*nij+9*nik+njk
end

# --- returns edge list from the motif id
function motifbackmap(motifnum)
    motifnum -= 1
    i = fld(motifnum,81)
    nij = fld(i,3)-1
    nji = mod(i,3)-1
    motifnum = mod(motifnum,81)
    i = fld(motifnum,9)
    nik = fld(i,3)-1
    nki = mod(i,3)-1
    i = mod(motifnum,9)
    njk = fld(i,3)-1
    nkj = mod(i,3)-1
    return (nij,nji,nik,nki,njk,nkj)
end

motifname = Dict();
motif2id = Dict();
id2motif = Dict();

id=0;
for ij=-1:1:1 for ji=-1:1:1 for ik=-1:1:1 for ki=-1:1:1 for jk=-1:1:1 for kj=-1:1:1
    mini = min(motifmap(ij,ji,ik,ki,jk,kj),motifmap(ji,ij,jk,kj,ik,ki),motifmap(kj,jk,ki,ik,ji,ij),motifmap(ik,ki,ij,ji,kj,jk),motifmap(jk,kj,ji,ij,ki,ik),motifmap(ki,ik,kj,jk,ij,ji))
    motifname[motifmap(ij,ji,ik,ki,jk,kj)] = mini
    if (!haskey(motif2id,mini))
        id += 1
        motif2id[mini] = id
        id2motif[id] = mini
    end
end end end end end end

# --- counts 3-node motif frequency in nw
function motifcount3(nw)
    n = length(nw[:,1])
#    ret = zeros(Float64,729)            # 3^6 edge combinations
    ret = zeros(Float64,138)            # 138 = # of distinct 3-motifs
    for i=1:n for j=1:n for k=1:n
        id = motif2id[motifname[motifmap(nw[i,j],nw[j,i],nw[i,k],nw[k,i],nw[j,k],nw[k,j])]]
        ret[id] += 1
    end end end
    ret /= n^3
end

max_stars_repo_name: akabakcioglu/network-dynamics

max_issues_repo_name: akabakcioglu/network-dynamics

max_forks_repo_name: akabakcioglu/network-dynamics


text:
using Ipopt
using MathOptInterface

const MOI = MathOptInterface

@testset "MOI wrapper" begin
    CASE57_SOLUTION = [
        1.0260825400262428,
        1.016218932360429,
        1.009467718490475,
        1.027857273911734,
        1.06,
        0.9962215167521776,
        0.9965415804936182,
        0.0,
        0.6,
        0.0,
        8.603379123083476,
        0.0,
        1.3982297290334753,
    ]

    datafile = joinpath(INSTANCES_DIR, "case57.m")
    nlp = ExaPF.ReducedSpaceEvaluator(datafile)

    @testset "ReducedSpaceEvaluator with L-BFGS" begin
        optimizer = Ipopt.Optimizer()
        MOI.set(optimizer, MOI.RawParameter("print_level"), 0)
        MOI.set(optimizer, MOI.RawParameter("limited_memory_max_history"), 50)
        MOI.set(optimizer, MOI.RawParameter("hessian_approximation"), "limited-memory")
        MOI.set(optimizer, MOI.RawParameter("tol"), 1e-4)

        solution = ExaPF.optimize!(optimizer, nlp)
        MOI.empty!(optimizer)
        @test solution.status ∈ [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
        @test solution.minimum ≈ 3.7589338e+04
        @test solution.minimizer ≈ CASE57_SOLUTION rtol=1e-5
    end

    # Test resolution with AugLagEvaluator and Ipopt, as used inside ProxAL
    @testset "AugLagEvaluator with Hessian" begin
        optimizer = Ipopt.Optimizer()
        MOI.set(optimizer, MOI.RawParameter("print_level"), 0)
        MOI.set(optimizer, MOI.RawParameter("tol"), 1e-4)
        MOI.set(optimizer, MOI.RawParameter("max_iter"), 30)
        aug = ExaPF.AugLagEvaluator(nlp, ExaPF.initial(nlp))

        solution = ExaPF.optimize!(optimizer, aug)
        MOI.empty!(optimizer)
        @test solution.status ∈ [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
        @test solution.minimum ≈ 25138.76310420395
    end
end


text:
using GR

const maps=["uniform", "temperature", "grayscale", "glowing", "rainbowlike", "geologic", "greenscale", "cyanscale", "bluescale", "magentascale", "redscale", "flame", "brownscale", "pilatus", "autumn", "bone", "cool", "copper", "gray", "hot", "hsv", "jet", "pink", "spectral", "spring", "summer", "winter", "gist_earth", "gist_heat", "gist_ncar", "gist_rainbow", "gist_stern", "afmhot", "brg", "bwr", "coolwarm", "CMRmap", "cubehelix", "gnuplot", "gnuplot2", "ocean", "rainbow", "seismic", "terrain", "viridis", "inferno", "plasma", "magma"]

function main()
    a = round.(Int32,LinRange(1000,1255,256))
    xl = 0
    setwsviewport(0, 0.25, 0, 0.125)
    setwswindow(0, 1, 0, 0.5)
    setviewport(0.05, 0.95, 0.025, 0.475)
    setcharheight(0.010)
    setcharup(-1,0)
    settextcolorind(255)
    for cmap in 0:47
        setcolormap(cmap)
        xr = (cmap + 1) / 48.0
        cellarray(xl+0.002, xr-0.002, 0.2, 1, 1, 256, a)
        settextalign(1,3)
        text(0.04 + xr * 0.9, 0.48, string(cmap))
        settextalign(3,3)
        text(0.04 + xr * 0.9, 0.1, string(maps[cmap+1]))
        xl = xr
    end
end

main()


text:
#############################################################
# A convenient macro to define a Composite Type, with
# default values for fields and a constructor that accept
# keyword parameters to initialize the fields. For example
#
# @defstruct MyStruct Any (
#   field1 :: Int = 0,
#   (field2 :: AbstractString = "", !isempty(field2))
# )
#
# where each field could be either
#
#   field_name :: field_type = default_value
#
# or put within a tuple, with the second element
# specifying a validation check on the field value.
# In the example above, the default value for
# field2 does not satisfy the assertion, this
# could be used to force user to provide a
# valid value when no meaningful default value
# is available.
#
# The macro will define a constructor that could accept
# the keyword arguments. Since the defined type is
# immutable, a "copy" function is also defined, which
# takes a prototype object, and accept extra keyword
# parameters that could be used to construct a new
# object with specified changes of fields.
#############################################################
import Base.copy
export      copy

macro defstruct(name, super_name, fields)
  @assert fields.head == :tuple
  fields = fields.args
  @assert length(fields) > 0
  name = esc(name)

  field_defs     = Array(Expr, length(fields))     # :(field2 :: Int)
  field_names    = Array(Symbol, length(fields))   # :field2
  field_defaults = Array(Expr, length(fields))     # :(field2 :: Int = 0)
  field_asserts  = Array(Expr, length(fields))     # :(field2 >= 0)

  for i = 1:length(fields)
    field = fields[i]
    if field.head == :tuple
      field_asserts[i] = field.args[2]
      field = field.args[1]
    end
    field_defs[i] = esc(field.args[1])
    field_names[i] = field.args[1].args[1]
    field_defaults[i] = Expr(:kw, field.args...)
  end

  # body of layer type, defining fields
  type_body = Expr(:block, field_defs...)

  # constructor
  asserts = map(filter(i -> isassigned(field_asserts,i), 1:length(fields))) do i
    :(@assert($(field_asserts[i])))
  end
  construct = Expr(:call, name, field_names...)
  ctor_body = Expr(:block, asserts..., construct)
  ctor_def = Expr(:call, name, Expr(:parameters, field_defaults...))
  ctor = Expr(:(=), ctor_def, ctor_body)

  # for copy constructor
  field_assigns = Expr(:block, [:(params[symbol($(esc(string(fname))))] = proto.$fname) for fname in field_names]...)
  field_expose = Expr(:block, [:($(esc(fname)) = params[symbol($(esc(string(fname))))]) for fname in field_names]...)
  assert_block = Expr(:block, asserts...)
  obj_construct = Expr(:call, name, field_names...)
  copy_fname = esc(:copy)

  quote
    immutable $(name) <: $super_name
      $type_body
    end

    $ctor

    function $copy_fname(proto::$name; kw...)
      params = Dict()
      $field_assigns

      for (k,v) in kw
        @assert haskey(params, k) "Unrecognized field " * string(k) * " for " * $(string(name.args[1]))
        params[k] = v
      end

      $field_expose
      $assert_block

      $obj_construct
    end
  end
end

#############################################################
# A macro used to characterize a layer. Example
#
# @characterize_layer(HDF5DataLayer,
#     is_source => true,
#     is_sink => false,
#     can_do_bp => false,
# )
#############################################################
macro characterize_layer(layer, properties...)
  defs = Array(Expr, length(properties))
  for (i,prop) in enumerate(properties)
    @assert(isa(prop, Expr) && prop.head == :(=>), "Property should be: property_name => value")

    prop_name = prop.args[1]
    prop_val  = prop.args[2]
    defs[i] = quote
      $(esc(prop_name))(::$(esc(layer))) = $prop_val
    end
  end

  Expr(:block, defs...)
end


text:
using Plots
rectangle(w, h, x, y) = Shape(x .+ [0, w, w, 0], y .+ [0, 0, h, h])

function plot_search_space(grid, com_grid, fname)

    plot(
        0:9,
        0:9,
        size = (500, 500),
        legend = false,
        xaxis = false,
        yaxis = false,
        aspect_ratio = :equal,
    )
    for i in 0:8, j in 0:8
        plot!(rectangle(1, 1, i, j), color = :white)
    end
    for i in [3, 6]
        plot!([i, i], [0, 9], color = :black, linewidth = 4)
    end
    for j in [3, 6]
        plot!([0, 9], [j, j], color = :black, linewidth = 4)
    end

    for vidx in keys(com_grid)
        if CS.isfixed(com_grid[vidx])
            x = vidx[2] - 0.5
            y = 10 - vidx[1] - 0.5
            if grid[vidx] != 0
                annotate!([x, y, text(string(CS.value(com_grid[vidx])), 20, :black)])
            else
                annotate!([x, y, text(string(CS.value(com_grid[vidx])), 20, :blue)])
            end
        else
            vals = CS.values(com_grid[vidx])
            sort!(vals)
            x = vidx[2] - 0.5
            y = 10 - vidx[1] - 0.2
            max_idx = min(3, length(vals))
            text_vals = join(vals[1:max_idx], ",")
            annotate!([x, y, text(text_vals, 7)])
            if length(vals) > 3
                y -= 0.25
                max_idx = min(6, length(vals))
                text_vals = join(vals[4:max_idx], ",")
                annotate!([x, y, text(text_vals, 7)])
            end
            if length(vals) > 6
                y -= 0.25
                max_idx = length(vals)
                text_vals = join(vals[7:max_idx], ",")
                annotate!([x, y, text(text_vals, 7)])
            end

        end
    end
    if fname != ""
        png("visualizations/images/$(fname)")
    end
end

function plot_killer(grid, sums, fname; fill = true, mark = nothing)
    plot(;
        size = (900, 900),
        legend = false,
        xaxis = false,
        yaxis = false,
        aspect_ratio = :equal,
    )
    for s in sums
        vidx = s.indices[1]
        x = vidx[2] - 0.75
        y = 10 - vidx[1] - 0.17
        annotate!(x, y, text(s.result, 15, :black))
        if mark === nothing
            for vidx in s.indices
                plot!(
                    rectangle(1, 1, vidx[2] - 1, 9 - vidx[1]),
                    color = s.color,
                    alpha = 0.4,
                )
            end
        else
            for vidx in s.indices
                if vidx in mark
                    plot!(
                        rectangle(1, 1, vidx[2] - 1, 9 - vidx[1]),
                        color = "red",
                        alpha = 0.4,
                    )
                else
                    plot!(
                        rectangle(1, 1, vidx[2] - 1, 9 - vidx[1]),
                        color = s.color,
                        alpha = 0.4,
                    )
                end
            end
        end
    end

    for i in [3, 6]
        plot!([i, i], [0, 9], color = :black, linewidth = 4)
    end
    for j in [3, 6]
        plot!([0, 9], [j, j], color = :black, linewidth = 4)
    end

    if fill
        for vidx in keys(grid)
            if CS.isfixed(grid[vidx])
                x = vidx[2] - 0.5
                y = 10 - vidx[1] - 0.5
                annotate!(x, y, text(value(grid[vidx]), 20, :black))
            else
                vals = values(grid[vidx])
                sort!(vals)
                x = vidx[2] - 0.3
                y = 10 - vidx[1] - 0.2
                max_idx = min(3, length(vals))
                text_vals = join(vals[1:max_idx], ",")
                annotate!(x, y, text(text_vals, 11))
                if length(vals) > 3
                    y -= 0.25
                    max_idx = min(6, length(vals))
                    text_vals = join(vals[4:max_idx], ",")
                    annotate!(x, y, text(text_vals, 11))
                end
                if length(vals) > 6
                    y -= 0.25
                    max_idx = length(vals)
                    text_vals = join(vals[7:max_idx], ",")
                    annotate!(x, y, text(text_vals, 11))
                end

            end
        end
    end

    png("/home/ole/Julia/ConstraintSolver/visualizations/images/$(fname)")
end

function plot_str8ts(grid, white, fname; search_space = nothing, only_fixed = true)
    plot(;
        size = (900, 900),
        legend = false,
        xaxis = false,
        yaxis = false,
        aspect_ratio = :equal,
    )

    for r in 1:9, c in 1:9
        plot!(rectangle(1, 1, c - 1, 9 - r), color = white[r, c] == 1 ? "white" : "black")
        x = c - 0.5
        y = 10 - r - 0.5
        if grid[r, c] != 0
            annotate!(
                x,
                y,
                text(string(grid[r, c]), 20, white[r, c] == 1 ? :black : :white),
            )
        end
        if search_space !== nothing
            if search_space[r, c] != 0 &&
               (search_space[r, c] isa Int || length(search_space[r, c]) == 1)
                @assert white[r, c] == 1
                annotate!(x, y, text(string(search_space[r, c]), 20, :blue))
            end
        end
    end

    png("/home/ole/Julia/ConstraintSolver/visualizations/images/$(fname)")
end


text:
#include("cr_sod_shock_main.jl")

using Roots

mutable struct SodCRParameters_noCRs

    rhol::Float64
    rhor::Float64
    Pl::Float64
    Pr::Float64
    Ul::Float64
    Ur::Float64
    cl::Float64
    cr::Float64
    M::Float64
    t::Float64
    x_contact::Float64
    Pe_ratio::Float64
    γ_th::Float64
    γ_cr::Float64
    γ_exp::Float64
    η2::Float64
    acc_function
    ξ::Float64
    first_guess::Float64

    function SodCRParameters_noCRs(;rhol::Float64=1.0,  rhor::Float64=0.125,
                                Pl::Float64=0.0,    Pr::Float64=0.0,
                                Ul::Float64=0.0,    Ur::Float64=0.0,
                                Mach::Float64=0.0,  t::Float64,
                                x_contact::Float64=70.0,
                                Pe_ratio::Float64=0.01,
                                γ_th::Float64=5.0/3.0,
                                γ_cr::Float64=4.0/3.0,
                                thetaB::Float64=0.0,
                                theta_crit::Float64=(π/4.0),
                                dsa_model::Int64=-1,
                                xs_first_guess::Float64=4.7)

        γ_exp    = ( γ_th - 1.0 )/( 2.0 * γ_th )
        η2       = (γ_th-1.0)/(γ_th+1.0)

        # calculate Ul and Pl depending on input
        if (Pl == 0.0) & (Ul != 0.0)
            Pl = ( γ_th - 1.0 ) * rhol * Ul
        elseif (Ul == 0.0) & (Pl != 0.0)
            Ul = Pl / ( (γ_th - 1.0) * rhol )
        else
            error("Both Ul and Pl are zero!")
        end

        if dsa_model == 0
            acc_function = KR07_acc
        elseif dsa_model == 1
            acc_function = KR13_acc
        elseif dsa_model == 2
            acc_function = Ryu19_acc
        elseif dsa_model == 3
            acc_function = CS14_acc
        elseif dsa_model == 4
            acc_function = P16_acc
        else
            error("Invalid DSA model selection!")

        end

        # calculate Ur and Pr depending on input
        if (Pr == 0.0) & (Ur != 0.0)
            Pr = ( γ_th - 1.0 ) * rhor * Ur
        elseif (Ur == 0.0) & (Pr != 0.0)
            Ur = Pr / ( (γ_th - 1.0) * rhor )
        elseif (Ur == 0.0) & (Pr == 0.0) & (Mach == 0.0)
            error("Ur, Pr and Mach are zero! Can't find solution!")
        else
            println("Both Ur and Pr are zero! Will calculate them depending on Machnumber.")
            Pr = solvePrfromMachCR(Pl, Pr,
                                   rhor, rhol,
                                   γ_th, γ_cr, γ_exp,
                                   Mach, acc_function)
            Ur = Pr / ( (γ_th - 1.0) * rhor )
        end

        # calculate B angle dependent efficiency following Pais+ 2018, MNRAS, 478, 5278
        delta_theta = π/18.0
        thetaB *= (π/180.0)
        etaB = 0.5*( tanh( (theta_crit - thetaB)/delta_theta ) + 1.0 )

        ξ = etaB*acc_function(Mach)/(1.0 - etaB*acc_function(Mach))

        # if dsa_model == 4
        #     ξ = etaB*acc_function(Mach)/(1.0 - etaB*acc_function(Mach))
        # else
        #     ξ = etaB*acc_function(Mach)
        # end
        #ξ = acc_function(Mach)/(1.0 - acc_function(Mach))

        cl = sqrt.( γ_th * Pl / rhol)
        cr = sqrt.( γ_th * Pr / rhor)

        new(rhol, rhor,
            Pl, Pr,
            Ul, Ur,
            cl, cr,
            Mach, t,
            x_contact,
            Pe_ratio,
            γ_th, γ_cr,
            γ_exp,
            η2,
            acc_function,
            ξ,
            xs_first_guess)

    end
end

# Riemann solver for standard hydro shocktube following Pfrommer et al 2006
# Use left and right state of initial condition to solve shocktube problems

"""
    Helper functions to set up Machnumber dependent IC
"""
function solvePrCR(Pr::Float64,
                 rhol::Float64, rhor::Float64,
                 Pl::Float64, Mach::Float64,
                 γ::Float64, γ_cr::Float64,
                 eff_function)

    γ_1 = γ - 1.0
    γ_pow = γ_1/(2.0*γ)
    η2 = (γ-1.0)/(γ+1.0)

    c_l = sqrt.(γ*Pl/rhol)
    c_r = sqrt.(γ*Pr/rhor)

    f(x) = ( x/Pr - 1.0 ) * sqrt.( (1.0 - η2) / (γ * ( x/Pr + η2 ) )) -
            2.0 / γ_1 * c_l / c_r * ( 1.0 - ( x/Pl )^γ_pow )

    P_m = find_zero(f, (Pr, Pl), Bisection())

    # CR part

    P_cr = eff_function(Mach) * P_m

    vm = 2.0 * c_l / γ_1 * ( 1.0 - (P_m/Pl)^γ_pow )

    C = 2.0*(γ - γ_cr)*P_cr/( γ_1* vm^2 * (γ_cr - 1.0) )

    ρ_mr0 = rhor * ( ( P_m + η2 * Pr ) / ( Pr + η2 * P_m))
    #ρ_ml = rhol * ( P_m / Pl )^(1.0/γ)

    find_rho1(rho1) = ρ_mr0/(1.0 - C/rho1) - rho1

    ρ_mr = find_zero( find_rho1, ρ_mr0 )

    vm = 2.0 * c_l / γ_1 * ( 1.0 - (P_m/Pl)^γ_pow )

    vs = vm / ( 1.0 - rhor/ρ_mr)

    M = vs / c_r

    return M - Mach
end

function MachSolver_HelperFunction(Pl::Float64, Pr::Float64,
                                   rhor::Float64, rhol::Float64,
                                   γ_th::Float64, γ_cr::Float64, γ_exp::Float64,
                                   M::Float64, acc_function, etaB::Float64)

    vs = 0.0

    cr = sqrt.( γ_th * Pr / rhor)
    cl = sqrt.( γ_th * Pl / rhol)

    ξ = etaB * acc_function(M)/( 1.0 - acc_function(M))

    if ξ > 0.0
        # solve Density
        f_z(xs) = find_xs(xs, rhor, Pr, Pl, cr, cl,
                            γ_th, γ_cr, γ_exp, ξ)

        xs = find_zero( f_z, 3.0 )

        rho4 = xs*rhor

        # solve pressure
        P34_tot = P4_f(xs, Pr, γ_th, γ_cr, ξ)

        # solve velocity
        v34 = 2.0 * cl / ( γ_th - 1.0 ) * ( 1.0 - (P34_tot / Pl)^γ_exp )

        vs = v34 / ( 1.0 - rhor / rho4 )
    end
    Mach = vs/cr

    return M - Mach
end

function solveMachfromPrCR(Pl::Float64, Pr::Float64,
                     rhor::Float64, rhol::Float64,
                     γ_th::Float64, γ_cr::Float64, γ_exp::Float64,
                     M::Float64, acc_function)

   f(M) = MachSolver_HelperFunction(Pl, Pr,
                                    rhor, rhol,
                                    γ_th, γ_cr, γ_exp,
                                    M, acc_function)

   Mach = find_zero(f, 5.0)

   return Mach
end

function solvePrfromMachCR(Pl::Float64, Pr::Float64,
                     rhor::Float64, rhol::Float64,
                     γ_th::Float64, γ_cr::Float64, γ_exp::Float64,
                     M::Float64, acc_function)

    f(Pr) = MachSolver_HelperFunction(Pl, Pr,
                                     rhor, rhol,
                                     γ_th, γ_cr, γ_exp,
                                     M, acc_function)

    Pr = find_zero(f, 0.03)

    return Pr
end

"""
    Shared functions
"""
function get_ξ(acc_function, M::Float64)
    return ξ = acc_function(M)/(1.0 - acc_function(M))
end

function xs_f(rho4::Float64, rhor::Float64)
    return rho4/rhor
end

function ys_f(xs::Float64, γ_th::Float64, γ_cr::Float64, ξ::Float64)
    return ( ξ * ( xs*(γ_cr - 1.0) - (γ_cr + 1.0) )*xs^γ_th - xs*(γ_th + 1.0) + (γ_th - 1.0) ) /
           ( ξ * ( xs*(γ_cr - 1.0) - (γ_cr + 1.0) )         + xs*(γ_th - 1.0) - (γ_th + 1.0) )
end

function A_f(xs::Float64, γ_th::Float64, γ_cr::Float64, ξ::Float64)
    #return (rho4 - rhor)/(rho4 + rhor)
    return ( ys_f(xs, γ_th, γ_cr, ξ)        + ξ * ( ys_f(xs, γ_th, γ_cr, ξ) - xs^γ_th)         - 1.0 ) /
           ( ys_f(xs, γ_th, γ_cr, ξ) * γ_th + ξ * ( ys_f(xs, γ_th, γ_cr, ξ) - xs^γ_th ) * γ_cr + γ_th )
end

function P_inj_f(xs::Float64, rhor::Float64, Pr::Float64, γ_th::Float64, γ_cr::Float64, ξ::Float64)
    return (( γ_cr - 1.0)/(γ_th - 1.0)) * ξ *
           ( ys_f(xs, γ_th, γ_cr, ξ) - xs^γ_th ) * Pr
end

function P4_f(xs::Float64, Pr::Float64, γ_th::Float64, γ_cr::Float64, ξ::Float64)
    ys(xs) = ys_f(xs, γ_th, γ_cr, ξ)

    return ( ys(xs) + ( γ_cr - 1.0)/(γ_th - 1.0) *
                ξ * ( ys(xs) - xs^γ_th ) )  * Pr
end

function find_xs(xs::Float64, rhor::Float64,
                 Pr::Float64,   Pl::Float64,
                 cr::Float64,   cl::Float64,
                 γ_th::Float64, γ_cr::Float64, γ_exp::Float64,
                 ξ::Float64)

    P4(xs) = P4_f(xs, Pr, γ_th, γ_cr, ξ)
    A(xs) = A_f(xs, γ_th, γ_cr, rhor)

    return ( P4(xs)/Pr - 1.0 ) * A(xs)/( 1.0 + A(xs))  -
              2.0 * γ_th / ( γ_th - 1.0 )^2 *
              cl^2 / cr^2 *
           ( 1.0 - (P4(xs)/Pl)^γ_exp )^2
end

"""
    Functions to solve individual segments of the shocktube
"""

"""
    Pressure
"""
function solveP34(par::SodCRParameters_noCRs, sol::SodCRSolution)
    # Solves the pressure of the middle state in a shocktube
    return P4_f(sol.xs, par.Pr, par.γ_th, par.γ_cr, par.ξ)
end

function solveP2(x::Float64, par::SodCRParameters_noCRs)
    # solves the pressure along the rarefaction wave
    return par.Pl * ( -par.η2 * x / ( par.cl * par.t ) + ( 1.0 - par.η2 ) )^( 1.0/par.γ_exp )
end

function solveP(x::Float64, par::SodCRParameters_noCRs, sol::SodCRSolution)
    # returns P values according to position in shocktube

    if x <= -par.cl * par.t
        return par.Pl, # P_tot
               par.Pl, # P_th
               0.0,    # P_cr
               0.0,    # P_cr_p
               0.0     # P_cr_e
    elseif -par.cl * par.t < x <= -sol.vt * par.t
        P2 = solveP2(x, par)
        return P2,
               P2,
               0.0,
               0.0,
               0.0
    elseif -sol.vt * par.t < x <= sol.v34 * par.t
        return sol.P34_tot,
               sol.P34_tot,
               0.0,
               0.0,
               0.0
    elseif sol.v34 * par.t < x <= sol.vs * par.t
        return sol.P34_tot,
               sol.P4_th,
               sol.P4_cr,
               sol.P4_cr * ( 1.0 - par.Pe_ratio),
               sol.P4_cr * par.Pe_ratio
    elseif sol.vs * par.t < x
        return par.Pr,
               par.Pr,
               0.0,
               0.0,
               0.0
    else
        println("Error!")
        return 0.0, 0.0, 0.0, 0.0, 0.0
    end
end

"""
    Density
"""
function solveRho2(x::Float64, par::SodCRParameters_noCRs)
    # solves density along the rarefaction wave
    return par.rhol * ( -par.η2 * x / ( par.cl * par.t ) +
                      ( 1.0 - par.η2 ) )^( 2.0 / (par.γ_th - 1.0 ) )
end

function solveRho3(par::SodCRParameters_noCRs, sol::SodCRSolution)
    # solves density left of contact discontiuity and right of refraction wave
    return par.rhol * ( sol.P34_tot/par.Pl )^( 1.0/par.γ_th )
end

function solveRho4(par::SodCRParameters_noCRs)
    # solves density left of shockfront following
    # Eq. B7-B11 , Pfrommer et. al. 2016

    f_z(xs) = find_xs(xs, par.rhor, par.Pr, par.Pl, par.cr, par.cl,
                        par.γ_th, par.γ_cr, par.γ_exp, par.ξ)

    xs = find_zero( f_z, par.first_guess )

    rho4 = xs*par.rhor

    return rho4

end

function solveRho(x::Float64, par::SodCRParameters_noCRs, sol::SodCRSolution)
    # returns P values according to position in shocktube

    if x <= -par.cl * par.t
        return par.rhol
    elseif -par.cl * par.t < x <= -sol.vt * par.t
        return solveRho2(x, par)
    elseif -sol.vt * par.t < x <= sol.v34 * par.t
        return sol.rho3
    elseif sol.v34 * par.t < x <= sol.vs * par.t
        return sol.rho4
    elseif sol.vs * par.t < x
        return par.rhor
    else
        println("Error!")
        return 0.0
    end
end

"""
    Velocity
"""
function solveV34(par::SodCRParameters_noCRs, sol::SodCRSolution)
    # solves velocity along the isopressure region 2-3
    return 2.0 * par.cl / ( par.γ_th - 1.0 ) * ( 1.0 - (sol.P34_tot / par.Pl)^par.γ_exp )
end

function solveV2(x::Float64, par::SodCRParameters_noCRs)
    # solves the velocity along the rarefaction wave
    return ( 1.0 - par.η2 ) * ( x / par.t + par.cl )
end

function solveVs(par::SodCRParameters_noCRs, sol::SodCRSolution)
    # solves the shock velocity
    return sol.v34 / ( 1.0 - par.rhor / sol.rho4 )
end

function solveVt(par::SodCRParameters_noCRs, sol::SodCRSolution)
    # solves the velocity of the tail of the rarefaction wave
    return par.cl - sol.v34/( 1.0 - par.η2 )
end

function solveV(x::Float64, par::SodCRParameters_noCRs, sol::SodCRSolution)
    # solves the velocity along the shocktube
    if x <= -par.cl * par.t
        return 0.0
    elseif -par.cl * par.t < x <= -sol.vt * par.t
        return solveV2(x, par)
    elseif -sol.vt * par.t < x <= sol.vs * par.t
        return sol.v34
    elseif sol.vs * par.t < x
        return 0.0
    else
        println("Error!")
        return 0.0
    end
end

"""
    Main function for solver
"""

function solveSodShockCR_noPrepopulation(x::Array{Float64,1}; par::SodCRParameters_noCRs)

    # set up datatype to store riemann solution
    sol = SodCRSolution(x)

    # transform into rest-frame of contact discontiuity
    x_in = sol.x .- par.x_contact

    # solve Density
    sol.rho4 = solveRho4(par)
    sol.xs  = sol.rho4/par.rhor

    sol.P34_tot = P4_f(sol.xs, par.Pr, par.γ_th, par.γ_cr, par.ξ)

    sol.rho3 = solveRho3(par, sol)

    sol.xr  = sol.rho3/par.rhol

    # solve pressure
    sol.P4_cr = P_inj_f(sol.xs, par.rhor, par.Pr, par.γ_th, par.γ_cr, par.ξ)

    sol.P4_th = sol.P34_tot - sol.P4_cr
    sol.P3_th = sol.P34_tot

    # solve velocity
    sol.v34 = solveV34(par, sol)

    sol.vs = solveVs(par, sol)
    sol.vt = solveVt(par, sol)

    for i=1:length(sol.x)
        sol.v[i] = solveV(x_in[i], par, sol)
    end

    for i = 1:length(sol.x)
        sol.P_tot[i], sol.P_th[i], sol.P_cr[i], sol.P_cr_p[i], sol.P_cr_e[i] = solveP(x_in[i], par, sol)
    end

    for i = 1:length(sol.x)
        sol.rho[i] = solveRho(x_in[i], par, sol)
    end

    sol.Mach = sol.vs/par.cr

    sol.U_tot = sol.P_tot ./ ( (par.γ_th - 1.0) .* sol.rho )
    sol.U_th = sol.P_th ./ ( (par.γ_th - 1.0) .* sol.rho )
    sol.E_cr_p = sol.P_cr_p ./ ( (par.γ_cr - 1.0) .* sol.rho )
    sol.E_cr_e = sol.P_cr_e ./ ( (par.γ_cr - 1.0) .* sol.rho )

    return sol
end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
export Shamanskii_λ

function Shamanskii_λ(
    nlp_stop,
    PData::PDataFact,
    dₙ::Vector,
    gt::Vector;
    λnull::Bool = false,
)
    nlp_at_x = nlp_stop.current_state
    x = nlp_at_x.x
    n = length(x)
    xdemi = x + dₙ
    gtemp = grad(nlp_stop.pb, xdemi)

    ϵ2 = sqrt(eps(eltype(gt)))
    Γ = max.(abs.(PData.Δ), ϵ2)

    if λnull
        d̃ = PData.L \ gtemp[PData.pp]
        d̂ = PData.L' \ (PData.Q * (PData.Q' * d̃ ./ Γ))
        d = -d̂[invperm(PData.pp)]
    else
        D2 = SymTridiagonal(Hermitian(PData.D + PData.λ * Matrix(1.0I, n, n)))
        Q2 = eigvecs(D2)
        d̃ = PData.L \ gtemp[PData.pp]
        d̂ = PData.L' \ (Q2 * (Q2' * d̃ ./ Γ))
        d = -d̂[invperm(PData.pp)]
    end

    dₕₒ = dₙ .+ d

    return dₕₒ
end


text:
find_llvm() = begin
    if haskey(ENV, "TRAVIS")
        return readchomp(`which llvm-config-3.3`)
    end
    tmp = joinpath(JULIA_HOME, "llvm-config")
    isfile(tmp) && return tmp
    tmp = readchomp(`which llvm-config`)
    isfile(tmp) && return tmp
end
if !haskey(ENV, "LLVM_CONFIG") ENV["LLVM_CONFIG"] = find_llvm() end

# on travis and some other systems, the llvm lib may not be in the
# default path.
push!(Libdl.DL_LOAD_PATH, readchomp(`$(ENV["LLVM_CONFIG"]) --libdir`))

cd(joinpath(Pkg.dir(), "Clang", "deps", "src") )
run(`make`)
if (!ispath("../usr"))
  run(`mkdir ../usr`)
end
if (!ispath("../usr/lib"))
  run(`mkdir ../usr/lib`)
end

run(`mv libwrapclang.$(Libdl.dlext) ../usr/lib`)


text:
using Profile
import Base.gc_num
export summarize_events

const Timestamp = UInt64

struct ProfilerResult
    samples::Vector{UInt64}
    lineinfo::AbstractDict
end

"""
identifies

space (category, id)
time (timeline, start, finish)

also tracks gc_num during this and profiling samples.
"""
struct Timespan
    category::Symbol
    id::Any
    timeline::Any
    start::Timestamp
    finish::Timestamp
    gc_diff::Base.GC_Diff
    profiler_samples::ProfilerResult
end

struct Event{phase}
    category::Symbol
    id::Any
    timeline::Any
    timestamp::Timestamp
    gc_num::Base.GC_Num
    profiler_samples::ProfilerResult
end

Event(phase::Symbol, category::Symbol,
      id, tl, time, gc_num, prof) =
    Event{phase}(category, id, tl, time, gc_num, prof)

"""
create a timespan given the strt and finish events
"""
function make_timespan(start::Event, finish::Event)
    @assert start.category == finish.category
    @assert start.id == finish.id
    @assert start.timeline == finish.timeline

    Timespan(start.category,
             start.id,
             start.timeline,
             start.timestamp,
             finish.timestamp,
             Base.GC_Diff(finish.gc_num,start.gc_num),
             mix_samples(start.profiler_samples, finish.profiler_samples))
end

"""
Various means of writing an event to something.
"""
struct NoOpLog end

function write_event(::NoOpLog, event::Event)
end

struct FilterLog
    f::Function
    inner_chan::Any
end

function write_event(c::FilterLog, event)
    if c.f(event)
        write_event(c.inner_chan, event)
    end
end

function write_event(io::IO, event::Event)
    serialize(io, event)
end

function write_event(chan::Union{RemoteChannel, Channel}, event::Event)
    put!(chan, event)
end

function write_event(arr::AbstractArray, event::Event)
    push!(arr, event)
end

#function write_event(sig::Signal, event::Event)
#    push!(sig, event)
#end
"""
represents a process local events array.

A context with log_sink set to LocalEventLog() will
cause events to be recorded into the
"""
struct LocalEventLog end

const _local_event_log = Any[]
clear_local_event_log() = empty!(_local_event_log)
function write_event(::LocalEventLog, event::Event)
    write_event(Dagger._local_event_log, event)
end

function raise_event(ctx, phase, category, id,tl, t, gc_num, prof, async)
    ev = Event(phase, category, id, tl, t, gc_num, prof)
    if async
        @async write_event(ctx, ev)
    else
        write_event(ctx, ev)
    end
end

empty_prof() = ProfilerResult(UInt[], Profile.getdict(UInt[]))

function timespan_start(ctx, category, id, tl, async=isasync(ctx.log_sink))
    isa(ctx.log_sink, NoOpLog) && return # don't go till raise
    if ctx.profile
        Profile.start_timer()
    end
    raise_event(ctx, :start, category, id, tl, time_ns(), gc_num(), empty_prof(), async)
    nothing
end

function timespan_end(ctx, category, id, tl, async=isasync(ctx.log_sink))
    isa(ctx.log_sink, NoOpLog) && return
    prof = UInt[]
    if ctx.profile
        Profile.stop_timer()
        prof = Profile.fetch()
        Profile.clear()
    end
    raise_event(ctx, :finish, category, id, tl,time_ns(), gc_num(), ProfilerResult(prof, Profile.getdict(prof)), async)
    nothing
end

isasync(x) = false
isasync(x::Union{Channel, RemoteChannel, IO}) = true
"""
Overall state used during visualization
"""
mutable struct State
    start_events::Dict # (category, id) => Event
    finish_events::Dict  # (category, id) => Event
    #completed::Dict      # timeline => category => Array
    completed::Vector
    start_time::Timestamp
    finish_time::Timestamp
end
State() = State(Dict(), Dict(), Any[], 0, 0)

"""
Add a Timespan to a given State under `tl` (timeline)
and `category`.
"""
function add_span(state, tl, category, span)
    push!(state.completed, span)
    if state.start_time == 0
        state.start_time = span.start
    else
        state.start_time = min(span.start, state.start_time)
    end
    if state.finish_time == 0
        state.finish_time = span.finish
    else
        state.finish_time = max(span.finish, state.finish_time)
    end
    state
end

"""When building state for real-time visualization,
   use next_state to progress gantt state."""
function next_state(state::State, event::Event{:start})
    key = (event.category, event.id)
    if haskey(state.finish_events, key) # finish event reached before start
        span = make_timespan(event, pop!(state.finish_events, key))
        add_span(state, event.timeline, event.category, span)
    else
        state.start_events[key] = event
    end
    state
end

function next_state(state::State, event::Event{:finish})
    key = (event.category, event.id)
    if haskey(state.start_events, key)
        span = make_timespan(pop!(state.start_events, key), event)
        add_span(state, event.timeline, event.category, span)
    else
        state.finish_events[key] = event
    end
    state
end
next_state(state::State, events::AbstractArray) =
    foldl(next_state, events, init=state)

# util

function pushkey(dict, key, thing)
    if haskey(dict, key)
        push!(dict[key],thing)
    else
        dict[key] = Any[thing]
    end
end

function pushkey(dict, key1, args...)
    if haskey(dict, key1)
        pushkey(dict[key1], args...)
    else
        dict[key1] = Dict()
        pushkey(dict[key1], args...)
    end
end

function mix_samples(a,b)
    ProfilerResult(vcat(a.samples, b.samples),
                   merge(a.lineinfo, b.lineinfo))
end

function build_timespans(events)
    next_state(State(), events)
end

macro logmsg(ex)
    #:(println($(esc(ex))))
end

"""
Get the logs from each process, clear it too
"""
function get_logs!(::LocalEventLog)
    logs = Dict()
    @sync for p in procs()
        @async logs[p] = remotecall_fetch(p) do
            log = copy(Dagger._local_event_log)
            clear_local_event_log()
            log
        end
    end
    build_timespans(vcat(values(logs)...)).completed
end

function add_gc_diff(x,y)
    Base.GC_Diff(
        x.allocd     + y.allocd,
        x.malloc     + y.malloc,
        x.realloc    + y.realloc,
        x.poolalloc  + y.poolalloc,
        x.bigalloc   + y.bigalloc,
        x.freecall   + y.freecall,
        x.total_time + y.total_time,
        x.pause      + y.pause,
        x.full_sweep + y.full_sweep
    )
end

function aggregate_events(xs)
    gc_diff = reduce(add_gc_diff, map(x -> x.gc_diff, xs))
    time_spent = sum(map(x -> x.finish - x.start, xs))
    profiler_samples = treereduce(mix_samples, map(x->x.profiler_samples, xs))
    time_spent, gc_diff, profiler_samples
end

function summarize_events(time_spent, gc_diff, profiler_samples)
    Base.time_print(time_spent, gc_diff.allocd, gc_diff.total_time, Base.gc_alloc_count(gc_diff))
    if !isempty(profiler_samples.samples)
        Profile.print(profiler_samples.samples, profiler_samples.lineinfo)
    end
end

summarize_events(xs) = summarize_events(aggregate_events(xs)...)


text:
module TestIndex
using Compat, Compat.Test, DataFrames
using DataFrames: Index

i = Index()
push!(i, :A)
push!(i, :B)

inds = Any[1,
           big(1),
           :A,
           [true, false],
           [1],
           [big(1)],
           big(1):big(1),
           [:A],
           Union{Bool, Missing}[true, false],
           Union{Int, Missing}[1],
           Union{BigInt, Missing}[big(1)],
           Union{Symbol, Missing}[:A],
           Any[1],
           Any[:A]]

for ind in inds
    if ind == :A || ndims(ind) == 0
        @test i[ind] == 1
    else
        @test (i[ind] == [1])
    end
end

# TODO: change this to throw an error after deprecation
# @test_throws MethodError i[1.0]
# @test_throws MethodError i[true]
# @test_throws MethodError i[false]
# @test_throws ArgumentError i[Any[1, missing]]
# @test_throws ArgumentError i[[1, missing]]
# @test_throws ArgumentError i[[true, missing]]
# @test_throws ArgumentError i[Any[true, missing]]
# @test_throws ArgumentError i[[:A, missing]]
# @test_throws ArgumentError i[Any[:A, missing]]
# @test_throws ArgumentError i[1.0:1.0]
# @test_throws ArgumentError i[[1.0]]
# @test_throws ArgumentError i[Any[1.0]]

@test i[1:1] == 1:1

@test_throws BoundsError i[[true]]
@test_throws BoundsError i[[true, false, true]]

@test_throws ArgumentError i[["a"]]
@test_throws ArgumentError i[Any["a"]]

@test i[[]] == Int[]
@test i[Int[]] == Int[]
@test i[Symbol[]] == Int[]

@test names(i) == [:A,:B]
@test names!(i, [:a,:a], makeunique=true) == Index([:a,:a_1])
@test_throws ArgumentError names!(i, [:a,:a])
@test names!(i, [:a,:b]) == Index([:a,:b])
@test rename(i, Dict(:a=>:A, :b=>:B)) == Index([:A,:B])
@test rename(i, :a => :A) == Index([:A,:b])
@test rename(i, :a => :a) == Index([:a,:b])
@test rename(i, [:a => :A]) == Index([:A,:b])
@test rename(i, [:a => :a]) == Index([:a,:b])
@test rename(x->Symbol(uppercase(string(x))), i) == Index([:A,:B])
@test rename(x->Symbol(lowercase(string(x))), i) == Index([:a,:b])

@test delete!(i, :a) == Index([:b])
push!(i, :C)
@test delete!(i, 1) == Index([:C])

i = Index([:A, :B, :C, :D, :E])
i2 = copy(i)
names!(i2, reverse(names(i2)))
names!(i2, reverse(names(i2)))
@test names(i2) == names(i)
for name in names(i)
  i2[name] # Issue #715
end

#= Aliasing & Mutation =#

# columns should not alias if scalar broadcasted
df = DataFrame(A=[0],B=[0])
df[1:end] = 0.0
df[1,:A] = 1.0
@test df[1,:B] === 0

end


text:
#= IM 3 sample that shows an image.

  Needs "im.lib", "iup.lib", "cd.lib" and "iupcd.lib".

  Usage: im_view <file_name>

    Example: im_view test.tif

  Click on image to open another file.

#include <iup.h>
#include <cd.h>
#include <cdiup.h>
#include <im.h>
#include <im_image.h>
=#

module im_view_

using Printf
using IUP

export
  im_view

# forward declaration to workaround '@cfunction'
function cbCanvasButton end
function cbCanvasRepaint end
function cbCanvasMap end
function cbDialogClose end

global disable_repaint = false      # used to optimize repaint, while opening a new file

# --------------------------------------------------------------------------------
function im_view()

    IupOpen()         #Initializes IUP
    dlg = CreateDialog()

     mainMenu = CreateMainMenu()
    # Associates handle "mainMenu" with mainMenu
    IupSetHandle("mainMenu", mainMenu)
    IupSetAttribute(dlg, "MENU","mainMenu")

    imlabCreateButtonImages()

    btOpen = IupButton("", "imImgWinOpen" )
    IupSetAttribute(btOpen,"TIP","Loads an image file from disk.")
    IupSetAttribute(btOpen,"IMAGE","imImgWinOpenButton")
    hbToolBar = IupHbox(
        IupSetAttributes(IupFill(),"SIZE=5, EXPAND=NO"),    # Blank space before first icon
        btOpen,
        IupSetAttributes(IupFill(),"SIZE=5, EXPAND=NO"),
        IupSetAttributes(IupFill(),"SIZE=200, EXPAND=NO"),
        IupSetAttributes(IupFill(),"EXPAND=YES"),
    )
    IupSetAttribute(hbToolBar, "FLAT", "YES");

    IupShow(dlg);

    file_name = "*.*"
    if (IupGetFile(file_name) == 0)
        ShowImage(file_name, dlg)
    end

    IupMainLoop()
    IupDestroy(dlg)
    IupClose()
    return 0
end

# --------------------------------------------------------------------------------
function PrintError(error::Integer)
    if (error == IM_ERR_OPEN)
        IupMessage("IM", "Error Opening File.")
    elseif (error == IM_ERR_MEM)
        IupMessage("IM", "Insuficient memory.")
    elseif (error == IM_ERR_ACCESS)
        IupMessage("IM", "Error Accessing File.")
    elseif (error == IM_ERR_DATA)
        IupMessage("IM", "Image type not Suported.")
    elseif (error == IM_ERR_FORMAT)
        IupMessage("IM", "Invalid Format.")
    elseif (error == IM_ERR_COMPRESS)
        IupMessage("IM", "Invalid or unsupported compression.")
    else
        IupMessage("IM", "Unknown Error.")
    end
end

# --------------------------------------------------------------------------------
function cbCanvasRepaint(iup_canvas::Ptr{Ihandle})
    global disable_repaint

    cd_canvas = convert(Ptr{cdCanvas}, IupGetAttribute(iup_canvas, "cdCanvas"))

    if (cd_canvas == C_NULL || disable_repaint)
        return IUP_DEFAULT
    end

    cdCanvasActivate(cd_canvas);
    cdCanvasClear(cd_canvas);

    image = convert(Ptr{imImage}, IupGetAttribute(iup_canvas, "imImage"))
    if (image == C_NULL)
        return IUP_DEFAULT
    end

    img = unsafe_load(image)        # Need to get access to the composite type, not to its pointer
    imcdCanvasPutImage(cd_canvas, img, 0, 0, img.width, img.height, 0, 0, 0, 0)
#    imcdCanvasPutImage(cd_canvas, img, 20,20,int64(img.width/2), int64(img.height/2),  0, 0, 0, 0)

    #=data = [convert(Ptr{UInt8}, unsafe_load(img.data,1)),
            convert(Ptr{UInt8}, unsafe_load(img.data,2)),
            convert(Ptr{UInt8}, unsafe_load(img.data,3))]
    cdCanvasPutImageRectRGB(cd_canvas, img.width, img.height,
                    data[1], data[2], data[3],
                    0, 0, img.width, img.height, 0, 0, 0, 0)=#

    cdCanvasFlush(cd_canvas)
    return IUP_DEFAULT
end

# --------------------------------------------------------------------------------
function ShowImage(file_name::String, iup_dialog::Ptr{Ihandle})
    image = IupGetAttribute(iup_dialog, "imImage")        # typeof(image) => Ptr{UInt8}
    image = convert(Ptr{imImage}, image)        # If I use Ptr{imImage} it Booms???
    if (image != C_NULL)
        imImageDestroy(image)
    end
    IupSetAttribute(iup_dialog, "imImage")

    error = pointer([Cint(0)])
    image = imFileImageLoadBitmap(file_name, Cint(0), error)
    error = unsafe_load(error)
    if (error != 0)    PrintError(error)    end
    if (image == C_NULL) return end

    IupSetAttribute(iup_dialog, "imImage", image);
    IupStoreAttribute(iup_dialog, "TITLE", file_name);

    cbCanvasRepaint(iup_dialog)   # we can do this because canvas inherit attributes from the dialog */
end

# --------------------------------------------------------------------------------
function cbCanvasButton(iup_canvas::Ptr{Ihandle}, but::Int, pressed::Int)
    global disable_repaint

    file_name = "*.*"
    if (but != IUP_BUTTON1 || pressed != 0)
        return IUP_DEFAULT
    end

    disable_repaint = true
    if (IupGetFile(file_name) != 0)
        disable_repaint = false
        return IUP_DEFAULT
    end

    disable_repaint = false
    ShowImage(file_name, IupGetDialog(iup_canvas))
    return IUP_DEFAULT
end

# --------------------------------------------------------------------------------
function cbCanvasMap(iup_canvas::Ptr{Ihandle})
    cd_canvas = cdCreateCanvas(cdContextIup(), iup_canvas)
    IupSetAttribute(IupGetDialog(iup_canvas), "cdCanvas", cd_canvas)
    return IUP_DEFAULT
end

# --------------------------------------------------------------------------------
function cbDialogClose(iup_dialog::Ptr{Ihandle})
    cd_canvas = convert(Ptr{cdCanvas}, IupGetAttribute(iup_dialog, "cdCanvas"))
    image = convert(Ptr{imImage}, IupGetAttribute(iup_dialog, "imImage"))

    if (cd_canvas != C_NULL) cdKillCanvas(cd_canvas)    end
    if (image != C_NULL) imImageDestroy(image)        end

    IupSetAttribute(iup_dialog, "cdCanvas")
    IupSetAttribute(iup_dialog, "imImage")

    return IUP_CLOSE
end

# --------------------------------------------------------------------------------
function CreateDialog()

    iup_canvas = IupCanvas()
    IupSetCallback(iup_canvas, "BUTTON_CB", @cfunction(cbCanvasButton, Int, (Ptr{Ihandle}, Int, Int)))
    IupSetCallback(iup_canvas, "ACTION", @cfunction(cbCanvasRepaint, Int, (Ptr{Ihandle},)))
    IupSetCallback(iup_canvas, "MAP_CB", @cfunction(cbCanvasMap, Int, (Ptr{Ihandle},)))

    iup_dialog = IupDialog(iup_canvas)
    IupSetCallback(iup_dialog, "CLOSE_CB", @cfunction(cbDialogClose, Int, (Ptr{Ihandle},)))
    IupSetAttribute(iup_dialog, "SIZE", "HALFxHALF")      # initial size

    return iup_dialog
end

# ---------------------------------------------------------------------------------------------
function CreateMainMenu()
    item_help = IupItem("Help", "item_help_act")
    menu_help  = IupMenu(item_help);
    submenu_help  = IupSubmenu("Help", menu_help)

    # Creates main menu with file menu
    mainMenu = IupMenu(submenu_help)

    return mainMenu

end

# ----------------------------------------------------------------------------------------------
function imlabCreateButtonImages()
    new_bits = UInt8[
        2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2
        ,3,2,2,1,3,2,2,1,2,2,2,2,2,2,2,2
        ,1,3,2,1,4,2,1,3,2,2,2,2,2,2,2,2
        ,2,1,3,1,3,1,3,0,0,0,0,0,0,0,0,2
        ,2,3,1,3,4,1,1,1,1,4,4,4,4,4,0,2
        ,1,1,1,4,3,3,4,4,4,4,4,4,4,4,0,2
        ,2,2,4,1,4,4,1,4,4,4,4,4,4,4,0,2
        ,1,2,2,1,3,4,4,4,4,4,4,4,4,4,0,2
        ,2,2,2,1,4,4,4,4,4,4,4,4,4,4,0,2
        ,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
        ,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
        ,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
        ,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
        ,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
        ,2,2,2,0,4,4,4,4,4,4,4,4,4,4,0,2
        ,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,2
    ]

    new_colors = [
        "0 0 0",
        "132 132 132",
        "BGCOLOR",
        "255 255 0",
        "255 255 255",
        ""
    ]

    open_bits = UInt8[
        1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
        1,1,1,1,1,1,1,1,1,1,0,0,0,1,1,1,
        1,1,1,1,1,1,1,1,1,0,1,1,1,0,1,0,
        1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,
        1,1,0,0,0,1,1,1,1,1,1,1,1,0,0,0,
        1,0,1,3,1,0,0,0,0,0,0,0,1,1,1,1,
        1,0,3,1,3,1,3,1,3,1,3,0,1,1,1,1,
        1,0,1,3,1,3,1,3,1,3,1,0,1,1,1,1,
        1,0,3,1,3,1,0,0,0,0,0,0,0,0,0,0,
        1,0,1,3,1,0,2,2,2,2,2,2,2,2,2,0,
        1,0,3,1,0,2,2,2,2,2,2,2,2,2,0,1,
        1,0,1,0,2,2,2,2,2,2,2,2,2,0,1,1,
        1,0,0,2,2,2,2,2,2,2,2,2,0,1,1,1,
        1,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,
        1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
        1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    ]

    open_colors = [
        "0 0 0",
        "BGCOLOR",
        "128 128 0",
        "255 255 0",
        ""
    ]

    CreateButtonImage(16,16, new_bits, new_colors, "imImgWinNewButton");
    CreateButtonImage(16,16, open_bits, open_colors, "imImgWinOpenButton");
end

# ----------------------------------------------------------------------------------------------
function CreateButtonImage(w::Int, h::Int, bits, colors, name::String)
    global total_images, images

    color_str = ["0","1","2","3","4","5","6","7","8","9","10","11","12","13","14","15"]
    i = 1

    iup_image = IupImage(w, h, bits)
    IupSetHandle(name, iup_image)

    while (colors[i] != "")
        if (i > 16)
            aux_color_str = @sprintf("%d", i)
            IupStoreAttribute(iup_image, aux_color_str, colors[i]);
        else
            IupStoreAttribute(iup_image, color_str[i], colors[i]);
        end
        i += 1
    end

end

end   # module


text:
## Here we show how to create the moments data bunnytime.mnts.jld and bunnytimetest.mnts, from the full simulated data (which is kept in a Zenodo repo).

Pkg.clone("https://github.com/jondea/MultipleScattering.jl.git")
using MultipleScattering

Pkg.add("JLD")
using JLD

model_to_label(f) = (f.volfrac,f.particles[1].r)
impulse = get_gaussian_freq_impulse(1.0,2.48)

# Download full simulated data. Note, can take a long time.
  #training set
  download("https://zenodo.org/record/1126642/files/bunny.mst.jld", "bunny.mst.jld")
  #test set
  download("https://zenodo.org/record/1126642/files/bunnytest.mst.jld", "bunnytest.mst.jld")

# There is also a copy of moments data on Zenodo, which can de download using:
  # # Download data from permenant repo
  # download("https://zenodo.org/record/1126642/files/train.mnts.jld", "train.mnts.jld")
  # download("https://zenodo.org/record/1126642/files/test.mnts.jld", "test.mnts.jld")

# test set: create the moments of the backscattered waves in time
  test = collect(values(load("bunnytest.mst.jld")))[1]
  test = sort(test, by = model_to_label)
  test_batched = collect(groupby(model_to_label, test));

  tmp_labels = [model_to_label(bs[1]) for bs in test_batched];
  if length(union(tmp_labels)) != length(tmp_labels)
    error("batched the test data inccorectly")
  end

  test_moments = map(test_batched) do bs
    time_arr = cut_interval(ω_to_t(bs[1].k_arr),[10.,96.])
    StatisticalMoments(
        [TimeSimulation(b;time_arr = time_arr, impulse = impulse) for b in bs],
        4;
        response_apply=real
    )
  end;
  save("train.mnts.jld", "Array{StatisticalMoments}", test_moments)

  test_batched = 0
  test_moments = 0
  gc() # garbage collection

# training set: create the moments of the backscattered waves in time
  train = collect(values(load("bunny.mst.jld")))[1]
  train = sort(train, by = model_to_label)
  train_batched = collect(groupby(model_to_label, train));
  train = 0
  gc()

  tmp_labels = [model_to_label(bs[1]) for bs in train_batched];
  if length(union(tmp_labels)) != length(tmp_labels)
    error("batched the test data inccorectly")
  end

  train_moments = map(train_batched) do bs
    time_arr = cut_interval(ω_to_t(bs[1].k_arr),[10.,96.])
    StatisticalMoments(
        [TimeSimulation(b;time_arr = time_arr, impulse = impulse) for b in bs],
        4;
        response_apply=real
    )
  end;
  save("test.mnts.jld", "Array{StatisticalMoments}", train_moments)


text:
using OrdinaryDiffEq, DiffEqBase, Base.Test
srand(100)

linear = (t,u) -> (1.01*u)
analytic_linear = (t,u0) -> u0*exp(1.01*t)
prob = ODETestProblem(linear,1/2,analytic_linear,(0.0,1.0))

sol =solve(prob,Tsit5(),dt=1//2^(6),tstops=[1/2])

@test 1//2 ∈ sol.t

sol =solve(prob,RK4(),dt=1//3,tstops=[1/2])

@test sol.t == [0,1/3,1/2,1/3+1/2,1]

sol =solve(prob,RK4(),tstops=[1/5,1/4,1/3,1/2,3/4])

@test sol.t == [0,1/5,1/4,1/3,1/2,3/4,1]

sol =solve(prob,RK4(),tstops=[0,1/5,1/4,1/3,1/2,3/4,1])

@test sol.t == [0,1/5,1/4,1/3,1/2,3/4,1]

sol = solve(prob,RK4(),tstops=0:1//16:1)

@test sol.t == collect(0:1//16:1)

sol = solve(prob,RK4(),tstops=linspace(0,1,100))

@test sol.t == collect(linspace(0,1,100))


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("../module/CustomStencil.jl")
include("../misc/misc.jl")
include("../misc/cpu_stencils.jl")

## Stencil Definition
coefs = Float64.([1;0.5;0.25])
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 2, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst = NewStencilInstance(st_def, m_step=false)

## Input Data size Definition
radius = 2
nx = 6
ny = 6
nz = 6
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 4
## GPU Implementation
data = CreateData(dx,dy,dz)
gpu_out = ApplyStencil(st_inst, data, t_steps)

CUDA.cuProfilerStart()
d2 = ApplyMultiGPU(2, st_inst, t_steps, data)


text:
export linear_svm

import ..LinearSVMModel

using mlpack._Internal.io

import mlpack_jll
const linear_svmLibrary = mlpack_jll.libmlpack_julia_linear_svm

# Call the C binding of the mlpack linear_svm binding.
function linear_svm_mlpackMain()
  success = ccall((:linear_svm, linear_svmLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module linear_svm_internal
  import ..linear_svmLibrary

import ...LinearSVMModel

# Get the value of a model pointer parameter of type LinearSVMModel.
function IOGetParamLinearSVMModel(paramName::String)::LinearSVMModel
  LinearSVMModel(ccall((:IO_GetParamLinearSVMModelPtr, linear_svmLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type LinearSVMModel.
function IOSetParamLinearSVMModel(paramName::String, model::LinearSVMModel)
  ccall((:IO_SetParamLinearSVMModelPtr, linear_svmLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeLinearSVMModel(stream::IO, model::LinearSVMModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeLinearSVMModelPtr, linear_svmLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeLinearSVMModel(stream::IO)::LinearSVMModel
  buffer = read(stream)
  LinearSVMModel(ccall((:DeserializeLinearSVMModelPtr, linear_svmLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    linear_svm(; [delta, epochs, input_model, labels, lambda, max_iterations, no_intercept, num_classes, optimizer, seed, shuffle, step_size, test, test_labels, tolerance, training, verbose])

An implementation of linear SVMs that uses either L-BFGS or parallel SGD
(stochastic gradient descent) to train the model.

This program allows loading a linear SVM model (via the `input_model` parameter)
or training a linear SVM model given training data (specified with the
`training` parameter), or both those things at once.  In addition, this program
allows classification on a test dataset (specified with the `test` parameter)
and the classification results may be saved with the `predictions` output
parameter. The trained linear SVM model may be saved using the `output_model`
output parameter.

The training data, if specified, may have class labels as its last dimension.
Alternately, the `labels` parameter may be used to specify a separate vector of
labels.

When a model is being trained, there are many options.  L2 regularization (to
prevent overfitting) can be specified with the `lambda` option, and the number
of classes can be manually specified with the `num_classes`and if an intercept
term is not desired in the model, the `no_intercept` parameter can be
specified.Margin of difference between correct class and other classes can be
specified with the `delta` option.The optimizer used to train the model can be
specified with the `optimizer` parameter.  Available options are 'psgd'
(parallel stochastic gradient descent) and 'lbfgs' (the L-BFGS optimizer).
There are also various parameters for the optimizer; the `max_iterations`
parameter specifies the maximum number of allowed iterations, and the
`tolerance` parameter specifies the tolerance for convergence.  For the parallel
SGD optimizer, the `step_size` parameter controls the step size taken at each
iteration by the optimizer and the maximum number of epochs (specified with
`epochs`). If the objective function for your data is oscillating between Inf
and 0, the step size is probably too large.  There are more parameters for the
optimizers, but the C++ interface must be used to access these.

Optionally, the model can be used to predict the labels for another matrix of
data points, if `test` is specified.  The `test` parameter can be specified
without the `training` parameter, so long as an existing linear SVM model is
given with the `input_model` parameter.  The output predictions from the linear
SVM model may be saved with the `predictions` parameter.

As an example, to train a LinaerSVM on the data '`data`' with labels '`labels`'
with L2 regularization of 0.1, saving the model to '`lsvm_model`', the following
command may be used:

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> labels = CSV.read("labels.csv"; type=Int)
julia> lsvm_model, _, _ = linear_svm(delta=1, labels=labels,
            lambda=0.1, num_classes=0, training=data)
```

Then, to use that model to predict classes for the dataset '`test`', storing the
output predictions in '`predictions`', the following command may be used:

```julia
julia> using CSV
julia> test = CSV.read("test.csv")
julia> _, predictions, _ = linear_svm(input_model=lsvm_model,
            test=test)
```

# Arguments

 - `delta::Float64`: Margin of difference between correct class and other
      classes.  Default value `1`.

 - `epochs::Int`: Maximum number of full epochs over dataset for psgd
      Default value `50`.

 - `input_model::LinearSVMModel`: Existing model (parameters).
 - `labels::Array{Int, 1}`: A matrix containing labels (0 or 1) for the
      points in the training set (y).
 - `lambda::Float64`: L2-regularization parameter for training.  Default
      value `0.0001`.

 - `max_iterations::Int`: Maximum iterations for optimizer (0 indicates no
      limit).  Default value `10000`.

 - `no_intercept::Bool`: Do not add the intercept term to the model.
      Default value `false`.

 - `num_classes::Int`: Number of classes for classification; if
      unspecified (or 0), the number of classes found in the labels will be
      used.  Default value `0`.

 - `optimizer::String`: Optimizer to use for training ('lbfgs' or 'psgd').
       Default value `lbfgs`.

 - `seed::Int`: Random seed.  If 0, 'std::time(NULL)' is used.  Default
      value `0`.

 - `shuffle::Bool`: Don't shuffle the order in which data points are
      visited for parallel SGD.  Default value `false`.

 - `step_size::Float64`: Step size for parallel SGD optimizer.  Default
      value `0.01`.

 - `test::Array{Float64, 2}`: Matrix containing test dataset.
 - `test_labels::Array{Int, 1}`: Matrix containing test labels.
 - `tolerance::Float64`: Convergence tolerance for optimizer.  Default
      value `1e-10`.

 - `training::Array{Float64, 2}`: A matrix containing the training set
      (the matrix of predictors, X).
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output_model::LinearSVMModel`: Output for trained linear svm model.
 - `predictions::Array{Int, 1}`: If test data is specified, this matrix is
      where the predictions for the test set will be saved.
 - `probabilities::Array{Float64, 2}`: If test data is specified, this
      matrix is where the class probabilities for the test set will be saved.

"""
function linear_svm(;
                    delta::Union{Float64, Missing} = missing,
                    epochs::Union{Int, Missing} = missing,
                    input_model::Union{LinearSVMModel, Missing} = missing,
                    labels = missing,
                    lambda::Union{Float64, Missing} = missing,
                    max_iterations::Union{Int, Missing} = missing,
                    no_intercept::Union{Bool, Missing} = missing,
                    num_classes::Union{Int, Missing} = missing,
                    optimizer::Union{String, Missing} = missing,
                    seed::Union{Int, Missing} = missing,
                    shuffle::Union{Bool, Missing} = missing,
                    step_size::Union{Float64, Missing} = missing,
                    test = missing,
                    test_labels = missing,
                    tolerance::Union{Float64, Missing} = missing,
                    training = missing,
                    verbose::Union{Bool, Missing} = missing,
                    points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, linear_svmLibrary), Nothing, ());

  IORestoreSettings("Linear SVM is an L2-regularized support vector machine.")

  # Process each input argument before calling mlpackMain().
  if !ismissing(delta)
    IOSetParam("delta", convert(Float64, delta))
  end
  if !ismissing(epochs)
    IOSetParam("epochs", convert(Int, epochs))
  end
  if !ismissing(input_model)
    linear_svm_internal.IOSetParamLinearSVMModel("input_model", convert(LinearSVMModel, input_model))
  end
  if !ismissing(labels)
    IOSetParamURow("labels", labels)
  end
  if !ismissing(lambda)
    IOSetParam("lambda", convert(Float64, lambda))
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(no_intercept)
    IOSetParam("no_intercept", convert(Bool, no_intercept))
  end
  if !ismissing(num_classes)
    IOSetParam("num_classes", convert(Int, num_classes))
  end
  if !ismissing(optimizer)
    IOSetParam("optimizer", convert(String, optimizer))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(shuffle)
    IOSetParam("shuffle", convert(Bool, shuffle))
  end
  if !ismissing(step_size)
    IOSetParam("step_size", convert(Float64, step_size))
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(test_labels)
    IOSetParamURow("test_labels", test_labels)
  end
  if !ismissing(tolerance)
    IOSetParam("tolerance", convert(Float64, tolerance))
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output_model")
  IOSetPassed("predictions")
  IOSetPassed("probabilities")
  # Call the program.
  linear_svm_mlpackMain()

  return linear_svm_internal.IOGetParamLinearSVMModel("output_model"),
         IOGetParamURow("predictions"),
         IOGetParamMat("probabilities", points_are_rows)
end


text:
function df_to_dict_with_id(df::DataFrame)
    Dict(col[1] => Dict(zip(df[!, 1], col[2])) for col in eachcol(df, true))
end

function df_to_dict(df::DataFrame)
    Dict(col[1] => col[2] for col in eachcol(df, true))
end

function dictzip(df::DataFrame, which::Tuple{Symbol, Symbol})
    return Dict(zip(df[which[1]], df[which[2]]))
end

function insert_default!(dict::Dict, default, keys)
	for k in keys
		haskey(dict, k) || (dict[k] = default)
	end
end

function has(x, y)
	x == y
end

function getKeyVector(df::Dict, value)
	key_value = [k for (k,v) in df if has(v, value)]
	return key_value
end