text:
# Script to run the complete Benders algorithm (divide and conquer, LO relaxation phase and integer phase) -- see
# Section 5.3 of Bertsimas and Mišić.
# This script produces the numbers needed for Table 7 of the paper.
#
# This script does not take any command line arguments. To run it,
# type
#
# > julia Toubia_Benders_DNCWarmStart.jl
#
# in the terminal.

using JuMP, Gurobi

include("../optimalPLD_code/optimalPLD_createPath.jl")
include("../optimalPLD_code/optimalPLD_solveSAO_Benders.jl")

include("../optimalPLD_code/optimalPLD_runDNCHeuristicObjFn.jl")

include("../optimalPLD_code/optimalPLD_solveSAO_LP_Benders.jl")
include("../optimalPLD_code/optimalPLD_createRanking_rfn.jl")

srand(1683)
machine = 2;

n_vec = [3584]
numPerm_vec = [330]

numReps = 1;

counter = 0;

data_prefix = "toubia2003"

expdirpath = optimalPLD_createPath(machine, "exp_dir", "Toubia_Benders_DNCWarmStart");
mkdir(expdirpath)

numInst = 8;
inst_vec = [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8]
b_feas_vec = [1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 10, -10, 20, -20, 50, -50]
n = 3584;
A_feas_mat = [ones(1,n);
			 -ones(1,n);
			 ones(1,n);
			 -ones(1,n);
			 ones(1,n);
			 -ones(1,n);
			 ones(1,n);
			 -ones(1,n);
			 ones(1,n);
			 -ones(1,n);
			 ones(1,n);
			 -ones(1,n);
			 ones(1,n);
			 -ones(1,n);
			 ones(1,n);
			 -ones(1,n)]

text_description_by_inst = ["sum x_i = 1", "sum x_i = 2", "sum x_i = 3",
							"sum x_i = 4", "sum x_i = 5", "sum x_i = 10",
							"sum x_i = 20", "sum x_i = 50"];

max_removals = Inf;
MIP_timeLimit = 6*3600;
isRelaxation = false;

outcsvfilepath = string(expdirpath, "out.csv");
outcsvhandle = open(outcsvfilepath, "w")
row_string = string("n", ",",
					"numPermutations", ",",
					"rep", ",",
					"instance", ",",
					"text_description", ",",
					"DNC_expRevenue", ",",
					"DNC_solutiontime", ",",
					"LPB_expRevenue", ",",
					"LPB_solutiontime", ",",
					"LPB_numIterations", ",",
					"LPB_numCutsGenerated", ",",
					"MIPB_expRevenue", ",",
					"MIPB_upperBound", ",",
					"MIPB_solutiontime", ",",
					"MIPB_gap", ",",
					"MIPB_total_num_lazy_checked", ",",
					"MIPB_total_num_lazy_generated")
row_string = string(row_string, "\n");

# print_unescaped(outcsvhandle, row_string);
print(outcsvhandle, row_string);
flush(outcsvhandle);

counter = 0;
for nind = 1:length(n_vec)
	n = n_vec[nind]
	for nPind = 1:length(numPerm_vec)
		numPermutations = numPerm_vec[nPind];

		println("**** n = $n, numPermutations = $numPermutations ****")

		println("**** Reading from data directory...")
		data_name = string(data_prefix,"_neq",n,"_Keq",numPermutations, "_v2")
		data_dir = optimalPLD_createPath(machine, "data_dir", data_name)
		revenues_mat = readcsv(string(data_dir, "revenues_mat.csv"))
		lambda_mat = readcsv(string(data_dir, "lambda_mat.csv"))
		orderings_mat = readcsv(string(data_dir, "orderings_mat.csv"))

		for instance=1:numInst
			rep = 1;
			# numPermutations = 10;
			revenues = revenues_mat[rep,:][:]
			lambda = lambda_mat[rep,:][1:numPermutations]
			# numPermutations = numPerm_vec[];
			# lambda = lambda_mat[rep,1:numPermutations];
			orderings = orderings_mat[numPermutations*(rep-1)+1:rep*numPermutations,:]

			relinds = find(inst_vec .== instance);
			@show relinds
			A_feas = A_feas_mat[relinds,:];
			b_feas = b_feas_vec[relinds];
			text_description = text_description_by_inst[instance];

			println("\t Creating ranking revenue function...")
			revenue_fn = optimalPLD_createRanking_rfn(n, [revenues;0], numPermutations, orderings, lambda)

			println("**** Running DNC with 10 random starting points ****")
			DNC_bestHeuristicRevenue, DNC_bestHeuristicSoln, DNC_heuristicTime = optimalPLD_runDNCHeuristicObjFn( b_feas[1], revenue_fn, 3584, 10, 1492 )
			# error("Stop here.")
			DNC_bestHeuristicSoln = convert(Array{Int64}, DNC_bestHeuristicSoln)

			x_initial = zeros(Float64, n);
			x_initial[DNC_bestHeuristicSoln[:]] = 1;

			println("**** Solving n = $n, nP = $numPermutations, rep = $rep ****")
			counter = counter + 1;

			# isRelaxation = true;
			# x_val, y_val, expRevenue, solutiontime, MIPgap = optimalPLD_solveSAO( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation )

			# x_initial = [];

			# A_feas = ones(1,n);
			# b_feas = [5];

			# x_initial_LP = [];
			best_x_val, LPB_expRevenue, LPB_solutiontime,
				rankingid_list, alpha_list, beta_list, gamma_list, LPB_numIterations = optimalPLD_solveSAO_LP_Benders( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, x_initial, A_feas, b_feas)

			# rankingid_list = [];
			# alpha_list = [];
			# beta_list =[];
			# gamma_list = [];

			writecsv( string(expdirpath, "rankingid_list_instance", instance, ".csv"), rankingid_list)
			writecsv( string(expdirpath, "alpha_list_instance", instance, ".csv"), alpha_list)
			writecsv( string(expdirpath, "beta_list_instance", instance, ".csv"), beta_list)
			writecsv( string(expdirpath, "gamma_list_instance", instance, ".csv"), gamma_list)

			LPB_numCutsGenerated = length(rankingid_list);

			x_val, MIPB_expRevenue, MIPB_upperBound, MIPB_solutiontime, MIPB_gap, MIPB_total_num_lazy_checked, MIPB_total_num_lazy_generated,  = optimalPLD_solveSAO_Benders( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial,
																			A_feas, b_feas,
																			rankingid_list, alpha_list, beta_list, gamma_list)

			@show MIPB_total_num_lazy_checked
			@show MIPB_total_num_lazy_generated
			# error("Stop here. Examine num lazy gen'd.")

			row_string = string(n, ",",
								numPermutations, ",",
								rep, ",",
								instance, ",",
								text_description, ",",
								DNC_bestHeuristicRevenue, ",",
								DNC_heuristicTime, ",",
								LPB_expRevenue, ",",
								LPB_solutiontime, ",",
								LPB_numIterations, ",",
								LPB_numCutsGenerated, ",",
								MIPB_expRevenue, ",",
								MIPB_upperBound, ",",
								MIPB_solutiontime, ",",
								MIPB_gap, ",",
								MIPB_total_num_lazy_checked, ",",
								MIPB_total_num_lazy_generated)
			row_string = string(row_string, "\n");
			# print_unescaped(outcsvhandle, row_string);
			print(outcsvhandle, row_string);
			flush(outcsvhandle);

		end
	end
end

close(outcsvhandle)

max_stars_repo_name: vvmisic/EFCPLO

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: vvmisic/EFCPLO

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: vvmisic/EFCPLO

max_forks_repo_licenses: MIT Unlicense


text:
@unix_only begin
require("testdefs.jl")

pipe_fds = Array(Cint,2)
@test 0 == ccall(:pipe, Cint, (Ptr{Cint},), pipe_fds)

function test_poll_fd(timeout_s)
    rc = poll_fd(RawFD(pipe_fds[1]), timeout_s; readable=true)
    produce(rc)
end

function test_timeout_fd(tval)
    tic()
    t = Task(()->test_poll_fd(tval))
    tr = consume(t)
    t_elapsed = toq()

    @test tr.timedout

    @test tval <= t_elapsed
end

function test_read(slval)
    tval = slval * 1.1
    tic()
    t = Task(()->test_poll_fd(tval))

    sleep(slval)
    @test 1 == ccall(:write, Csize_t, (Cint, Ptr{Uint8},Csize_t), pipe_fds[2], bytestring("A"), 1)

    tr = consume(t)
    t_elapsed = toq()

    @test isreadable(tr) || iswritable(tr)

    dout = Array(Uint8, 1)
    @test 1 == ccall(:read, Csize_t, (Cint, Ptr{Uint8},Csize_t), pipe_fds[1], dout, 1)
    @test dout[1] == int8('A')

    @test slval <= t_elapsed
end

test_timeout_fd(.1)
test_timeout_fd(1)
test_read(.1)
test_read(1)

ccall(:close, Cint, (Cint,), pipe_fds[1])
ccall(:close, Cint, (Cint,), pipe_fds[2])

end

ext: jl

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
#! /usr/bin/julia

# Rosetta Code, Linear congruential generator

function lcg_maker{T<:Integer}(r::T, a::T, c::T, m::T, sh::T)
    state = r
    function lcg_rand()
        state = mod(a*state + c, m)
        return state >> sh
    end
    return lcg_rand
end

snum = 10
seed = 0
bsd_rand = lcg_maker(seed, 1103515245, 12345, 2^31, 0)

print("The first ", snum, " results for a BSD rand() seeded with ")
println(seed, ":")

for i in 1:snum
    println(@sprintf "%14d" bsd_rand())
end

ms_rand = lcg_maker(seed, 214013, 2531011, 2^31, 16)

println()
print("The first ", snum, " results for a M\$ rand() seeded with ")
println(seed, ":")

