text:
using Base

export saveImage, append!, poseIndex, odoDiff
export RoverPose

struct RoverImage
    timestamp::Float64
    camJpeg::AbstractString
end

mutable struct RoverPose
    poseId::String
    poseIndex::Int64
    timestamp::Float64
    x::Float64
    y::Float64
    theta::Float64
    camImages::Vector{RoverImage}
    prevPose::Nullable{RoverPose}
    # Other stuff like AprilTags here

    # RoverPose(timestamp::Float64, locX::Float64, locY::Float64, theta::Float64, camJpeg::AbstractString) = new(timestamp, locX, locY, theta, camJpeg)
    RoverPose(prevPose::RoverPose) = new(string(Base.Random.uuid4()), prevPose.poseIndex+1, prevPose.timestamp, prevPose.x, prevPose.y, prevPose.theta, Vector{RoverImage}(), prevPose);
    RoverPose() = new(string(Base.Random.uuid4()), 1, 0, 0, 0, 0, Vector{RoverImage}(0), Nullable{RoverPose}());
end

Base.show(io::IO, roverPose::RoverPose) = @printf "[%0.3f] (%s) %s: At (%0.2f, %0.2f) with heading %0.2f, delta %s, and %d images" roverPose.timestamp roverPose.poseId string(poseIndex(roverPose)) roverPose.x roverPose.y roverPose.theta string(odoDiff(roverPose)) length(roverPose.camImages);

"""
Save an image to a file.
"""
function saveImage(roverImage :: RoverImage, imFile :: AbstractString)
    out = open(imFile,"w")
    write(out,roverImage.camJpeg)
    close(out)
end

"""
Returns the pose index symbol, e.g. x1.
"""
function poseIndex(roverPose::RoverPose)
    return Symbol("x"*string(roverPose.poseIndex))
end

"""
Returns the diff (dx, dy, dtheta) between roverPose and roverPose.prevPose.
If no prevPose, assumes origin is (0,0,0).
"""
function odoDiff(roverPose::RoverPose)
    if isnull(roverPose.prevPose)
        return [roverPose.x, roverPose.y, roverPose.theta]
    end
    prevPose = get(roverPose.prevPose)

    return [roverPose.x - prevPose.x, roverPose.y - prevPose.y, roverPose.theta - prevPose.theta]
end

"""
Append the state to this pose, increasing the movement by the given update and adding in the sensory data.
"""
function append!(roverPose::RoverPose, newRoverState::PyCall.PyObject, maxImages = 100)
    duration = newRoverState[:getDuration]()
    endTime = newRoverState[:getEndTime]()
    jpegBytes = newRoverState[:getImage]()
    push!(roverPose.camImages, RoverImage(endTime, jpegBytes))
    while(length(roverPose.camImages) > maxImages)
        pop!(roverPose.camImages)
    end
    # Timestamp
    roverPose.timestamp = endTime
    # Should be either rotation or translation but not both at same time
    roverPose.theta += duration * newRoverState[:getRotSpeedNorm]()
    roverPose.x -= duration * sin(roverPose.theta) * newRoverState[:getTransSpeedNorm]()
    roverPose.y += duration * cos(roverPose.theta) * newRoverState[:getTransSpeedNorm]()
    #println("Duration = $(duration), updated position = $(roverPose.x), $(roverPose.y)")
end

function isPoseWorthy(roverPose::RoverPose)
    # If no previous pose, always true
    if isnull(roverPose.prevPose)
        return true
    end
    prevPose = get(roverPose.prevPose)

    deltaT = roverPose.timestamp - prevPose.timestamp
    deltaDist = sqrt((roverPose.x - prevPose.x)^2 + (roverPose.y - prevPose.y)^2)
    deltaAbsAngRad = abs(roverPose.theta - prevPose.theta)
    while deltaAbsAngRad < 0
        deltaAbsAngRad += 2.0*pi
    end
    while deltaAbsAngRad > 2.0*pi
        deltaAbsAngRad -= 2.0*pi
    end
    return deltaT > 20 || deltaDist > 0.5 || deltaAbsAngRad > pi / 8
end

println("Thanks for importing RoverPose :)")


text:
"""
    ModelCallbacks

The function callbacks for the dynamics and constraints.
"""
mutable struct ModelCallbacks

    ffct_cb::Callback
    dfdx_vec_cb::Callback
    dfdu_vec_cb::Callback
    dfdp_vec_cb::Callback

    lfct_cb::Callback
    dldx_cb::Callback
    dldu_cb::Callback
    dldp_cb::Callback

    Vfct_cb::Callback
    dVdx_cb::Callback
    dVdp_cb::Callback
    dVdT_cb::Callback

    gfct_cb::Callback
    dgdx_vec_cb::Callback
    dgdu_vec_cb::Callback
    dgdp_vec_cb::Callback

    hfct_cb::Callback
    dhdx_vec_cb::Callback
    dhdu_vec_cb::Callback
    dhdp_vec_cb::Callback

    gTfct_cb::Callback
    dgTdx_vec_cb::Callback
    dgTdp_vec_cb::Callback
    dgTdT_vec_cb::Callback

    hTfct_cb::Callback
    dhTdx_vec_cb::Callback
    dhTdp_vec_cb::Callback
    dhTdT_vec_cb::Callback

    dfdx_cb::Callback
    dfdxtrans_cb::Callback
    dfdt_cb::Callback
    dHdxdt_cb::Callback

    function ModelCallbacks()
        return new( nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing,
                    nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing,
                    nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing )
    end
end

"""
    Model( model::Model; states=0, inputs=0, parameters=0, eqcon=0, ineqcon=0, terminal_eqcon=0, terminal_ineqcon=0 )

The model describing the problem to solve using GRAMPC. The arguments set the dimensions of the created model.

* `states` is the number of states the dynamics have
* `inputs` is the number of inputs the dynamics have
* `parameters` is the number of parameters in the problem
* `eqcon` is the number of equality constraints in the problem
* `ineqcon` is the number of inequality constraints in the problem
* `terminal_eqcon` is the number of terminal equality constraints in the problem
* `terminal_ineqcon` is the number of terminal inequality constraints in the problem
"""
mutable struct Model{T}
    # C pointer to the underlying GRAMPC solver problem
    problem::Ptr{CGRAMPC}
    Ccallbacks::Cgrampcjl_callbacks

    # Problem size
    Nx::Int     # Number of state variables
    Nu::Int     # Number of control inputs
    Np::Int     # Number of problem parameters
    Ng::Int     # Number of equality constraints
    Nh::Int     # Number of inequality constraints
    NgT::Int    # Number of terminal equality constraints
    NhT::Int    # Number of terminal inequality constraints

    # Which optimization parts are enabled
    optimizeu::Bool
    optimizep::Bool
    optimizeT::Bool
    eqCon::Bool
    ineqCon::Bool
    termEqCon::Bool
    termIneqCon::Bool
    integralCost::Bool
    terminalCost::Bool

    # Dynamics information, cost and constraints
    callbacks::ModelCallbacks
    massMatrix::Matrix{Float64}

    # User-supplied parameters that are passed to the callbacks
    param::T

    function Model(; states=0, inputs=0, parameters=0, eqcon=0, ineqcon=0, terminal_eqcon=0, terminal_ineqcon=0 )
        model = new{Vector{Float64}}( C_NULL, Cgrampcjl_callbacks(),
                                      states, inputs, parameters, eqcon, ineqcon, terminal_eqcon, terminal_ineqcon,
                                      true, false, false, false, false, false, false, false, false,
                                      ModelCallbacks(),
                                      Matrix{Float64}( undef, 0, 0 ),
                                      Vector{Float64}( undef, 0 ) )

        finalizer( GRAMPC.clean!, model )
        return model
    end

end

"""
    is_initialized( model::Model )

Return if the underlying C structure for the GRAMPC model has been initialized.
"""
function _is_initialized( model::Model )
    return model.problem != C_NULL
end

function clean!( model::Model )
    if _is_initialized( model )
        ccall( (:grampcjl_free, libgrampc_julia), Cvoid, (Ptr{CGRAMPC},), model.problem )
    end
end

# Julia wrapper for typeGRAMPCsol, which contains the (public) solution data
struct Solution
    x::Vector{Float64}
    u::Vector{Float64}
    p::Vector{Float64}

    T::Float64
    J::Vector{Float64} # A 2-element array
    constraint_norm::Float64
    penalty_norm::Float64

    grad_iter::Vector{Int32}

    status::GRAMPCStatus

    function Solution( model::Model )
        _is_initialized( model ) || return

        x = Float64[]
        u = Float64[]
        p = Float64[]
        J = Float64[]
        iter = Int32[]

        problem = unsafe_load( model.problem )
        csol    = unsafe_load( problem.sol )

        # The C array will be unallocated if there are no uses of that variable
        ( model.Nx > 0 ) && ( x = unsafe_wrap( Array, csol.xnext, model.Nx ) )
        ( model.Nu > 0 ) && ( u = unsafe_wrap( Array, csol.unext, model.Nu ) )
        ( model.Np > 0 ) && ( p = unsafe_wrap( Array, csol.pnext, model.Np ) )

        J = unsafe_wrap( Array, csol.J , 2 )

        # Number of gradient iterations every
        numIter = _get_setting( model, :MaxMultIter )
        iter = unsafe_wrap( Array, csol.iter, numIter )

        return new( x, u, p, csol.Tnext, J, csol.cfct, csol.pen, iter, GRAMPCStatus( csol.status ) )
    end
end

ext: jl


text:
"""
    GetIndexWithKMP(string::String, sub_string::String, ignore_case::Bool)::Int

Find if a string contain a given sub string with KMP Search, as an explanation might be too lengthy
As such, a detailed explanation can be found at https://towardsdatascience.com/pattern-search-with-the-knuth-morris-pratt-kmp-algorithm-8562407dba5b

Questions and answers
1. Why KMP instead of naive search?
- Given n = length of string , m = length of string P
- Then, the naive search will have time complexity O(mn) to find all occurrences of pattern P in S
- However, knuth-morris-pratt (KMP) algorithm will have tim complexity of O(m+n) to find all occurrences of pattern P in S
- Hence, KMP is more efficient than the naive approach

Example

```julia
GetIndexWithKMP("Hello I am Gervin", "Hello I am Gervin", false) # returns 1
GetIndexWithKMP("ABABDABACDABABCABAB", "ABABCABAB", false) # returns 11
GetIndexWithKMP("SeEms Like IgNOrE CaSe Work", "seems like ignore case work", true) # returns 1
```

Contributed By:- [Gervin Fung](https://github.com/GervinFung)
Note: This function will also allow ignoring cases also

"""

# by convention -1 will be returned if the string does not contain the sub-string given
# in this case, 0 will be returned for Julia with start index of 1
# reason is that -1 is for language with start index of 0
# Should be convenient for others as most programming language use -1 if the string does not contain the sub-string given
const NO_SUBSTRING_INDEX = 0
# 1 is the first index, unlike others, where 0 is the first index
# Should be convenient for others as most programming language use 0 as first index
const JULIA_FIRST_INDEX = 1

function create_suffic_array(
    pattern_length::Int,
    sub_string::String,
)::Vector{Int}

    # Longest Proper Prefix which is Suffix array
    lps::Vector{Int} = ones(Int, pattern_length)

    index::Int = JULIA_FIRST_INDEX
    i::Int = 1 + JULIA_FIRST_INDEX

    while (i < pattern_length)
        if sub_string[i] == sub_string[index]
            lps[i] = index + 1
            index += 1
            i += 1
        else
            if index != 1
                index = lps[index-1]
            else
                lps[i] = 1
                i += 1
            end
        end
    end

    return lps
end

# this function will be used to obtain the index which the substring was found
function get_index_with_kmp(
    string::String,
    sub_string::String,
    ignore_case::Bool,
)::Int
    string = ignore_case ? lowercase(string) : string
    sub_string = ignore_case ? lowercase(sub_string) : sub_string

    string_length::Int = length(string)
    substring_length::Int = length(sub_string)

    lps::Vector{Int} = create_suffic_array(substring_length, sub_string)

    # likewise, 1 is the first index, unlike others, where 0 is the first index
    i::Int = JULIA_FIRST_INDEX
    j::Int = JULIA_FIRST_INDEX

    while i < string_length + 1 && j < substring_length + 1
        if string[i] == sub_string[j]
            i += 1
            j += 1
        else
            if j != 1
                j = lps[j-1]
            else
                i += 1
            end
        end
    end

    return j == substring_length + 1 ? i - j + 1 : NO_SUBSTRING_INDEX
end

# optional function that returns boolean if string does contain the sub-string given
function contain_substring_with_kmp(
    string::String,
    sub_string::String,
    ignore_case::Bool,
)::Bool
    return get_index_with_kmp(string, sub_string, ignore_case) !=
           NO_SUBSTRING_INDEX
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
export letsjam, play

function letsjam(f::Function,
                 samprate::Integer = 44100,
                 channels::Integer = 1,
                 dtype::DataType = Float32)
    global jamstream
    initialize()
    jamstream = openStream(samprate, channels, dtype)
    f()
    closeStream(jamstream)
    jamstream = Nothing()
    terminate()
end

function play(data::Array)
    global jamstream
    if typeof(jamstream) != Nothing
        writeStream(jamstream, data)
    else
        error("play can only be called within a letsjam statement")
    end
end


text:
using AtomsBase
using Test

@testset "Chemical formula" begin
    @test chemical_formula([:H])                  == "H"
    @test chemical_formula([:H, :H])              == "H₂"
    @test chemical_formula([:H, :O, :H])          == "H₂O"
    @test chemical_formula([:O, :H, :O, :H])      == "H₂O₂"
    @test chemical_formula([:Ga, :N, :O, :H, :H]) == "GaH₂NO"
end


text:
###################################
## Bayesian Deep Learning
## via Gaussian Process Regression

using Distributions
using PDMats
using PyPlot

function calc_kernel_matrix(X, kernel)
    N = length(X)
    K = zeros(N, N)
    for n1 in 1 : N
        for n2 in n1 : N
            K[n1, n2] = kernel(X[n1], X[n2])
        end
    end
    return Symmetric(K)
end

function calc_kernel_sequence(X, x, kernel)
    N = length(X)
    seq = [kernel(X[n], x) for n in 1 : N]
    return seq
end