for i in 1:snum
    println(@sprintf "%14d" ms_rand())
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module LinsolveAD
    using KrylovKit, LinearAlgebra
    using Random, Test
    using ChainRulesCore, ChainRulesTestUtils, Zygote, FiniteDifferences

    fdm = ChainRulesTestUtils._fdm
    precision(T::Type{<:Number}) = eps(real(T))^(2/3)
    n = 10
    N = 30

    function build_mat_example(A, b; tol = precision(eltype(A)), kwargs...)
        Avec, A_fromvec = to_vec(A)
        bvec, b_fromvec = to_vec(b)
        T = eltype(A)

        function mat_example(Av, bv)
            A′ = A_fromvec(Av)
            b′ = b_fromvec(bv)
            x, info = linsolve(A′, b′, zero(b′), GMRES(; tol = tol, kwargs...))
            info.converged == 0 && @warn "linsolve did not converge"
            xv, = to_vec(x)
            return xv
        end
        return mat_example, Avec, bvec
    end

    function build_fun_example(A, b, c, d, e, f; tol = precision(eltype(A)), kwargs...)
        Avec, matfromvec = to_vec(A)
        bvec, vecfromvec = to_vec(b)
        cvec, = to_vec(c)
        dvec, = to_vec(d)
        evec, scalarfromvec = to_vec(e)
        fvec, = to_vec(f)

        function fun_example(Av, bv, cv, dv, ev, fv)
            A′ = matfromvec(Av)
            b′ = vecfromvec(bv)
            c′ = vecfromvec(cv)
            d′ = vecfromvec(dv)
            e′ = scalarfromvec(ev)
            f′ = scalarfromvec(fv)

            x, info = linsolve(b′, zero(b′), GMRES(; tol = tol, kwargs...), e′, f′) do y
                A′*y + c′ * dot(d′, y)
            end
            # info.converged > 0 || @warn "not converged"
            xv, = to_vec(x)
            return xv
        end
        return fun_example, Avec, bvec, cvec, dvec, evec, fvec
    end

    @testset "Small linsolve AD test" begin
        @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
            A = 2*(rand(T, (n,n)) .- one(T)/2)
            b = 2*(rand(T, n) .- one(T)/2)
            b /= norm(b)

            mat_example, Avec, bvec = build_mat_example(A, b; tol = cond(A)*eps(real(T)), krylovdim = n, maxiter = 1)

            (JA, Jb) = FiniteDifferences.jacobian(fdm, mat_example, Avec, bvec)
            (JA′, Jb′) = Zygote.jacobian(mat_example, Avec, bvec)
            @test JA ≈ JA′ rtol=cond(A)*precision(T)
            @test Jb ≈ Jb′ rtol=cond(A)*precision(T)
        end
    end

    @testset "Large linsolve AD test" begin
        for T in (Float64, ComplexF64)
            A = rand(T,(N,N)) .- one(T)/2
            A = I-(9//10)*A/maximum(abs, eigvals(A))
            b = 2*(rand(T, N) .- one(T)/2)
            c = 2*(rand(T, N) .- one(T)/2)
            d = 2*(rand(T, N) .- one(T)/2)
            e = rand(T)
            f = rand(T)

            fun_example, Avec, bvec, cvec, dvec, evec, fvec = build_fun_example(A, b, c, d, e, f; tol = precision(T), krylovdim = 20)

            (JA, Jb, Jc, Jd, Je, Jf) = FiniteDifferences.jacobian(fdm, fun_example,
                Avec, bvec, cvec, dvec, evec, fvec)
            (JA′, Jb′, Jc′, Jd′, Je′, Jf′) = Zygote.jacobian(fun_example, Avec, bvec, cvec, dvec, evec, fvec)
            @test JA ≈ JA′
            @test Jb ≈ Jb′
            @test Jc ≈ Jc′
            @test Jd ≈ Jd′
            @test Je ≈ Je′
            @test Jf ≈ Jf′
        end
    end
end


text:
#=
The main module for Femshop.
=#
module Femshop

# Public macros and functions
export @generateFor, @domain, @mesh, @solver, @stepper, @setSteps, @functionSpace, @trialFunction, @matrixFree,
        @testFunction, @nodes, @order, @boundary, @addBoundaryID, @referencePoint, @variable, @coefficient, @parameter, @testSymbol, @initial,
        @timeInterval, @weakForm, @fluxAndSource, @LHS, @RHS, @exportCode, @importCode,
        @customOperator, @customOperatorFile,
        @outputMesh, @useLog, @finalize
export init_femshop, set_language, set_custom_gen_target, dendro, set_solver, set_stepper, set_specified_steps, set_matrix_free,
        reformat_for_stepper, reformat_for_stepper_fv,
        add_mesh, output_mesh, add_boundary_ID, add_test_function,
        add_initial_condition, add_boundary_condition, add_reference_point, set_rhs, set_lhs, set_lhs_surface, set_rhs_surface,
        set_symexpressions, set_assembly_loops
export build_cache_level, build_cache, build_cache_auto
export sp_parse
export generate_code_layer, generate_code_layer_surface, generate_code_layer_fv
export Variable, add_variable
export Coefficient, add_coefficient
export Parameter, add_parameter
export FVInfo

### Module's global variables ###
# configuration, output files, etc.
config = nothing;
prob = nothing;
project_name = "unnamedProject";
output_dir = pwd();
language = 0;
gen_framework = 0;
generate_external = false;
solver = nothing;
codegen_params = nothing;
#log
use_log = false;
log_file = "";
log_line_index = 1;
#mesh
mesh_data = nothing;    # The basic element information as read from a MSH file or generated here.
grid_data = nothing;    # The full collection of nodes(including internal nodes) and other mesh info in the actual DOF ordering.
geo_factors = nothing;  # Geometric factors
fv_info = nothing;      # Finite volume info
refel = nothing;        # Reference element (will eventually be an array of refels)
fv_refel = nothing;        # Reference element for FV
elemental_order = [];   # Determines the order of the elemental loops
parent_maps = nothing;  # For high order FV
fv_grid = nothing;      # For high order FV
fv_geo_factors = nothing;# Geometric factors for FV
#problem variables
var_count = 0;
variables = [];
coefficients = [];
parameters = [];
test_functions = [];
indexers = [];
#generated functions
genfunc_count = 0;
genfunctions = [];
#rhs
linears = [];
face_linears = [];
#lhs
bilinears = [];
face_bilinears = [];
#assembly loop functions
assembly_loops = [];
#symbolic layer
symexpressions = [[],[],[],[]];
# string versions of generated code for printing
code_strings = [[],[],[],[]];
#time stepper
time_stepper = nothing;
specified_dt = 0;
specified_Nsteps = 0;
use_specified_steps = false;

use_cachesim = false;

#handles for custom code gen functions
custom_gen_funcs = [];

###############################################################
include("femshop_includes.jl");
include("macros.jl");
include("femshop_interface.jl"); # included here after globals are defined
###############################################################

config = Femshop_config(); # These need to be initialized here
prob = Femshop_prob();

# This has already been initialized on loading, so this is more of a reset function
function init_femshop(name="unnamedProject")
    global config = Femshop_config();
    global prob = Femshop_prob();
    global project_name = name;
    global language = JULIA;
    global framework = 0;
    global generate_external = false;
    global codegen_params = nothing;
    global log_file = "";
    global use_log = false;
    global log_line_index = 1;
    global mesh_data = nothing;
    global grid_data = nothing;
    global geo_factors = nothing;
    global fv_info = nothing;
    global refel = nothing;
    global fv_refel = nothing;
    global elemental_order = [];
    global parent_maps = nothing;
    global fv_grid = nothing;
    global fv_geo_factors = nothing;
    global var_count = 0;
    global variables = [];
    global coefficients = [];
    global parameters = [];
    global test_functions = [];
    global indexers = [];
    global genfunc_count = 0;
    global genfunctions = [];
    global linears = [];
    global bilinears = [];
    global face_linears = [];
    global face_bilinears = [];
    global assembly_loops = [];
    global symexpressions = [[],[],[],[]];
    global code_strings = [[],[],[],[]];
    global time_stepper = nothing;
    global specified_dt = 0;
    global specified_Nsteps = 0;
    global use_specified_steps = false;
    global use_cachesim = false;
end

# Sets the code generation target
function set_included_gen_target(lang, framework, dirpath, name; head="")
    #TODO
    println("Premade generation targets currently under construction. Sorry. Reverting to Julia target.");
    return;

    global language = lang;
    global gen_framework = framework;
    global output_dir = dirpath;
    global project_name = name;
    global generate_external = true;
    init_codegenerator(dirpath, name, head);

    # Need to set these three functions

    set_generation_target(get_external_language_elements, generate_external_code_layer, generate_external_files);
end

# Setting a custom target requires three functions
# 1. get_external_language_elements() - file extensions, comment chars etc.
# 2. generate_external_code_layer(var, entities, terms, lorr, vors) - Turns symbolic expressions into code
# 3. generate_external_files(var, lhs_vol, lhs_surf, rhs_vol, rhs_surf) - Writes all files based on generated code
function set_custom_gen_target(lang_elements, code_layer, file_maker, dirpath, name; head="", params=nothing)
    global language = CUSTOM_GEN_TARGET;
    global gen_framework = CUSTOM_GEN_TARGET;
    global output_dir = dirpath;
    global project_name = name;
    global generate_external = true;
    init_code_generator(dirpath, name, head);
    set_generation_target(lang_elements, code_layer, file_maker);
end

# Set parameters used by code generation. This is specific to the target.
function set_codegen_parameters(params)
    global codegen_params = params;
end

# Set the solver module and type
function set_solver(stype)
    config.solver_type = stype;
    if stype == DG
        global solver = DGSolver;
    elseif stype == CG
        global solver = CGSolver;
    elseif stype == FV
        global solver = FVSolver;
        sourcedir = @__DIR__
        add_custom_op_file(string(sourcedir)*"/fv_ops.jl");
    end
end

# Sets the time stepper
function set_stepper(type, cfl)
    global time_stepper = Stepper(type, cfl);
    global config.stepper = type;
    log_entry("Set time stepper to "*type, 1);
end

# Time steps can be chosen automatically or specified. This is for manual specifying.
function set_specified_steps(dt, steps)
    global specified_dt = dt;
    global specified_Nsteps = steps;
    global use_specified_steps = true;
    prob.time_dependent = true;
    prob.end_time = dt*steps;
    log_entry("Set time stepper values to dt="*string(dt)*", Nsteps="*string(steps), 1);
end

# Uses matrix free iterative solver with the given max iterations and error tolerance.
function set_matrix_free(max, tol)
    config.linalg_matrixfree = true;
    config.linalg_matfree_max = max;
    config.linalg_matfree_tol = tol;
end

# Adds a mesh and builds the full grid and reference elements.
# Actually, the mesh object passed here could be either just a MeshData struct,
# or a tuple of MeshData, Refel and Grid already built.
function add_mesh(mesh)
    if typeof(mesh) <: Tuple
        global mesh_data = mesh[1];
        global refel = mesh[2];
        global grid_data = mesh[3];

    else
        global mesh_data = mesh;
        global refel;
        global grid_data;
        (refel, grid_data) = grid_from_mesh(mesh_data);

    end

    constantJ = true;
    do_faces = false;
    do_vol = false;
    if config.solver_type == DG || config.solver_type == FV
        do_faces = true;
    end
    if config.solver_type == FV
        do_vol = true;
    end

    global geo_factors = build_geometric_factors(refel, grid_data, do_face_detj=do_faces, do_vol_area=do_vol, constant_jacobian=constantJ);

    log_entry("Added mesh with "*string(mesh_data.nx)*" vertices and "*string(mesh_data.nel)*" elements.", 1);
    log_entry("Full grid has "*string(size(grid_data.allnodes,2))*" nodes.", 2);

    # If using FV, set up the FV_info
    if config.solver_type == FV
        global fv_info = build_FV_info(grid_data);
    end

    # Set elemental loop ordering to match the order from mesh for now.
    global elemental_order = 1:mesh_data.nel;
end

# Write the mesh to a MSH file
function output_mesh(file, format)
    write_mesh(file, format, mesh_data);
    log_entry("Wrote mesh data to file.", 1);
end

function set_parent_and_child(p_maps, c_grid)
    global parent_maps = p_maps;
    global fv_grid = c_grid;
    dim = config.dimension;
    nfaces = size(fv_grid.element2face,1);
    global fv_refel = refel = build_refel(dim, 1, nfaces, config.elemental_nodes);
    global fv_geo_factors = build_geometric_factors(fv_refel, fv_grid, do_face_detj=true, do_vol_area=true, constant_jacobian=true);
    global fv_info = build_FV_info(fv_grid);
    log_entry("FCreated parent-child grid with "*string(size(c_grid.allnodes,2))*" nodes and "*string(size(c_grid.loc2glb,2))*" elements.", 2);

    # If CELL variables exist, resize their values
    N = size(fv_grid.loc2glb, 2);
    for i=1:length(variables)
        if variables[i].location == CELL
            if var.type == SCALAR
                val_size = (1, N);
            elseif var.type == VECTOR
                val_size = (config.dimension, N);
            elseif var.type == TENSOR
                val_size = (config.dimension*config.dimension, N);
            elseif var.type == SYM_TENSOR
                val_size = (Int((config.dimension*(config.dimension+1))/2), N);
            elseif var.type == VAR_ARRAY
                if typeof(var.indexer) <: Array
                    comps = 1;
                    for j=1:length(var.indexer)
                        comps *= length(var.indexer[j].range);
                    end
                elseif typeof(var.indexer) == Indexer
                    comps = length(var.indexer.range);
                else
                    comps = 1;
                end
                val_size = (comps, N);
            end

            variables[i].values = zeros(val_size);
        end
    end
end

# Maybe remove. This is in grid.jl
function add_boundary_ID(bid, on_bdry)
    add_boundary_ID_to_grid(bid, on_bdry, grid_data);
end

# Defines a test function symbol by creating a special coefficient object.
function add_test_function(v, type)
    varind = length(test_functions) + 1;
    # make SymType
    dim = config.dimension;
    components = 1;
    if type == SCALAR
        components = 1;
    elseif type == VECTOR
        components = dim;
    elseif type == TENSOR
        components = dim*dim;
    elseif type == SYM_TENSOR
        if dim == 1
            components = 1;
        elseif dim == 2
            components = 3;
        elseif dim == 3
            components = 6;;
        elseif dim == 4
            components = 10;
        end
    elseif type == VAR_ARRAY
        components = 1;
    end
    symvar = sym_var(string(v), type, components);

    push!(test_functions, Femshop.Coefficient(v, symvar, varind, type, NODAL, []););
    log_entry("Set test function symbol: "*string(v)*" of type: "*type, 2);
end

# Adds a variable and allocates everything associated with it.
function add_variable(var)
    global var_count += 1;

    # adjust values arrays
    if grid_data === nothing
        N = 1;
    else
        if var.location == CELL
            if fv_grid === nothing
                N = size(grid_data.loc2glb, 2);
            else
                N = size(fv_grid.loc2glb, 2);
            end
        else
            N = size(grid_data.allnodes,2);
        end
    end

    if var.type == SCALAR
        val_size = (1, N);
        comps = 1;
    elseif var.type == VECTOR
        val_size = (config.dimension, N);
        comps = config.dimension;
    elseif var.type == TENSOR
        val_size = (config.dimension*config.dimension, N);
        comps = config.dimension*config.dimension;
    elseif var.type == SYM_TENSOR
        val_size = (Int((config.dimension*(config.dimension+1))/2), N);
        comps = val_size[1];
    elseif var.type == VAR_ARRAY
        if typeof(var.indexer) <: Array
            comps = 1;
            for i=1:length(var.indexer)
                comps *= length(var.indexer[i].range);
            end
        elseif typeof(var.indexer) == Indexer
            comps = length(var.indexer.range);
        else
            comps = 1;
        end
        val_size = (comps, N);
    end
    var.total_components = comps;

    if language == JULIA || language == 0
        var.values = zeros(val_size);
    end

    # make symbolic layer variable symbols
    var.symvar = sym_var(string(var.symbol), var.type, var.total_components);

    global variables = [variables; var];

    global linears = [linears; nothing];
    global bilinears = [bilinears; nothing];
    global face_linears = [face_linears; nothing];
    global face_bilinears = [face_bilinears; nothing];
    global assembly_loops = [assembly_loops; nothing];
    global symexpressions[1] = [symexpressions[1]; nothing];
    global symexpressions[2] = [symexpressions[2]; nothing];
    global symexpressions[3] = [symexpressions[3]; nothing];
    global symexpressions[4] = [symexpressions[4]; nothing];
    global code_strings[1] = [code_strings[1]; ""];
    global code_strings[2] = [code_strings[2]; ""];
    global code_strings[3] = [code_strings[3]; ""];
    global code_strings[4] = [code_strings[4]; ""];

    log_entry("Added variable: "*string(var.symbol)*" of type: "*var.type*", location: "*var.location, 2);
end

# Adds a coefficient with either constant value or some generated function of (x,y,z,t)
function add_coefficient(c, type, location, val, nfuns)
    global coefficients;
    vals = [];
    if nfuns == 0 # constant values
        vals = val;
        if length(vals) == 1 && !(typeof(vals) <: Array)
            vals = [val];
        end

    else # genfunction values
        if typeof(val) <: Array
            ind = length(genfunctions) - nfuns + 1;
            for i=1:length(val)
                if typeof(val[i]) == String
                    push!(vals, genfunctions[ind]);
                    ind += 1;
                else
                    push!(vals, val[i]);
                end
            end
        else
            push!(vals, genfunctions[end]);
        end
    end
    components = length(vals);

    symvar = sym_var(string(c), type, components);

    index = length(coefficients);
    push!(coefficients, Coefficient(c, symvar, index, type, location, vals));

    log_entry("Added coefficient "*string(c)*" : "*string(val), 2);

    return coefficients[end];
end

# Adds a parameter entity.
function add_parameter(p, type, val)
    index = length(parameters);
    push!(parameters, Parameter(p, index, type, val));

    log_entry("Added parameter "*string(p)*" : "*string(val), 2);

    return parameters[end];
end

# Adds an Indexer entity
function add_indexer(indexer)
    push!(indexers, indexer);
    if length(indexer.range) < 3
        log_entry("Added indexer "*string(indexer.symbol)*" : "*string(indexer.range));
    else
        log_entry("Added indexer "*string(indexer.symbol)*" : ["*string(indexer.range[1])*", ... "*string(indexer.range[end])*"]");
    end

end

# Needed to insert parameter expressions into the weak form expressions.
function swap_parameter_xyzt(ex)
    if typeof(ex) == Symbol
        if ex === :x
            return :parameterCoefficientForx ;
        elseif ex === :y
            return :parameterCoefficientFory ;
        elseif ex === :z
            return :parameterCoefficientForz ;
        elseif ex === :t
            return :parameterCoefficientFort ;
        end
    elseif typeof(ex) == Expr
        for i=1:length(ex.args)
            ex.args[i] = swap_parameter_xyzt(ex.args[i]); # Recursively swap
        end
    end
    return ex;
end

# Sets initial condition for the given variable, but does not initialize.
function add_initial_condition(varindex, ex, nfuns)
    global prob;
    while length(prob.initial) < varindex
        prob.initial = [prob.initial; nothing];
    end
    if typeof(ex) <: Array
        vals = [];
        ind = length(genfunctions) - nfuns + 1;
        for i=1:length(ex)
            if typeof(ex[i]) == String
                push!(vals, genfunctions[ind]);
                ind += 1;
            else
                push!(vals, ex[i]);
            end
        end
        prob.initial[varindex] = vals;
    else
        if typeof(ex) == String
            prob.initial[varindex] = genfunctions[end];
        else
            prob.initial[varindex] = ex;
        end
    end

    log_entry("Initial condition for "*string(variables[varindex].symbol)*" : "*string(prob.initial[varindex]), 2);
    # hold off on initializing till solve or generate is determined.
end

# Sets boundary condition for a variable and BID.
function add_boundary_condition(var, bid, type, ex, nfuns)
    global prob;
    # make sure the arrays are big enough
    if size(prob.bc_func)[1] < var_count || size(prob.bc_func,2) < bid
        if !(grid_data===nothing)
            nbid = length(grid_data.bids);
        else
            # For some targets no grid is created
            nbid = 1;
        end

        tmp1 = Array{String,2}(undef, (var_count, nbid));
        tmp2 = Array{Any,2}(undef, (var_count, nbid));
        tmp3 = zeros(Int, (var_count, nbid));
        fill!(tmp1, "");
        fill!(tmp2, GenFunction("","","",0,0));
        tmp1[1:size(prob.bc_func,1), 1:size(prob.bc_func,2)] = Base.deepcopy(prob.bc_type);
        tmp2[1:size(prob.bc_func,1), 1:size(prob.bc_func,2)] = Base.deepcopy(prob.bc_func);
        tmp3[1:size(prob.bc_func,1), 1:size(prob.bc_func,2)] = prob.bid;
        prob.bc_type = tmp1;
        prob.bc_func = tmp2;
        prob.bid = tmp3;
    end
    if typeof(ex) <: Array
        vals = [];
        ind = length(genfunctions) - nfuns + 1;
        for i=1:length(ex)
            if typeof(ex[i]) == String
                push!(vals, genfunctions[ind]);
                ind += 1;
            else
                push!(vals, ex[i]);
            end
        end
        prob.bc_func[var.index, bid] = vals;
    else
        if typeof(ex) == String
            prob.bc_func[var.index, bid] = [genfunctions[end]];
        else
            prob.bc_func[var.index, bid] = [ex];
        end
    end
    prob.bc_type[var.index, bid] = type;
    prob.bid[var.index, bid] = bid;

    log_entry("Boundary condition: var="*string(var.symbol)*" bid="*string(bid)*" type="*type*" val="*string(ex), 2);
end

# A reference point is a single point with a defined value.
# It's treated like a Dirichlet boundary with one boundary node.
function add_reference_point(var, pos, val)
    global prob;
    # make sure the array is big enough
    if size(prob.ref_point,1) < var_count
        tmp = Array{Any,2}(undef, (var_count, 3));
        for i=1:size(prob.ref_point,1)
            tmp[i,:] = prob.ref_point[i,:];
        end
        for i=(size(prob.ref_point,1)+1):var_count
            tmp[i,1] = false;
            tmp[i,2] = [0,0];
            tmp[i,3] = [0];
        end
        prob.ref_point = tmp;
    end
    if typeof(pos) <: Number
        pos = pos*ones(config.dimension);
    end
    if typeof(val) <: Number
        val = val*ones(length(var.symvar));
    end

    # Find the closest vertex to pos
    # The stored pos is actually the index into glbvertex pointing to the closest vertex
    ind = [1,1];
    mindist = 12345;
    for ei=1:mesh_data.nel
        for i=1:size(grid_data.glbvertex,1)
            d = 0;
            for comp=1:length(pos)
                d = d + abs(grid_data.allnodes[comp, grid_data.glbvertex[i, ei]] - pos[comp]);
            end
            if d<mindist
                ind = [i,ei];
                mindist = d;
            end
        end
    end

    prob.ref_point[var.index, 1] = true;
    prob.ref_point[var.index, 2] = ind;
    prob.ref_point[var.index, 3] = val;

    log_entry("Reference point: var="*string(var.symbol)*" position="*string(pos)*" value="*string(val), 2);
end

# Sets the RHS(linear) function for the given variable(s).
function set_rhs(var, code="")
    global linears;
    if language == 0 || language == JULIA
        if typeof(var) <:Array
            for i=1:length(var)
                linears[var[i].index] = genfunctions[end];
                code_strings[2][var[i].index] = code;
            end
        else
            linears[var.index] = genfunctions[end];
            code_strings[2][var.index] = code;
        end

    else # external generation
        if typeof(var) <:Array
            for i=1:length(var)
                linears[var[i].index] = code;
            end
        else
            linears[var.index] = code;
        end
    end
end

# Sets the LHS(bilinear) function for the given variable(s).
function set_lhs(var, code="")
    global bilinears;
    if language == 0 || language == JULIA
        if typeof(var) <:Array
            for i=1:length(var)
                bilinears[var[i].index] = genfunctions[end];
                code_strings[1][var[i].index] = code;
            end
        else
            bilinears[var.index] = genfunctions[end];
            code_strings[1][var.index] = code;
        end

    else # external generation
        if typeof(var) <:Array
            for i=1:length(var)
                bilinears[var[i].index] = code;
            end
        else
            bilinears[var.index] = code;
        end
    end
end

# Sets the surface RHS(linear) function for the given variable(s).
function set_rhs_surface(var, code="")
    global face_linears;
    if language == 0 || language == JULIA
        if typeof(var) <:Array
            for i=1:length(var)
                face_linears[var[i].index] = genfunctions[end];
                code_strings[4][var[i].index] = code;
            end
        else
            face_linears[var.index] = genfunctions[end];
            code_strings[4][var.index] = code;
        end

    else # external generation
        if typeof(var) <:Array
            for i=1:length(var)
                face_linears[var[i].index] = code;
            end
        else
            face_linears[var.index] = code;
        end
    end
end

# Sets the surface LHS(bilinear) function for the given variable(s).
function set_lhs_surface(var, code="")
    global face_bilinears;
    if language == 0 || language == JULIA
        if typeof(var) <:Array
            for i=1:length(var)
                face_bilinears[var[i].index] = genfunctions[end];
                code_strings[3][var[i].index] = code;
            end
        else
            face_bilinears[var.index] = genfunctions[end];
            code_strings[3][var.index] = code;
        end

    else # external generation
        if typeof(var) <:Array
            for i=1:length(var)
                face_bilinears[var[i].index] = code;
            end
        else
            face_bilinears[var.index] = code;
        end
    end
end

function set_symexpressions(var, ex, lorr, vors)
    # If the ex is empty, don't set anything
    if ex == [] || ex == [[]]
        return;
    end

    if typeof(var) <:Array
        for i=1:length(var)
            set_symexpressions(var[i], ex, lorr, vors);
        end

    else
        global symexpressions;
        if lorr == LHS
            if vors == "volume"
                ind = 1;
            else
                ind = 2;
            end
        else # rhs
            if vors == "volume"
                ind = 3;
            else
                ind = 4;
            end
        end
        symexpressions[ind][var.index] = ex;
    end
end

function set_assembly_loops(var, code="")
    global assembly_loops;
    if language == 0 || language == JULIA
        if typeof(var) <:Array
            for i=1:length(var)
                assembly_loops[var[i].index] = genfunctions[end];
            end
        else
            assembly_loops[var.index] = genfunctions[end];
        end

    else # external generation
        if typeof(var) <:Array
            for i=1:length(var)
                assembly_loops[var[i].index] = code;
            end
        else
            assembly_loops[var.index] = code;
        end
    end
end

function finalize()
    finalize_femshop();
end

end # module


text:
function double_it(outs::Vector{Ptr{Void}}, ins::Vector{Ptr{Void}})
  try
    mex_return(outs, [ 2*v for v in mex_args(ins) ]...)
  catch e
    mex_showerror(e)
  end
end

max_stars_repo_name: bermanmaxim/jlcall

max_issues_repo_name: bermanmaxim/jlcall

max_forks_repo_name: bermanmaxim/jlcall


text:
module LibOSXUVC

using Compat

export
    UVCControlInfo,
    UVCRange,
    UVCCameraControl,
    UVCControlType,
    setInterfaceIndex,
    getInterfaceIndex,
    setProcessingUnitId,
    getProcessingUnitId,
    setupByLocationId,
    setBoolValue,
    getBoolValue,
    setNormalizedValue,
    getNormalizedValue,
    getControlInfo,
    getRangeForControl

deps = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(deps)
    include(deps)
else
    error("LibOSXUVC not properly installed. Please run Pkg.build(\"LibOSXUVC\")")
end

const version = convert(VersionNumber,
    unsafe_string(ccall((:OSXUVCVersion, libOSXUVC), Ptr{Cchar}, ())))

macro uvccall(f, rettype, argtypes, args...)
    args = map(esc, args)
    quote
        r = ccall(($f, LibOSXUVC.libOSXUVC),
              $rettype, $argtypes, $(args...))
        if r != Void && r != 0
            error("[OSXUVC]: failed in $($f), code: $r")
        end
        r
    end
end

type UVCControlInfo
    size::UInt16
    selector::UInt16
    unit::UInt16
end

type UVCRange
    min::Cint
    max::Cint
end

type UVCCameraControl
    handle::Ptr{Void}

    function UVCCameraControl()
        handle = ccall((:OSXUVCUVCCameraControlCreate, libOSXUVC),
            Ptr{Void}, ())
        p = new(handle)
        finalizer(p, obj -> ccall((:OSXUVCUVCCameraControlDestroy, libOSXUVC),
            Void,
            (Ptr{Void},), obj.handle))
        p
    end
end

module UVCControlType

import ..LibOSXUVC: libOSXUVC

min_control_type = ccall((:UVCGetMinControlType, libOSXUVC), Cint,())
max_control_type = ccall((:UVCGetMaxControlType, libOSXUVC), Cint,())
for v in min_control_type:max_control_type
    name = ccall((:UVCGetControlTypeShortString, libOSXUVC),
            Ptr{Cchar}, (Int32,), v) |> unsafe_string |> Symbol
    @eval const $name = $v
end

end

function setupByLocationId(ucc::UVCCameraControl, locationId::AbstractString)
    intLocatiionId = ccall(
        (:OSXUVCConvertLocationIdStringToUInt32, libOSXUVC),
        UInt32, (Ptr{Cchar},), pointer(locationId))
    setupByLocationId(ucc, intLocatiionId)
end

function setupByLocationId(ucc::UVCCameraControl, locationId::UInt32)
    @uvccall(:OSXUVCUVCCameraControlSetupByLocationId,
        Cint, (Ptr{Void}, UInt32), ucc.handle, locationId)
    ucc
end

function setInterfaceIndex(ucc::UVCCameraControl, idx)
    @uvccall(:OSXUVCUVCCameraControlSetInterfaceIndex, Cint,
        (Ptr{Void}, UInt16), ucc.handle, idx)
end

function getInterfaceIndex(ucc::UVCCameraControl)
    idx = UInt16[1]
    @uvccall(:OSXUVCUVCCameraControlGetInterfaceIndex, Cint,
        (Ptr{Void}, Ptr{UInt16}), ucc.handle, pointer(idx))
    idx[1]
end

function setProcessingUnitId(ucc::UVCCameraControl, id)
    @uvccall(:OSXUVCUVCCameraControlSetProcessingUnitId, Cint,
        (Ptr{Void}, UInt16), ucc.handle, id)
end

function getProcessingUnitId(ucc::UVCCameraControl)
    id = UInt16[1]
    @uvccall(:OSXUVCUVCCameraControlGetProcessingUnitId, Cint,
        (Ptr{Void}, Ptr{UInt16}), ucc.handle, pointer(id))
    id[1]
end

# convenient constructor
function UVCCameraControl(idx, id,
        locationId::Union{UInt32,AbstractString})
    ucc = UVCCameraControl()
    setInterfaceIndex(ucc, idx)
    setProcessingUnitId(ucc, id)
    setupByLocationId(ucc, locationId)
    ucc
end

### set/get ###

function setBoolValue(ucc::UVCCameraControl, typ, val)
    @uvccall(:OSXUVCUVCCameraControlSetBoolValue,
        Cint, (Ptr{Void}, Cint, Cint), ucc.handle, typ, val)
end

function getBoolValue(ucc::UVCCameraControl, typ)
    val = Cint[1]
    @uvccall(:OSXUVCUVCCameraControlGetBoolValue,
        Cint, (Ptr{Void}, Cint, Ptr{Cint}), ucc.handle, typ, pointer(val))
    val[1]
end

function setNormalizedValue(ucc::UVCCameraControl, typ, val)
    @uvccall(:OSXUVCUVCCameraControlSetNormalizedValue,
        Cint, (Ptr{Void}, Cint, Cfloat), ucc.handle, typ, val)
end

function getNormalizedValue(ucc::UVCCameraControl, typ)
    val = Cfloat[1]
    @uvccall(:OSXUVCUVCCameraControlGetNormalizedValue,
        Cint, (Ptr{Void}, Cint, Ptr{Cfloat}), ucc.handle, typ, pointer(val))
    val[1]
end

# low-level interface

function getControlInfo(ucc::UVCCameraControl, typ)
    info = UVCControlInfo(0,0,0)
    @uvccall(:OSXUVCUVCCameraControlGetControlInfo,
        Cint, (Ptr{Void}, Cint, Ptr{UVCControlInfo}),
        ucc.handle, typ, &info)
    info
end

function getRangeForControl(ucc::UVCCameraControl, control::UVCControlInfo)
    range = UVCRange(0,0)
    @uvccall(:OSXUVCUVCCameraControlGetRangeForControl,
        Cint, (Ptr{Void}, Ref{UVCControlInfo}, Ptr{UVCRange}),
        ucc.handle, control, &range)
    range
end

### utils ###

function Base.dump(ucc::UVCCameraControl)
    ccall((:OSXUVCUVCCameraControlDump, libOSXUVC),
        Void, (Ptr{Void},), ucc.handle)
end

end # module


text:
using DataFrames
using CSV

cd("/Users/kiristern/Documents/GitHub/SDM/")

predators = CSV.read("data/preds/originals/Puma_concolor2.csv")

# select only ursus arctos from imported all predators csv file
#df = predators[predators.species .== "Ursus arctos", :]

prey1 = CSV.read("data/prey_sp/Sus_scrofa2.csv")
prey2 = CSV.read("data/prey_sp/Odocoileus_hemionus2.csv")

#create one dataframe with predators/prey dataframes
pred_prey = vcat(predators, prey1, prey2)

pred_prey = select(pred_prey, [:species, :decimalLatitude, :decimalLongitude, :year])

pred_prey = dropmissing!(pred_prey, [:decimalLatitude, :decimalLongitude, :year])

# rename!(pred_prey, :decimalLatitude => :latitude)
# rename!(pred_prey, :decimalLongitude => :longitude)

CSV.write("data/pred_prey/Puma_concolor.csv", pred_prey)

#select for one species from dataframe
species = [pred_prey[pred_prey.species .== u,:] for u in unique(pred_prey.species)]
occ = species[1]

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

mutable struct CONFIGStruct

    dim::Int64;
    ncx::Int64;
    ncv::Int64;
    ncq::Int64;
    nce::Int64;
    ncf::Int64;
    ncs::Int64;
    nci::Int64;
    nct::Int64;
    ncg::Int64;
    ncz::Int64;
    ncm::Int64;
    nco::Int64;
    ncl::Int64;
    ncp::Int64;
    nconfigs::Int64;
    natom::Array{Int64,2};   # number of atoms per each configuration
    natomall::Int64;# number of atoms for all configurations

    training::Array{Int64,1};
    validation::Array{Int64,1};
    testing::Array{Int64,1};

    # simulation box for each configuration
    a::Array{Float64,2}; # the 1st principal vector of the simulation box
    b::Array{Float64,2}; # the 2nd principal vector of the simulation box
    c::Array{Float64,2}; # the 3rd principal vector of the simulation box
    lattice::Array{Float64,2};  # energies for all configurations
    pbc::Array{Int64,2};  # energies for all configurations
    e::Array{Float64,2};  # energies for all configurations
    stress::Array{Float64,2};  # energies for all configurations

    we::Array{Float64,2};  # energies weights for all configurations
    wf::Array{Float64,2};  # energies weights for all configurations
    ws::Array{Float64,2};  # energies weights for all configurations

    t::Array{Int64,2};   # atom types for all configurations
    Z::Array{Int64,2};   # atom types for all configurations
    move::Array{Int64,2};   # atom types for all configurations
    tags::Array{Int64,2};   # atom types for all configurations
    group::Array{Int64,2};   # atom types for all configurations
    mass::Array{Float64,2}; # atom positions for all configurations
    x::Array{Float64,2}; # atom positions for all configurations
    q::Array{Float64,2}; # atom charges for all configurations
    v::Array{Float64,2}; # atom velocities for all configurations
    f::Array{Float64,2}; # atom forces for all configurations

    CONFIGStruct() = new();
end

function initializeconfig(app)

    config = CONFIGStruct();
    nconfigs = 0;            # number of configurations
    dim = 0
    nci = 0;
    nct = 0;
    ncg = 0;
    ncx = 0;
    ncv = 0;
    ncf = 0;
    ncq = 0;
    nce = 0;
    ncs = 0;
    ncz = 0;
    ncm = 0;
    nco = 0;
    ncl = 0;
    ncp = 0;

    config.nconfigs = nconfigs;
    config.dim = dim;
    config.ncx = ncx;
    config.ncv = ncv;
    config.ncq = ncq;
    config.nce = nce;
    config.ncf = ncf;
    config.ncs = ncs;
    config.nci = nci;
    config.nct = nct;
    config.ncg = ncg;
    config.ncz = ncz;
    config.ncm = ncm;
    config.nco = nco;
    config.ncl = ncl;
    config.ncp = ncp;

    config.natom = ones(1, nconfigs);   # number of atoms per each configuration
    config.natomall = sum(config.natom);# number of atoms for all configurations

    config.training = Array(1:nconfigs)  # training configurations
    config.validation = zeros(Int64,nconfigs)   # validation configurations
    config.testing = zeros(Int64,nconfigs)      # testing configurations

    # simulation box for each configuration
    config.a = zeros(dim, nconfigs); # the 1st principal vector of the simulation box
    config.b = zeros(dim, nconfigs); # the 2nd principal vector of the simulation box
    config.c = zeros(dim, nconfigs); # the 3rd principal vector of the simulation box
    config.lattice = zeros(ncl, nconfigs); # stresses for all configurations
    config.pbc = zeros(ncp, nconfigs); # periodic boundary conditions
    config.e = zeros(nce, nconfigs);  # energies for all configurations
    config.stress = zeros(ncs, nconfigs); # stresses for all configurations

    config.we = ones(1, nconfigs);      # energy weight per each configuration
    config.wf = ones(1, nconfigs);      # force weight per each configuration
    config.ws = ones(1, nconfigs);      # stress weight per each configuration

    config.Z = zeros(ncz, config.natomall);    # atom numbers for all configurations
    config.mass = zeros(ncm, config.natomall); # atom masses for all configurations
    config.move = zeros(nco, config.natomall);   # atom move masks for all configurations
    config.tags = zeros(nci, config.natomall); # atom tags for all configurations
    config.t = zeros(nct, config.natomall);   # atom types for all configurations
    config.group = zeros(ncg, config.natomall);   # atom groups for all configurations
    config.x = zeros(ncx, config.natomall); # atom positions for all configurations
    config.q = zeros(ncq, config.natomall); # atom charges for all configurations
    config.v = zeros(ncv, config.natomall); # atom velocities for all configurations
    config.f = zeros(ncf, config.natomall); # atom forces for all configurations

    return config;
end

function catconfig(config1, config2)

    config = initializeconfig(0)
    config.nconfigs = config1.nconfigs + config2.nconfigs
    config.dim = config1.dim;
    config.ncx = config1.ncx;
    config.ncv = config1.ncv;
    config.ncq = config1.ncq;
    config.nce = config1.nce;
    config.ncf = config1.ncf;
    config.ncs = config1.ncs;
    config.nci = config1.nci;
    config.nct = config1.nct;
    config.ncg = config1.ncg;
    config.ncz = config1.ncz;
    config.ncm = config1.ncm;
    config.nco = config1.nco;
    config.ncl = config1.ncl;
    config.ncp = config1.ncp;

    config.natom = cat(config1.natom, config2.natom, dims=2)
    config.natomall = sum(config.natom);# number of atoms for all configurations

    config.a = cat(config1.a, config2.a, dims=2)
    config.b = cat(config1.b, config2.b, dims=2)
    config.c = cat(config1.c, config2.c, dims=2)
    config.lattice = cat(config1.lattice, config2.lattice, dims=2)
    config.pbc = cat(config1.pbc, config2.pbc, dims=2)
    config.e = cat(config1.e, config2.e, dims=2)
    config.stress = cat(config1.stress, config2.stress, dims=2)

    config.we = cat(config1.we, config2.we, dims=2)
    config.wf = cat(config1.wf, config2.wf, dims=2)
    config.ws = cat(config1.ws, config2.ws, dims=2)

    config.Z = cat(config1.Z, config2.Z, dims=2)
    config.mass = cat(config1.mass, config2.mass, dims=2)
    config.move = cat(config1.move, config2.move, dims=2)
    config.tags = cat(config1.tags, config2.tags, dims=2)
    config.t = cat(config1.t, config2.t, dims=2)
    config.group = cat(config1.group, config2.group, dims=2)
    config.x = cat(config1.x, config2.x, dims=2)
    config.q = cat(config1.q, config2.q, dims=2)
    config.v = cat(config1.v, config2.v, dims=2)
    config.f = cat(config1.f, config2.f, dims=2)

    return config
end

function extractarray(a, indices)

    if size(a,2) > 1
        a = a[:,indices]
    end

    return a
end

function extractconfig(config, indices)

    natom = config.natom[:]
    config.natom = config.natom[:,indices]
    config.natomall = sum(config.natom);
    config.nconfigs = length(config.natom)

    config.a = extractarray(config.a, indices)
    config.b = extractarray(config.b, indices)
    config.c = extractarray(config.c, indices)
    config.e = extractarray(config.e, indices)
    config.pbc = extractarray(config.pbc, indices)
    config.lattice = extractarray(config.lattice, indices)
    config.stress = extractarray(config.stress, indices)

    config.we = extractarray(config.we, indices)
    config.wf = extractarray(config.wf, indices)
    config.ws = extractarray(config.ws, indices)

    natom = [0; cumsum(natom)]
    ind = []
    for i = 1:length(indices)
        j = indices[i]
        if i == 1
            ind = Array((natom[j]+1):natom[j+1]);
        else
            ind = [ind; Array((natom[j]+1):natom[j+1])];
        end
    end

    config.Z = extractarray(config.Z, ind)
    config.mass = extractarray(config.mass, ind)
    config.move = extractarray(config.move, ind)
    config.tags = extractarray(config.tags, ind)
    config.t = extractarray(config.t, ind)
    config.group = extractarray(config.group, ind)
    config.x = extractarray(config.x, ind)
    config.q = extractarray(config.q, ind)
    config.v = extractarray(config.v, ind)
    config.f = extractarray(config.f, ind)

    return config
end

using Random

function traingconfigs(nconfigs, training_percentage, validation_percentage, testing_percentage, randomize)

    if randomize==1
        ind = randperm(nconfigs)
    else
        ind = Array(1:nconfigs)
    end

    ntrain = Int64(round(training_percentage*nconfigs/100.0))
    nvalid = Int64(round(validation_percentage*nconfigs/100.0))
    ntest = Int64(round(testing_percentage*nconfigs/100.0))

    training = ind[1:ntrain]
    validation = ind[(ntrain+1):(ntrain+nvalid)]
    testing = ind[(ntrain+nvalid+1):min(ntrain+nvalid+ntest, nconfigs)]

    return training, validation, testing
end

function setconfig(config, training_percentage, validation_percentage, testing_percentage, randomize)

    config.training, config.validation, config.testing =
        traingconfigs(config.nconfigs, training_percentage, validation_percentage, testing_percentage, randomize)

    return config
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
module Hadoop

using Blocks
using HDFS
using HDFS.MapReduce

importall Blocks

export Block

function Block(f::HdfsURL)
    worker_ids = workers()
    worker_ips = map(x->getaddrinfo(string(isa(x, LocalProcess)?getipaddr():x.host)), map(x->Base.worker_from_id(x), worker_ids))

    block_dist = hdfs_blocks(f, 1, 0, true)
    block_wrkr_ids = map(ips->worker_ids[findin(worker_ips, ips)], block_dist)
    filestat = stat(f)
    block_sz = filestat.block_sz
    file_sz = filestat.size

    data = [(f, ((x-1)*block_sz+1):(min(file_sz,x*block_sz))) for x in 1:length(block_dist)]
    Block(f, data, block_wrkr_ids, as_it_is, as_it_is)
end

end # module Hadoop


text:
include("SubpassDescription.jl")
include("SubpassDependency.jl")
include("AttachmentDescription.jl")
include("ExtensionStructure.jl")

mutable struct RenderPassCreateInfo
    flags::vk.VkRenderPassCreateFlags
    mAttachments::Vector{AttachmentDescription}
    mSubpasses::Vector{SubpassDescription}
    mNext::Ref{ExtensionStructure}

    mVkDependencies::Vector{vk.VkSubpassDependency}

    mHandleRef::Ref{vk.VkRenderPassCreateInfo}

    function RenderPassCreateInfo()
        this = new()
        this.mAttachments = Vector{AttachmentDescription}()
        this.mSubpasses = Vector{SubpassDescription}()
        this.mVkDependencies = Vector{vk.VkSubpassDependency}()
        return this
    end
end
@class RenderPassCreateInfo

function createSimpleForward(::Type{RenderPassCreateInfo}, colorFormat::vk.VkFormat)::RenderPassCreateInfo
    info = RenderPassCreateInfo()

    depth = createWithDepth32float(AttachmentDescription)
    clear(depth)
    discard(depth)
    commit(depth)
    addAttachment(info, depth)

    color = createWithColor(AttachmentDescription, colorFormat)
    clear(color)
    finalLayout(color, vk.VK_IMAGE_LAYOUT_PRESENT_SRC_KHR)
    commit(color)
    addAttachment(info, color)

    addSubpass(info, createFullSubpass(info))
    autoAddDependencies(info)
    commit(info)

    return info
end

function addDependency(this::RenderPassCreateInfo, dep::vk.VkSubpassDependency)
    push!(this.mVkDependencies, dep);
end

function addSubpass(this::RenderPassCreateInfo, sub::SubpassDescription)
    push!(this.mSubpasses, sub);
end

function addAttachment(this::RenderPassCreateInfo, att::AttachmentDescription)
    push!(this.mAttachments, att);
end

function getAttachments(this::RenderPassCreateInfo)::Vector{AttachmentDescription}
    return this.mAttachments
end

function createFullSubpass(this::RenderPassCreateInfo)::SubpassDescription
    sub = SubpassDescription()
    attCount = length(this.mAttachments)
    for i = 1 : attCount
        attachment = UInt32(i - 1)
        if (aspectsOf(getFormat(this.mAttachments[i])) & vk.VK_IMAGE_ASPECT_COLOR_BIT != 0)
            push!(sub.colorAttachments, vk.VkAttachmentReference(
                                            attachment, #attachment::UInt32
                                            vk.VK_IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL#layout::VkImageLayout
                                        ))
        else
            sub.depthStencilAttachment = vk.VkAttachmentReference(
                                            attachment, #attachment::UInt32
                                            vk.VK_IMAGE_LAYOUT_DEPTH_STENCIL_ATTACHMENT_OPTIMAL#layout::VkImageLayout
                                        )
        end
    end
    commit(sub)
    return sub;
end

function autoAddDependencies(this::RenderPassCreateInfo)
    n = length(this.mSubpasses)
    for i = 0 : n
        srcSubpass = (i == 0) ? vk.VK_SUBPASS_EXTERNAL : UInt32(i - 1)
        dstSubpass = (i == n) ? vk.VK_SUBPASS_EXTERNAL : UInt32(i)
        srcStageMask = dstStageMask = 0
        srcAccessMask = dstAccessMask = 0
        if i > 0
            srcSubpass = UInt32(i - 1)
            srcStageMask = vk.VK_PIPELINE_STAGE_COLOR_ATTACHMENT_OUTPUT_BIT
            srcAccessMask = vk.VK_ACCESS_COLOR_ATTACHMENT_WRITE_BIT
        else
            srcSubpass = vk.VK_SUBPASS_EXTERNAL;
            srcStageMask = vk.VK_PIPELINE_STAGE_BOTTOM_OF_PIPE_BIT
        end

        if i < n
            dstSubpass = UInt32(i)
            dstStageMask = vk.VK_PIPELINE_STAGE_FRAGMENT_SHADER_BIT
            dstAccessMask = vk.VK_ACCESS_SHADER_READ_BIT
        else
            dstSubpass = vk.VK_SUBPASS_EXTERNAL;
            dstStageMask = vk.VK_PIPELINE_STAGE_BOTTOM_OF_PIPE_BIT
        end
        dep = vk.VkSubpassDependency(
            srcSubpass, #srcSubpass::UInt32
            dstSubpass, #dstSubpass::UInt32
            srcStageMask, #srcStageMask::VkPipelineStageFlags
            dstStageMask, #dstStageMask::VkPipelineStageFlags
            srcAccessMask, #srcAccessMask::VkAccessFlags
            dstAccessMask, #dstAccessMask::VkAccessFlags
            vk.VK_DEPENDENCY_BY_REGION_BIT #dependencyFlags::VkDependencyFlags
        )

        push!(this.mVkDependencies, dep)
    end
end

function commit(this::RenderPassCreateInfo)
    subpassCount = length(this.mSubpasses)
    vkSubpasses = Vector{vk.VkSubpassDescription}()
    for subpass in this.mSubpasses
        push!(vkSubpasses, handleRef(subpass)[])
    end

    attachCount = length(this.mAttachments)
    vkAttachs = Vector{vk.VkAttachmentDescription}()
    for attach in this.mAttachments
        push!(vkAttachs, handleRef(attach)[])
    end

    this.mHandleRef = Ref(vk.VkRenderPassCreateInfo(
                            vk.VK_STRUCTURE_TYPE_RENDER_PASS_CREATE_INFO, #sType::VkStructureType
                            isdefined(this, :mNext) ? pointer(this.mNext) : C_NULL, #pNext::Ptr{Cvoid}
                            0, #flags::VkRenderPassCreateFlags
                            attachCount, #attachmentCount::UInt32
                            pointer(vkAttachs), #pAttachments::Ptr{VkAttachmentDescription}
                            subpassCount, #subpassCount::UInt32
                            pointer(vkSubpasses), #pSubpasses::Ptr{VkSubpassDescription}
                            length(this.mVkDependencies), #dependencyCount::UInt32
                            pointer(this.mVkDependencies) #pDependencies::Ptr{VkSubpassDependency}
                        ))
end

function handleRef(this::RenderPassCreateInfo)::Ref{vk.VkRenderPassCreateInfo}
    return this.mHandleRef
end


text:
using Logging

logger = SimpleLogger(stdout, Logging.Info)
global_logger(logger)

println("HELLO from ", ENV["SERVER_ID"])

server_id = parse(Int64, ENV["SERVER_ID"])

using Pkg
Pkg.activate(".")
Pkg.offline()

using Pidoh
# using SparseArrays

experiment = loadexperiment("./data.jld2")
run(experiment, server_id)

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
# This script reads a file in gen format containing polygons and stores them
# in a Julia array. The sequence of polygons can be plotted using a plotting
# backend, such as GR, PyPlot, or Gadfly.
#
# EXAMPLES:
#
# We load the results for the BEAM model, which contain 20000 polygons:
#
#     julia> include("plot_gen.jl")
#     julia> using Reachability
#     julia> import GR
#     julia> Xk = read_gen("SpaceEx_plots/beam_t_x89_20_spaceex.txt");
#     julia> R = ReachSolution(Xk, Options(:plot_vars => [1,2]));
#     julia> plot(R)
#     alternative, more spohisticated plotting command
#     julia> plot(R, color="red", xlabel="t", ylabel = "x89", ylims=(-55000, 55000))

# choose a plotting backend, eg. PyPlot or Gadfly
using Plots, LazySets

"""
    read_gen(filename)

Read a sequence of polygons stored in vertex representation (gen format).

### Input

- `filename` -- path of the file containing the polygons

The x and y coordinates of each vertex are separated by an empty space, and
polygons are separated by empty lines. For example:

```
1.01 1.01
0.99 1.01
0.99 0.99
1.01 0.99

0.908463 1.31047
0.873089 1.31047
0.873089 1.28452
0.908463 1.28452
```

### Output

A sequence of matrices, where the first row is for the x coordinates, and the
second row for the y coordinates.

In the example of above:

```julia
    julia> P = read_gen("test.gen")
    2-element Array{Array{Float64,2},1}:
     [1.01 0.99 0.99 1.01; 1.01 1.01 0.99 0.99]
     [0.908463 0.873089 0.873089 0.908463; 1.31047 1.31047 1.28452 1.28452]
```

### Warning

The input file should end with at least one empty line before the end of files
this is to detect the last polygon.
"""
function read_gen(filename::String)
    Mi = Vector{Vector{Float64}}()
    P = Vector{VPolygon{Float64}}()
    # detects when we finished reading a new polygon, needed because polygons
    # may be separated by more than one end-of-line
    new_polygon = true
    open(filename) do f
        for line in eachline(f)
          if !isempty(line)
              push!(Mi, map(x -> parse(Float64, x), split(line)))
              new_polygon = true
          elseif isempty(line) && new_polygon
              push!(P, VPolygon(Mi))
              Mi = Vector{Vector{Float64}}()
              new_polygon = false
          end
        end
    end
    return P
end

function read_gen_old(filename::String)::Array{Matrix{Float64}, 1}
    Mi = Array{Vector{Float64}, 1}()
    P = Array{Matrix{Float64}, 1}()
    # detects when we finished reading a new polygon, needed because polygons
    # may be separated by more than one end-of-line
    new_polygon = true
    open(filename) do f
        for line in eachline(f)
          if !isempty(line)
              push!(Mi, map(x -> parse(Float64, x), split(line)))
              new_polygon = true
          elseif isempty(line) && new_polygon
              push!(P, hcat(Mi...))
              Mi = Array{Vector{Float64}, 1}()
              new_polygon = false
          end
        end
    end
    return P
end

"""
    plot_polygon(P; backend, [name], [gridlines], [plot_polygon])

Plot one or more polygons given as an array of matrices.

INPUT:

- ``P`` -- a polygon, given as a list of matrices, where the first row is for
           the x coordinates and the second row for the y coordinates
- ``backend``     -- (optional, default: ``'pyplot'``): select the plot backend; valid
  options are:
                 -  ``'pyplot_savefig'`` -- use PyPlot package, save to a file
                 -  ``'pyplot_inline'`` -- use PyPlot package, showing in external program
                 - ``'gadfly'`` -- use Gadfly package, showing in browser
                 - ``''`` -- (empty string), return nothing, without plotting
- ``name``        -- (optional, default: ``'plot.png'``) the filename of the plot
  (if it is saved to disk)
- ``gridlines``   -- (optional, default: false) to display or not gridlines in
                   the output plot
- ``color``       -- (optional, default: ``'red'``) the color of the plotting lines
- ``plot_labels`` -- (optional, default: ``['', '']``) the labels of the axes
"""
function plot_polygon(P::Array{Matrix{Float64}, 1};
                      backend="pyplot_savefig",
                      name="plot.png",
                      gridlines=false,
                      color="red",
                      plot_labels::Vector{String}=["", ""])
    if backend in ["pyplot_inline", "pyplot_savefig"]
        if !isdefined(:PyPlot)
            error("this backend requires that your script loads the PyPlot module")
        end
        eval(Expr(:using, :PyPlot))
        if backend == "pyplot_savefig"
            PyPlot.ioff()  # turn off interactive plotting
            fig = PyPlot.figure()
        end
        gridlines ? grid("on") : grid("off")
        # check if P is iterable
        applicable(start, P) ? nothing : P = [P]
        PyPlot.xlabel(plot_labels[1])
        PyPlot.ylabel(plot_labels[2])
        for pi in P
            xcoords = pi[1, :]
            ycoords = pi[2, :]
            PyPlot.plot(xcoords, ycoords, color=color, linewidth=0.6)
        end
        if backend == "pyplot_savefig"
            PyPlot.savefig(name, bbox_inches="tight")
            PyPlot.close(fig)
        end

    elseif backend in ["gadfly"]
        if !isdefined(:Gadfly)
            error("this backend requires that your script loads the Gadfly module")
        end
        eval(Expr(:using, :Gadfly))
        layers_array = []
        # check if P is iterable
        applicable(start, P) ? nothing : P = [P]
        for pi in P
            xcoords = pi[1, :]
            ycoords = pi[2, :]
            push!(layers_array, Gadfly.layer(x=xcoords, y=ycoords, Geom.polygon(preserve_order=false, fill=true)))
        end
        Gadfly.plot(layers_array...)

    elseif backend == ""
        return nothing
    else
        error("plot backend not valid")
    end
end

"""
    plot_all(path)

Plot all txt files (in gen format) of a specified folder.

INPUT:

- ``path`` -- (optional, default: ``'.'``, the current directory) this is the
              directory where the data files are located

OUTUPT:

One or more files with the extension ``.txt.png``.

To give an idea, for 6 files where each file contains 20000 polygons,
it takes around 30 minutes.
"""
function plot_all(path::String=".")
    files = filter(x -> contains(x, ".txt"), readdir(path))
    for f in files
        println("Plotting ", f)
        plot_polygon(read_gen(f), name=f * ".png", gridlines=true)
    end
end

return nothing

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
module simgym

using LinearAlgebra, Distributions

export Cell, Simulation
export integrate!, len_jones, makeCellArray, pixelate

include("Cell.jl")
include("Simulation.jl")
include("potentials.jl")
include("pixelate.jl")
include("cellDivide.jl")

end # module

max_stars_repo_name: andi-bergeson/simgym

max_issues_repo_name: andi-bergeson/simgym

max_forks_repo_name: andi-bergeson/simgym


text:
import CUDA
using ModernGL

# 6 calculate posterior
# get likelihood, CPU version
function depth_ll(depth_predict::Matrix, depth_real::Matrix)

    n_rows, n_cols = size(depth_predict)

    ll = 0.0
    vaild_count = 0
    for c_idx = 1:n_cols
        for r_idx = 1:n_rows
            real_d_value = depth_real[r_idx, c_idx]
            predict_d_value = depth_predict[r_idx, c_idx]
            if (predict_d_value != 0.0) && (real_d_value != 0.0)
                diff = predict_d_value - real_d_value
                # @show diff
                ll += diff^2
                vaild_count += 1
            end
        end
    end
    return ll / max(1, vaild_count)
end

"""
for best performance, set tmp from a outside array
"""
function depth_ll_cuda(frame_buffer::Surfels.GLFrameBuffer,
                        depth_real_f32_cu::CUDA.CuArray{Float32, 2};
                        tmp::CUDA.CuArray{Float32, 2}=similar(depth_real_f32_cu))
    w = frame_buffer.cam.width
    h = frame_buffer.cam.height
    Surfels.glCheckError()
    glCopyImageSubData(frame_buffer.depth_render_buf, GL_TEXTURE_2D, 0, 0, 0, 0,
                       frame_buffer.depth_2nd_tex, GL_TEXTURE_2D, 0, 0, 0, 0,
                       w, h, 1);
    Surfels.glCheckError()

    # # ccall to cuda calculate
    # using CUDA
    # depth_real_f32 = convert(Matrix{Float32}, depth_real_uint) ./ 1000.0
    # depth_real_f32_cu = cu(depth_real_f32)
    # out_log_likelihood_cu = cu(depth_real_f32)
    # a = proj_mat[3, 3]
    # b = proj_mat[3, 4]
    # c = proj_mat[4, 3]
    # void cu_depth_log_likelihood(float* depth_real_ptr,
    #                                void** depth_predict_tex_cuda_res_ptr,
    #                                int width,
    #                                int height,
    #                                double a,
    #                                double b,
    #                                double c,
    #                                float max_depth,
    #                                float* out_log_likelihood);
    out_log_likelihood_cu = tmp
    ccall((:cu_depth_log_likelihood, :libsurfels),
          Cvoid,
    	  (CUDA.CuPtr{Cfloat}, Ptr{Ptr{Cvoid}}, Cint, Cint, Cdouble, Cdouble,
    	  Cdouble, Cfloat, Cfloat, Cint, Cint, CUDA.CuPtr{Cfloat}),
    	  depth_real_f32_cu.storage.buffer.ptr,
    	  Base.unsafe_convert(Ptr{Ptr{Cvoid}}, frame_buffer.depth_cuda_res),
    	  w,
    	  h,
    	  frame_buffer.a,
    	  frame_buffer.b,
    	  frame_buffer.c,
    	  frame_buffer.min_depth,
    	  frame_buffer.max_depth,
    	  16, # thread block size: 16 x 16
    	  Int32(0),
    	  out_log_likelihood_cu.storage.buffer.ptr)

    # out = collect(out_log_likelihood_cu)
    # imshow(Gray.(out))

    cu_ll = sum(out_log_likelihood_cu) #/ sum(out_log_likelihood_cu .!= 0.0)
    cu_count = sum(out_log_likelihood_cu .!= 0.0)
    return cu_ll / max(1, cu_count)
end

function DMP_RGBD(
    rgbd_reader::T,
    model_fpath::String,
    id2Twc, # traj
    cam::CameraModel.RgbdCamParams;
    min_depth::Float64 = 0.3,
    max_depth::Float64 = 15.0,
    confidence_thres::Float64=0.1,
    max_id::Int64=100000000000000,
    tmp::CUDA.CuArray{Float32, 2}=CUDA.CuArray{Float32, 2}(
        undef, Int64(cam.height), Int64(cam.width)
    )
) where T <: SLAMData.RgbdReader

    hidden_window = Surfels.create_hidden_gl_window(cam.width, cam.height)

    program::ModernGL.GLuint = Surfels.create_gl_render_program()

    frame_buffer::Surfels.GLFrameBuffer = Surfels.GLFrameBuffer(
        cam.width, cam.height, cam, min_depth, max_depth, use_texture = true
    )

    # load GlModel
    model::Surfels.GlSurfelsData = Surfels.read_surfels_file2gl(model_fpath)

    ll = 0.0
    count = 0
    Surfels.glCheckError()

    ids = sort(collect(keys(id2Twc)))

    for frame_id = ids
        if frame_id > max_id
            continue
        end
        #render to get predicted observation
        Surfels.render(
            [model],
            id2Twc[frame_id], # Twc
            confidence_thres,
            1,
            200000,
            program,
            frame_buffer
        )
        Surfels.glCheckError()

        # real observation
        # depth_real_f32_cu = frames_d_cu[frame_id]
        img, depth_real_uint = SLAMData.readRGBD(rgbd_reader, id=frame_id)
        depth_real_f = convert(Matrix{Float32}, depth_real_uint) ./ 1000.0
        depth_real_f32_cu = CUDA.cu(depth_real_f)

        ll += depth_ll_cuda(frame_buffer, depth_real_f32_cu, tmp=tmp)

        Surfels.glCheckError()
        count += 1
    end
    # TODO destroy program
    finalize(frame_buffer)
    finalize(model)
    finalize(hidden_window)

    return ll #/ Float64(count)
end


text:
# This is part of Quiz.jl. License is MIT.

"""
`struct qT <: qElem`\n
Quiz generic text element.
"""
struct qT <: qElem
    val::AbstractString
end

# Type functor: full extraction and assembly
"""
`(x::qT)()`\n
`qT` functor that returns its value as a `string`.
"""
(x::qT)() = string(x.val)

# the: convenient value extracting function
the(x::qT) = x()

# export
export qT

#----------------------------------------------------------------------------------------------#
#                                       qC -- Constants                                        #
#----------------------------------------------------------------------------------------------#

"""
`struct qC <: qElem`\n
Quiz constant element.
"""
struct qC <: qElem
    nam::Symbol
    val::N64
    uni::Unitful.Units
    qC(n::Symbol, x::Number) = new(n, val_unit(x)...)
    qC(n::Symbol, v::Number, u::Unitful.Units) = begin
        val, uni = val_unit(v)
        nou = unit(1) # a non-unit
        new(n, val, uni == nou ? u : (u == nou ? uni : u*uni))
    end
end

# Type functor: full extraction and assembly
"""
`(x::qC)()`\n
`qC` functor that returns its value as a `Pair{Symbol, Quantity}`.
"""
(x::qC)() = x.nam => Quantity{typeof(x.val), dimension(x.uni), typeof(x.uni)}(x.val * x.uni)

# the: convenient value extracting function
the(x::qC) = x().second

# iE: convenient constant instantiation expression
iE(x::qC) = begin a, b = x(); :($a = $b); end

# export
export qC

#----------------------------------------------------------------------------------------------#
#                                       qV -- Variables                                        #
#----------------------------------------------------------------------------------------------#

"""
`struct qV <: qElem`\n
Quiz variant (or variable) element.
"""
struct qV <: qElem
    nam::Symbol
    val::Array{<:N64, 1}
    uni::Unitful.Units
    qV(n::Symbol, v::AbstractVector{<:Number}) = begin
        vSep = map(x -> val_unit(x), v)
        uSet = [Set(i[2] for i in vSep)...]
        if length(uSet) == 1
            new(n, [i[1] for i in vSep], uSet[1])
        else
            throw(ArgumentError("non-uniform units on input"))
        end
    end
    qV(n::Symbol, v::AbstractVector{<:Number}, u::Unitful.Units) = begin
        vSep = map(x -> val_unit(x), v)
        uSet = [Set(i[2] for i in vSep)...]
        if length(uSet) == 1
            Uni, nou = uSet[1], unit(1)
            UNI = Uni == nou ? u : (u == nou ? Uni : u*Uni)
            new(n, [i[1] for i in vSep], UNI)
        else
            throw(ArgumentError("non-uniform units on input"))
        end
    end
end

# Type functor: full extraction and assembly
"""
`(x::qV)()`\n
`qV` functor that returns its value as a `Pair{Symbol, Vector{Quantity}}`.
"""
function (x::qV)()
    TYP, DIM, UNI = typeof(x.val[1]), dimension(x.uni), typeof(x.uni)
    x.nam => [ Quantity{TYP, DIM, UNI}(i * x.uni) for i in x.val ]
end

"""
`(x::qV)(i::Integer)`\n
`qV` functor that returns its `i`th value as a `Pair{Symbol, Quantity}`.
"""
function (x::qV)(i::Integer)
    X = x()
    X.first => X.second[i]
end

# the: convenient value extracting function
the(x::qV) = x().second

# iE: convenient constant instantiation expression
iE(x::qV) = begin a, b = x(); :($a = rand($b)); end

# Extending rand(::qV)
import Random.rand
rand(x::qV) = rand(the(x))

# qC(::qV)
qC(v::qV) = qC(v.nam, rand(v))

# export
export qV

#----------------------------------------------------------------------------------------------#
#                                        qA -- Answers                                         #
#----------------------------------------------------------------------------------------------#

"""
`struct qA <: qElem`\n
Quiz answer element.
"""
struct qA <: qElem
    wgt::N64
    ask::AbstractString
    nam::Symbol
    uni::Unitful.Units
end

# export
export qA


text:
function igyp01()
#
#  generalization of Igarash and Ypma
#
    p = poly([2.35*[1;1;1];2.56])
    z = [2.35 3;2.56 1];
    p, PolyZeros(z)
end


text:
# LSL.jl: Julia interface for Lab Streaming Layer
# Copyright (C) 2019 Samuel Powell

# Resolver.jl: functions and types for resolving streams

export resolve_streams, resolve_byprop, resolve_bypred

const MAX_RESOLVE = 1024

"""
    Continuous Resolver

A convenience class resolving streams continuously in the background. This object can be
queried at any time for the set of streams that are currently visible on the network.
"""
mutable struct ContinuousResolver
  handle::lib.lsl_continuous_resolver

  function ContinuousResolver(handle)
    resolver = new(handle)
    finalizer(_destroy, resolver)
  end

end

# Define conversion to pointer
unsafe_convert(::Type{lib.lsl_continuous_resolver}, resolver::ContinuousResolver) = resolver.handle

function _destroy(resolver::ContinuousResolver)
  if resolver.handle != C_NULL
    lib.lsl_destroy_continuous_resolver(resolver)
    resolver.handle = C_NULL
  end
  return nothing
end

"""
    ContinuousResolver(; forget_after = 5.0)

Construct a continuous resolver which resolves all streams on the network.

# Keyword arguments
`forget_after::Number`: When a stream is no longer visible on the network (e.g., because it
                        was shut down), this is the time in seconds after which it is no
                        longer reported by the resolver.
"""
function ContinuousResolver(; forget_after = 5.0)
  handle = lib.lsl_create_continuous_resolver(forget_after)
  if handle == C_NULL
    error("liblsl library returned NULL pointer during ContinuousResolver creation")
  end
  return ContinuousResolver(handle)
end

"""
  ContinuousResolver(prop = "", value = ""; forget_after = 5.0)

Construct a continuous resolver that resolves all streams with a specific value for a
given property.

This is analogous to the functionality provided by the function resolve_stream(prop,value).

# Arguments
-`prop::String`: The StreamInfo property that should have a specific value (e.g., "name",
"type", "source_id", or "desc/manufaturer").
-`value::String`: The string value that the property should have (e.g., "EEG" as the type
                         property).

# Keyword arguments
-`forget_after::Number`: When a stream is no longer visible on the network (e.g., because it
                        was shut down), this is the time in seconds after which it is no
                        longer reported by the resolver.
"""
function ContinuousResolver(prop::String, value::String; forget_after = 5.0)
  handle = lib.lsl_create_continuous_resolver_byprop(prop, value, forget_after)
  if handle == C_NULL
    error("liblsl library returned NULL pointer during ContinuousResolver creation")
  end
  return ContinuousResolver(handle)
end

"""
  ContinuousResolver(predicate; forget_after = 5.0)

Construct a new continuous resolver that resolves all streams that match a given XPath
1.0 predicate.

This is analogous to the functionality provided by the free function resolve_bypred

# Arugments
-`predicate::String`: Predicate string, e.g. "name='BioSemi'" or "type='EEG' and
starts-with(name,'BioSemi') and count(description/desc/channels/channel)=32"

# Keyword arguments
-`forget_after::Number`: When a stream is no longer visible on the network (e.g., because it
                        was shut down), this is the time in seconds after which it is no
                        longer reported by the resolver.
"""
function ContinuousResolver(predicate::String; forget_after = 5.0)
  handle = lib.lsl_create_continuous_resolver_bypred(predicate, forget_after)
  if handle == C_NULL
    error("liblsl library returned NULL pointer during ContinuousResolver creation")
  end
  return ContinuousResolver(handle)
end

"""
    resolve_streams(;timeout = 1.0)

Returns all currently available streams from any outlet on the network.

The network is usually the subnet specified at the local router, but may also include a
group of machines visible to each other via multicast packets (given that the network
supports it), or list of hostnames. These details may optionally be customized by the
experimenter in a configuration file (see Network Connectivity in the LSL wiki).

Returns a vector of StreamInfo objects (with empty desc field), any of which can
subsequently be used to open an inlet. The full description can be retrieved from the inlet.

# Keyword arguments:
-`timeout`: The waiting time for the operation, in seconds, to search for streams. Warning:
            If this is too short (<0.5s) only a subset (or none) of the outlets that are
            present on the network may be returned.
"""
function resolve_streams(;timeout = 1.0)
  infoptrs = Vector{lib.lsl_streaminfo}(undef, MAX_RESOLVE)
  nstreams = handle_error(lib.lsl_resolve_all(infoptrs, MAX_RESOLVE, timeout))
  infos = Vector{StreamInfo}(undef, nstreams)
  for i in 1:nstreams
    format_type = _jl_channel_format(lib.lsl_get_channel_format(infoptrs[i]))
    infos[i] = StreamInfo(infoptrs[i], format_type)
  end
  return infos
end

"""
    resolve_byprop(prop, value; minimum = 1, timeout = LSL_FOREVER)

Resolve all streams with a specific value for a given property.

If the goal is to resolve a specific stream, this method is preferred over resolving all
streams and then selecting the desired one.

Returns a vector of matching StreamInfo objects (with empty desc field), any of which can
subsequently be used to open an inlet.

# Arguments
-`prop::String`: The StreamInfo property that should have a specific value (e.g., "name",
                "type", "source_id", or "desc/manufaturer").
-`value::String`: The string value that the property should have (e.g., "EEG" as the type
                 property).

# Keyword arguments
-`minimum::Integer`: Return at least this many streams.
-`timeout::Number`: A timeout of the operation, in seconds. If the timeout expires, less than
                   the desired number of streams (possibly none) will be returned.

# Example
results = resolve_byprop("type", "EEG")
"""
function resolve_byprop(prop::String, value::String; minimum = 1, timeout = FOREVER)
  infoptrs = Vector{lib.lsl_streaminfo}(undef, MAX_RESOLVE)
  nstreams = handle_error(lib.lsl_resolve_byprop(infoptrs, MAX_RESOLVE,
                                                 prop, value, minimum, timeout))
  infos = Vector{StreamInfo}(undef, nstreams)
  for i in 1:nstreams
    format_type = _jl_channel_format(lib.lsl_get_channel_format(infoptrs[i]))
    infos[i] = StreamInfo(infoptrs[i], format_type)
  end
  return infos
end

"""
    resolve_bypred(predicate; minimum = 1, timeout = LSL_FOREVER)

Resolve all streams that match a given predicate.

Advanced query that allows to impose more conditions on the retrieved streams; the given
string is an XPath 1.0 predicate for the <description> node (omitting the surrounding []'s),
see also http://en.wikipedia.org/w/index.php?title=XPath_1.0&oldid=474981951.

Returns a vector of matching StreamInfo objects (with empty desc field), any of which can
subsequently be used to open an inlet.

# Argumennts:
-`predicate::String`: Predicate string, e.g. "name='BioSemi'" or "type='EEG' and
                     starts-with(name,'BioSemi') and count(description/desc/channels/channel)=32"

# Keyword arguments
-`minimum::Integer`: Return at least this many streams.
-`timeout::Number`: A timeout of the operation, in seconds. If the timeout expires, less than
                  the desired number of streams (possibly none) will be returned.
"""
function resolve_bypred(predicate::String; minimum = 1, timeout = FOREVER)
  infoptrs = Vector{lib.lsl_streaminfo}(undef, MAX_RESOLVE)
  nstreams = handle_error(lib.lsl_resolve_bypred(infoptrs, MAX_RESOLVE,
                                                 predicate, minimum, timeout))
  infos = Vector{StreamInfo}(undef, nstreams)
  for i in 1:nstreams
    format_type = _jl_channel_format(lib.lsl_get_channel_format(infoptrs[i]))
    infos[i] = StreamInfo(infoptrs[i], format_type)
  end
  return infos
end


text:
"""
    PathGraph{T<:Integer}

A path graph, with each node linked to the previous and next one.
"""
struct PathGraph{T<:Integer} <: LG.AbstractGraph{T}
    nv::T
    function PathGraph{T}(nv::T) where {T<:Integer}
        _nv = nv >= 0 ? T(nv) : zero(T)
        new{T}(_nv)
    end
end

PathGraph(nv::T) where {T<:Integer} = PathGraph{T}(nv)

LG.eltype(::PathGraph{T}) where {T} = T
LG.edgetype(::PathGraph{T}) where {T} = LG.Edge{T}
LG.is_directed(::Type{<:PathGraph}) = false
LG.nv(g::PathGraph) = g.nv
LG.ne(g::PathGraph) = max(0, Int(LG.nv(g)) - 1)
LG.vertices(g::PathGraph) = 1:LG.nv(g)

LightGraphs.edges(g::PathGraph) = [LG.Edge(i, i+1) for i in 1:LG.nv(g)-1]

LightGraphs.has_vertex(g::PathGraph, v) = 1 <= v <= LG.nv(g)

function LightGraphs.outneighbors(g::PathGraph, v)
    LG.has_vertex(g, v) || return Int[]
    nv(g) > 1 || return Int[]
    if v == 1
        return [2]
    end
    if v == nv(g)
        return [nv(g)-1]
    end
    return [v-1, v+1]
end

LightGraphs.inneighbors(g::PathGraph, v) = outneighbors(g, v)

function LightGraphs.has_edge(g::PathGraph, v1, v2)
    if !has_vertex(g, v1) || !has_vertex(g, v2)
        return false
    end
    return abs(v1-v2) == 1
end

# =======================================================
#         overrides
# =======================================================

# we use this check so that we have the same convention as in LightGraphs
LG.is_connected(g::PathGraph) = nv(g) > 0

function LG.connected_components(g::PathGraph)

    nv(g) == 0 && return typeof(vertices(g))[]
    return [vertices(g)]
end

# has_self_loops is defined in terms of this
LG.num_self_loops(::PathGraph) = 0

LG.is_bipartite(g::PathGraph) = true

function LG.squash(g::PathGraph)
    nvg = nv(g)
    for T ∈ (UInt8, UInt16, UInt32, UInt64)
        nv(g) < typemax(T) && return PathGraph(T(nvg))
    end
end

# ---- degree -----------------------------------

function LG.Δ(g::PathGraph)
    nvg = nv(g)
    return typemin(Int) * (nvg == 0) + (nvg >= 2) + (nvg >= 3)
end

LG.Δout(g::PathGraph) = Δ(g)
LG.Δin(g::PathGraph) = Δ(g)

function LG.δ(g::PathGraph)
    nvg = nv(g)
    return ifelse(nvg == 0, typemax(Int), ifelse(nvg == 1, 0, 1))
end

LG.δout(g::PathGraph) = δ(g)
LG.δin(g::PathGraph) = δ(g)

# =======================================================
#         matrices
# =======================================================

# TODO If there is a library that provides this, we could use a Toeplitz matrix instead
function LG.adjacency_matrix(g::PathGraph, T::DataType=Int; dir=:out)

    nvg = nv(g)
    dv = zeros(T, nv(g))
    ev = ones(T, max(0, Int(nv(g)) - 1))
    return SymTridiagonal(dv, ev)
end


text:
"""Vector geometry."""

mutable struct Vector
    x
    y
end

Vector() = Vector(0.0, 0.0)

import Base.+, Base.-, Base.*, Base./, Base.==
add(a::Vector, b::Vector) = Vector(a.x + b.x, a.y + b.y)
sub(a::Vector, b::Vector) = Vector(a.x - b.x, a.y - b.y)
mult(a::Vector, scalar) = Vector(a.x * scalar, a.y * scalar)
div(a::Vector, scalar) = Vector(a.x / scalar, a.y / scalar)
neg(v::Vector) = mult(v, -1)
(+)(a::Vector, b::Vector) = add(a, b)
(-)(a::Vector, b::Vector) = sub(a, b)
(*)(a::Vector, scalar) = mult(a, scalar)
(/)(a::Vector, scalar) = div(a, scalar)
(-)(v::Vector) = neg(v)
function (==)(a::Vector, b::Vector)
    a.x == b.x && a.y == b.y && return true
    return false
end

function add!(a::Vector, b::Vector)
    a.x += b.x
    a.y += b.y
    return a
end

function sub!(a::Vector, b::Vector)
    a.x -= b.x
    a.y -= b.y
    return a
end

function mult!(v::Vector, scalar)
    v.x *= scalar
    v.y *= scalar
    return v
end

function div!(v::Vector, scalar)
    v.x /= scalar
    v.y /= scalar
    return v
end

magnitude(v::Vector) = sqrt(v.x * v.x + v.y * v.y)

function normalize!(v::Vector)
    m = magnitude(v)
    if m == 0
        v.x = 0.0
        v.y = 0.0
    else
        v.x /= m
        v.y /= m
    end
    return v
end

dot(a::Vector, b::Vector) = a.x * b.x + a.y * b.y
det(a::Vector, b::Vector) = a.x * b.y - a.y * b.x
cross(v::Vector) = Vector(-v.x, v.y)

function rotate!(v::Vector, angle, point)
    angle == 0. && return
    dx = v.x - point.x
    dy = v.y - point.y
    v.x = point.x + (dx * cos(angle) - dy * sin(angle))
    v.y = point.y + (dx * sin(angle) + dy * cos(angle))
    return v
end

rotate!(v::Vector, angle) = rotate!(v, angle, Vector(0., 0.))

import Base.isapprox
function isapprox(a::Vector, b::Vector; atol = 1e-16)
    isapprox(a.x, b.x; atol = atol) && isapprox(a.y, b.y; atol = atol) && return true
    return false
end


text:
const JNI_VERSION_1_1 = convert(Cint, 0x00010001)
const JNI_VERSION_1_2 = convert(Cint, 0x00010002)
const JNI_VERSION_1_4 = convert(Cint, 0x00010004)
const JNI_VERSION_1_6 = convert(Cint, 0x00010006)
const JNI_VERSION_1_8 = convert(Cint, 0x00010008)

const JNI_TRUE = convert(Cchar, 1)
const JNI_FALSE = convert(Cchar, 0)

# Return Values
const JNI_OK           = convert(Cint, 0)               #/* success */
const JNI_ERR          = convert(Cint, -1)              #/* unknown error */
const JNI_EDETACHED    = convert(Cint, -2)              #/* thread detached from the VM */
const JNI_EVERSION     = convert(Cint, -3)              #/* JNI version error */
const JNI_ENOMEM       = convert(Cint, -4)              #/* not enough memory */
const JNI_EEXIST       = convert(Cint, -5)              #/* VM already created */
const JNI_EINVAL       = convert(Cint, -6)              #/* invalid arguments */

const JAVA_HOME_CANDIDATES = ["/usr/lib/jvm/default-java/",
                              "/usr/lib/jvm/default/"]

function javahome_winreg()
    keys = ["SOFTWARE\\JavaSoft\\Java Runtime Environment", "SOFTWARE\\JavaSoft\\Java Development Kit", "SOFTWARE\\JavaSoft\\JDK"]

    for key in keys
        try
            value = querykey(WinReg.HKEY_LOCAL_MACHINE, key, "CurrentVersion")
            key *= "\\" * value
            return querykey(WinReg.HKEY_LOCAL_MACHINE, key, "JavaHome")
        catch
            # Try the next value in the loop.
        end
    end

    error("Cannot find an installation of Java in the Windows Registry. Please install a JRE/JDK, or set the JAVA_HOME environment variable if one is already installed.")
end

@static Sys.isunix() ? (global const libname = "libjvm") : (global const libname = "jvm")

function findjvm()
    javahomes = Any[]
    libpaths = Any[]

    if haskey(ENV, "JAVA_HOME")
        push!(javahomes, ENV["JAVA_HOME"])
    else
        @static if Sys.iswindows()
            ENV["JAVA_HOME"] = javahome_winreg()
            push!(javahomes, ENV["JAVA_HOME"])
        end
        @static if Sys.isunix()
            # Find default javahome by checking location of the java command
            try
                javapath = chomp(read(`which java`,String))
                while(islink(javapath))
                    javapath = readlink(javapath)
                end
                javapath = dirname(javapath)
                javapath = match(r"(.*)(/jre)+(/bin)+",javapath)[1]
                push!(javahomes,javapath)
            catch err
                @info "JavaCall could not determine javapath from `which java`"
            end
        end
    end
    isfile("/usr/libexec/java_home") && push!(javahomes, chomp(read(`/usr/libexec/java_home`, String)))

    for fname ∈ JAVA_HOME_CANDIDATES
        isdir(fname) && push!(javahomes, fname)
    end

    push!(libpaths, pwd())
    for n in javahomes
        @static if Sys.iswindows()
            push!(libpaths, joinpath(n, "bin", "server"))
            push!(libpaths, joinpath(n, "jre", "bin", "server"))
            push!(libpaths, joinpath(n, "bin", "client"))
        end
        @static if Sys.islinux()
            if Sys.WORD_SIZE==64
                push!(libpaths, joinpath(n, "jre", "lib", "amd64", "server"))
                push!(libpaths, joinpath(n, "lib", "amd64", "server"))
            elseif Sys.WORD_SIZE==32
                push!(libpaths, joinpath(n, "jre", "lib", "i386", "server"))

                push!(libpaths, joinpath(n, "lib", "i386", "server"))
            end
        end
        push!(libpaths, joinpath(n, "jre", "lib", "server"))
        push!(libpaths, joinpath(n, "lib", "server"))
    end

    ext = @static Sys.iswindows() ? "dll" : (@static Sys.isapple() ? "dylib" : "so")
    ext = "."*ext

    try
        for n in libpaths
            libpath = joinpath(n,libname*ext);
            if isfile(libpath)
                if Sys.iswindows()
                    bindir = dirname(dirname(libpath))
                    m = filter(x -> occursin(r"msvcr(?:.*).dll",x), readdir(bindir))
                    if !isempty(m)
                        Libdl.dlopen(joinpath(bindir,m[1]))
                    end
                end
                global libjvm = Libdl.dlopen(libpath)
                @debug("Loaded $libpath")
                return
            end
        end
    catch
    end

    errorMsg =
        [
         "Cannot find java library $libname$ext\n",
         "Search Path:"
         ];
    for path in libpaths
        push!(errorMsg,"\n   $path")
    end
    throw(JavaCallError(reduce(*,errorMsg)))
end

struct JavaVMOption
    optionString::Ptr{UInt8}
    extraInfo::Ptr{Nothing}
end

struct JavaVMInitArgs
    version::Cint
    nOptions::Cint
    options::Ptr{JavaVMOption}
    ignoreUnrecognized::Cchar
end

@static Sys.isunix() ? (const sep = ":") : nothing
@static Sys.iswindows() ? (const sep = ";") : nothing
cp = OrderedSet{String}()
opts = OrderedSet{String}()

"""
    JavaCall.getClassPath()

    Obtains the Java classpath.

    Before JavaCall.init(), this reports the classpath constructed using
    JavaCall.addClasspath.

    After JavaCall.init(), this reports System.getProperty("java.class.path")
"""
function getClassPath()
    if isloaded()
        # jls = @jimport java.lang.System
        jls = JavaObject{Symbol("java.lang.System")}
        return jcall(jls, "getProperty", JString, (JString,), "java.class.path")::String
    else
        ccp = collect(cp)
        return join(ccp,sep)
    end
end

"""
    JavaCall.addClassPath(s::String)

    Add a string to the classpath. Must be called before JavaCall.init

    foo/*.jar will add all the jars in the directory foo to the classpath
    foo/*     will add all the jars and directories recursively including foo

    See also addClassPathRecursive and addJarsToClassPath
"""
function addClassPath(s::String)
    if isloaded()
        @warn("JVM already initialised. This call has no effect")
        return
    end
    if s==""; return; end
    dirname, pattern = splitdir(s)
    if pattern == "*.jar" && isdir(dirname)
        _addJarsToClassPath(dirname)
    elseif pattern == "*" && isdir(dirname)
        _addClassPathRecursive(dirname)
    else
        push!(cp, s)
    end
    return
end

function _addClassPathRecursive(dirname::String)
    push!(cp,dirname)
    for (root,dirs,filenames) in walkdir(dirname)
        addJarsToClassPath(filenames,root)
        union!(cp, [joinpath(root,dir) for dir in dirs])
    end
end
"""
    JavaCall.addClassPathRecursive(dirname::String)

    Adds dirname and all jars and directories in dirname recursively to the classpath
"""
addClassPathRecursive(dirname) = isloaded() ?
    @warn("JVM already initialized. This call has no effect") :
    _addClassPathRecursive(dirname)

function _addJarsToClassPath(files::Array{String,1},path::String="")
    for filename in files
        if endswith(filename,".jar")
            push!(cp, joinpath(path,filename) )
        end
    end
end
_addJarsToClassPath(dirname::String) = addJarsToClassPath(readdir(dirname),dirname)
"""
    JavaCall.addJarsToClassPath(dirname::String)

    Add jars in dirname to the classpath.
    Equivalent to JavaCall.addJarsToClasspath(readdir(dirname),dirname)

    JavaCall.addJarsToClassPath(files::Array{String,1}, [path::String])

    Add files which end in ".jar" to the classpath prefixed by an optional path
"""
addJarsToClassPath(args...) = isloaded() ?
    @warn("JVM already initialized. This call has no effect") :
    _addJarsToClassPath(args...)

addOpts(s::String) = isloaded() ? @warn("JVM already initialised. This call has no effect") : push!(opts, s)

function init()
    if isempty(cp)
        init(opts)
    else
        ccp = collect(cp)
        options = collect(opts)
        classpath = "-Djava.class.path="*join(ccp,sep)
        init(vcat(options, classpath))
    end
end

const ROOT_TASK_ERROR = JavaCallError(
    "Either the environmental variable JULIA_COPY_STACKS must be 1 " *
    "OR JavaCall must be used on the root Task.")

const JULIA_COPY_STACKS_ON_WINDOWS_ERROR = JavaCallError(
    "JULIA_COPY_STACKS should not be set on Windows.")

# JavaCall must run on the root Task or JULIA_COPY_STACKS is enabled
isroottask() = Base.roottask === Base.current_task()
@static if Sys.iswindows()
    isgoodenv() = ! JULIA_COPY_STACKS
    assertroottask_or_goodenv() = isgoodenv() ? nothing : throw(JULIA_COPY_STACKS_ON_WINDOWS_ERROR)
else
    isgoodenv() = JULIA_COPY_STACKS || isroottask()
    assertroottask_or_goodenv() = isgoodenv() ? nothing : throw(ROOT_TASK_ERROR)
end

isloaded() = isdefined(JavaCall, :jnifunc) && isdefined(JavaCall, :penv) && penv != C_NULL

assertloaded() = isloaded() ? nothing : throw(JavaCallError("JVM not initialised. Please run init()"))
assertnotloaded() = isloaded() ? throw(JavaCallError("JVM already initialised")) : nothing

# Pointer to pointer to pointer to pointer alert! Hurrah for unsafe load
function init(opts)
    assertnotloaded()
    assertroottask_or_goodenv()
    opt = [JavaVMOption(pointer(x), C_NULL) for x in opts]
    ppjvm = Array{Ptr{JavaVM}}(undef, 1)
    ppenv = Array{Ptr{JNIEnv}}(undef, 1)
    vm_args = JavaVMInitArgs(JNI_VERSION_1_6, convert(Cint, length(opts)),
                             convert(Ptr{JavaVMOption}, pointer(opt)), JNI_TRUE)
    res = ccall(create, Cint, (Ptr{Ptr{JavaVM}}, Ptr{Ptr{JNIEnv}}, Ptr{JavaVMInitArgs}), ppjvm, ppenv,
                Ref(vm_args))
    res < 0 && throw(JavaCallError("Unable to initialise Java VM: $(res)"))
    global penv = ppenv[1]
    global pjvm = ppjvm[1]
    jnienv = unsafe_load(penv)
    jvm = unsafe_load(pjvm)
    global jvmfunc = unsafe_load(jvm.JNIInvokeInterface_)
    global jnifunc = unsafe_load(jnienv.JNINativeInterface_) #The JNI Function table
    return
end

"""
    init_current_vm()

Allow initialization from running VM. Uses the first VM it finds.

# Example using JNA

```java
package zot.julia;

import com.sun.jna.Native;

public class Julia {
    static {
        Native.register("julia");
        jl_init__threading();
    }

    public static double bubba = Math.random();

    public static native void jl_init__threading();
    public static native void jl_eval_string(String code);
    public static native void jl_atexit_hook(int status);

    public static void main(String args[]) {
        System.out.println("test");
        jl_eval_string("println(\"test from Julia\")");
        jl_eval_string("using JavaCall");
        jl_eval_string("JavaCall.init_current_vm()");
        jl_eval_string("println(\"initialized VM\")");
        jl_eval_string("jlm = @jimport java.lang.Math");
        jl_eval_string("println(jcall(jlm, \"sin\", jdouble, (jdouble,), pi/2))");
        jl_eval_string("jl = @jimport zot.julia.Julia");
        System.out.println("Bubba should be " + bubba);
        jl_eval_string("println(\"bubba: \", jfield(jl, \"bubba\", jdouble))");
        jl_eval_string("println(\"Done with tests\")");
        jl_atexit_hook(0);
    }
}
```
"""
function init_current_vm()
    ppjvm = Array{Ptr{JavaVM}}(undef, 1)
    ppenv = Array{Ptr{JNIEnv}}(undef, 1)
    pnum = Array{Cint}(undef, 1)
    ccall(Libdl.dlsym(libjvm, :JNI_GetCreatedJavaVMs), Cint, (Ptr{Ptr{JavaVM}}, Cint, Ptr{Cint}), ppjvm, 1, pnum)
    global pjvm = ppjvm[1]
    jvm = unsafe_load(pjvm)
    global jvmfunc = unsafe_load(jvm.JNIInvokeInterface_)
    ccall(jvmfunc.GetEnv, Cint, (Ptr{Nothing}, Ptr{Ptr{JNIEnv}}, Cint), pjvm, ppenv, JNI_VERSION_1_8)
    global penv = ppenv[1]
    jnienv = unsafe_load(penv)
    global jnifunc = unsafe_load(jnienv.JNINativeInterface_) #The JNI Function table
end

function destroy()
    if (!isdefined(JavaCall, :penv) || penv == C_NULL)
        throw(JavaCallError("Called destroy without initialising Java VM"))
    end
    assertroottask_or_goodenv()
    res = ccall(jvmfunc.DestroyJavaVM, Cint, (Ptr{Nothing},), pjvm)
    res < 0 && throw(JavaCallError("Unable to destroy Java VM"))
    global penv=C_NULL; global pjvm=C_NULL;
end


text:
module NER

module Convolution

using Merlin
using ProgressMeter
using HDF5

const BACKEND = CPUBackend()

include("BIOES.jl")
include("ner_conv/decoder.jl")
include("ner_conv/nn.jl")

end

module Lstm

using Merlin
using ProgressMeter
using HDF5

const BACKEND = CPUBackend()

include("BIOES.jl")
include("ner_lstm/decoder.jl")
include("ner_lstm/nn.jl")

end

end