function predict(X_train, Y_train, X, kernel, sigma2_y)
    N = length(X)
    mu = zeros(N)
    sigma2 = zeros(N)
    N_train = length(Y_train)
    K = calc_kernel_matrix(X_train, kernel)
    invmat = inv(sigma2_y * eye(N_train) + K)
    for n in 1 : N
        seq = calc_kernel_sequence(X_train, X[n], kernel) # N dim
        mu[n] = seq' * invmat * Y_train
        sigma2[n] = sigma2_y + kernel(X[n], X[n]) - seq'*invmat*seq
    end

    return mu, sigma2
end

function plot_predict(X, mu, sigma2)
    fill_between(X, mu + sqrt.(sigma2), mu - sqrt.(sigma2), color="c", alpha=0.5)
    plot(X, mu, "-b")
end

###################
# common setting

# input
x_min = - 5.0
x_max = + 5.0
N = 100
X = linspace(x_min, x_max, N)

# noise parameter
sigma2_y = 0.1

###################
# kernels

# multinomial covariance function
M = 3
Sigma_W = diagm([10.0, 1.0, 0.1, 0.01])
multi(x, M) = [x^m for m in 0:M]
kernel_m(x1, x2) = trace(multi(x1, M)*multi(x2, M)'*Sigma_W)
K_m = calc_kernel_matrix(X, kernel_m)

# RBF covariance function
alpha = 1.0
beta = 1.0
kernel_g(x1, x2) = alpha * exp(-0.5*inv(beta^2)*(x1 - x2)^2)
K_g = calc_kernel_matrix(X, kernel_g)

# neural network (erf) covariance function
Sigma = 10.0*eye(2)
aug(x) = Array([1, x])
kernel_e(x1, x2) = (2/pi)*asin((2*trace(aug(x2)*aug(x1)'*Sigma))/sqrt((1+2*trace(aug(x1)*aug(x1)'*Sigma))*(1+2*trace(aug(x2)*aug(x2)'*Sigma))))
K_e = calc_kernel_matrix(X, kernel_e)

# neural network (ReLU) covariance function
rad(x1, x2) = acos(max(min(sum(x1.*x2)/(norm(x1)*norm(x2)), 1.0), -1.0))
kernel_r(x1, x2) = (1.0 / pi) * norm(aug(x1)) * norm(aug(x2)) * (sin(rad(aug(x1), aug(x2))) + (pi - rad(aug(x1), aug(x2)))*cos(rad(aug(x1), aug(x2))))
K_r = calc_kernel_matrix(X, kernel_r)

# deep neural network (ReLU) covariance function
L = 8
sigma2_b = 1.0
sigma2_w = 2.0
kernel_tmp(x1, x2, L) = L > 0 ? (k_11=kernel_tmp(x1, x1, L-1);k_22=kernel_tmp(x2, x2, L-1);theta=acos(kernel_tmp(x1, x2, L-1) / sqrt(k_11 * k_22));sigma2_b + (sigma2_w/(2*pi))*sqrt(k_11*k_22)*(sin(theta) + (pi - theta)*cos(theta)) ) : sigma2_b + sigma2_w*(x1'*x2 / length(x1))
kernel_d(x1, x2) = kernel_tmp(x1, x2, L)
K_d = calc_kernel_matrix(X, kernel_d)

###################
# training

# data
N_train_all = 40
X_train_all = linspace(-3.0, 3.0, N_train_all)
Y_train_all = sin.(X_train_all .* 2*pi / (maximum(X_train_all) - minimum(X_train_all)))

# sample from priors
Y_m = []
Y_g = []
Y_e = []
Y_r = []
Y_d = []
tiny = 1.0e-5
num_sample = 10
for _ in 1 : num_sample
    push!(Y_m, rand(MvNormal(zeros(N), tiny*eye(N) + K_m)))
    push!(Y_g, rand(MvNormal(zeros(N), tiny*eye(N) + K_g)))
    push!(Y_e, rand(MvNormal(zeros(N), tiny*eye(N) + K_e)))
    push!(Y_r, rand(MvNormal(zeros(N), tiny*eye(N) + K_r)))
    push!(Y_d, rand(MvNormal(zeros(N), tiny*eye(N) + K_d)))
end

dir_name = "fig"
if !isdir(dir_name)
    mkdir(dir_name)
end

figure("function sample", figsize=(16,8))
for n in 1 : N_train_all
    clf()

    N_train = n
    X_train = X_train_all[1:n]
    Y_train = Y_train_all[1:n]

    # prediction
    mu_m, sigma2_m = predict(X_train, Y_train, X, kernel_m, sigma2_y)
    mu_g, sigma2_g = predict(X_train, Y_train, X, kernel_g, sigma2_y)
    mu_e, sigma2_e = predict(X_train, Y_train, X, kernel_e, sigma2_y)
    mu_r, sigma2_r = predict(X_train, Y_train, X, kernel_r, sigma2_y)
    mu_d, sigma2_d = predict(X_train, Y_train, X, kernel_d, sigma2_y)

    ###################
    # visualization

    # prior
    for i in 1 : num_sample
        # sample Y
        subplot(2,5,1)
        title("cubic")
        plot(X, Y_m[i], "-")

        subplot(2,5,2)
        title("RBF")
        plot(X, Y_g[i], "-")

        subplot(2,5,3)
        title("NN (erf)")
        plot(X, Y_e[i], "-")

        subplot(2,5,4)
        title("NN (ReLU)")
        plot(X, Y_r[i], "-")

        subplot(2,5,5)
        title("DNN (ReLU)")
        plot(X, Y_d[i], "-")
    end

    # prediction
    y_min = -3
    y_max = +3

    subplot(2,5,6)
    plot_predict(X, mu_m, sigma2_m)
    plot(X_train, Y_train, "xk")
    xlim([x_min, x_max])
    ylim([y_min, y_max])

    subplot(2,5,7)
    plot_predict(X, mu_g, sigma2_g)
    plot(X_train, Y_train, "xk")
    xlim([x_min, x_max])
    ylim([y_min, y_max])

    subplot(2,5,8)
    plot_predict(X, mu_e, sigma2_e)
    plot(X_train, Y_train, "xk")
    xlim([x_min, x_max])
    ylim([y_min, y_max])

    subplot(2,5,9)
    plot_predict(X, mu_r, sigma2_r)
    plot(X_train, Y_train, "xk")
    xlim([x_min, x_max])
    ylim([y_min, y_max])

    subplot(2,5,10)
    plot_predict(X, mu_d, sigma2_d)
    plot(X_train, Y_train, "xk")
    xlim([x_min, x_max])
    ylim([y_min, y_max])

    savefig(@sprintf("fig/%03d.png", n))
end

max_stars_repo_name: sammy-suyama/mlblog

max_issues_repo_name: sammy-suyama/mlblog

max_forks_repo_name: sammy-suyama/mlblog


text:
# Julia wrapper for header: /usr/include/scip/pub_matrix.h
# Automatically generated using Clang.jl wrap_c

function SCIPmatrixGetColValPtr(matrix, col)
    ccall((:SCIPmatrixGetColValPtr, libscip), Ptr{Cdouble}, (Ptr{SCIP_MATRIX}, Cint), matrix, col)
end

function SCIPmatrixGetColIdxPtr(matrix, col)
    ccall((:SCIPmatrixGetColIdxPtr, libscip), Ptr{Cint}, (Ptr{SCIP_MATRIX}, Cint), matrix, col)
end

function SCIPmatrixGetColNNonzs(matrix, col)
    ccall((:SCIPmatrixGetColNNonzs, libscip), Cint, (Ptr{SCIP_MATRIX}, Cint), matrix, col)
end

function SCIPmatrixGetNColumns(matrix)
    ccall((:SCIPmatrixGetNColumns, libscip), Cint, (Ptr{SCIP_MATRIX},), matrix)
end

function SCIPmatrixGetColUb(matrix, col)
    ccall((:SCIPmatrixGetColUb, libscip), Cdouble, (Ptr{SCIP_MATRIX}, Cint), matrix, col)
end

function SCIPmatrixGetColLb(matrix, col)
    ccall((:SCIPmatrixGetColLb, libscip), Cdouble, (Ptr{SCIP_MATRIX}, Cint), matrix, col)
end

function SCIPmatrixGetColNUplocks(matrix, col)
    ccall((:SCIPmatrixGetColNUplocks, libscip), Cint, (Ptr{SCIP_MATRIX}, Cint), matrix, col)
end

function SCIPmatrixGetColNDownlocks(matrix, col)
    ccall((:SCIPmatrixGetColNDownlocks, libscip), Cint, (Ptr{SCIP_MATRIX}, Cint), matrix, col)
end

function SCIPmatrixGetVar(matrix, col)
    ccall((:SCIPmatrixGetVar, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_MATRIX}, Cint), matrix, col)
end

function SCIPmatrixGetColName(matrix, col)
    ccall((:SCIPmatrixGetColName, libscip), Cstring, (Ptr{SCIP_MATRIX}, Cint), matrix, col)
end

function SCIPmatrixGetRowValPtr(matrix, row)
    ccall((:SCIPmatrixGetRowValPtr, libscip), Ptr{Cdouble}, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetRowIdxPtr(matrix, row)
    ccall((:SCIPmatrixGetRowIdxPtr, libscip), Ptr{Cint}, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetRowNNonzs(matrix, row)
    ccall((:SCIPmatrixGetRowNNonzs, libscip), Cint, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetRowName(matrix, row)
    ccall((:SCIPmatrixGetRowName, libscip), Cstring, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetNRows(matrix)
    ccall((:SCIPmatrixGetNRows, libscip), Cint, (Ptr{SCIP_MATRIX},), matrix)
end

function SCIPmatrixGetRowLhs(matrix, row)
    ccall((:SCIPmatrixGetRowLhs, libscip), Cdouble, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetRowRhs(matrix, row)
    ccall((:SCIPmatrixGetRowRhs, libscip), Cdouble, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixIsRowRhsInfinity(matrix, row)
    ccall((:SCIPmatrixIsRowRhsInfinity, libscip), UInt32, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetNNonzs(matrix)
    ccall((:SCIPmatrixGetNNonzs, libscip), Cint, (Ptr{SCIP_MATRIX},), matrix)
end

function SCIPmatrixGetRowMinActivity(matrix, row)
    ccall((:SCIPmatrixGetRowMinActivity, libscip), Cdouble, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetRowMaxActivity(matrix, row)
    ccall((:SCIPmatrixGetRowMaxActivity, libscip), Cdouble, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetRowNMinActNegInf(matrix, row)
    ccall((:SCIPmatrixGetRowNMinActNegInf, libscip), Cint, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetRowNMinActPosInf(matrix, row)
    ccall((:SCIPmatrixGetRowNMinActPosInf, libscip), Cint, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetRowNMaxActNegInf(matrix, row)
    ccall((:SCIPmatrixGetRowNMaxActNegInf, libscip), Cint, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetRowNMaxActPosInf(matrix, row)
    ccall((:SCIPmatrixGetRowNMaxActPosInf, libscip), Cint, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixGetCons(matrix, row)
    ccall((:SCIPmatrixGetCons, libscip), Ptr{SCIP_CONS}, (Ptr{SCIP_MATRIX}, Cint), matrix, row)
end

function SCIPmatrixUplockConflict(matrix, col)
    ccall((:SCIPmatrixUplockConflict, libscip), UInt32, (Ptr{SCIP_MATRIX}, Cint), matrix, col)
end

function SCIPmatrixDownlockConflict(matrix, col)
    ccall((:SCIPmatrixDownlockConflict, libscip), UInt32, (Ptr{SCIP_MATRIX}, Cint), matrix, col)
end

function SCIPmatrixCreate(scip, matrixptr, initialized, complete)
    ccall((:SCIPmatrixCreate, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_MATRIX}}, Ptr{UInt32}, Ptr{UInt32}), scip, matrixptr, initialized, complete)
end

function SCIPmatrixFree(scip, matrix)
    ccall((:SCIPmatrixFree, libscip), Cvoid, (Ptr{SCIP_}, Ptr{Ptr{SCIP_MATRIX}}), scip, matrix)
end

function SCIPmatrixPrintRow(scip, matrix, row)
    ccall((:SCIPmatrixPrintRow, libscip), Cvoid, (Ptr{SCIP_}, Ptr{SCIP_MATRIX}, Cint), scip, matrix, row)
end

function SCIPmatrixGetParallelRows(scip, matrix, scale, pclass)
    ccall((:SCIPmatrixGetParallelRows, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_MATRIX}, Ptr{Cdouble}, Ptr{Cint}), scip, matrix, scale, pclass)
end

function SCIPmatrixGetParallelCols(scip, matrix, scale, pclass, varineq)
    ccall((:SCIPmatrixGetParallelCols, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_MATRIX}, Ptr{Cdouble}, Ptr{Cint}, Ptr{UInt32}), scip, matrix, scale, pclass, varineq)
end


text:
# This file was generated, do not modify it. # hide
cv3 = CV(; nfolds=3)
res = evaluate(rf_mdl, X, y, resampling=CV(shuffle=true),
               measure=rms, verbosity=0)


text:
Pkg.add("DataFrames")
Pkg.add("Gadfly")
Pkg.add("Cairo")
Pkg.add("Fontconfig")
using DataFrames
using Gadfly
using Cairo
using Fontconfig

# Plot all validities

style = Theme(
    major_label_font_size=35pt,
    minor_label_font_size=30pt,
    background_color="white",
    major_label_font="CMU Serif",
    minor_label_font="CMU Serif"
)

println("All Validities...")
table = readtable("output/allValidities.csv");
labels = ["0", "200", "400", "600", "800", "1000"];
yticks=log2.([1,10,100,1000]);
output = plot(
    table,
    x=:validity,
    Geom.histogram(bincount=50),
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Coord.Cartesian(xmin=0, xmax=1, ymin=0,ymax=11),
    Guide.yticks(ticks=yticks),
    Guide.xlabel("Validity"),
    Guide.ylabel("# of Sessions"),
    style
);
draw(PDF("graphs/allValidities.pdf", 1600px, 800px), output);

println("Anonymous Validities...")
table=readtable("output/anonymousValidities.csv")
output = plot(
    table,
    x=:validity,
    Geom.histogram(bincount=50),
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Coord.Cartesian(xmin=0, xmax=1, ymin=0,ymax=11),
    Guide.yticks(ticks=yticks),
    Guide.xlabel("Validity"),
    Guide.ylabel("# of Sessions"),
    #Guide.title("Validities of anonymous users"),
    Theme(
        major_label_font_size=35pt,
        minor_label_font_size=30pt,
        background_color="white",
        major_label_font="CMU Serif",
        minor_label_font="CMU Serif",
        default_color="red"
    )
);
draw(PDF("graphs/anonymousValidities.pdf", 1600px, 800px), output);

println("Microworker Validities...")
table = readtable("output/microworkerValidities.csv")
output = plot(
    table,
    x=:validity,
    Geom.histogram(bincount=50),
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Coord.Cartesian(xmin=0, xmax=1, ymin=0,ymax=11),
    Guide.yticks(ticks=yticks),
    Guide.xlabel("Validity"),
    Guide.ylabel("# of Sessions"),
    #Guide.title("Validities of Microworkers"),
    Theme(
            major_label_font_size=35pt,
            minor_label_font_size=30pt,
            background_color="white",
            major_label_font="CMU Serif",
            minor_label_font="CMU Serif",
            default_color="green"
        )
);
draw(PDF("graphs/microworkerValidities.pdf", 1600px, 800px), output);

println("Player Validities...")
table = readtable("output/playerValidities.csv")
output = plot(
    table,
    x=:validity,
    Geom.histogram(bincount=50),
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Coord.Cartesian(xmin=0, xmax=1, ymin=0,ymax=11),
    Guide.yticks(ticks=yticks),
    Guide.xlabel("Validity"),
    Guide.ylabel("# of Sessions"),
    #Guide.title("Validities of Players"),
    Theme(
            major_label_font_size=35pt,
            minor_label_font_size=30pt,
            background_color="white",
            major_label_font="CMU Serif",
            minor_label_font="CMU Serif",
            default_color="orange"
        )
);
draw(PDF("graphs/playerValidities.pdf", 1600px, 800px), output);

println("All Answers...")
table = readtable("output/allAnswers.csv")
output = plot(
    table,
    x=:answer,
    y=:count,
    Guide.yticks(ticks=[0;50000;100000;150000]),
    Scale.y_continuous(labels=d-> @sprintf("%d",d/10000)),
    Geom.bar(position=:dodge),
    Coord.Cartesian(ymax=150000),
    Guide.xlabel("Answer"),
    Guide.ylabel("# of answers x10'000"),
    Theme(
        major_label_font_size=35pt,
        minor_label_font_size=30pt,
        background_color="white",
        major_label_font="CMU Serif",
        minor_label_font="CMU Serif",
        bar_spacing=150px
    )
);
draw(PDF("graphs/allAnswers.pdf", 1600px, 800px), output);

println("Microworkers Answers...")
table = readtable("output/microworkerAnswers.csv")
output = plot(
    table,
    x=:answer,
    y=:count,
    Guide.yticks(ticks=[0;50000;100000;150000]),
    Scale.y_continuous(labels=d-> @sprintf("%d",d/10000)),
    Geom.bar(position=:dodge),
    Coord.Cartesian(ymax=150000),
    Guide.xlabel("Answer"),
    Guide.ylabel("# of answers x10'000"),
    Theme(
            major_label_font_size=35pt,
            minor_label_font_size=30pt,
            background_color="white",
            major_label_font="CMU Serif",
            minor_label_font="CMU Serif",
            default_color="green",
            bar_spacing=150px
        )
);
draw(PDF("graphs/microworkerAnswers.pdf", 1600px, 800px), output);

println("Player Answers...")
table = readtable("output/playerAnswers.csv")
output = plot(
    table,
    x=:answer,
    y=:count,
    Guide.yticks(ticks=[0;50000;100000;150000]),
    Scale.y_continuous(labels=d-> @sprintf("%d",d/10000)),
    Geom.bar(position=:dodge),
    Coord.Cartesian(ymax=150000),
    Guide.xlabel("Answer"),
    Guide.ylabel("# of answers x10'000"),
    Theme(
            major_label_font_size=35pt,
            minor_label_font_size=30pt,
            background_color="white",
            major_label_font="CMU Serif",
            minor_label_font="CMU Serif",
            default_color="orange",
            bar_spacing=150px
        )
);
draw(PDF("graphs/playerAnswers.pdf", 1600px, 800px), output);

println("Anonymous Answers...")
table = readtable("output/anonymousAnswers.csv")
output = plot(
    table,
    x=:answer,
    y=:count,
    Guide.yticks(ticks=[0;50000;100000;150000]),
    Scale.y_continuous(labels=d-> @sprintf("%d",d/10000)),
    Geom.bar(position=:dodge),
    Coord.Cartesian(ymax=150000),
    Guide.xlabel("Answer"),
    Guide.ylabel("# of answers x10'000"),
    Theme(
            major_label_font_size=35pt,
            minor_label_font_size=30pt,
            background_color="white",
            major_label_font="CMU Serif",
            minor_label_font="CMU Serif",
            default_color="red",
            bar_spacing=150px
        )
);
draw(PDF("graphs/anonymousAnswers.pdf", 1600px, 800px), output);

labels = ["0","10 min","20 min","30 min","40 min", "50 min"];
xticks = [0,600,1200,1800,2400,3000];

println("All Durations...")
table = readtable("output/allDurations.csv")
output = plot(
    table,
    x=:duration,
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Geom.histogram(),
    Coord.Cartesian(xmin=0, xmax=3100, ymin=0,ymax=10),
    Guide.yticks(ticks=yticks),
    Scale.x_continuous(labels = x -> get(labels, round(Int, x/600)+1, "")),
    Guide.xticks(ticks=xticks),
    Guide.xlabel("Time"),
    Guide.ylabel("# of Sessions"),
    style
);
draw(PDF("graphs/allDurations.pdf", 1600px, 800px), output);

println("Microworker Durations...")
table = readtable("output/microworkerDurations.csv")
output = plot(
    table,
    x=:duration,
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Geom.histogram(),
    Coord.Cartesian(xmin=0, xmax=3100, ymin=0,ymax=10),
    Guide.yticks(ticks=yticks),
    Scale.x_continuous(labels = x -> get(labels, round(Int, x/600)+1, "")),
    Guide.xticks(ticks=xticks),
    Guide.xlabel("Time"),
    Guide.ylabel("# of Sessions"),
    Theme(
            major_label_font_size=35pt,
            minor_label_font_size=30pt,
            background_color="white",
            major_label_font="CMU Serif",
            minor_label_font="CMU Serif",
            default_color="green"
        )
);
draw(PDF("graphs/microworkerDurations.pdf", 1600px, 800px), output);

println("Player Durations...")
table = readtable("output/playerDurations.csv")
output = plot(
    table,
    x=:duration,
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Geom.histogram(),
    Coord.Cartesian(xmin=0, xmax=3100, ymin=0,ymax=10),
    Guide.yticks(ticks=yticks),
    Scale.x_continuous(labels = x -> get(labels, round(Int, x/600)+1, "")),
    Guide.xticks(ticks=xticks),
    Guide.xlabel("Time"),
    Guide.ylabel("# of Sessions"),
    Theme(
            major_label_font_size=35pt,
            minor_label_font_size=30pt,
            background_color="white",
            major_label_font="CMU Serif",
            minor_label_font="CMU Serif",
            default_color="orange"
        )
);
draw(PDF("graphs/playerDurations.pdf", 1600px, 800px), output);

println("Anonymous Durations...")
table = readtable("output/anonymousDurations.csv")
output = plot(
    table,
    x=:duration,
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Geom.histogram(),
    Coord.Cartesian(xmin=0, xmax=3100, ymin=0,ymax=10),
    Guide.yticks(ticks=yticks),
    Scale.x_continuous(labels = x -> get(labels, round(Int, x/600)+1, "")),
    Guide.xticks(ticks=xticks),
    Guide.xlabel("Time"),
    Guide.ylabel("# of Sessions"),
    Theme(
            major_label_font_size=35pt,
            minor_label_font_size=30pt,
            background_color="white",
            major_label_font="CMU Serif",
            minor_label_font="CMU Serif",
            default_color="red"
        )
);
draw(PDF("graphs/anonymousDurations.pdf", 1600px, 800px), output);

labels = ["0","0.5 min","1 min","1.5 min","2 min"];
xticks = [0,30,60,90,120];

println("All Durations Normalized...")
table = readtable("output/allDurationsNormalized.csv")
output = plot(
    table,
    x=:durationNormalized,
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Geom.histogram(),
    Coord.Cartesian(xmin=0, xmax=125, ymin=0,ymax=12),
    Guide.yticks(ticks=yticks),
    Scale.x_continuous(labels = x -> get(labels, round(Int, x/30)+1, "")),
    Guide.xticks(ticks=xticks),
    Guide.xlabel("Time"),
    Guide.ylabel("# of Sessions"),
    style
);
draw(PDF("graphs/allDurationsNormalized.pdf", 1600px, 800px), output);

println("Microworker Durations Normalized...")
table = readtable("output/microworkerDurationsNormalized.csv")
output = plot(
    table,
    x=:durationNormalized,
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Geom.histogram(),
    Coord.Cartesian(xmin=0, xmax=125, ymin=0,ymax=12),
    Guide.yticks(ticks=yticks),
    Scale.x_continuous(labels = x -> get(labels, round(Int, x/30)+1, "")),
    Guide.xticks(ticks=xticks),
    Guide.xlabel("Time"),
    Guide.ylabel("# of Sessions"),
    Theme(
            major_label_font_size=35pt,
            minor_label_font_size=30pt,
            background_color="white",
            major_label_font="CMU Serif",
            minor_label_font="CMU Serif",
            default_color="green"
        )
);
draw(PDF("graphs/microworkerDurationsNormalized.pdf", 1600px, 800px), output);

println("Player Durations Normalized...")
table = readtable("output/playerDurationsNormalized.csv")
output = plot(
    table,
    x=:durationNormalized,
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Geom.histogram(),
    Coord.Cartesian(xmin=0, xmax=125, ymin=0,ymax=12),
    Guide.yticks(ticks=yticks),
    Scale.x_continuous(labels = x -> get(labels, round(Int, x/30)+1, "")),
    Guide.xticks(ticks=xticks),
    Guide.xlabel("Time"),
    Guide.ylabel("# of Sessions"),
    Theme(
            major_label_font_size=35pt,
            minor_label_font_size=30pt,
            background_color="white",
            major_label_font="CMU Serif",
            minor_label_font="CMU Serif",
            default_color="orange"
        )
);
draw(PDF("graphs/playerDurationsNormalized.pdf", 1600px, 800px), output);

println("Anonymous Durations Normalized...")
table = readtable("output/anonymousDurationsNormalized.csv")
output = plot(
    table,
    x=:durationNormalized,
    Scale.y_log2(labels=d-> @sprintf("%d",2^d)),
    Geom.histogram(),
    Coord.Cartesian(xmin=0, xmax=125, ymin=0,ymax=12),
    Guide.yticks(ticks=yticks),
    Scale.x_continuous(labels = x -> get(labels, round(Int, x/30)+1, "")),
    Guide.xticks(ticks=xticks),
    Guide.xlabel("Time"),
    Guide.ylabel("# of Sessions"),
    Theme(
            major_label_font_size=35pt,
            minor_label_font_size=30pt,
            background_color="white",
            major_label_font="CMU Serif",
            minor_label_font="CMU Serif",
            default_color="red"
        )
);
draw(PDF("graphs/anonymousDurationsNormalized.pdf", 1600px, 800px), output);

yticks=[0,1,2,3,4,5,6,7,8,9,10];

println("Unique workers batch 1...")
table = readtable("output/Batch_2823797_batch_resultsWorkers.csv")
output = plot(
    table,
    x=:count,
    Scale.y_sqrt(labels=d-> @sprintf("%d",d^2)),
    Coord.Cartesian(xmin=0, xmax=90, ymin=0,ymax=10),
    Guide.yticks(ticks=yticks),
    Geom.histogram(),
    Guide.xlabel("# of HITs"),
    Guide.ylabel("# of Workers"),
    style
);
draw(PDF("graphs/Batch-2823797-batch-resultsWorkers.pdf", 1600px, 800px), output);

println("Unique workers batch 2...")
table = readtable("output/Batch_2823956_batch_resultsWorkers.csv")
output = plot(
    table,
    x=:count,
    Scale.y_sqrt(labels=d-> @sprintf("%d",d^2)),
    Coord.Cartesian(xmin=0, xmax=90, ymin=0,ymax=10),
    Guide.yticks(ticks=yticks),
    Geom.histogram(),
    Guide.xlabel("# of HITs"),
    Guide.ylabel("# of Workers"),
    style
);
draw(PDF("graphs/Batch-2823956-batch-resultsWorkers.pdf", 1600px, 800px), output);

println("Unique workers batch 3...")
table = readtable("output/Batch_2848581_batch_resultsWorkers.csv")
output = plot(
    table,
    x=:count,
    Scale.y_sqrt(labels=d-> @sprintf("%d",d^2)),
    Coord.Cartesian(xmin=0, xmax=90, ymin=0,ymax=10),
    Guide.yticks(ticks=yticks),
    Geom.histogram(),
    Guide.xlabel("# of HITs"),
    Guide.ylabel("# of Workers"),
    style
);
draw(PDF("graphs/Batch-2848581-batch-resultsWorkers.pdf", 1600px, 800px), output);

println("Unique workers batch 4...")
table = readtable("output/Batch_2853164_batch_resultsWorkers.csv")
output = plot(
    table,
    x=:count,
    Scale.y_sqrt(labels=d-> @sprintf("%d",d^2)),
    Coord.Cartesian(xmin=0, xmax=90, ymin=0,ymax=10),
    Guide.yticks(ticks=yticks),
    Geom.histogram(),
    Guide.xlabel("# of HITs"),
    Guide.ylabel("# of Workers"),
    style
);
draw(PDF("graphs/Batch-2853164-batch-resultsWorkers.pdf", 1600px, 800px), output);

yticks=[0,2,4,6,8,10,12,14,16];

println("Unique workers of all batches...")
table = readtable("output/allWorkers.csv")
output = plot(
    table,
    x=:count,
    Scale.y_sqrt(labels=d-> @sprintf("%d",d^2)),
    Coord.Cartesian(xmin=0, xmax=90, ymin=0,ymax=16),
    Guide.yticks(ticks=yticks),
    Geom.histogram(),
    Guide.xlabel("# of HITs"),
    Guide.ylabel("# of Workers"),
    style
);
draw(PDF("graphs/allWorkers.pdf", 1600px, 800px), output);

println("Played games...")
table = readtable("output/gamesPlayed.csv")
output = plot(
    table,
    x=:gamesPlayed,
    Geom.histogram(),
    Coord.Cartesian(xmax=320),
    Guide.xlabel("# of players"),
    Guide.ylabel("# of games played"),
    style
);
draw(PDF("graphs/gamesPlayed.pdf", 1600px, 800px), output);

println("Number of days playing...")
table = readtable("output/daysOfPlaying.csv")
output = plot(
    table,
    x=:days,
    Geom.histogram(),
    Coord.Cartesian(ymin=0),
    Guide.xlabel("days"),
    Guide.ylabel("# of players"),
    style
);
draw(PDF("graphs/daysOfPlaying.pdf", 1600px, 800px), output);

println("Games per day...")
table = readtable("output/gamesPerDay.csv")
output = plot(
    table,
    x=:games,
    Geom.histogram(),
    Guide.xlabel("Games per day"),
    Guide.ylabel("# of players"),
    style
);
draw(PDF("graphs/gamesPerDay.pdf", 1600px, 800px), output);

println("Games per day overall...")
table = readtable("output/gamesPerDayOverall.csv")
output = plot(
    table,
    x=:games,
    Geom.histogram(),
    Guide.xlabel("Games per day"),
    Guide.ylabel("# of players"),
    style
);
draw(PDF("graphs/gamesPerDayOverall.pdf", 1600px, 800px), output);

println("Random Validities...")
table = readtable("output/randomValidities.csv")
output = plot(
    table,
    x=:validity,
    Geom.histogram(),
    Guide.xlabel("Validity"),
    Guide.ylabel("# of sessions"),
    style
);
draw(PDF("graphs/randomValidities.pdf", 1600px, 800px), output);


text:
"""
    basicExperiment(; genj::GeneticJuliaStructure = GenJ)

Runs a single experiment with the tools and methods specified in `genj`.

# Keyword Arguments
- `genj::GeneticJuliaStructure = GenJ`: the main structure. For `code user`: don't modify the default value unless multiple experiments are going to be run.

See also: [`run`](@ref)
"""
function basicExperiment(; genj::GeneticJuliaStructure = GenJ)

    initTime(genj)

    genPopulation!(genj)
    evaluate!(genj, genj._population)

    initBestFitness(genj)

    while !(reached(genj))

        saveResults(genj)

        selectedParents = selectParents(genj)

        offspring = cross(genj, selectedParents)
        mutatedOffspring = mutate(genj, offspring)

        evaluate!(genj, mutatedOffspring)

        replacePopulation!(genj, mutatedOffspring)
    end

    saveResults(genj)
end # function

"""
    SPEA(archiveSize::Integer, distanceFunction::Function=euclideanDistance; genj::GeneticJuliaStructure = GenJ)

documentation
"""
function SPEA(archiveSize::Integer, distanceFunction::Function=euclideanDistance; genj::GeneticJuliaStructure = GenJ)

    getDominant(population::Array{Individual}) = filter(x->getGlobalFitness(x)==1, population)

    function truncateArchive(archive::Array{Individual})

        k = round(Integer, sqrt(archiveSize + genj._generator._popSize))
        rest = length(archive) - archiveSize
        fitnesses = getFitness(archive)
        distanceMatrix = pairwiseDistance(fitnesses, fitnesses, distanceFunction)

        for i=1:rest
            nearestDistances = sum(distanceMatrix[1, getKnearest(distanceMatrix[1,:], k)])
            bestIndex = 1
            for j=2:size(distanceMatrix)[1]
                current = sum(distanceMatrix[j, getKnearest(distanceMatrix[j,:], k)])
                if current < nearestDistances
                    nearestDistances = current
                    bestIndex = j
                end
            end
            distanceMatrix = distanceMatrix[1:end .!= bestIndex, 1:end .!= bestIndex]
            archive = archive[1:end .!= bestIndex]
        end

        return archive
    end

    initTime(genj)
    genPopulation!(genj)
    evaluate!(genj, genj._population)

    initBestFitness(genj)

    archive = getDominant(genj._population)

    while !reached(genj)

        saveResults(genj)

        if length(archive) > archiveSize
            archive=truncateArchive(archive)
        end
        selectedParents = selectParents(genj, union(genj._population, archive))

        offspring = cross(genj, selectedParents)
        mutatedOffspring = mutate(genj, offspring)

        evaluate!(genj, mutatedOffspring)

        replacePopulation!(genj, mutatedOffspring)

        archive = getDominant(union(genj._population, archive))
    end

    saveResults(genj)
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
#https://discourse.julialang.org/t/meshgrid-function-in-julia/48679/4
function meshgrid(x,y)
    lx = length(x)
    ly = length(y)
    X = x'.*ones(lx)
    Y = ones(ly)'.*y
    return X,Y
end

#Unrolls the dictionary to form name-value pairs of the form
function namevaluepairs(options)
    optionsString=""
    for option in options
        optionName = option[1]
        optionValue = option[2]
        if optionValue isa String
            optionsString *= ",\"$optionName\",\"$optionValue\""
        else
            optionsString *= ",\"$optionName\",$optionValue"
        end
    end
    return optionsString
end


text:
function standardDeviationSD(x)

    sqrt(varianceSD(x))

end


text:
using Test
using MiniCheetahRobot
using RigidBodyDynamics

@testset "MiniCheetahRobot" begin
    mechanism = MiniCheetahRobot.mechanism()
    @test num_velocities(mechanism) == 21
    meshdir = joinpath(MiniCheetahRobot.packagepath(), "cheetah_description", "meshes")
    @test isdir(meshdir)
    @test count(x -> endswith(x, ".obj"), readdir(meshdir)) == 4
end


text:
function Base.show(io::IO, g::Grouped)::Void
    groupbys = g.groupbys
    # nargs = length(groupbys)
    predicates = g.metadata[:predicates]
    @printf(io, "Grouped %s\n", name(typeof(g).parameters[1]))
    println(io, "Groupings by:")
    for groupby in groupbys
        print(io, " "^4)
        print_groupby(io, groupby, predicates)
    end
    println()
    print(io, "Source: "); show(io, g.src)
    return
end

function print_groupby(io, groupby, predicates)::Void
    if haskey(predicates, groupby)
        predicate = predicates[groupby]
        @printf(
            io, "%s (with alias :%s) \n", string(predicate), string(groupby)
        )
    else
        @printf(io, "%s \n", string(groupby))
    end
    return
end


text:
#==============================================================================

          Solving the HJB with a state constraint using
                the explicit method

               Translated from matlab code on Ben Moll's website:
                              http://www.princeton.edu/~moll/HACTproject.htm

==============================================================================#

using LinearAlgebra, SparseArrays, Plots

γ = 2 # parameter from CRRA utility
r = 0.03 # the interest rate
ρ = 0.05 # the discount rate

z1=.1
z2 =.2
z= [z1 z2]

λ1 = 0.02
λ2 = 0.03

λ = [λ1 λ2]

H = 500 # number of points in the grid space
amin=-0.02
amax=2

a = LinRange(amin,amax,H)
da = (amax-amin)/(H-1)

aa = [a a]
zz = ones(H,1)*z

maxit = 100 # the maximum number of iterations we allow the finite differencing algorithm
crit = 10^(-6) # our critical value
Δ = 1000

dVf=zeros(H,2)
dVb = zeros(H,2)
c= zeros(H,2)

# Create the matrix that captures the systems dynamic
Aswitch = [-sparse(I,H,H)*λ[1] sparse(I,H,H)*λ[1]; sparse(I,H,H)*λ[2] -sparse(I,H,H)*λ[2]]

# Inital Guess
v0=zeros(H,2)
v0[:,1] = (z[1] .+ r.*a).^(1-γ)/(1-γ)/ρ
v0[:,2] = (z[2] .+ r.*a).^(1-γ)/(1-γ)/ρ
global v=v0

dist=[]
V_n=[]

# The finite differeing loop

for n in 1:maxit
    global V=v
    push!(V_n,V)
    # Forward differencing
    dVf[1:H-1,:]=(V[2:H,:]-V[1:H-1,:])/da
    dVf[H,:] = (z .+ r.*amax).^(-γ) # impose state constraint at the max, just in case
    # Backward differencing
    dVb[2:H,:]=(V[2:H,:]-V[1:H-1,:])/da
    dVb[1,:]= (z .+r.*amin).^(-γ)

    # Find consumption and savings with the forward difference
    cf = dVf.^(-1/γ)
    sf = zz +r.*aa -cf

    # Find consumption and savings with the backward difference
    cb= dVb.^(-1/γ)
    sb = zz +r.*aa -cb

    # Find consumption and savings at steady state
    c0 = zz + r.*aa
    dV0 = c0.^(-γ)

    # Now implement the upwind scheme in order to select the best differencing method
    If = sf .> 0
    Ib = sb .< 0
    I0 = (1 .- If - Ib)

    # State constraint at amin is automatically implemented

    global dV_Upwind = dVf.*If + dVb.*Ib + dV0.*I0

    global c = dV_Upwind.^(-1/γ)
    global u = c.^(1-γ)/(1-γ)
    # Construct matrix for the evolution of the system
    X = -min.(sb,0)/da
    Y= -max.(sf,0)/da + min.(sb,0)/da
    Z = max.(sf,0)/da

    A1 = spdiagm(0=>Y[:,1], -1 => X[2:H,1], 1=> Z[1:H-1,1])
    A2 = spdiagm(0=>Y[:,2], -1 => X[2:H,2], 1=> Z[1:H-1,2])
    AA = [A1 spzeros(H,H); spzeros(H,H) A2]

    A = AA + Aswitch

    B = (ρ +1/Δ)*sparse(I,2*H,2*H) -A
    u_stacked = [u[:,1]; u[:,2]]
    V_stacked = [V[:,1]; V[:,2]]

    b = u_stacked + V_stacked/Δ
    V_stacked = B\b # Solves the system of equations

    V = [V_stacked[1:H] V_stacked[H+1:2*H]]

    V_change = V - v

    global v = V

    push!(dist, findmax(abs.(V_change))[1])

    if dist[n] < crit
        println("Value Function converged, Iteration=")
        println(n)
        break
    end
end

# Graphs
adot = zz + r.*aa -c

plot(dist[:], xlabel="Iteration", ylabel="\$\\lvert V^{n+1}-V^{n}\\rvert\$")

plot(a,V, xlabel="a", ylabel="\$V_{i}(a)\$", xlims=(amin,amax),legend=false)

plot(a,c,xlabel="a", ylabel="\$c_{i}(a)\$", xlims=(amin,amax), legend=false)

plot(a, adot, xlabel="a", ylabel="\$s_{i}(a)\$", xlims=(amin,amax))


text:
using Gadfly
using Colors
using DataFrames
using Query
using CSV
using Statistics
import Cairo, Fontconfig
import ColorSchemes

Gadfly.push_theme(Theme(major_label_font="Helvetica", major_label_font_size=20pt,
                        minor_label_font="Helvetica", minor_label_font_size=12pt,
                        key_title_font="Helvetica", key_title_font_size=16pt,
                        key_label_font="Helvetica", key_label_font_size=14pt,
                        line_width=0.8mm, point_size=1.0mm,
                        highlight_width=0.1mm,
                        lowlight_color=c->RGBA{Float32}(c.r, c.g, c.b, 0.2),
                        default_color=colorant"#000000"))

colors = [colorant"#a6cee3", colorant"#1f78b4", colorant"#b2df8a",
          colorant"#33a02c", colorant"#fb9a99", colorant"#e31a1c"]

function op_color(p::Int64)
    # eval(Meta.parse(string("ColorSchemes.Set1_", max(p, 3))))
    # reshape(ColorSchemes.Paired_8, (2, 4))'[:]
    ColorSchemes.Set1_4
end

function get_evolution(logfile::String; id::String="", goalgen::Int64=10000)
    res = CSV.File(logfile; delim=" ",
                   header=[:timestamp, :darwin, :info, :seed, :generation,
                           :maxfit, :meanfit, :stdfit]) |> DataFrame
    if length(id) > 0
        res[:id] = id
    else
        res[:id] = string(res[:seed][1])
    end
    fits = ["coverage_fitness", "water_search_fitness",
            "water_memorize_fitness", "cross_search_fitness",
            "cross_memorize_fitness", "cross_strategy_fitness"]
    res[:fit] = mod.(floor.(Int64, (res[:generation] .- 1) / goalgen), 4) .+ 1
    res[:fit_id] = string.(res[:id], ", ", res[:fit])
    res
end

function plot_evolution(res::DataFrame; title::String="Evolution",
                        filename="evolution.pdf")
    plt = plot(res, x=:generation, y=:maxfit, color=:id,
               # shape=:id,
               ygroup=:fit,
               Geom.subplot_grid(free_y_axis=true,
               # Geom.point,
               Geom.smooth(smoothing=0.2)),
               Guide.xlabel("Generation"),
               Guide.ylabel("Fitness"),
               # Scale.color_discrete_hue(op_color),
               Guide.title(title))
    draw(PDF(filename, 8inch, 8inch), plt)
end

function plot_all()
    res = get_evolution("results/88778/evolution.log", id="CGP",
                        goalgen=100)
    res = vcat(res, get_evolution("results/88779/evolution.log", id="Neurodevo",
                                  goalgen=100))
    plot_evolution(res; filename="evolution_100.pdf")

    res = get_evolution("results/88794/evolution.log",
                        id="CGP", goalgen=10)
    res = vcat(res, get_evolution("results/88796/evolution.log",
                                  id="Neurodevo", goalgen=10))
    plot_evolution(res; filename="evolution_10.pdf")
end


text:
using Distributions, Random
using Test, LinearAlgebra, PDMats

v = 7.0
S = Matrix(1.0I, 2, 2)
S[1, 2] = S[2, 1] = 0.5

W = Wishart(v,S)
IW = InverseWishart(v,S)

rng = MersenneTwister(123)

@testset "Testing matrix-variates with $key" for (key, func) in
    Dict("rand(...)" => [rand, rand],
         "rand(rng, ...)" => [dist -> rand(rng, dist), (dist, n) -> rand(rng, dist, n)])

    for d in [W,IW]
        local d
        @test size(d) == size(func[1](d))
        @test length(d) == length(func[1](d))
        @test rank(d) == rank(func[1](d))
        @test typeof(d)(params(d)...) == d
        @test partype(d) == Float64
    end

    @test partype(Wishart(7, Matrix{Float32}(I, 2, 2))) == Float32
    @test partype(InverseWishart(7, Matrix{Float32}(I, 2, 2))) == Float32

    @test isapprox(mean(func[2](W,100000)) , mean(W) , atol=0.1)
    @test isapprox(mean(func[2](IW,100000)), mean(IW), atol=0.1)
    @test isapprox(cov(hcat(vec.(func[2](W, 1000000))...)'), cov(W) , atol=0.1)
    @test isapprox(cov(hcat(vec.(func[2](IW, 1000000))...)'), cov(IW), atol=0.1)

    v = 3.0

    @test isapprox(pdf(Wishart(v,S), S)         , 0.04507168, atol=1e-8)
    @test isapprox(pdf(Wishart(v,S), inv(S))    , 0.01327698, atol=1e-8)
    @test isapprox(pdf(Wishart(v,inv(S)),S)     , 0.0148086 , atol=1e-8)
    @test isapprox(pdf(Wishart(v,inv(S)),inv(S)), 0.01901462, atol=1e-8)

    @test isapprox(pdf(Wishart(v,S), [S, S])               , [0.04507168, 0.04507168], atol=1e-6)
    @test isapprox(pdf(Wishart(v,S), [inv(S), inv(S)])     , [0.01327698, 0.01327698], atol=1e-6)
    @test isapprox(pdf(Wishart(v,inv(S)), [S, S])          , [0.0148086, 0.0148086]  , atol=1e-6)
    @test isapprox(pdf(Wishart(v,inv(S)), [inv(S), inv(S)]), [0.01901462, 0.01901462], atol=1e-6)

    @test logpdf(Wishart(v,S), S)           ≈ log.(pdf(Wishart(v,S), S))
    @test logpdf(Wishart(v,S), inv(S))      ≈ log.(pdf(Wishart(v,S), inv(S)))
    @test logpdf(Wishart(v,inv(S)), S)      ≈ log.(pdf(Wishart(v,inv(S)), S))
    @test logpdf(Wishart(v,inv(S)), inv(S)) ≈ log.(pdf(Wishart(v,inv(S)), inv(S)))

    @test logpdf(Wishart(v,S), [S, S])                ≈ log.(pdf(Wishart(v,S), [S, S]))
    @test logpdf(Wishart(v,S), [inv(S), inv(S)])      ≈ log.(pdf(Wishart(v,S), [inv(S), inv(S)]))
    @test logpdf(Wishart(v,inv(S)), [S, S])           ≈ log.(pdf(Wishart(v,inv(S)), [S, S]))
    @test logpdf(Wishart(v,inv(S)), [inv(S), inv(S)]) ≈ log.(pdf(Wishart(v,inv(S)), [inv(S), inv(S)]))

    @test isapprox(pdf(InverseWishart(v,S), S)         , 0.04507168 , atol=1e-8)
    @test isapprox(pdf(InverseWishart(v,S), inv(S))    , 0.006247377, atol=1e-8)
    @test isapprox(pdf(InverseWishart(v,inv(S)),S)     , 0.03147137 , atol=1e-8)
    @test isapprox(pdf(InverseWishart(v,inv(S)),inv(S)), 0.01901462 , atol=1e-8)

    @test isapprox(pdf(InverseWishart(v,S), [S, S])               , [0.04507168,  0.04507168] , atol=1e-6)
    @test isapprox(pdf(InverseWishart(v,S), [inv(S), inv(S)])     , [0.006247377, 0.006247377], atol=1e-6)
    @test isapprox(pdf(InverseWishart(v,inv(S)), [S, S])          , [0.03147137,  0.03147137] , atol=1e-6)
    @test isapprox(pdf(InverseWishart(v,inv(S)), [inv(S), inv(S)]), [0.01901462,  0.01901462] , atol=1e-6)

    @test logpdf(InverseWishart(v,S), S)           ≈ log(pdf(InverseWishart(v,S), S))
    @test logpdf(InverseWishart(v,S), inv(S))      ≈ log(pdf(InverseWishart(v,S), inv(S)))
    @test logpdf(InverseWishart(v,inv(S)), S)      ≈ log(pdf(InverseWishart(v,inv(S)), S))
    @test logpdf(InverseWishart(v,inv(S)), inv(S)) ≈ log(pdf(InverseWishart(v,inv(S)), inv(S)))

    @test logpdf(InverseWishart(v,S), [S, S])                ≈ log.(pdf(InverseWishart(v,S), [S, S]))
    @test logpdf(InverseWishart(v,S), [inv(S), inv(S)])      ≈ log.(pdf(InverseWishart(v,S), [inv(S), inv(S)]))
    @test logpdf(InverseWishart(v,inv(S)), [S, S])           ≈ log.(pdf(InverseWishart(v,inv(S)), [S, S]))
    @test logpdf(InverseWishart(v,inv(S)), [inv(S), inv(S)]) ≈ log.(pdf(InverseWishart(v,inv(S)), [inv(S), inv(S)]))
end

@testset "Testing matrix-variates with $key" for (key, func) in
    Dict("rand!(..., false)" => [(dist, X) -> rand!(dist, X, false), rand!],
         "rand!(rng, ..., false)" => [(dist, X) -> rand!(rng, dist, X, false), (dist, X) -> rand!(rng, dist, X)])

    for d in [W,IW]
        local d
        local m = Matrix{Float64}(undef, size(d))
        local x = func[2](d, m)
        @test x ≡ m
        @test size(d) == size(x)
        @test length(d) == length(func[2](d, m))
        @test typeof(d)(params(d)...) == d
        @test partype(d) == Float64
    end

    m = [Matrix{partype(W)}(undef, size(W)) for _ in Base.OneTo(100000)]
    x = func[1](W,m)
    @test x ≡ m
    @test isapprox(mean(x) , mean(W) , atol=0.1)
    x = func[1](IW,m)
    @test x ≡ m
    @test isapprox(mean(x), mean(IW), atol=0.1)

    m3 = Array{partype(W), 3}(undef, size(W)..., 100000)
    x = func[2](W,m3)
    @test x ≡ m3
    @test isapprox(mean(x) , mean(mean(W)) , atol=0.1)
    x = func[2](IW,m3)
    @test x ≡ m3
    @test isapprox(mean(x), mean(mean(IW)), atol=0.1)
end

@testset "Testing matrix-variates with $key" for (key, func) in
    Dict("rand!(..., true)" => (dist, X) -> rand!(dist, X, true),
         "rand!(rng, true)" => (dist, X) -> rand!(rng, dist, X, true))

    m = Vector{Matrix{partype(W)}}(undef, 100000)
    x = func(W,m)
    @test x ≡ m
    @test isapprox(mean(x), mean(W) , atol=0.1)
    x = func(IW,m)
    @test x ≡ m
    @test isapprox(mean(x), mean(IW), atol=0.1)

end

repeats = 10
m = Vector{Matrix{partype(W)}}(undef, repeats)
rand!(W, m)
@test isassigned(m, 1)
m1=m[1]
rand!(W, m)
@test m1 ≡ m[1]
rand!(W, m, true)
@test m1 ≢ m[1]
m1 = m[1]
rand!(W, m, false)
@test m1 ≡ m[1]

@testset "Wishart conversion" for elty in (Float32, Float64, BigFloat)

    Del1 = convert(Wishart{elty}, W)
    Del2 = convert(Wishart{elty}, v, PDMat(S), W.c0)

    @test partype(Del1) == elty
    @test partype(Del2) == elty

end

@testset "InverseWishart conversion" for elty in (Float32, Float64, BigFloat)

    Del1 = convert(InverseWishart{elty}, IW)
    Del2 = convert(InverseWishart{elty}, v, PDMat(S), IW.c0)

    @test partype(Del1) == elty
    @test partype(Del2) == elty

end


text:
using SEC, Test, Statistics, Distributions, Random
import SEC.Sample: markov_link, markov_chain

@testset "markov link correctness" begin
    Random.seed!(1234)
    nll(X) = 10
    proposal(X) = 2
    current_X = 1
    current_nll = nll(current_X)
    (; new_X, new_nll, proposal_X, proposal_nll) = markov_link(
        nll = nll,
        proposal = proposal,
        current_nll = current_nll,
        current_X = current_X
    )
    @test new_X == 2
    @test new_nll == 10
    @test proposal_X == 2
    @test proposal_nll == 10
end

@testset "markov chain reproducibility" begin
    nll(X) = X^2 / 2
    proposal(X) = X + randn()
    seed_X = 0.0
    chain_length = 10000
    Random.seed!(1234)
    (; chain_X, chain_nll) = markov_chain(nll, proposal, seed_X, chain_length)
    original_chain_X = copy(chain_X)
    original_chain_nll = copy(chain_nll)
    Random.seed!(1234)
    (; chain_X, chain_nll) = markov_chain(nll, proposal, seed_X, chain_length)
    @test all((original_chain_X - chain_X) .== 0.0)
    @test all((original_chain_nll - chain_nll) .== 0.0)
    # hand unrolling loop
    Random.seed!(1234)
    current_X = 0.0
    current_nll = nll(current_X)
    proposal_X = proposal(seed_X)
    proposal_nll = nll(proposal_X)
    Δ = current_nll - proposal_nll
    accept_bool = log(rand(Uniform(0, 1))) < Δ
    current_X = accept_bool ? proposal_X : current_X
    @test chain_X[2] == current_X
end

@testset "markov chain convergence" begin
    Random.seed!(1234)
    nll(X) = X^2 / 2
    proposal(X) = X + randn()
    seed_X = 0.15
    chain_length = 10000
    (; chain_X, chain_nll) = markov_chain(nll, proposal, seed_X, chain_length)
    @test abs(mean(chain_X)) < 1 / sqrt(chain_length)
    @test abs(std(chain_X) - 1.0) < 1 / sqrt(chain_length)
end

max_stars_repo_name: bischtob/SEC

max_issues_repo_name: bischtob/SEC

max_forks_repo_name: bischtob/SEC


text:
"""
    ILP(optimizer, max_iter)

ILP iteratively solves a linearized primal optimization to compute maximum allowable disturbance.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: hyperrectangle
3. Output: halfspace

# Return
`AdversarialResult`

# Method
Iteratively solve a linear encoding of the problem.
Default `optimizer` is `GLPKSolverMIP()`. Default `max_iter` is `10`.

# Property
Sound but not complete.

# Reference
[O. Bastani, Y. Ioannou, L. Lampropoulos, D. Vytiniotis, A. Nori, and A. Criminisi,
"Measuring Neural Net Robustness with Constraints,"
in *Advances in Neural Information Processing Systems*, 2016.](https://arxiv.org/abs/1605.07262)
"""
@with_kw struct ILP{O<:AbstractMathProgSolver}
    optimizer::O    = GLPKSolverMIP()
    max_iter::Int64 = 10
end

function solve(solver::ILP, problem::Problem)
    x = problem.input.center
    i = 0
    while i < solver.max_iter
        model = Model(solver = solver.optimizer)
        act_pattern = get_activation(problem.network, x)

        neurons = init_neurons(model, problem.network)
        add_complementary_set_constraint!(model, problem.output, last(neurons))
        encode_relaxed_lp!(model, problem.network, neurons, act_pattern)
        o = max_disturbance!(model, first(neurons) - problem.input.center)

        status = solve(model, suppress_warnings = true)
        if status != :Optimal
            return AdversarialResult(:Unknown)
        end
        x = getvalue(first(neurons))
        if match_activation(problem.network, x, act_pattern)
            radius = getvalue(o)
            if radius >= minimum(problem.input.radius)
                return AdversarialResult(:SAT, radius)
            else
                return AdversarialResult(:UNSAT, radius)
            end
        end
        i += 1
    end
    return AdversarialResult(:Unknown)
end

function match_activation(nnet::Network, x::Vector{Float64}, act_pattern::Vector{Vector{Bool}})
    curr_value = x
    for (i, layer) in enumerate(nnet.layers)
        curr_value = layer.weights * curr_value + layer.bias
        for (j, val) in enumerate(curr_value)
            act = act_pattern[i][j]
            act  && val < 0.0 && return false
            !act && val > 0.0 && return false
        end
        curr_value = layer.activation(curr_value)
    end
    return true
end


text:
export CUSPARSEError

struct CUSPARSEError <: Exception
    code::cusparseStatus_t
    msg::AbstractString
end
Base.show(io::IO, err::CUSPARSEError) = print(io, "CUSPARSEError(code $(err.code), $(err.msg))")

function CUSPARSEError(code::cusparseStatus_t)
    msg = status_message(code)
    return CUSPARSEError(code, msg)
end

function status_message( status )
    if status == CUSPARSE_STATUS_SUCCESS
        return "success"
    end
    if status == CUSPARSE_STATUS_NOT_INITIALIZED
        return "not initialized"
    end
    if status == CUSPARSE_STATUS_ALLOC_FAILED
        return "allocation failed"
    end
    if status == CUSPARSE_STATUS_INVALID_VALUE
        return "invalid value"
    end
    if status == CUSPARSE_STATUS_ARCH_MISMATCH
        return "architecture mismatch"
    end
    if status == CUSPARSE_STATUS_MAPPING_ERROR
        return "mapping error"
    end
    if status == CUSPARSE_STATUS_EXECUTION_FAILED
        return "execution failed"
    end
    if status == CUSPARSE_STATUS_INTERNAL_ERROR
        return "internal error"
    end
    if status == CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED
        return "matrix type not supported"
    end
end

macro check(sparse_func)
    quote
        local err = $(esc(sparse_func::Expr))
        if err != CUSPARSE_STATUS_SUCCESS
            throw(CUSPARSEError(cusparseStatus_t(err)))
        end
        err
    end
end


text:
@testset "get_volumes.jl" begin

end


text:
#############################################################################
# relative_entropy.jl
# relative entropy (ie, sum_i( x_i log (x_i/y_i) ) of expressions x and y
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################

export relative_entropy, log_perspective
export sign, curvature, monotonicity, evaluate

# TODO: make this work for a *list* of inputs, rather than just for scalar/vector/matrix inputs

struct RelativeEntropyAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr,AbstractExpr}
    size::Tuple{Int, Int}

    function RelativeEntropyAtom(x::AbstractExpr, y::AbstractExpr)
        if sign(x) == ComplexSign() || sign(y) == ComplexSign()
            error("Both the arguments should be real but these are instead $(sign(x)) and $(sign(y))")
        else
            children = (x, y)
            return new(:entropy, hash(children), children, size(x))
        end
    end
end

function sign(x::RelativeEntropyAtom)
    return NoSign()
end

function monotonicity(x::RelativeEntropyAtom)
    return (NoMonotonicity(),NoMonotonicity())
end

function curvature(x::RelativeEntropyAtom)
    return ConvexVexity()
end

function evaluate(e::RelativeEntropyAtom)
    x = evaluate(e.children[1])
    y = evaluate(e.children[2])
    if any(isnan, y) return Inf end

    out = x.*log.(x./y)
    # fix value when x=0:
    # out will only be NaN if x=0, in which case the correct value is 0
    out[isnan.(out)] = 0
    return out
end

function conic_form!(e::RelativeEntropyAtom, unique_conic_forms::UniqueConicForms=UniqueConicForms())
    if !has_conic_form(unique_conic_forms, e)
        # transform to conic form:
        # x log x/y <= z
        # x log y/x >= -z
        # log y/x >= -z/x
        # y/x >= exp(-z/x)
        # y >= x exp(-z/x)
        # and cf the standard form for the exponential cone {(x,y,z): y*exp(x/y) <= z}
        z = Variable(e.size)
        x = e.children[1]
        y = e.children[2]
        objective = conic_form!(z, unique_conic_forms)
        for i=1:size(x,1)
            for j=1:size(x,2)
                conic_form!(ExpConstraint(-z[i,j], x[i,j], y[i,j]), unique_conic_forms)
            end
        end
        # need to constrain x>=0 and y>0.
        # x>=0 we get for free from the form of the exponential cone, so just add
        conic_form!(y>=0, unique_conic_forms) # nb we don't know how to ask for strict inequality
        cache_conic_form!(unique_conic_forms, e, objective)
    end
    return get_conic_form(unique_conic_forms, e)
end

relative_entropy(x::AbstractExpr, y::AbstractExpr) = sum(RelativeEntropyAtom(x, y))
# y*log(x/y)
log_perspective(x::AbstractExpr, y::AbstractExpr) = -relative_entropy(y, x)


text:
module ToyStencilsTests

include("test_trapezoid.jl")

end  # module ToyStencilsTests


text:
function test_NfOrd()
  print("Testing NfOrd ...")

  Qx, x = PolynomialRing(FlintQQ, "x")

  K1, a1 = NumberField(x^3 - 2, "a")
  O1 = EquationOrder(K1)

  K2, a2 = NumberField(x - 2, "a")
  O2 = EquationOrder(K2)

  f3 = x^64 - 64*x^62 + 1952*x^60 - 37760*x^58 + 520144*x^56 - 5430656*x^54 + 44662464*x^52 - 296854272*x^50 + 1623421800*x^48 - 7398867840*x^46 + 28362326720*x^44 - 92043777280*x^42 + 254005423840*x^40 - 597659820800*x^38 + 1200442440064*x^36 - 2057901325824*x^34 + 3006465218196*x^32 - 3732682723968*x^30 + 3922021702720*x^28 - 3467892873984*x^26 + 2561511781920*x^24 - 1565841089280*x^22 + 782920544640*x^20 - 315492902400*x^18 + 100563362640*x^16 - 24754058496*x^14 + 4559958144*x^12 - 602516992*x^10 + 53796160*x^8 - 2968064*x^6 + 87296*x^4 - 1024*x^2 + 2

  K3, a3 = NumberField(f3, "a")
  O3 = Order(nf_elem[ a3^i for i in 0:63])

  K4, a4 = NumberField(x^2 - 5, "a")
  O4 = Order(K4, Hecke.FakeFmpqMat(ZZ[1 0; 0 2], fmpz(1)))

  K6, a6 = NumberField(x^2 - 180, "a")
  O6 = EquationOrder(K6)

  O7 = Order(K6, Hecke.FakeFmpqMat(ZZ[6 0; 0 1], ZZ(6)))

  # deepcopy

  O5 = @inferred deepcopy(O2)

  @test nf(O2) == nf(O5)

  # field access

  b = @inferred parent(O2)
  @test b == @inferred parent(O5)

  @test K1 == @inferred nf(O1)
  @test K2 == @inferred nf(O2)
  @test K3 == @inferred nf(O3)

  @test @inferred isequationorder(O1)
  @test @inferred isequationorder(O2)
  @test @inferred !isequationorder(O3)
  @test @inferred !isequationorder(O4)
  @test @inferred isequationorder(O5)

  b = @inferred basis(O1)
  @test b == [ O1(1), O1(a1), O1(a1^2) ]

  b = @inferred basis(O2)
  @test b == [ O2(1) ]

  b = @inferred basis(O3)
  @test b == [ O3(a3^i) for i in 0:63]

  b = @inferred basis(O4)
  @test b == [ O4(1), O4(2*a4) ]

  @test O1.basis_nf == [ a1^0, a1, a1^2]

  @test O2.basis_nf == [ K2(1) ]

  @test O3.basis_nf == [ a3^i for i in 0:63]

  @test O4.basis_nf == [ a4^0, 2*a4 ]

  b = @inferred basis_mat(O1)
  @test b == Hecke.FakeFmpqMat(one(MatrixSpace(FlintZZ, 3, 3)), one(FlintZZ))

  b = @inferred basis_mat(O2)
  @test b == Hecke.FakeFmpqMat(one(MatrixSpace(FlintZZ, 1, 1)), one(FlintZZ))

  b = @inferred basis_mat(O3)
  @test b == Hecke.FakeFmpqMat(one(MatrixSpace(FlintZZ, 64, 64)), one(FlintZZ))

  b = @inferred basis_mat(O4)
  @test b == Hecke.FakeFmpqMat(ZZ[1 0; 0 2], one(FlintZZ))

  b = @inferred basis_mat_inv(O1)
  @test b == Hecke.FakeFmpqMat(one(MatrixSpace(FlintZZ, 3, 3)), one(FlintZZ))

  b = @inferred basis_mat_inv(O2)
  @test b == Hecke.FakeFmpqMat(one(MatrixSpace(FlintZZ, 1, 1)), one(FlintZZ))

  b = @inferred basis_mat_inv(O3)
  @test b == Hecke.FakeFmpqMat(one(MatrixSpace(FlintZZ, 64, 64)), one(FlintZZ))

  b = @inferred basis_mat_inv(O4)
  @test b == Hecke.FakeFmpqMat(ZZ[2 0; 0 1], ZZ(2))

  # index

  b = @inferred gen_index(O1)
  @test b == 1
  b = @inferred index(O1)
  @test b == 1

  b = @inferred gen_index(O2)
  @test b == 1
  b = @inferred index(O2)
  @test b == 1

  b = @inferred gen_index(O3)
  @test b == 1
  b = @inferred index(O3)
  @test b == 1

  b = @inferred gen_index(O4)
  @test b == QQ(1, 2)
  @test_throws ErrorException index(O4)

  b = @inferred gen_index(O7)
  @test b == QQ(6)
  b = @inferred index(O7)
  @test b == 6

  @test !@inferred is_index_divisor(O1, 2)
  @test !@inferred is_index_divisor(O1, 3)
  @test @inferred is_index_divisor(O7, 2)
  @test @inferred is_index_divisor(O7, fmpz(3))
  @test !@inferred is_index_divisor(O7, 5)

  # discriminant

  b = @inferred discriminant(O1)
  @test b == -108

  b = @inferred discriminant(O2)
  @test b == 1

  b = @inferred discriminant(O3)
  @test b == fmpz(2)^447

  b = @inferred discriminant(O4)
  @test b == fmpz(80)

  @test 3 == @inferred degree(O1)
  @test 1 == @inferred degree(O2)
  @test 64 == @inferred degree(O3)
  @test 2 == @inferred degree(O4)

  @test (1, 1) == @inferred signature(O1)
  @test (1, 0) == @inferred signature(O2)
  @test (64, 0) == @inferred signature(O3)
  @test (2, 0) == @inferred signature(O4)

  # minkowski mat

  RR = ArbField(64)

  b = RR[ RR(1) sqrt(RR(2)) RR(0); (exp(1//RR(3) * log(RR(2)))) (-exp(-1//RR(6) * log(RR(2)))) (sqrt(RR(3)) * exp(-1//RR(6) * log(RR(2)))); (exp(1//RR(3) * log(RR(4)))) (-exp(1//RR(6) * log(RR(2)))) (-exp(1//RR(6) * log(RR(54)))) ]
  bb = @inferred minkowski_mat(O1, 256)

  @test overlaps(b, bb)
  for i in 1:3
    for j in 1:3
      @test Hecke.radiuslttwopower(bb[i, j], -256)
    end
  end

  b = one(MatrixSpace(RR, 1, 1))

  bb = @inferred minkowski_mat(O2, 1024)

  @test overlaps(b, bb)
  for i in 1:1
    for j in 1:1
      @test Hecke.radiuslttwopower(bb[i, j], -1024)
    end
  end

  bb = @inferred minkowski_mat(O3, 1024)

  for i in 1:64
    for j in 1:64
      @test Hecke.radiuslttwopower(bb[i, j], -1024)
    end
  end

  @test contains(RR("19063561108769878656033240617946635072004849200892084525390959717509 +/- 1"), det(bb))

  b = RR[ RR(1) RR(1); -2*sqrt(RR(5)) 2*sqrt(RR(5))]

  bb = @inferred minkowski_mat(O4, 1024)

  @test overlaps(b, bb) ||
  for i in 1:2
    for j in 1:2
      @test Hecke.radiuslttwopower(bb[i, j], -1024)
    end
  end

  # element inclusion

  b = @inferred in(a1, O1)
  @test b

  b = @inferred in(a1//2, O1)
  @test !b

  b = @inferred in(a2, O2)
  @test b

  b = @inferred in(a2//3, O2)
  @test !b

  b = @inferred in(a3^4, O3)
  @test b

  b = @inferred in(a3^3//2, O3)
  @test !b

  b = @inferred in(2*a4, O4)
  @test b

  b = @inferred in(a4, O4)
  @test !b

  # denominator

  b = @inferred den(a1, O1)
  @test b == 1
  b = @inferred den(a1//7, O1)
  @test b == 7

  b = @inferred den(a2, O2)
  @test b == 1
  b = @inferred den(a2//4, O2)
  @test b == 2

  b = @inferred den(a3, O3)
  @test b == 1
  b = @inferred den(a3^3//2, O3)
  @test b == 2

  b = @inferred den(a4, O4)
  @test b == 2
  b = @inferred den(a4//2, O4)
  @test b == 4

  # addition

  O6_2 = Order(K6, Hecke.FakeFmpqMat(ZZ[2 0; 0 1], ZZ(2)))
  O6_3 = Order(K6, Hecke.FakeFmpqMat(ZZ[3 0; 0 1], ZZ(3)))

  b = @inferred O6_2 + O6_3
  @test basis_mat(b) == Hecke.FakeFmpqMat(ZZ[6 0; 0 1], ZZ(6))

  @test discriminant(b) == 20

  @test_throws ErrorException O4 + O4
  @test_throws ErrorException O6_2 + O6_2

  println("PASS")
end


text:
@info "Loading test packages..."
using GMC_NS, Test, Distributions, Serialization
import GMC_NS:box_move, box_reflect, boundary_norm, reflect, r_perturb
import LinearAlgebra:normalize

@info "Beginning tests..."
using Random
Random.seed!()

include("galilean_unit_tests.jl")
include("ensemble_tests.jl")

@info "Tests complete!"


text:
# Load Turing.
using Turing
# Load other dependencies
using Distributions, RDatasets

# Example data set - Fisher's Iris data set
data = dataset("datasets", "iris")
species = data[!, "Species"]
dat = data[!, 1:4]

@model pPCA(x, ::Type{T} = Float64) where {T} = begin

  # Dimensionality of the problem.
  N, D = size(x)

  # latent variable z
  z = Matrix{T}(undef, D, N)
  for d in 1:D
    z[d, :] ~ MvNormal(N, 1.)
  end

  # weights/loadings w
  w = Matrix{T}(undef, D, D)
  for d in 1:D
    w[d, :] ~ MvNormal(D, 1.)
  end

  mu = w * z

  for d in 1:D
    x[:,d] ~ MvNormal(mu[d,:], 1.)
  end

end

ppca = pPCA(dat)

# Hamiltonian Monte Carlo (HMC) sampler parameters
n_iterations = 10000
ϵ = 0.05
τ = 10

chain = sample(ppca, HMC(ϵ, τ), n_iterations)

# Extract paramter estimates for plotting.
w = reshape(mean(group(chain, :w))[:,2], (4,4))
z = reshape(mean(group(chain, :z))[:,2], size(dat))

#  plot(z[1, :], z[2,:])

## Extend to ARD
## Original paper
#
#  #  for (d in 1:D){
    #  w[d] ~ normal(0, sigma * alpha);
  #  }
  #  sigma ~ lognormal(0, 1);
  #  alpha ~ inv_gamma(1, 1);

max_stars_repo_name: leachim/TuringTutorials

max_issues_repo_name: leachim/TuringTutorials

max_forks_repo_name: leachim/TuringTutorials


text:
using Plots
using LinearAlgebra
using MatrixEquations
using CPUTime
using ArraysOfArrays
using Parameters

const M = 5.0  # 車両質量[kg]
const m = 1.0  # 振子質量[kg]
const L = 1.5  # 振子長さ[m]
const l = 0.75
const D = 0.01  # 車両と地面の動摩擦係数
const d = 0.01  # 振子と車両の動摩擦係数
const g = 9.80665  # 重力加速度[m/s²]

function split_vec_of_arrays(u)
    vec.(u) |>
    x ->
    VectorOfSimilarVectors(x).data |>
    transpose |>
    VectorOfSimilarVectors
end

"""ルンゲクッタ法（4次）"""
function solve_RungeKutta(dx, x₀, t_span, Δt, args...)

    t = range(t_span..., step = Δt)  # 時間軸
    x = Vector{typeof(x₀)}(undef, length(t))  # 解を格納する1次元配列

    x[1] = x₀  # 初期値
    for i in 1:length(x)-1
        k₁ = dx(t[i], x[i], args)
        k₂ = dx(t[i]+Δt/2, x[i]+k₁*Δt/2, args)
        k₃ = dx(t[i]+Δt/2, x[i]+k₂*Δt/2, args)
        k₄ = dx(t[i]+Δt, x[i]+k₃*Δt, args)
        x[i+1] = x[i] .+ (k₁ .+ 2k₂ .+ 2k₃ .+ k₄) .* Δt/6
    end

    t, x
end

"""パラメータ"""
@with_kw struct LQRParam{T}
    Q::Matrix{T}
    R::Matrix{T}
end

function dx(t, x, K::Tuple{Matrix{T}}) where T
    """LQRのときの微分方程式"""

    multi = [
        1.0 0.0 0.0 0.0
        0.0 0.0 1.0 0.0
        0.0 M+m 0.0 m*l*cos(x[3])
        0.0 m*l*cos(x[3]) 0.0 4/3*m*l^2
    ]
    offset_x = [
        x[2]
        x[4]
        m*l*x[4]^2*sin(x[3])-D*x[2]
        m*g*l*sin(x[3])-d*x[4]
    ]
    offset_u = [
        0.0
        0.0
        1.0
        0.0
    ]

    Fa = inv(multi) * offset_x
    Fb = inv(multi) * offset_u

    u = -K[1] * x

    dx = Fa .+ Fb.*u
end

function run_simulation(p::LQRParam{T}, Δt::T=0.01, TIME_SPAN::T=5.0) where T
    """LQRで制御"""

    # 原点近傍での線形化システム
    multi = [
        1.0 0.0 0.0 0.0
        0.0 0.0 1.0 0.0
        0.0 M+m 0.0 m*L
        0.0 m*L 0.0 4/3*m*l^2
    ]
    offset_x = [
        0.0 1.0 0.0 0.0
        0.0 0.0 0.0 1.0
        0.0 -D 0.0 0.0
        0.0 0.0 m*g*l -d
    ]
    offset_u = [
        0.0
        0.0
        1.0
        0.0
    ]

    A = inv(multi) * offset_x
    B = inv(multi) * offset_u

    # リカッチ方程式を解く
    S = zero(B)
    P, _, _ = arec(A, B, p.R, p.Q, S)

    K = inv(p.R) * B' * P  # 最適フィードバックゲイン
    println("K = ", K)

    ### 数値シミュレーション ###
    x0 = [0.0, 0.0, π/6, 0.0]
    tspan = (0.0, TIME_SPAN)  # 範囲を設定
    t, x = solve_RungeKutta(dx, x0, tspan, Δt, K)
    x, v, θ, ω = split_vec_of_arrays(x)

    # plot
    #plot([x, v, θ, ω], label=["x" "v" "θ" "ω"])

    anim = @animate for i in 1:2:length(t)
        now = round(i * Δt, digits=1)
        s = "t = " * string(now) * " [s]"
        plot(
            [x[i], x[i]+L*cos(-θ[i]+π/2)],
            [0, L*sin(-θ[i]+π/2)],
            #size=(400,400),
            linewidth=5,
            xlims = (-3L, 3L),
            ylims = (-1.1L, 1.1L),
            aspect_ratio = 1,
            legend=false,
            annotate=(0, -1, s),
        )
    end
    gif(anim, "example_8_jl.gif", fps = 30)
    #mp4(anim, "example_8_jl.mp4", fps = 30)

end

p = LQRParam(
    diagm(0 => [10.0, 10.0, 100.0, 100.0]),  # 重み行列
    ones((1,1)) .* 0.01  # 重み行列
)

@time run_simulation(p)


text:
# Evaluations

@doc """
### evalCurve

Determine the value of curve at certain points

**Inputs**

*  `u` : Array of coordinates where the curve needs to be computed
*  `U` : Knot vector
*  `order` : order of B-spline basis functions, (order = p+1)
*  `P` : Array of control points
*  `C` : Resulting curve values

**Outputs**

*  None

SOURCE: The NURBS book 2nd Edition, Algorithm A3.1
      : Gaetan's pyspline/src/eval_curve
"""->

function evalCurve(u, U, order, P, C)

  @assert length(P) + order == length(U)
  p = order - 1  # Degree of B-spline basis function
  nctl = length(P)
  N = Array{Float64}(order) # Array of basis functions 1D

  for i = 1:length(u)
    span = findSpan(u[i], U, order, nctl)
    basisFunctions(U, order, u[i], span, N)
    C[i] = 0.0
    for j = 1:order
      C[i] += N[j]*P[span-order+j]
    end  # End for j = 1:p+1
  end    # End for i = 1:length(u)

  return nothing
end

@doc """
### evalVolume

Determine the the coordinates of all points in a 3D volume using B-splines. The
symbol convention used in the function is from the book

"The NURBS book 2nd Edition"

**Arguments**

*  `map` : Object of mapping type
*  `Vol` : (x,y,z) coordinates of the embedded volume within the contol points

####Method 2

Applies only to DG meshes so fae

**Input**

* `map`
* `mesh` : Abstract Pumi DG Mesh

**Output**

* `vertices` : Array of updated vertices same shape as mesh.vert_coords

"""->

function evalVolume(map::Mapping, Vol::AbstractArray{Tmsh,4}) where Tmsh

  fill!(Vol, 0.0) # Zero out all entries of Vol

  for k = 1:map.numnodes[3]
    for j = 1:map.numnodes[2]
      for i = 1:map.numnodes[1]
        xyz =sview(Vol, i,j,k,:)
        evalVolumePoint(map, map.xi[i,j,k,:], xyz)
      end  # End for i = 1:map.numnodes[3]
    end    # End for j = 1:map.numnodes[2]
  end      # End for k = 1:map.numnodes[1]

  return nothing
end
#=
function evalVolume{Tffd}(map::PumiMapping{Tffd}, mesh::AbstractDGMesh)

  vertices = zeros(Tffd, size(mesh.vert_coords))
  arr = zeros(Tffd, 3)
  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      fill!(arr,0.0)
      evalVolumePoint(map, map.xi[:,j,i], arr)
      for k = 1:map.ndim
        vertices[k,j,i] = arr[k]
      end
    end
  end

  return vertices
end
=#
#=
@doc """
evalSurface

Evaluate surface points (edge points for 2D) of a a pumi mesh object

**Arguments**

*  `map`  : PumiMapping object
*  `mesh` : Pumi mesh object
*  `sbp`  : Summation-By-Parts object

"""->

function evalSurface{Tffd}(map::PumiMapping{Tffd}, mesh::AbstractCGMesh,
                           sbp::AbstractSBP)

  if map.ndim == 2
    x = zeros(Tffd, 3)
    for itr = 1:length(map.geom_faces)
      geom_face_number = map.geom_faces[itr]
      # get the boundary array associated with the geometric edge
      itr2 = 0
      for itr2 = 1:mesh.numBC
        if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
          break
        end
      end
      start_index = mesh.bndry_offsets[itr2]
      end_index = mesh.bndry_offsets[itr2+1]
      idx_range = start_index:end_index
      bndry_facenums =sview(mesh.bndryfaces, start_index:(end_index - 1))
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices on the boundary face (local face number)
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          fill!(x, 0.0)
          evalVolumePoint(map, map.xi[itr][:,j,i], x)
          mesh.coords[:,vtx_arr[j],bndry_i.element] = x[1:2]
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end  # End for itr = 1:length(map.geom_faces)
  else
    for itr = 1:length(map.geom_faces)
      geom_face_number = map.geom_faces[itr]
      # get the boundary array associated with the geometric edge
      itr2 = 0
      for itr2 = 1:mesh.numBC
        if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
          break
        end
      end
      start_index = mesh.bndry_offsets[itr2]
      end_index = mesh.bndry_offsets[itr2+1]
      idx_range = start_index:end_index
      bndry_facenums = sview(mesh.bndryfaces, start_index:(end_index - 1))
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices on the boundary face (local face number)
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          xyz =sview(mesh.coords,:,vtx_arr[j],bndry_i.element)
          evalVolumePoint(map, map.xi[itr][:,j,i], arr)
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end  # End for itr = 1:length(map.geom_faces)
  end    # End if map.ndim == 2

  return nothing
end
=#

"""
  Evalutes the new locations of the surface points.  Users should update
  `map.cp_xyz` and then call this function.

  **Inputs**

   * map: a PumiMapping object

  **Inputs/Outputs**

   * pts: a mesh.dim x map.numFacePts array to be overwritten with the new
          surface points (ordered according to the numbering supplied to
          the `PumiMapping` constructor.
"""
function evalSurface(map::PumiMapping{Tffd}, pts::Array{Tffd, 2}) where Tffd

  @assert size(pts, 1) == map.ndim
  @assert size(pts, 2) == map.numFacePts

  if map.ndim == 2
    check2DSymmetry(map.cp_xyz)
  end

  x = zeros(Tffd, 3)  # compatability between 2D and 3D

  for i=1:map.numFacePts

    xi_i = sview(map.xi, :, i)
    evalVolumePoint(map, xi_i, x)

    for j=1:map.ndim
      pts[j, i] = x[j]
    end
  end

  return nothing
end

#=
"""
  Computes the surface node coordinates based on the current control point
  positions.  Users should update the values in `map.cp_xyz` and then call this
  function.

  **Inputs**

   * map: the PumiMappping object
   * mesh: the Pumi mesh object

  **Outputs**

   * vertices: An array of arrays holding the new surface coordinates.
               The length of the outer array is the number of geometric
               entities the FFD box encapsulates.  Each inner array is
               `mesh.dim` x `mesh.coord_numNodesPerFace` x the number of faces
               on the geometric entity

  Note this function is currently limited to linear coordinate fields.
"""
function evalSurface{Tffd}(map::PumiMapping{Tffd}, mesh::AbstractDGMesh)

  @assert mesh.coord_order == 1

  nwall_faces = getnWallFaces(mesh, map.bc_nums)
  vertices = Array{Array{Tffd,3}}(length(map.bc_nums))
  defineVertices(mesh, map.bc_nums, vertices)

  # vertices = zeros(Tffd, size(mesh.vert_coords))
  x = zeros(Tffd, 3)
  for (idx, itr) in enumerate(map.bc_nums)
    start_index = mesh.bndry_offsets[itr]
    end_index = mesh.bndry_offsets[itr+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
    nfaces = length(bndry_facenums)
    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      # get the local index of the vertices on the boundary face (local face number)
      # map.xi has the xi coordinates of the face nodes only
      for j = 1:mesh.coord_numNodesPerFace
        fill!(x, 0.0)
        evalVolumePoint(map, map.xi[idx][:,j,i], x)
        for k = 1:map.ndim
          vertices[idx][k, j, i] = x[k]
        end
      end  # End for j = 1:length(vtx_arr)
    end    # End for i = 1:nfaces
  end  # End for itr = 1:length(map.bc_nums)

  return vertices
end
=#
@doc """
### evalVolumePoint

Computes a point in the FFD volume. The symbol convention used in this function
is from "The NURBS book 2nd Edition".

**Arguments**

*  `map` : Object of Mapping type
*  `xi`  : Parametric FFD coordinates (referred to as the (s,t,u) coordinate
           systems within FFD functions and as (u,v,w) in this function)
*  `xyz` : (x,y,z) coordinates of the parametric point in the FFD volume
"""->

function evalVolumePoint(map::AbstractMappingType, xi::AbstractArray{Tffd,1}, xyz) where Tffd

  fill!(xyz, 0.0)

  Nu = zeros(Tffd, map.order[1])
  Nv = zeros(Tffd, map.order[2])
  Nw = zeros(Tffd, map.order[3])

  # Work with u
  span = findSpan(xi[1], map.edge_knot[1], map.order[1], map.nctl[1])
  basisFunctions(map.edge_knot[1], map.order[1], xi[1], span, Nu)
  startu = span - map.order[1]

  # Work with v
  span = findSpan(xi[2], map.edge_knot[2], map.order[2], map.nctl[2])
  basisFunctions(map.edge_knot[2], map.order[2], xi[2], span, Nv)
  startv = span - map.order[2]

  # Work with w
  span = findSpan(xi[3], map.edge_knot[3], map.order[3], map.nctl[3])
  basisFunctions(map.edge_knot[3], map.order[3], xi[3], span, Nw)
  startw = span - map.order[3]

  for ii = 1:map.order[1]
    for jj = 1:map.order[2]
      for kk = 1:map.order[3]
        for idim = 1:3
          xyz[idim] += Nu[ii]*Nv[jj]*Nw[kk]*
                       map.cp_xyz[idim, startu+ii, startv+jj, startw+kk]
        end
      end  # End for kk = 1:map.order[3]
    end    # End for jj = 1:map.order[2]
  end      # End for ii = 1:map.order[1]

  return nothing
end

@doc """
### contractWithdGdB

Used to contract the matrix dG/dB, the derivative of the mesh node coordinates
with respect to the mapping control points, with a given array, here labelled
dJdGrid and stored in (j,k,m,:) format. This can be used to find
dJ/dB = (dJ/dG)(dG/dB) for some objective J. In practice, the objective used is
actually

             J_practice = J_obj + psi^T dot R

where J_obj is the actual objective, psi are the flow adjoints, and R is the
flow residual.

**Arguments**

*  `map`     : Object of mapping type
*  `dJdGrid` : The derivative of the objective w.r.t the grid coordintes

"""->

function contractWithdGdB(map::AbstractMappingType{Tffd}, xi, dJdG) where Tffd

  # Evaluate the knot vectors and basis values
  Nu = zeros(Tffd, map.order[1])
  Nv = zeros(Tffd, map.order[2])
  Nw = zeros(Tffd, map.order[3])
  span = zeros(Int, 3)

  # Work with u
  span[1] = findSpan(xi[1], map.edge_knot[1], map.order[1], map.nctl[1])
  basisFunctions(map.edge_knot[1], map.order[1], xi[1], span[1], Nu)
  startu = span[1] - map.order[1]

  # Work with v
  span[2] = findSpan(xi[2], map.edge_knot[2], map.order[2], map.nctl[2])
  basisFunctions(map.edge_knot[2], map.order[2], xi[2], span[2], Nv)
  startv = span[2] - map.order[2]

  # Work with w
  span[3] = findSpan(xi[3], map.edge_knot[3], map.order[3], map.nctl[3])
  basisFunctions(map.edge_knot[3], map.order[3], xi[3], span[3], Nw)
  startw = span[3] - map.order[3]

  for ii = 1:map.order[1]
    for jj = 1:map.order[2]
      for kk = 1:map.order[3]
        coeff = Nu[ii]*Nv[jj]*Nw[kk]
        for idim = 1:3
          map.work[idim, startu+ii, startv+jj, startw+kk] += coeff*dJdG[idim]
        end
      end  # End for kk = 1:map.order[3]
    end    # End for jj = 1:map.order[2]
  end      # End for ii = 1:map.order[1]

  return nothing
end  # End function contractWithdGdB(map, dJdGrid)

"""
  Evalutes a transposed Jacobian-vector product

  Xcp_bar =  (dXs/dXcp)^T Xs_bar

  where Xs are the surface coordinates, and Xcp are the control point
  coordinates.

  **Inputs**

   * map: a PumiMapping object
   * Xs_bar: the array to multiply the transposed jacobian against,
             size 3 x `map.numFacePts`

  **Inputs/Outputs**

   * Xcp_bar: arrayto be overwritten with results, same size as `map.cp_xyz`
"""
function evaldXdControlPointTransposeProduct(map::PumiMapping{Tffd},
                          Xs_bar::AbstractMatrix, Xcp_bar::AbstractArray{T, N}) where {T, Tffd, N}

  if map.ndim == 2
    @assert N == 3
    @assert size(Xcp_bar, 1) == 2
    @assert size(Xcp_bar, 2) == map.nctl[1]
    @assert size(Xcp_bar, 3) == map.nctl[2]
    check2DSymmetry(map.cp_xyz)

  else
    @assert N == 4
    @assert size(Xcp_bar, 1) == 3  # all meshes are 3 dimensional to FFD
    @assert size(Xcp_bar, 2) == map.nctl[1]
    @assert size(Xcp_bar, 3) == map.nctl[2]
    @assert size(Xcp_bar, 4) == map.nctl[3]

  end
  @assert size(Xs_bar, 1) == map.ndim
  @assert size(Xs_bar, 2) == map.numFacePts

  Xs_bar_i = zeros(Tffd, 3)  # compatability with 2D and 3D
  fill!(map.work, 0.0)
  for i=1:map.numFacePts
    xi_i = sview(map.xi, :, i)
    for j=1:map.ndim
      Xs_bar_i[j] = Xs_bar[j, i]
    end
    contractWithdGdB(map, xi_i, Xs_bar_i)
  end

  work_slice = sview(map.work, 1:3, :, :, :)

  # in 2D, we also need to do the reverse mode of setControlPoints
  if map.ndim == 2
    for i=1:size(Xcp_bar, 3)
      for j=1:size(Xcp_bar, 2)
        for k=1:size(Xcp_bar, 1)
          Xcp_bar[k, j, i] = work_slice[k, j, i, 1] + work_slice[k, j, i, 2]
        end
      end
    end
  else
    copy!(Xcp_bar, work_slice)
  end

  return nothing
end

"""
  Evaluates the Jacobian-vector product

  Xs_dot = (dXs/dXcp) Xcp_dot

  where Xs are the surface coordinates, and Xcp are the control point
  coordinates.

  It is not safe to complex step this function.

  **Inputs**

   * map: a PumiMapping object
   * Xcp_dot: array to multiply the jacobian against, same size as `map.cp_xyz`

  **Inputs/Outputs**

   * Xs_dot: the array to be overwritten with the results
             size 3 x `map.numFacePts`

  **Implementation Notes**

  This function uses complex step internally to do the calculation, using
  `map.map_cs`.  `Tffd` does not have to be a complex type for the `map`
  passed to this function.

"""
function evaldXdControlPointProduct(map::PumiMapping{Tffd},
               Xcp_dot::AbstractArray{T, N}, Xs_dot::AbstractMatrix ) where {T, Tffd, N}

  @assert size(Xs_dot, 1) == map.ndim
  @assert size(Xs_dot, 2) == map.numFacePts
  if map.ndim == 2
    check2DSymmetry(map.cp_xyz)
    @assert N == 3
    @assert size(Xcp_dot, 1) == 2
    for i=2:3
      @assert size(Xcp_dot, i) == size(map.cp_xyz, i)
    end

  else
    @assert N == 4
    for i=1:4
      @assert size(Xcp_dot, i) == size(map.cp_xyz, i)
    end
  end

  # copy things into map.map_cs
  map_cs = map.map_cs
  copy!(map_cs._cp_xyz, map._cp_xyz)
#  copy!(map_cs.xi, map.xi)

  h = 1e-20
  pert = Complex128(0, h)
  if map.ndim == 2
    # update both planes of control points
    for i=1:size(Xcp_dot, 3)
      for j=1:size(Xcp_dot, 2)
        for k=1:size(Xcp_dot, 1)
          map_cs._cp_xyz[k, j, i, 1] += Xcp_dot[k, j, i]*pert
          map_cs._cp_xyz[k, j, i, 2] += Xcp_dot[k, j, i]*pert
        end
      end
    end
  else
    # update all points
    @simd for i=1:length(Xcp_dot)
      map_cs._cp_xyz[i] += Xcp_dot[i]*pert
    end
  end

  vertices = map_cs.vertices
  evalSurface(map_cs, vertices)

  @simd for i=1:length(Xs_dot)
    Xs_dot[i] = imag(vertices[i])/h
  end

  # remove the perturbation
  @simd for i=1:length(Xcp_dot)
    map_cs._cp_xyz[i] = real(map_cs._cp_xyz[i])
  end

  return nothing
end

#=
"""
  evaluates (dXs/dXcp)^T v, where v is a user-supplied vector, Xs are the
  surface grid points and Xcp are the control points (note that Xs is called
  G and Xcp is called B in other parts of the code).

  Specifically, Xs is the non-unqiue representation of the surface vertices
  (it is per-geometric entity).

  **Inputs**

   * map: the PumiMapping object
   * mesh: the Pumi mesh object
   * Xs_bar: the v vector that the transposed jacobian is multiplied against.
             It is an array of arrays, where the outer array has length equal to
             the number of geometric entities the FFD box encapsulates.  Each
             inner array has dimensions `mesh.dim` x `mesh.coord_numNodesPerFace`
             x the number of faces on the geometric entity.

  **Inputs/Outputs**

   * Xcp_bar: array, same shape as map.cp_xyz.  This array will be overwritten
              with the result

  Note: Xs_bar allows duplicate entries (if a vertex is on the intersection of
        two geometric entites), but only one entry can be non-zero
        in order to get the correct result
"""
function evaldXdControlPointTransposeProduct{Tmsh}(map::PumiMapping, mesh::AbstractDGMesh{Tmsh}, Xs_bar::Array{Array{Complex128, 3}, 1}, Xcp_bar::AbstractArray{Complex128, 4})

  @assert size(Xcp_bar, 1) == 3  # all meshes are 3 dimensional to FFD
  @assert size(Xcp_bar, 2) == map.nctl[1]
  @assert size(Xcp_bar, 3) == map.nctl[2]
  @assert size(Xcp_bar, 4) == map.nctl[3]

  @assert length(Xs_bar) == length(map.bc_nums)
  for i=1:length(Xs_bar)
    @assert size(Xs_bar[i], 1) == mesh.dim
    @assert size(Xs_bar[i], 2) == mesh.coord_numNodesPerFace
    @assert size(Xs_bar[i], 3) == size(map.xi[i], 3)
  end

  fill!(map.work, 0.0)  # prepare for accumulation
  for (idx, itr) in enumerate(map.bc_nums)
    start_index = mesh.bndry_offsets[itr]
    end_index = mesh.bndry_offsets[itr+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
    nfaces = length(bndry_facenums)
    Xs_bar_itr = Xs_bar[idx]
    Xs_bar_j = zeros(Tmsh, 3)  # vector of length 3 gives compatability between
                               # 2D and 3D

    for i=1:nfaces
      for j=1:mesh.coord_numNodesPerFace
        # copy value into Xs_bar_j
        for k=1:mesh.dim
          Xs_bar_j[k] = Xs_bar_itr[k, j, i]
        end
        contractWithdGdB(map, map.xi[idx][:, j, i], Xs_bar_j)

        # contractWithdGdB(map, map.xi[itr][:,j,i], dJdVert_arr[:,ctr])
      end  # end loop j
    end  # end loop i
  end  # end loop itr

  # copy result into user-provided array
  # was 1:mesh.dim
  work_slice = sview(map.work, 1:3, :, :, :)
  copy!(Xcp_bar, work_slice)

  return nothing
end
=#

#=
"""
### evaldXdControlPointProduct

Higher level function that is used for computing the derivative

  ∂J/∂B = (∂J/∂G)(∂G/∂B)

for pumi meshes. Here J denotes the functional, G denotes the grid, and B denotes
the control point.

**Arguments**

* `map` : Object of PumiMapping type
* `mesh` : Object of AbstractDGMesh type
* `dJdVert` : derivative of the a functional w.r.t. the pumi mesh vertices.

"""
function evaldXdControlPointProduct(map::PumiMapping, mesh::AbstractDGMesh,
                                    dJdVert::AbstractArray{Float64,1})

  my_rank = MPI.Comm_rank(MPI.COMM_WORLD)

  fill!(map.work, 0.0)
  ctr = 1
  local_vertnum_history = Int[]
  dJdVert_arr = reshape(dJdVert, 3, convert(Int,length(dJdVert)/3))

  # Gather nface_verts from all ranks
  nlocal_face_verts = getLocalNumFaceVerts_unique(mesh, map.geom_faces)
  nlocal_face_verts_arr = MPI.Allgather(nlocal_face_verts, MPI.COMM_WORLD)
  ranks_for_bndry_faces = Int[]
  for i = 1:length(nlocal_face_verts_arr)
    if nlocal_face_verts_arr[i] > 0
      push!(ranks_for_bndry_faces, i-1)
    end # End if nface_verts_arr[i] > 0
  end

  for itr = 1:length(map.geom_faces)
    geom_face_number = map.geom_faces[itr]
    # get the boundary array associated with the geometric edge
    itr2 = 0
    for itr2 = 1:mesh.numBC
      if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
        break
      end
    end
    start_index = mesh.bndry_offsets[itr2]
    end_index = mesh.bndry_offsets[itr2+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums =sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
    nfaces = length(bndry_facenums)

    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
      for j = 1:length(vtx_arr)
        local_vertnum = mesh.element_vertnums[vtx_arr[j],bndry_i.element]
        if findfirst(local_vertnum_history, local_vertnum) == 0
          # Check if this local vertex exists on other ranks
          if haskey(mesh.vert_sharing.rev_mapping, local_vertnum)
            # Check which all MPI ranks have this vertex and its position in
            # that rank's mesh.vert_sharing.
            rank_arr = first(mesh.vert_sharing.rev_mapping[local_vertnum])
            localIdx_arr = last(mesh.vert_sharing.rev_mapping[local_vertnum])
            # Lower rank always has ownership rights
            # if my_rank < minimum(rank_arr)
            #   contractWithdGdB(map, map.xi[itr][:,j,i], dJdVert_arr[:,ctr])
            #   ctr += 1
            # else
            #   mpi_neighbor_rank = getNeighborRank(mesh, bndry_i)
            #   if my_rank < mpi_neighbor_rank
            #     contractWithdGdB(map, map.xi[itr][:,j,i], dJdVert_arr[:,ctr])
            #     ctr += 1
            #   end
            # end # End if my_rank < minimum(rank_arr)
            intersect_arr = intersect(rank_arr, ranks_for_bndry_faces)
            if intersect_arr == []
              contractWithdGdB(map, map.xi[itr][:,j,i], dJdVert_arr[:,ctr])
              ctr += 1
            elseif my_rank < minimum(intersect_arr)
              contractWithdGdB(map, map.xi[itr][:,j,i], dJdVert_arr[:,ctr])
              ctr += 1
            end # End if intersect_arr == []
          else # The vertex exist only on one MPI rank
            contractWithdGdB(map, map.xi[itr][:,j,i], dJdVert_arr[:,ctr])
            ctr += 1
          end # End if haskey
        end
        push!(local_vertnum_history, local_vertnum)
      end  # End for j = 1:length(vtx_arr)
    end    # End for i = 1:nfaces

    # Now do an All reduce operation on map.work
    recv_arr = zeros(map.work)
    MPI.Allreduce!(map.work, recv_arr, MPI.SUM, MPI.COMM_WORLD)
    for i = 1:length(map.work)
      map.work[i] = recv_arr[i]
    end

  end  # End for itr = 1:length(map.geom_faces)

  return nothing
end
=#


text:
using MultipathChannel
using Plots
using Random

Random.seed!(1)

# TDL-C channel
samplingrate = 15.36e6
rmsdelay = 300e-9
initialtime = 0.0
maxdoppler = 111.1
channelmodel = "TDLC"
rmsdelay = 300e-9
chan = initChan(maxdoppler, samplingrate, channelmodel, rmsdelay=rmsdelay)

# Impulse
impulse = zeros(Float64, 10000)
impulse[1] = 1.0
response = applyChan(impulse, chan)
plot(1:50, impulse[1:50], line=:stem, legend=false)
plot!(1:50, response[1:50], line=:stem, legend=false)
xlabel!("Delay (in sample)"); ylabel!("Envelope")
savefig("TDLC_impulse.png")

# Square wave
sqwave = repeat([1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0]', 20)
response = applyChan(sqwave, chan)
plot(1:50, sqwave[1:50], legend=false)
plot!(1:50, real(response[1:50]), legend=false)
xlabel!("Delay (in sample)"); ylabel!("Envelope")
savefig("TDLC_square_wave.png")

# Sine wave
t = 0:0.05:50
sinewave = sin.(2π*t)
response = applyChan(sinewave, chan)
plot(1:50, real(sinewave[1:50]), legend=false)
plot!(1:50, real(response[1:50]), legend=false)
xlabel!("Delay (in sample)"); ylabel!("Envelope")
savefig("TDLC_sine_wave.png")