text:
# This example requires Cloudy to be installed (it's best to install the master
# branch), which can be done by:
#] add Cloudy#master
using Cloudy
using Cloudy.ParticleDistributions
using Cloudy.KernelTensors

# Import modules
using Distributions  # probability distributions and associated functions
using StatsBase
using LinearAlgebra
using StatsPlots
using Plots
using JLD2
using Random

# Import Calibrate-Emulate-Sample modules
using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.Observations
using EnsembleKalmanProcesses.ParameterDistributionStorage
using EnsembleKalmanProcesses.DataStorage

# Import the module that runs Cloudy
include("DynamicalModel.jl")
using .DynamicalModel

rng_seed = 41
Random.seed!(rng_seed)

homedir = pwd()
figure_save_directory = homedir * "/output/"
data_save_directory = homedir * "/output/"
if ~isdir(figure_save_directory)
    mkdir(figure_save_directory)
end
if ~isdir(data_save_directory)
    mkdir(data_save_directory)
end

###
###  Define the (true) parameters and their priors
###

# Define the parameters that we want to learn
# We assume that the true particle mass distribution is a Gamma
# distribution with parameters N0_true, θ_true, k_true
par_names = ["N0", "θ", "k"]
n_par = length(par_names)
N0_true = 300.0  # number of particles (scaling factor for Gamma distribution)
θ_true = 1.5597  # scale parameter of Gamma distribution
k_true = 0.0817  # shape parameter of Gamma distribution
# Note that dist_true is a Cloudy distribution, not a Distributions.jl
# distribution
ϕ_true = [N0_true, θ_true, k_true]
dist_true = ParticleDistributions.GammaPrimitiveParticleDistribution(ϕ_true...)

###
###  Define priors for the parameters we want to learn
###

# Define constraints
lbound_N0 = 0.4 * N0_true
lbound_θ = 1.0e-1
lbound_k = 1.0e-4
c1 = bounded_below(lbound_N0)
c2 = bounded_below(lbound_θ)
c3 = bounded_below(lbound_k)
constraints = [[c1], [c2], [c3]]

# We choose to use normal distributions to represent the prior distributions of
# the parameters in the transformed (unconstrained) space. i.e log coordinates
d1 = Parameterized(Normal(4.5, 1.0))  #truth is 5.19
d2 = Parameterized(Normal(0.0, 2.0))  #truth is 0.378
d3 = Parameterized(Normal(-1.0, 1.0)) #truth is -2.51
distributions = [d1, d2, d3]

priors = ParameterDistribution(distributions, constraints, par_names)

###
###  Define the data from which we want to learn the parameters
###

data_names = ["M0", "M1", "M2"]
moments = [0.0, 1.0, 2.0]
n_moments = length(moments)

###
###  Model settings
###

# Collision-coalescence kernel to be used in Cloudy
coalescence_coeff = 1 / 3.14 / 4 / 100
kernel_func = x -> coalescence_coeff
kernel = CoalescenceTensor(kernel_func, 0, 100.0)

# Time period over which to run Cloudy
tspan = (0.0, 1.0)

###
###  Generate (artificial) truth samples
###

model_settings_true = ModelSettings(kernel, dist_true, moments, tspan)
G_t = run_dyn_model(ϕ_true, model_settings_true)
n_samples = 100
y_t = zeros(length(G_t), n_samples)
# In a perfect model setting, the "observational noise" represents the
# internal model variability. Since Cloudy is a purely deterministic model,
# there is no straightforward way of coming up with a covariance structure
# for this internal model variability. We decide to use a diagonal
# covariance, with entries (variances) largely proportional to their
# corresponding data values, G_t
Γy = convert(Array, Diagonal([100.0, 5.0, 30.0]))
μ = zeros(length(G_t))

# Add noise
for i in 1:n_samples
    y_t[:, i] = G_t .+ rand(MvNormal(μ, Γy))
end

truth = Observations.Obs(y_t, Γy, data_names)
truth_sample = truth.mean

###
###  Calibrate: Ensemble Kalman Inversion
###

N_ens = 50 # number of ensemble members
N_iter = 8 # number of EKI iterations
# initial parameters: N_par x N_ens
initial_par = construct_initial_ensemble(priors, N_ens; rng_seed)
ekiobj = EnsembleKalmanProcess(initial_par, truth_sample, truth.obs_noise_cov, Inversion(), Δt = 0.1)

# Initialize a ParticleDistribution with dummy parameters. The parameters
# will then be set within `run_dyn_model`
dummy = ones(n_par)
dist_type = ParticleDistributions.GammaPrimitiveParticleDistribution(dummy...)
model_settings = DynamicalModel.ModelSettings(kernel, dist_type, moments, tspan)
# EKI iterations
for n in 1:N_iter
    θ_n = get_u_final(ekiobj)
    # Transform parameters to physical/constrained space
    ϕ_n = mapslices(x -> transform_unconstrained_to_constrained(priors, x), θ_n; dims = 1)
    # Evaluate forward map
    G_n = [run_dyn_model(ϕ_n[:, i], model_settings) for i in 1:N_ens]
    G_ens = hcat(G_n...)  # reformat
    EnsembleKalmanProcessModule.update_ensemble!(ekiobj, G_ens)
end

# EKI results: Has the ensemble collapsed toward the truth?
θ_true = transform_constrained_to_unconstrained(priors, ϕ_true)
println("True parameters (unconstrained): ")
println(θ_true)

println("\nEKI results:")
println(mean(get_u_final(ekiobj), dims = 2))

u_stored = get_u(ekiobj, return_array = false)
g_stored = get_g(ekiobj, return_array = false)
@save data_save_directory * "parameter_storage_eki.jld2" u_stored
@save data_save_directory * "data_storage_eki.jld2" g_stored

#plots
gr(size = (1800, 600))

u_init = get_u_prior(ekiobj)
for i in 1:N_iter
    u_i = get_u(ekiobj, i)

    p1 = plot(u_i[1, :], u_i[2, :], seriestype = :scatter, xlims = extrema(u_init[1, :]), ylims = extrema(u_init[2, :]))
    plot!(
        p1,
        [θ_true[1]],
        xaxis = "u1",
        yaxis = "u2",
        seriestype = "vline",
        linestyle = :dash,
        linecolor = :red,
        label = false,
        title = "EKI iteration = " * string(i),
    )
    plot!(p1, [θ_true[2]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum")

    p2 = plot(u_i[2, :], u_i[3, :], seriestype = :scatter, xlims = extrema(u_init[2, :]), ylims = extrema(u_init[3, :]))
    plot!(
        p2,
        [θ_true[2]],
        xaxis = "u1",
        yaxis = "u2",
        seriestype = "vline",
        linestyle = :dash,
        linecolor = :red,
        label = false,
        title = "EKI iteration = " * string(i),
    )
    plot!(p2, [θ_true[3]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum")

    p3 = plot(u_i[3, :], u_i[1, :], seriestype = :scatter, xlims = extrema(u_init[3, :]), ylims = extrema(u_init[1, :]))
    plot!(
        p3,
        [θ_true[3]],
        xaxis = "u1",
        yaxis = "u2",
        seriestype = "vline",
        linestyle = :dash,
        linecolor = :red,
        label = false,
        title = "EKI iteration = " * string(i),
    )
    plot!(p3, [θ_true[1]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = "optimum")

    p = plot(p1, p2, p3, layout = (1, 3))
    display(p)
    sleep(0.5)
end


text:
using ModelingToolkit, LinearAlgebra, OrdinaryDiffEq, Test
MT = ModelingToolkit

# Repressilator model
@parameters t α₀ α K n δ β μ
@variables m(t) P(t) R(t)
rxs = [
       Reaction(α₀, nothing, [m]),
       Reaction(α / (1 + (R/K)^n), nothing, [m]),
       Reaction(δ, [m], nothing),
       Reaction(β, [m], [m,P]),
       Reaction(μ, [P], nothing)
      ]

specs = [m,P,R]
pars  = [α₀,α,K,n,δ,β,μ]
@named rs = ReactionSystem(rxs, t, specs, pars)

# using ODESystem components
@named os₁ = convert(ODESystem, rs)
@named os₂ = convert(ODESystem, rs)
@named os₃ = convert(ODESystem, rs)
connections = [os₁.R ~ os₃.P,
               os₂.R ~ os₁.P,
               os₃.R ~ os₂.P]
@named connected = ODESystem(connections, t, [], [], systems=[os₁,os₂,os₃])
oderepressilator = structural_simplify(connected)

pvals = [os₁.α₀ => 5e-4,
         os₁.α => .5,
         os₁.K => 40.0,
         os₁.n => 2,
         os₁.δ => (log(2)/120),
         os₁.β => (20*log(2)/120),
         os₁.μ => (log(2)/600),
         os₂.α₀ => 5e-4,
         os₂.α => .5,
         os₂.K => 40.0,
         os₂.n => 2,
         os₂.δ => (log(2)/120),
         os₂.β => (20*log(2)/120),
         os₂.μ => (log(2)/600),
         os₃.α₀ => 5e-4,
         os₃.α => .5,
         os₃.K => 40.0,
         os₃.n => 2,
         os₃.δ => (log(2)/120),
         os₃.β => (20*log(2)/120),
         os₃.μ => (log(2)/600)]
u₀    = [os₁.m => 0.0, os₁.P => 20.0, os₂.m => 0.0, os₂.P => 0.0, os₃.m => 0.0, os₃.P => 0.0]
tspan = (0.0, 100000.0)
oprob = ODEProblem(oderepressilator, u₀, tspan, pvals)
sol = solve(oprob, Tsit5())

# hardcoded network
function repress!(f, y, p, t)
    α = p.α; α₀ = p.α₀; β = p.β; δ = p.δ; μ = p.μ; K = p.K; n = p.n
    f[1] = α / (1 + (y[6] / K)^n) - δ * y[1] + α₀
    f[2] = α / (1 + (y[4] / K)^n) - δ * y[2] + α₀
    f[3] = α / (1 + (y[5] / K)^n) - δ * y[3] + α₀
    f[4] = β * y[1] - μ * y[4]
    f[5] = β * y[2] - μ * y[5]
    f[6] = β * y[3] - μ * y[6]
    nothing
end
ps = (α₀=5e-4, α=.5, K=40.0, n=2, δ=(log(2)/120), β=(20*log(2)/120), μ=(log(2)/600))
u0 = [0.0,0.0,0.0,20.0,0.0,0.0]
oprob2 = ODEProblem(repress!, u0, tspan, ps)
sol2 = solve(oprob2, Tsit5())
tvs = 0:1:tspan[end]

indexof(sym,syms) = findfirst(isequal(sym),syms)
i = indexof(os₁.P, states(oderepressilator))
@test all(isapprox(u[1],u[2],atol=1e-4) for u in zip(sol(tvs, idxs=2), sol2(tvs, idxs=4)))

# using ReactionSystem components

# @named rs₁ = ReactionSystem(rxs, t, specs, pars)
# @named rs₂ = ReactionSystem(rxs, t, specs, pars)
# @named rs₃ = ReactionSystem(rxs, t, specs, pars)
# connections = [rs₁.R ~ rs₃.P,
#                rs₂.R ~ rs₁.P,
#                rs₃.R ~ rs₂.P]
# @named csys = ODESystem(connections, t, [], [])
# @named repressilator = ReactionSystem(t; systems=[csys,rs₁,rs₂,rs₃])
# @named oderepressilator2 = convert(ODESystem, repressilator)
# sys2 = structural_simplify(oderepressilator2)  # FAILS currently

ext: jl


text:
include("valueiteration.jl")
include("mcts.jl")


text:
mutable struct Input
    A::Matrix{}
    b::Vector{}
    c::Vector{}
    standard::Bool
    n::Int
    m::Int
    rho::Float64
    alpha::Float64
    tol::Float64
    max_iter::Int
    verbose::Bool
    solver::Int
    crossover::Bool
    crossover_tol::Float64
end

Base.@kwdef mutable struct MidTerm
    termination_status::Int = 0
    iter::Int = 0
    d::Vector{} = []
    base::Vector{Int} = []
    nbase::Vector{Int} = []
    i::Int = 0
    j::Int = 0
    z::Float64 = 0.0
    x::Vector{Float64} = []
    red_cost::Vector{Float64} = []
end

struct OutputIP
    x::Vector{}
    s::Vector{}
    p::Vector{}
    X::Matrix{}
    S::Matrix{}
    P::Matrix{}
    gap::Vector{}
    mu::Float64
    z::Float64
    termination_status::Int
    iter::Integer
end

struct OutputSimplex
    x::Vector{Float64}
    z::Float64
    termination_status::Int
    base::Vector{Int}
    nbase::Vector{Int}
end


text:
using Fakerjl

function phone()
    return Fakerjl.genericgenerator("person_phone")
end

function email()
    return Fakerjl.genericgenerator("person_email")
end


text:
using Tabulars
using Base.Test

include("DictTabular.jl")
include("ArrayTabular.jl")
include("TupleTabular.jl")
include("StructTabular.jl")
include("PermutedDimsTabular.jl")
include("SubTabular.jl")
include("show.jl")


text:
export
		get_number_of_dimensions,
		get_size,
		get_size_of_dimension

function get_at{T <: Integer}( dmat::DMatrix, coord::Array{T,1} )
	retval = ccall( (:dmatrix_GetAtCoord, LIB_SMILE), Float64, (Ptr{Void},Ptr{Int32},UInt32), dmat.ptr, Int32(coord), length(coord) )
	return retval
end

# the order of the probabilities is given by considering the state of
# the first parent as "most significant", then the second parent, then
# the third (and so on), and finally considering the coordinate of the
# node itself as the least significant one.

# Example: consider a node C with two parents, A and B
# A has 2 instances: A1 & A2
# B has 3 instances: B1 & B2 & B3
# C has 2 instances: C1 & C2
# the first index, 0, corresponds to [A1 B1 C1]
# the next index,  1, corresponds to [A1 B1 C2]
# then             2,       ->       [A1 B2 C1]
#                  3,       ->       [A1 B2 C2]
#                  4,       ->       [A1 B3 C1]
#     ...
#                  11       ->       [A2 B3 C2]

Base.getindex(dmat::DMatrix, index::Integer) =
    ccall( (:dmatrix_GetAtInd, LIB_SMILE), Float64, (Ptr{Void},Int32), dmat.ptr, index )
Base.setindex!(dmat::DMatrix, value::Real, ind::Integer) =
    ccall( (:dmatrix_SetAtInd, LIB_SMILE), Void, (Ptr{Void},Int32,Float64), dmat.ptr, ind, value )

get_number_of_dimensions( dmat::DMatrix ) =
	ccall( (:dmatrix_GetNumberOfDimensions, LIB_SMILE), Int32, (Ptr{Void},), dmat.ptr )

function get_size( dmat::DMatrix )
	retval = ccall( (:dmatrix_GetSize, LIB_SMILE), Int32, (Ptr{Void},), dmat.ptr )
	return retval
end

function get_size_of_dimension( dmat::DMatrix, aDimension::Integer )
	# indexing starts at 0
	# returns 0 if you pick a bad value for aDimension
	retval = ccall( (:dmatrix_GetSizeOfDimension, LIB_SMILE), Int32, (Ptr{Void},Int32), dmat.ptr, aDimension )
	return retval
end

coordinates_to_index(dmat::DMatrix, arr::IntArray) =
	ccall( (:dmatrix_CoordinatesToIndex, LIB_SMILE), Int32, (Ptr{Void},Ptr{Void}), dmat.ptr, arr.ptr )


text:
@testset "DataF1 tests" begin show_testset_section() #Scope for test data

using MDDatasets
using Test

#==Input data
===============================================================================#
xingsall = CrossType(:all)
xingsrise = CrossType(:rise)
xingsfall = CrossType(:fall)

d1 = DataF1(1:10.0)
d2 = xshift(d1, 4.5) + 12
d3 = d1 + 12
d4 = DataF1(d1.x, d1.y[end:-1:1])
d9 = xshift(d1, 100)

#Required mostly to test cross function:
y = [0,0,0,0,1,-3,4,5,6,7,-10,-5,0,0,0,-5,-10,10,-3,1,-4,0,0,0,0,0,0]
d10=DataF1(collect(1:length(y)), y)

#==Basic Tests
===============================================================================#
@testset "Basic operations on DataF1" begin show_testset_description()
	r1 = d1+d2
	r2 = d1+d3
	#r3 = d1+d4
	@test r1.x == [1.0, 2.0, 3.0, 4.0, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.5, 12.5, 13.5, 14.5]
	@test r1.y == [1.0, 2.0, 3.0, 4.0, 5.0, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 18.0, 19.0, 20.0, 21.0, 22.0]
	@test r2.x == d1.x
	@test r2.y == [14.0, 16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0]
	@test length(d1) == 10
	@test length(d2) == 10
	@test length(d3) == 10
	@test length(r1) == 20 #Shifted x values.
	@test length(r2) == 10
end

@testset "Functionnality of clip()" begin show_testset_description()
	_d1 = clip(d1, xmin=1, xmax=9)
	@test _d1.x == collect(1:9.0)
	@test _d1.y == _d1.x

	_d1 = clip(d1, xmin=2, xmax=5)
	@test _d1.x == collect(2:5.0)
	@test _d1.y == _d1.x

	_d1 = clip(d1, 3:10)
	@test _d1.x == collect(3:10.0)
	@test _d1.y == _d1.x

	_d1 = clip(d1, xmin=3)
	@test _d1.x == collect(3:d1.x[end])
	@test _d1.y == _d1.x

	_d1 = clip(d1, xmax=8.5)
	@test _d1.x == [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 8.5]
	@test _d1.y == _d1.x

	_d1 = clip(d1, xmin=1.5, xmax=5.5)
	@test _d1.x == [1.5, 2.0, 3.0, 4.0, 5.0, 5.5]
	@test _d1.y == _d1.x

	_d3 = clip(d3, xmin=1.5, xmax=5.5)
	@test _d3.x == [1.5, 2.0, 3.0, 4.0, 5.0, 5.5]
	@test _d3.y == [13.5, 14.0, 15.0, 16.0, 17.0, 17.5]
end

@testset "Functionnality of x/ycross()" begin show_testset_description()
	rtol = 1e-4

	#Test ycross:
	_x = xcross(d10, allow=xingsall)
	@test _x.y ≈ [4.0, 5.25, 6.42857, 10.4118, 14.0, 17.5, 18.7692, 19.75, 20.2, 22.0] rtol=rtol
	@test _x.x == _x.y
	_x = xcross(d10)
	@test _x.y ≈ [5.25, 6.42857, 10.4118, 17.5, 18.7692, 19.75, 20.2] rtol=rtol
	@test _x.x == _x.y
	_x = xcross(d10, allow=xingsrise)
	@test _x.y ≈ [6.42857, 17.5, 19.75] rtol=rtol
	@test _x.x == _x.y
	_x = xcross(d10, allow=xingsfall)
	@test _x.y ≈ [5.25, 10.4118, 18.7692, 20.2] rtol=rtol
	@test _x.x == _x.y
	_x = xcross(d10-.5)
	@test _x.y ≈ [4.5, 5.125, 6.5, 10.3824, 17.525, 18.7308, 19.875, 20.1] rtol=rtol
	@test _x.x == _x.y

	#Test ycross:
	_x = ycross(d10, .5, allow=xingsall)
	@test _x.x ≈ [4.5, 5.125, 6.5, 10.3824, 17.525, 18.7308, 19.875, 20.1] rtol=rtol
	@test _x.y ≈ zeros(length(_x.y)) .+ 0.5 rtol=rtol
	_x2 = ycross(d10, .5)
	@test _x.x ≈ _x2.x rtol=rtol
	@test _x.y ≈ _x2.y rtol=rtol
	_x = ycross(d10, .5, allow=xingsrise)
	@test _x.x ≈ [4.5, 6.5, 17.525, 19.875] rtol=rtol
	@test _x.y ≈ zeros(length(_x.y)) .+ 0.5 rtol=rtol
	_x = ycross(d10, .5, allow=xingsfall)
	@test _x.x ≈ [5.125, 10.3824, 18.7308, 20.1] rtol=rtol
	@test _x.y ≈ zeros(length(_x.y)) .+ 0.5 rtol=rtol

	#Test xcross1/ycross1:
	@test xcross1(d10, n=1) == 5.25
	@test xcross1(d10-.5, n=2) == 5.125
	@test ycross1(d10, .5, n=3) ≈ 0.5 rtol=rtol
end

@testset "Functionnality of meas() interface" begin show_testset_description()
	_xref = xcross(d10, allow=xingsrise)
	_x = meas(:xcross, d10, allow=xingsrise)
	@test _x.x == _xref.x
	@test _x.y == _xref.y

	_x = meas(:xcross, Event, d10, allow=xingsrise)
	@test _x.x == collect(1:length(_x.x))
	@test _x.y == _xref.y
end

@testset "Functionnality value()" begin show_testset_description()
	for i in 1:length(d3)
		@test d3.y[i] == value(d3, x=d3.x[i])
	end
end

@testset "Functionnality of ensure()" begin show_testset_description()
	function fail_to_ensure()
		@warn("TODO: fix bad ensure-do syntax - sounds like body should execute when predicate is true")

		#Test ensure-do syntax:
		ensure(false) do
			#Allows user to build a more complex error object...
			ArgumentError("Predicate is false.")
		end
		return
	end

	#Simple syntax:
	@test ensure(true, "No error") == nothing
	@test_throws SystemError ensure(false, SystemError("Some system error"))
	@test_throws ArgumentError fail_to_ensure()
end

end


text:
export getService, shutdown

"""
    getService(server)

Returns the service that this ServiceServer will advertise on.
"""
getService(s::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::ServiceServer")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::ServiceServer")},(false, false, false)},N} where N}) = unsafe_string(@cxx s->getService())

"""
    shutdown(server)

Unadvertise the service associated with this ServiceServer.
"""
shutdown(s::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::ServiceServer")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::ServiceServer")},(false, false, false)},N} where N}) = @cxx s->shutdown()


text:
include("../src/spline_module.jl")

module stest
using ..spline
using PyPlot
PyPlot.pygui(true)

function test_discrete()
    x = [1.0, 1.5, 2.1, 3., 3.1, 3.4, 4.,  5.5, 6, 7, 8]
    y = [10.0, 10.2, 11,  8,  5,   9,   13,  12,  10, 10, 10]

    natural_spline = spline.fit_natural(x, y)

    xq = collect( range(x[1], x[end], 10000) )
    yq = spline.interpolate(natural_spline, xq)

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot(x, y, marker = :o)
    PyPlot.plot(xq, yq)
end

function test_sin()

    x = collect( range(0, 2*pi, 15) )
    y = sin.(x)

    natural_spline = spline.fit_natural(x, y)

    xq = collect( range(x[1], x[end], 1000) )
    yq = spline.interpolate(natural_spline, xq)

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot(x, y, marker = :o)
    PyPlot.plot(xq, yq)
end

function test_lim_sin()

    x = collect( range(0, 2*pi, 15) )
    y = sin.(x)
    y[y .<= 0] .= 0.0

    natural_spline = spline.fit_natural(x, y)

    xq = collect( range(x[1], x[end], 1000) )
    yq = spline.interpolate(natural_spline, xq)

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot(x, y, marker = :o)
    PyPlot.plot(xq, yq)
end

function test_clamped()
    x = [1.0, 1.5, 2.1, 3., 3.1, 3.4, 4.,  5.5, 6]
    y = [10.0, 10.2, 11,  8,  5,   9,   13,  12,  10]
    g1 = -10
    g2 = +20
    clamped_spline = spline.fit_clamped(x, y, g1, g2)

    xq = collect( range(x[1], x[end], 1000) )
    yq = spline.interpolate(clamped_spline, xq)
    dyq = spline.diff_interpolate(clamped_spline, xq)

    PyPlot.figure()

    PyPlot.subplot(2,1,1)
    PyPlot.grid()
    PyPlot.plot(x, y, marker = :o)
    PyPlot.plot(xq, yq)

    PyPlot.subplot(2,1,2)
    PyPlot.grid()
    PyPlot.plot(xq, dyq)
end

x = collect( range(0, 2*2*pi, 15) )
y = sin.(x)
# y[y .>= 0.5]  .= 0.5
# y[y .<= -0.5] .= -0.5

akima_spline = spline.fit_modified_akima(x, y)
xq = collect( range(x[1], x[end], 100) )
yq = spline.interpolate(akima_spline, xq)
dyq = spline.diff_interpolate(akima_spline, xq)

PyPlot.figure()
PyPlot.grid()
PyPlot.plot(x, y, marker = :o)
PyPlot.plot(xq, yq)

PyPlot.figure()
PyPlot.grid()
PyPlot.plot(xq, dyq)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
export viz_policy

function viz_policy(nnet_folder::AbstractString="",table_folder::AbstractString="",vnum::Int64=4, hu::Int64=45, epoch::Int64=200)

    # Make sure one of the folder paths is given
    if nnet_folder=="" && table_folder==""
        println("Need at least either a valid nnet or table folder path!")
        return
    end

    # Important variables
    currentPra = -1
    nnet = []
    policy = []
    grid = RectangleGrid(relhs,dh0s,dh1s,taus)

    # Colors
    ra_1 = RGB(1.,1.,1.) # white
    ra_2 = RGB(.0,1.0,1.0) # cyan
    ra_3 = RGB(144.0/255.0,238.0/255.0,144.0/255.0) # lightgreen
    ra_4 = RGB(30.0/255.0,144.0/255.0,1.0) # dodgerblue
    ra_5 = RGB(0.0,1.0,.0) # lime
    ra_6 = RGB(0.0,0.0,1.0) # blue
    ra_7 = RGB(34.0/255.0,139.0/255.0,34.0/255.0) # forestgreen
    ra_8 = RGB(0.0,0.0,128.0/255.0) # navy
    ra_9 = RGB(0.0,100.0/255.0,0.0) # darkgreen
    colors = [ra_1;ra_2;ra_3;ra_4;ra_5;ra_6;ra_7;ra_8;ra_9]
    bg_colors = [RGB(1.0,1.0,1.0)]

    # Create scatter plot classes for color key
    sc_string = "{"
    for i=1:9
        define_color("ra_$i",  colors[i])
        if i==1
            sc_string *= "ra_$i={mark=square, style={black, mark options={fill=ra_$i}, mark size=6}},"
        else
            sc_string *= "ra_$i={style={ra_$i, mark size=6}},"
        end
    end

    # Color key as a scatter plot
    sc_string=sc_string[1:end-1]*"}"
    xx = [-1.5,-1.5,-1.5, -1.5, -1.5, -1.5, -1.5, 0.4 ,.4,]
    yy = [1.65,1.15,0.65, 0.15, -0.35, -0.85, -1.35, 1.65, 1.15]
    zz = ["ra_1","ra_2","ra_3","ra_4","ra_5","ra_6","ra_7","ra_8","ra_9"]
    sc = string(sc_string)

    # Set up interactive display
    @manipulate for nbin = 100,
        savePlot = [false,true],
        xmin = 0.0,
        xmax = 40.0,
        ymin = -1000.0,
        ymax = 1000.0,
        dh0 = 0.0,
        dh1 = 0.0,
        pra = action_names

        # Get previous RA index
        pra = find(pra.==action_names)[1]
        if pra != currentPra
            if table_folder !=""
                alpha = h5open(@sprintf("%s/VertCAS_TrainingData_v2_%02d.h5",table_folder,pra), "r") do file
                    read(file, "y")
                end
                alpha = alpha';
                policy = read_policy(reshape(actions,(1,9)), alpha)
            end

            if nnet_folder != ""
                if epoch>0
                    nnet = NNet(@sprintf("%s/VertCAS_pra%02d_v%d_%dHU_%03d.nnet",nnet_folder,pra,vnum,hu,epoch))
                else
                    nnet = NNet(@sprintf("%s/VertCAS_pra%02d_v%d_%dHU.nnet",nnet_folder,pra,vnum, hu))
                end
            end
            currentPra = pra
        end

        # Evaluate network at pixel points
        if nnet_folder != ""
            inputsNet = hcat([[relh,dh0,dh1,tau] for tau=linspace(xmin, xmax,nbin) for relh=linspace(ymin,ymax,nbin)]...)
            q_nnet = evaluate_network_multiple(nnet,inputsNet)
        end

        # Q Table Heat Map
        function get_heat1(x::Float64, y::Float64)
            tau = x
            relh = y
            qvals = evaluate(policy, get_belief([relh,dh0,dh1,tau],grid,false))
            return actions[indmax(qvals)]
        end # function get_heat1

        ind = 1
        #Neural Net Heat Map
        function get_heat2(x::Float64, y::Float64)
            qvals = q_nnet[:,ind]
            ind +=1
            return actions[indmax(qvals)]
        end # function get_heat2

        #Plot table or neural network policies if possible
        g = GroupPlot(2, 2, groupStyle = "horizontal sep=4cm, vertical sep=2.5cm")
        if table_folder!=""
            push!(g, Axis([
                Plots.Image(get_heat1, (xmin, xmax), (ymin, ymax), zmin = 1, zmax = 9,
                xbins = nbin, ybins = nbin, colormap = ColorMaps.RGBArrayMap(colors), colorbar=false),
                ], xmin=xmin, xmax=xmax, ymin=ymin,ymax=ymax, width="10cm", height="8cm",
                   xlabel="Tau (s)", ylabel="Relative Alt (ft)", title="Original Table Advisories"))
        end
        if nnet_folder!=""
            push!(g, Axis([
                Plots.Image(get_heat2, (xmin, xmax), (ymin, ymax), zmin = 1, zmax = 9,
                            xbins = nbin, ybins = nbin, colormap = ColorMaps.RGBArrayMap(colors), colorbar=false),
                ], xmin=xmin, xmax=xmax, ymin=ymin,ymax=ymax,width="10cm", height="8cm",
                   xlabel="Tau (s)", ylabel="Relative Alt (ft)", title="Neural Network Advisories"))
        end

        # Save policy to a tex file to be used in papers
        if savePlot
            PGFPlots.save("PolicyPlot.tex", g, include_preamble=true)
        else

            # Create Color Key
            f = (x,y)->x # Dummy function for background white image
            push!(g, Axis([
                Plots.Image(f, (-2,2), (-2,2),colormap = ColorMaps.RGBArrayMap(bg_colors),colorbar=false),
                Plots.Scatter(xx, yy, zz, scatterClasses=sc),
                Plots.Node("RA 1: COC ",0.15,0.915,style="black,anchor=west", axis="axis description cs"),
                Plots.Node("RA 2: DNC ",0.15,0.790,style="black,anchor=west", axis="axis description cs"),
                Plots.Node("RA 3: DND",0.15,0.665,style="black,anchor=west", axis="axis description cs"),
                Plots.Node("RA 4: DES15000",0.15,0.540,style="black,anchor=west", axis="axis description cs"),
                Plots.Node("RA 5: CL1500 ",0.15,0.415,style="black,anchor=west", axis="axis description cs"),
                Plots.Node("RA 6: SDES1500",0.15,0.290,style="black,anchor=west", axis="axis description cs"),
                Plots.Node("RA 7: SCL1500",0.15,0.165,style="black,anchor=west", axis="axis description cs"),
                Plots.Node("RA 8:  SDES2500",0.63,0.915,style="black,anchor=west", axis="axis description cs"),
                Plots.Node("RA 9:  SCL2500",0.63,0.790,style="black,anchor=west", axis="axis description cs"),
                ],width="10cm",height="8cm", hideAxis =true, title="KEY"))
        end
    end
end

max_stars_repo_name: alisamiii/VerticalCAS

max_issues_repo_name: alisamiii/VerticalCAS

max_forks_repo_name: bethelhall/CoCoCAS


text:
using MetabolicFluxAnalysis
using CoordinateTransformations
using OpenCL

    # current tests..
    fluxi_1_f = open("testdata/test_01.fluxi")
    fluxi_1_s = readstring(fluxi_1_f)
    fluxi1 = MetabolicFluxAnalysis.parseFluxi( fluxi_1_s )

    fluxi1.net.fluxes_db[fluxi1.net.fluxes[1]].transition
    emudecomp_1 = MetabolicFluxAnalysis.compute_emu_decomp( fluxi1.net , [ fluxi1.measurements["MassDistribution_01"][1].sass ] )
    emu_graph_1 = MetabolicFluxAnalysis.create_emu_graph(emudecomp_1)

    configured_sas = MetabolicFluxAnalysis.get_configured_reactands(fluxi1,"SC_010")
    #fluxi_ast_1 = MetabolicFluxAnalysis.create_ast_from_emu_graph( emu_graph_1 , configured_sas )
    fluxi_ast_1 = MetabolicFluxAnalysis.create_ast_from_emu_graph( emu_graph_1 )

    unrolled_ast = Vector{Any}()
    MetabolicFluxAnalysis.ast_traverse_df!( fluxi_ast_1 , unrolled_ast)

    # get the target semus from the fluxi measurements..
    target_smrs = MetabolicFluxAnalysis.construct_all_mass_ratios_for_measurement(fluxi1,"MassDistribution_01")
    julia_fw_sim = MetabolicFluxAnalysis.ast_to_juliacode_01(fluxi1.net,fluxi_ast_1,target_smrs)
    #julia_fw_sim_fast = MetabolicFluxAnalysis.ast_to_juliacode_01(fluxi1.net,fluxi_ast_1,target_smrs;fast_lu=true)

    ocl_fw_sim = MetabolicFluxAnalysis.ast_to_opencl_01(fluxi1.net,fluxi_ast_1,target_smrs)

print("\n\n\nDA KERNEL..\n\n")
print(ocl_fw_sim.code)

kernel_str = ocl_fw_sim.code

#semu_data = convert( Array{Float32,1} , repmat([0.00 ; 1.00; 1.00; 0.00; 0.00; 1.0 ],1))
semu_data = convert( Array{Float32,1} , repmat([ 0.0 ; 1.0 ; 1.0 ; 0.0 ; 0.0 ; 1.0 ; 0.0 ; 0.0 ; 1.0 ; 0.0 ; 0.0 ],1))
flux_data = convert( Array{Float32,1} , repmat([100.; 110 ; 50 ; 20 ; 20 ; 80 ],4) )

a         = Array{Float32,1}(100)

ocl_devices = OpenCL.cl.devices()
ocl_device  = ocl_devices[2]
ocl_context =  MetabolicFluxAnalysis.OCLContext(ocl_device)

#device, ctx, queue = OpenCL.cl.create_compute_context()
device = ocl_context.ocl_device
ctx    = ocl_context.ocl_ctx
queue  = ocl_context.ocl_queue

semu_buff = OpenCL.cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=semu_data)
flux_buff = OpenCL.cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=flux_data)
c_buff = OpenCL.cl.Buffer(Float32, ctx, :w, length(a))
p = OpenCL.cl.Program(ctx, source=kernel_str)
OpenCL.cl.build!(p)

#p = cl.Program(ctx, source=kernel_02) |> cl.build!
print(kernel_str)
#cl.CL_device_info(1)

k = cl.Kernel(p, "fwsim")
#queue(k, [8;8], [8;8], semu_buff, flux_buff, c_buff)
queue(k, [8,8], [8,8], semu_buff, flux_buff, c_buff)

queue(k, [32,8], [8,8], semu_buff, flux_buff, c_buff)

#queue(k, [8;32], [8;8], semu_buff, flux_buff, c_buff)

#cl.api.clGetKernelWorkGroupInfo(cl.CL_KERNEL_WORK_GROUP_SIZE)
r = cl.read(queue, c_buff)

#cl.work_group_info( k , cl.CL_KERNEL_WORK_GROUP_SIZE ,device )
#cl.work_group_info( k , cl.CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE  ,device )

julia_fw_sim

##

test_A = [ -150.000   0.000   0.000  50.000   0.000
  0.000  -150.000   0.000   0.000  50.000
 20.000   0.000  -20.000   0.000   0.000
 110.000  20.000   0.000  -130.000   0.000
  0.000  110.000  20.000   0.000  -130.000]
test_B = [  0.000  -100.000
 -100.000   0.000
  0.000   0.000
  0.000   0.000
  0.000   0.000]
test_Y = [  1.000   0.000
  0.000   1.000]

test_X = test_A \ (test_B*test_Y)

lu_result = lu(test_A)

 (test_B*test_Y)

after_l_solve = lu_result[1] \ (test_B*test_Y)

after_u_solve = lu_result[2] \ (after_l_solve)

##

mA = [ -150.000  50.000
 110.000  -130.000]

mB = [ -100.000   0.000
  0.000  -20.000]

mY = [  0.000   1.000   0.000
  0.000   0.000   0.000]

mX2 = mA \ (mB * mY)


text:
using KrakenAPI
using Base.Test

@test 1+1 == 2


text:
type TensorDesc
  ptr::Ptr{Void}

  function TensorDesc{T}(a::CuArray{T})
    csize = Cint[size(a,i) for i=ndims(a):-1:1]
    cstrides = Cint[stride(a,i) for i=ndims(a):-1:1]
    p = Ptr{Void}[0]
    cudnnCreateTensorDescriptor(p)
    desc = new(p[1])
    finalizer(desc, cudnnDestroyTensorDescriptor)
    cudnnSetTensorNdDescriptor(desc, datatype(T), ndims(a), csize, cstrides)
    desc
  end
end

Base.unsafe_convert(::Type{Ptr{Void}}, td::TensorDesc) = td.ptr


text:
using Optim
export distxy!, demo_distxy

# Grab the 2n coordinates as a single vector
function vector_out(X::GraphEmbedding)
    vv = vlist(X.G)
    n = length(vv)

    x = zeros(2*n)
    for k=1:n
        v = X.xy[vv[k]]
        x[2*k-1] = v[1]
        x[2*k]   = v[2]
    end

    return x
end

# Inverse of vector_out
function vector_in!(X::GraphEmbedding, x::Vector)
    vv = vlist(X.G)
    n = length(vv)

    for k=1:n
        v = vv[k]
        X.xy[v] = x[2*k-1 : 2*k]
    end

    return
end

function distxy!(X::GraphEmbedding,
                 nits::Integer=0, verbose::Bool=true)
    D = dist_matrix(X.G)
    x0 = vector_out(X)

    function score(x::Vector)
        nn = length(x)
        n = NV(X.G)
        s = 0.

        for u=1:(n-1)
            pu = x[2*u-1 : 2*u]
            for v = (u+1):n
                pv = x[2*v-1 : 2*v]
                duv = D[u,v]
                term = (duv - norm(pu-pv))^2 / duv^(1.75)
                s += term
            end
        end
        return s
    end

    if nits > 0
        res = optimize(score,x0,iterations=nits)
    else
        res = optimize(score,x0)
    end

    x1 = res.minimizer
    vector_in!(X,x1)

    if verbose
        msg = "$(res.f_calls) function calls\t" *
          "score = $(score(x1))"
        @info msg
    end

    return score(x1)
end

function distxy!(X::GraphEmbedding,
                 tol::Float64=0.001, verb::Bool=true)
    x0 = Inf
    x1 = 0.0
    # if verb
    #     tic()
    # end
    while true
        x1 = distxy!(X,0,verb)
        if abs(x1-x0)/x0 < tol
            break
        end
        x0 = x1
    end
    # if verb
    #     toc()
    # end
    return x1
end

"""
`demo_distxy(G,tol=1e-3)` presents an animation showing the evolving
drawing found by `distxy!`.
"""
function demo_distxy(G::SimpleGraph=BuckyBall(), tol=1e-3)
    return demo_distxy(GraphEmbedding(G), tol)
end

function demo_distxy(X::GraphEmbedding,tol::Real=1e-3)
    tic()
    x0 = 1.e10
    figure(1)
    while true
        clf()
        draw(X)
        title("Embedding with distxy!")
        x1 = distxy!(X,0)
        if abs(x1-x0)/x0 < tol
            break
        end
        x0 = x1
    end
    title("Finished")
    toc()
    return X
end


text:
using SafeTestsets

@safetestset "Maths" begin
    include("./mathsTests.jl")
end

@safetestset "Structs" begin
    include("./structTests.jl")
end

@safetestset "IO" begin
    include("./ioTests.jl")
end

max_stars_repo_name: dcelisgarza/Lifting


text:
# Julia wrapper for header: /usr/include/libavutil/pixelutils.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_pixelutils_get_sad_fn

function av_pixelutils_get_sad_fn(w_bits::Integer, h_bits::Integer, aligned::Integer, log_ctx)
    ccall((:av_pixelutils_get_sad_fn, libavutil), av_pixelutils_sad_fn, (Cint, Cint, Cint, Ptr{Cvoid}), w_bits, h_bits, aligned, log_ctx)
end


text:
using StructJuMP, JuMP
using StructJuMPSolverInterface

include("select_solver.jl")

#############
# A sample model
#############
scen = 2
m = StructuredModel(num_scenarios=scen)
@variable(m, x[1:2])
@NLobjective(m, Min, x[1]^2+x[2]^2)
@NLconstraint(m, x[1] * x[2] <= 100)

for i in 1:scen
    bl = StructuredModel(parent=m)
    @variable(bl, y[1:2])
    @NLconstraint(bl, x[2]*x[1] + x[1]*y[1]*y[2] <= 10)
    @NLobjective(bl, Min, (y[1]+y[2])^2)
end

structJuMPSolve(m)

getVarValue(m)


text:
abstract type WrapperNode end

struct Definition{S<:Spec} <: WrapperNode
    spec::S
    p::Dict
end

const EnumDefinition = Definition{SpecEnum}
const BitmaskDefinition = Definition{SpecBitmask}
const ConstantDefinition = Definition{SpecConstant}

struct StructDefinition{HL,S<:Union{SpecStruct,SpecUnion}} <: WrapperNode
    spec::S
    p::Dict
end

StructDefinition{T}(spec, p) where {T} = StructDefinition{T,typeof(spec)}(spec, p)

struct HandleDefinition <: WrapperNode
    spec::SpecHandle
    p::Dict
end

struct Constructor{T,X} <: WrapperNode
    p::Dict
    to::T
    from::X
end

struct Documented{W<:WrapperNode} <: WrapperNode
    def::W
    p::Dict
end

Documented(def::WrapperNode) = Documented(def, "")
function Documented(def::WrapperNode, doc::AbstractString)
    doc = lstrip(doc, '\n')
    Documented(def, docstring(to_expr(def), doc))
end

abstract type MethodDefinition <: WrapperNode end

struct Convert{A,B} <: MethodDefinition
    T::A
    x::B
    p::Dict
end

struct GetProperty{D} <: MethodDefinition
    def::D
    p::Dict
end

struct APIFunction{S} <: MethodDefinition
    spec::S
    with_func_ptr::Bool
    p::Dict
end

struct Parent <: MethodDefinition
    def::HandleDefinition
    p::Dict
end

VulkanSpec.has_parent(def::HandleDefinition) = has_parent(def.spec)

struct StructureType <: MethodDefinition
    spec::SpecStruct
    ex::Expr
end

abstract type TypeMapping <: MethodDefinition end

struct HLTypeMapping <: TypeMapping
    spec::Union{SpecStruct,SpecUnion}
    ex::Expr
end

struct CoreTypeMapping <: TypeMapping
    spec::Union{SpecStruct,SpecUnion}
    ex::Expr
end

struct IntermediateTypeMapping <: TypeMapping
    spec::Union{SpecStruct,SpecUnion}
    ex::Expr
end

struct AliasDeclaration <: WrapperNode
    source::Symbol
    target::Symbol
end

function AliasDeclaration((source, target))
    f = startswith(string(source), "vk") ? Base.Fix2(function_name, true) : remove_vk_prefix
    AliasDeclaration(f(source), f(target))
end

to_expr(def::Union{StructureType, TypeMapping}) = def.ex
to_expr((; source, target)::AliasDeclaration) = :(const $source = $target)
to_expr(def::WrapperNode) = resolve_types(to_expr(def.p))
to_expr(def::Union{Documented, ConstantDefinition, EnumDefinition, BitmaskDefinition}) = to_expr(def.p)
to_expr(def::StructDefinition{true,<:SpecStruct}) = :(@auto_hash_equals $(resolve_types(to_expr(def.p))))

function documented(def::WrapperNode)
    doc = Documented(def)
    # Avoid empty docstrings.
    isempty(strip(doc.p[:docstring], '\n')) && return to_expr(def)
    to_expr(doc)
end

name(def::WrapperNode) = name(def.p)
name(alias::AliasDeclaration) = alias.source

function exports(def::WrapperNode)
    @match n = name(def) begin
        ::Symbol => n
        :($_.$n) => n
    end
end

exports(def::Union{EnumDefinition, BitmaskDefinition}) = [name(def); name.(def.p[:values])]

struct VulkanWrapper
    independent::Vector{Expr}
    interdependent::Vector{Expr}
    dependent::Vector{Expr}
    exports::Vector{Symbol}
    function VulkanWrapper(independent::Vector{Expr}, interdependent::Vector{Expr}, dependent::Vector{Expr}, exports::Vector{Symbol})
        new(independent, interdependent, dependent, unique(exports))
    end
end

include("wrap/classification.jl")
include("wrap/identifiers.jl")
include("wrap/defaults.jl")
include("wrap/call.jl")
include("wrap/return.jl")
include("wrap/decl.jl")
include("wrap/initialization.jl")
include("wrap/parent.jl")

include("wrap/constants.jl")
include("wrap/enums.jl")
include("wrap/bitmasks.jl")
include("wrap/handles.jl")
include("wrap/structs.jl")
include("wrap/unions.jl")
include("wrap/functions.jl")
include("wrap/reflection.jl")
include("wrap/docs.jl")

function VulkanWrapper(config::WrapperConfig)
    f = filter_specs(config)

    constants = ConstantDefinition.(filter(include_constant, f(spec_constants)))
    enums = EnumDefinition.(f(spec_enums))
    bitmasks = BitmaskDefinition.(f(spec_bitmasks))
    handles = HandleDefinition.(f(spec_handles))

    # Structures.
    structs = StructDefinition{false}.(f(spec_structs))
    structs_hl = StructDefinition{true}.(structs)

    struct_constructors = Constructor.(structs)
    struct_constructors_from_hl = [Constructor(T, x) for (T, x) in zip(structs, structs_hl)]
    struct_constructors_from_ll = [Constructor(T, x) for (T, x) in zip(structs_hl, structs)]
    struct_constructors_from_core = [Constructor(T, x) for (T, x) in zip(structs_hl, f(spec_structs))]
    struct_constructors_hl = Constructor.(structs_hl)

    ## Do not overwrite the default constructor (leads to infinite recursion).
    filter!(struct_constructors_hl) do def
        !isempty(def.p[:kwargs])
    end

    # Unions.
    unions = StructDefinition{false}.(f(spec_unions))
    unions_hl = StructDefinition{true}.(unions)

    union_constructors = [constructors.(unions)...;]
    union_constructors_from_hl = [Constructor(T, x) for (T, x) in zip(unions, unions_hl)]
    union_constructors_hl = [constructors.(unions_hl)...;]

    union_getproperty_hl = GetProperty.(unions_hl)

    enum_converts_to_integer = [Convert(enum, enum_val_type(enum)) for enum in enums]
    enum_converts_to_enum = [Convert(enum_val_type(enum), enum) for enum in enums]
    enum_converts_from_spec = [Convert(enum, spec_enum.name) for (enum, spec_enum) in zip(enums, f(spec_enums))]
    enum_converts_to_spec = [Convert(spec_enum.name, enum) for (enum, spec_enum) in zip(enums, f(spec_enums))]
    struct_converts_to_ll = [Convert(T, x) for (T, x) in zip(structs, structs_hl)]
    union_converts_to_ll = [Convert(T, x) for (T, x) in zip(unions, unions_hl)]

    funcs = APIFunction.(f(spec_funcs), false)
    funcs_fptr = APIFunction.(f(spec_funcs), true)
    funcs_hl = promote_hl.(funcs)
    funcs_hl_fptr = promote_hl.(funcs_fptr)

    create_func_wrappers = APIFunction{CreateFunc}[]
    create_func_wrappers_fptr = APIFunction{CreateFunc}[]
    create_func_wrappers_hl = APIFunction{APIFunction{CreateFunc}}[]
    create_func_wrappers_hl_fptr = APIFunction{APIFunction{CreateFunc}}[]
    handle_constructors = Constructor{HandleDefinition,APIFunction{CreateFunc}}[]
    handle_constructors_fptr = Constructor{HandleDefinition,APIFunction{CreateFunc}}[]
    handle_constructors_hl = Constructor{HandleDefinition,APIFunction{APIFunction{CreateFunc}}}[]
    handle_constructors_hl_fptr = Constructor{HandleDefinition,APIFunction{APIFunction{CreateFunc}}}[]
    handle_constructors_api = Constructor{HandleDefinition,APIFunction{SpecFunc}}[]
    handle_constructors_api_fptr = Constructor{HandleDefinition,APIFunction{SpecFunc}}[]
    handle_constructors_api_hl = Constructor{HandleDefinition,APIFunction{APIFunction{SpecFunc}}}[]
    handle_constructors_api_hl_fptr = Constructor{HandleDefinition,APIFunction{APIFunction{SpecFunc}}}[]

    for handle in handles
        cs = f(filter(x -> x.handle == handle.spec && !x.batch, spec_create_funcs))
        for api_constructor in cs
            (; func) = api_constructor
            f1 = APIFunction(func, false)
            f2 = APIFunction(func, true)
            f1_p = promote_hl(f1)
            f2_p = promote_hl(f2)
            if !isnothing(api_constructor.create_info_param)
                cf1 = APIFunction(api_constructor, false)
                cf2 = APIFunction(api_constructor, true)
                push!(create_func_wrappers, cf1)
                push!(create_func_wrappers_fptr, cf2)
                cf1_p = promote_hl(cf1)
                cf2_p = promote_hl(cf2)
                push!(create_func_wrappers_hl, cf1_p)
                push!(create_func_wrappers_hl_fptr, cf2_p)
                if can_wrap(handle.spec, cs, api_constructor)
                    if contains_api_structs(cf1)
                        push!(handle_constructors, Constructor(handle, cf1))
                        push!(handle_constructors_fptr, Constructor(handle, cf2))
                    end
                    push!(handle_constructors_hl, Constructor(handle, cf1_p))
                    push!(handle_constructors_hl_fptr, Constructor(handle, cf2_p))
                end
            end
            can_wrap(handle.spec, cs, api_constructor) || continue
            push!(handle_constructors_api_hl, Constructor(handle, f1_p))
            push!(handle_constructors_api_hl_fptr, Constructor(handle, f2_p))
            if contains_api_structs(f1)
                push!(handle_constructors_api, Constructor(handle, f1))
                push!(handle_constructors_api_fptr, Constructor(handle, f2))
            end
        end
    end

    parent_overloads = Parent.(filter(has_parent, handles))

    stypes = StructureType.(filter(x -> haskey(structure_types, x.name), f(spec_structs)))
    hl_type_mappings = [HLTypeMapping.(f(spec_structs)); HLTypeMapping.(f(spec_unions))]
    core_type_mappings = [CoreTypeMapping.(f(spec_structs)); CoreTypeMapping.(f(spec_unions))]
    intermediate_type_mappings = IntermediateTypeMapping.(filter(has_intermediate_type, f(spec_structs)))

    # For SPIR-V, there is no platform-dependent behavior, so no need to call `f`.
    spirv_exts = spec_spirv_extensions
    spirv_caps = map(spec_spirv_capabilities) do spec
        feats = map(spec.enabling_features) do feat
            FeatureCondition(struct_name(follow_alias(feat.type), true), nc_convert(SnakeCaseLower, feat.member), feat.core_version, feat.extension)
        end
        props = map(spec.enabling_properties) do prop
            bit = isnothing(prop.bit) ? nothing : remove_vk_prefix(prop.bit)
            PropertyCondition(struct_name(prop.type, true), nc_convert(SnakeCaseLower, prop.member), prop.core_version, prop.extension, prop.is_bool, bit)
        end
        SpecCapabilitySPIRV(spec.name, spec.promoted_to, spec.enabling_extensions, feats, props)
    end

    exported_symbols = Symbol[
        exports.(constants); exports.(enums)...; exports.(bitmasks)...; exports.(handles); exports.(structs); exports.(unions); exports.(structs_hl); exports.(unions_hl); exports.(funcs); exports.(funcs_hl); :SPIRV_EXTENSIONS; :SPIRV_CAPABILITIES; :CORE_FUNCTIONS; :INSTANCE_FUNCTIONS; :DEVICE_FUNCTIONS;
    ]

    aliases = AliasDeclaration[]
    for (source, target) in collect(alias_dict)
        startswith(string(target), "vk") && continue
        al = AliasDeclaration(source => target)
        al.target in exported_symbols && push!(aliases, al)
    end
    function_aliases = Expr[]
    _spec_funcs = f(spec_funcs)
    for (source, target) in collect(alias_dict)
        al = AliasDeclaration(source => target)
        startswith(string(source), "vk") || continue
        f = func_by_name(target)
        handle = dispatch_handle(f)
        param = @match handle begin
            :(device($x)) || :(instance($x)) || x::Symbol => x
            nothing => nothing
        end
        args = Any[:(args...)]
        !isnothing(param) && pushfirst!(args, param)
        if f in _spec_funcs
            push!(function_aliases,
                # :($source(args..., fptr::FunctionPtr; kwargs...) = $target(args..., fptr; kwargs...)),
                :($(al.source)($(args...); kwargs...) = @dispatch $source $handle $(al.target)($(args...); kwargs...))
            )
            push!(exported_symbols, al.source)
        end
    end

    append!(exported_symbols, exports.(aliases))

    VulkanWrapper(
        Expr[
            documented.(constants); documented.(enums); documented.(enum_converts_to_enum); documented.(enum_converts_to_integer); documented.(enum_converts_from_spec); documented.(enum_converts_to_spec); documented.(bitmasks);
        ],
        Expr[
            documented.(handles); documented.(structs); documented.(unions); documented.(structs_hl); documented.(unions_hl);
        ],
        Expr[
            documented.(parent_overloads);
            documented.(union_constructors);
            documented.(union_constructors_hl);
            documented.(union_constructors_from_hl);
            documented.(union_converts_to_ll);
            documented.(union_getproperty_hl);
            documented.(struct_constructors);
            documented.(struct_constructors_hl);
            to_expr.(struct_constructors_from_hl);
            to_expr.(struct_constructors_from_ll);
            to_expr.(struct_constructors_from_core);
            documented.(struct_converts_to_ll);

            documented.(funcs);
            to_expr.(funcs_fptr);
            documented.(funcs_hl);
            to_expr.(funcs_hl_fptr);

            documented.(create_func_wrappers);
            to_expr.(create_func_wrappers_fptr);
            documented.(create_func_wrappers_hl);
            to_expr.(create_func_wrappers_hl_fptr);

            documented.(handle_constructors);
            to_expr.(handle_constructors_fptr);
            documented.(handle_constructors_hl);
            to_expr.(handle_constructors_hl_fptr);
            documented.(handle_constructors_api);
            to_expr.(handle_constructors_api_fptr);
            documented.(handle_constructors_api_hl);
            to_expr.(handle_constructors_api_hl_fptr);

            to_expr.(stypes);
            to_expr.(hl_type_mappings);
            to_expr.(core_type_mappings);
            to_expr.(intermediate_type_mappings);

            to_expr.(aliases);
            function_aliases;

            :(const SPIRV_EXTENSIONS = [$(spirv_exts...)]);
            :(const SPIRV_CAPABILITIES = [$(spirv_caps...)]);
            :(const CORE_FUNCTIONS = $core_functions);
            :(const INSTANCE_FUNCTIONS = $instance_functions);
            :(const DEVICE_FUNCTIONS = $device_functions);
        ],
        exported_symbols,
    )
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    UniversalKriging(γ, degree, dim)
    UniversalKriging(X, z, γ, degree)

Universal Kriging with variogram model `γ` and polynomial
`degree` on a spatial domain of dimension `dim`.

Optionally, pass the coordinates `X` and values `z`
to the [`fit`](@ref) function.

### Notes

* [`OrdinaryKriging`](@ref) is recovered for 0th degree polynomial
* For non-polynomial mean, see [`ExternalDriftKriging`](@ref)
"""
struct UniversalKriging{G<:Variogram} <: KrigingEstimator
  γ::G
  degree::Int
  dim::Int
  exponents::Matrix{Int}

  function UniversalKriging{G}(γ, degree, dim) where {G<:Variogram}
    @assert degree ≥ 0 "degree must be nonnegative"
    @assert dim > 0 "dimension must be positive"
    exponents = UKexps(degree, dim)
    new(γ, degree, dim, exponents)
  end
end

UniversalKriging(γ, degree, dim) = UniversalKriging{typeof(γ)}(γ, degree, dim)

UniversalKriging(X, z, γ, degree) = GeoStatsBase.fit(UniversalKriging(γ, degree, size(X,1)), X, z)

function UKexps(degree::Int, dim::Int)
  # multinomial expansion
  expmats = [hcat(collect(multiexponents(dim, d))...) for d in 0:degree]
  exponents = hcat(expmats...)

  # sort expansion for better conditioned Kriging matrices
  sorted = sortperm(vec(maximum(exponents, dims=1)), rev=true)

  exponents[:,sorted]
end

nconstraints(estimator::UniversalKriging) = size(estimator.exponents, 2)

function set_constraints_lhs!(estimator::UniversalKriging, LHS::AbstractMatrix, X::AbstractMatrix)
  exponents = estimator.exponents
  nobs = size(X, 2)
  nterms = size(exponents, 2)
  T = eltype(LHS)

  # set polynomial drift blocks
  for i=1:nobs, j=1:nterms
    LHS[nobs+j,i] = prod(X[:,i].^exponents[:,j])
    LHS[i,nobs+j] = LHS[nobs+j,i]
  end

  # set zero block
  LHS[nobs+1:end,nobs+1:end] .= zero(T)

  nothing
end

factorize(estimator::UniversalKriging, LHS::AbstractMatrix) = lu(LHS, check=false)

function set_constraints_rhs!(estimator::FittedKriging{E,S},
                              xₒ::AbstractVector) where {E<:UniversalKriging,S<:KrigingState}
  exponents = estimator.estimator.exponents
  RHS = estimator.state.RHS
  nobs = size(estimator.state.X, 2)
  nterms = size(exponents, 2)

  for j in 1:nterms
    RHS[nobs+j] = prod(xₒ.^exponents[:,j])
  end

  nothing
end


text:
module Unregistered
import Example
export plusfive

plusfive(x) = Example.domath(x)

end # module


text:
using DataFrames, CSV, Dates, ProgressBars

include("initialise.jl")
include("forward_propagation.jl")
include("cost.jl")
include("check_accuracy.jl")
include("backpropagation.jl")
include("gradient_descent.jl")

"""
Train the network using the desired architecture by matching the training inputs (DMatrix) to their
corresponding outputs (Y) over some number of iterations (epochs) and a learning rate (η).
"""
function train_network(layer_dims::Vector{Int64}, DMatrix::Vector{Vector{Float64}}, Y::Vector{Vector{Float64}}, mapping::String; η = 0.001, epochs = 500, seed = 758907)

    η_first = η
    costs = BigFloat[BigFloat(0.0) for x in 1:length(DMatrix) * epochs]
    iters = Int[trunc((x-1) / length(DMatrix))+1 for x in 1:length(DMatrix) * epochs]
    confidence = BigFloat[BigFloat(0.0) for x in 1:length(DMatrix) * epochs]
    YvŶ = Vector{Char}[Char[' ' for x in 1:length(DMatrix) * epochs], Char[' ' for x in 1:length(DMatrix) * epochs]]
    mapping *= "#"

    # initialise
    params = initialise_model(layer_dims, seed)

    # LOGGING
    d = string(Dates.now())[1:13] * "-" * string(Dates.now())[15:16]

    # train
    for i = 1:epochs

        Ŷ = Any[[] for x in 1:length(DMatrix)]
        caches = [[] for x in 1:length(DMatrix)]

        iter = ProgressBar(1:length(DMatrix))
        set_description(iter, "Epoch $i 1/6 Forward Propagating")

        Threads.@threads for datum in iter
            Ŷ[datum], caches[datum] = forward_propagate_model_weights(DMatrix[datum], params)
        end # for

        set_description(iter, "Epoch $i 2/6 Calculating Costs") # I sorta modified this step...

        Threads.@threads for datum in iter
            costs[(i-1) * length(DMatrix) + datum] = calculate_cost(Ŷ[datum], Y[datum])
        end # for

        set_description(iter, "Epoch $i 3/6 Calculating Confidence")

        Threads.@threads for datum in iter
            confidence[(i-1) * length(DMatrix) + datum] = assess_accuracy(Ŷ[datum], Y[datum])
        end # for

        set_description(iter, "Epoch $i 4/6 Calculating Gradients")

        ∇ = Any[missing for x in 1:length(DMatrix)]
        Threads.@threads for datum in iter
            ∇[datum] = back_propagate_model_weights(Ŷ[datum], Y[datum], caches[datum])
        end # for

        iter = ProgressBar(1:length(collect(keys(params))))
        set_description(iter, "Epoch $i 5/6 Updating Model Params")
        ∇_avg = Dict()

        for x in iter
            ∇_avg[string("δ", collect(keys(params))[x])] = sum([δ[string("δ", collect(keys(params))[x])] for δ in ∇]) / length(∇)
        end # for

        params = update_model_weights(params, ∇_avg, η)

        iter = ProgressBar(1:length(DMatrix))
        set_description(iter, "Epoch $i 6/6 Logging calcs")

        for datum in iter
            YvŶ[1][(i-1) * length(DMatrix) + datum] = mapping[findmax(Y[datum])[2]]
            YvŶ[2][(i-1) * length(DMatrix) + datum] = mapping[findmax(Ŷ[datum])[2][1]]
        end # for

        open("case_names/logs/$d LOG.csv", "w") do io
            CSV.write(io, DataFrame([iters, costs, confidence, [(YvŶ[1][x] == YvŶ[2][x]) for x in 1:length(YvŶ[1])], YvŶ[1], YvŶ[2]], :auto), header=["Iteration", "Cost", "Accuracy", "Correct", "Answer", "Guess"])
        end # do
        open("case_names/logs/$d PAR.csv", "w") do io
            CSV.write(io, DataFrame([collect(values(params))], :auto), header=collect(keys(params)))
        end # do

    end # for

    return (cost = costs, iterations = iters, confidence = confidence, parameters = params)

end # function

max_stars_repo_name: Squalm/NeuralNetworks

max_issues_repo_name: Squalm/NeuralNetworks

max_forks_repo_name: Squalm/NeuralNetworks


text:
#Modify the previous program such that only multiples of three or five are considered in the sum, e.g. 3, 5, 6, 9, 10, 12, 15 for n=17

println("Please provide a integer to sum")
num = parse(Int64,readline())

function n_sum_3(x::Int64)
    se = 1:x
    sum( [i % 3 == 0 || i % 5 == 0  ? i : 0 for i in se] )
end

println(n_sum_3(num))

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
# # Conway's game of life

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../game of life.mp4" type="video/mp4">
# </video>
# ```

# [Game of life on wikipedia](https://en.wikipedia.org/wiki/Conway%27s_Game_of_Life).

# It is also available from the `Models` module as [`Models.game_of_life`](@ref).

using Agents, Random

# ## 1. Define the rules

# Conway's game of life is a *cellular automaton*, where each cell of the discrete space
# contains one agent only.

# The rules of Conway's game of life are defined based on four numbers:
# Death, Survival, Reproduction, Overpopulation, grouped as (D, S, R, O)
# Cells die if the number of their living neighbors is <D or >O,
# survive if the number of their living neighbors is ≤S,
# come to life if their living neighbors are  ≥R and ≤O.
rules = (2, 3, 3, 3) # (D, S, R, O)
nothing # hide

# ## 2. Build the model

# First, define an agent type. It needs to have the compulsary `id` and `pos` fields,
# as well as a `status` field that is `true` for cells that are alive and `false` otherwise.

mutable struct Cell <: AbstractAgent
    id::Int
    pos::Dims{2}
    status::Bool
end

# The following function builds a 2D cellular automaton given some `rules`.
# `dims` is a tuple of integers determining the width and height of the grid environment.
# `metric` specifies how to measure distances in the space, and in our example
# it actually decides whether cells connect to their diagonal neighbors.

# This function creates a model where all cells are dead.
function build_model(; rules::Tuple, dims = (100, 100), metric = :chebyshev, seed = 120)
    space = GridSpace(dims; metric)
    properties = Dict(:rules => rules)
    model = ABM(Cell, space; properties, rng = MersenneTwister(seed))
    idx = 1
    for x in 1:dims[1]
        for y in 1:dims[2]
            add_agent_pos!(Cell(idx, (x, y), false), model)
            idx += 1
        end
    end
    return model
end
nothing # hide

# Now we define a stepping function for the model to apply the rules to agents.
# We will also perform a synchronous agent update (meaning that the value of all
# agents changes after we have decided the new value for each agent individually).
function ca_step!(model)
    new_status = fill(false, nagents(model))
    for agent in allagents(model)
        n = alive_neighbors(agent, model)
        if agent.status == true && (n ≤ model.rules[4] && n ≥ model.rules[1])
            new_status[agent.id] = true
        elseif agent.status == false && (n ≥ model.rules[3] && n ≤ model.rules[4])
            new_status[agent.id] = true
        end
    end

    for id in allids(model)
        model[id].status = new_status[id]
    end
end

function alive_neighbors(agent, model) # count alive neighboring cells
    c = 0
    for n in nearby_agents(agent, model)
        if n.status == true
            c += 1
        end
    end
    return c
end
nothing # hide

# now we can instantiate the model:
model = build_model(rules = rules, dims = (50, 50))

# Let's make some random cells on
for i in 1:nagents(model)
    if rand(model.rng) < 0.2
        model.agents[i].status = true
    end
end

# ## 3. Animate the model

# We use the [`InteractiveDynamics.abm_video`](@ref) for creating an animation and saving it to an mp4

using InteractiveDynamics
import CairoMakie

ac(x) = x.status == true ? :black : :white
am(x) = x.status == true ? '■' : '□'
abm_video(
    "game of life.mp4",
    model,
    dummystep,
    ca_step!;
    title = "Game of Life",
    ac = :black,
    as = 12,
    am,
    framerate = 5,
    scatterkwargs = (strokewidth = 0,),
)
nothing # hide
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../game of life.mp4" type="video/mp4">
# </video>
# ```


text:
MathTypes{T,C} = Union{AbstractRGB{T},TransparentRGB{C,T},AbstractGray{T},TransparentGray{C,T}}

# provided by https://github.com/JuliaLang/julia/pull/35094
function register_hints()
    if isdefined(Base, :Experimental) && isdefined(Base.Experimental, :register_error_hint)
        Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs
            if exc.f in (one,) && argtypes[1] <: Union{Type{<:AbstractRGB}, AbstractRGB}
                print(io, "\nYou may need to `using ColorVectorSpace`.")
            end
            if exc.f in (ones,) && argtypes[1] <: Type{<:AbstractRGB}
                print(io, "\nYou may need to `using ColorVectorSpace`.")
            end
            # In theory we could list every function supported by ColorVectorSpace.
            # This list of functions is far from comprehensive but may be enough to catch many users.
            if exc.f in (+, -, *, /) && any(T->T<:MathTypes || T<:Type{<:MathTypes}, argtypes)
                print(io, "\nMath on colors is deliberately undefined in ColorTypes, but see the ColorVectorSpace package.")
            end
            if (exc.f === *) && all(T->T<:MathTypes || T<:Type{<:MathTypes}, argtypes)
                print(io, "\nYou may also need `⋅`, `⊙`, or `⊗`.")
            end
        end
    end
end


text:
module Strategy

import SearchBall, GeometryTypes
const sb = SearchBall
const gt = GeometryTypes

function init_strategies(config; kwargs...)
    # sb.init_text!(config, "distance", 0.0)
    return
end

function calc_myPlayers(config::sb.Config, dt::Number, set_type::Symbol, moved_object::Union{sb.Object{<:Number}, Nothing}=nothing)
    ws = config.ws

    # do stuff

    return
end

function calc_opponend(config::sb.Config, dt::Number, set_type::Symbol, moved_object::Union{sb.Object{<:Number}, Nothing}=nothing)
    ws = config.ws

    # do stuff

    return
end

function run_strategies(config::sb.Config, dt::Number=0, set_type::Symbol=:default, moved_object::Union{sb.Object{<:Number}, Nothing}=nothing)
    calc_myPlayers(config, dt, set_type, moved_object),
    calc_opponend(config, dt, set_type, moved_object)
    return
end

end


text:
using Oceananigans.Architectures
using Oceananigans.Architectures: device
import Oceananigans.Architectures: architecture
using CUDA, CUDA.CUSPARSE
using KernelAbstractions: @kernel, @index

using LinearAlgebra, SparseArrays, IncompleteLU
using SparseArrays: fkeep!

# Utils for sparse matrix manipulation

@inline constructors(::CPU, A::SparseMatrixCSC) = (A.n, A.n, A.colptr, A.rowval, A.nzval)
@inline constructors(::GPU, A::SparseMatrixCSC) = (CuArray(A.colptr), CuArray(A.rowval), CuArray(A.nzval),  (A.n, A.n))
@inline constructors(::CPU, A::CuSparseMatrixCSC) = (A.dims[1], A.dims[2], Int64.(Array(A.colPtr)), Int64.(Array(A.rowVal)), Array(A.nzVal))
@inline constructors(::GPU, A::CuSparseMatrixCSC) = (A.colPtr, A.rowVal, A.nzVal,  A.dims)
@inline constructors(::CPU, n::Number, constr::Tuple) = (n, n, constr...)
@inline constructors(::GPU, n::Number, constr::Tuple) = (constr..., (n, n))

@inline unpack_constructors(::CPU, constr::Tuple) = (constr[3], constr[4], constr[5])
@inline unpack_constructors(::GPU, constr::Tuple) = (constr[1], constr[2], constr[3])
@inline copy_unpack_constructors(::CPU, constr::Tuple) = deepcopy((constr[3], constr[4], constr[5]))
@inline copy_unpack_constructors(::GPU, constr::Tuple) = deepcopy((constr[1], constr[2], constr[3]))

@inline arch_sparse_matrix(::CPU, constr::Tuple) = SparseMatrixCSC(constr...)
@inline arch_sparse_matrix(::GPU, constr::Tuple) = CuSparseMatrixCSC(constr...)
@inline arch_sparse_matrix(::CPU, A::CuSparseMatrixCSC) = SparseMatrixCSC(constructors(CPU(), A)...)
@inline arch_sparse_matrix(::GPU, A::SparseMatrixCSC)   = CuSparseMatrixCSC(constructors(GPU(), A)...)
@inline arch_sparse_matrix(::CPU, A::SparseMatrixCSC)   = A
@inline arch_sparse_matrix(::GPU, A::CuSparseMatrixCSC) = A

# We need to update the diagonal element each time the time step changes!!
function update_diag!(constr, arch, problem_size, diag, Δt)
    M = prod(problem_size)
    colptr, rowval, nzval = unpack_constructors(arch, constr)

    loop! = _update_diag!(device(arch), 256, M)
    event = loop!(nzval, colptr, rowval, diag, Δt; dependencies=Event(device(arch)))
    wait(event)

    constr = constructors(arch, M, (colptr, rowval, nzval))
end

@kernel function _update_diag!(nzval, colptr, rowval, diag, Δt)
    col = @index(Global, Linear)
    map = 1
    for idx in colptr[col]:colptr[col+1] - 1
       if rowval[idx] == col
           map = idx
            break
        end
    end
    nzval[map] += diag[col] / Δt^2
end

@kernel function _get_inv_diag!(invdiag, colptr, rowval, nzval)
    col = @index(Global, Linear)
    map = 1
    for idx in colptr[col]:colptr[col+1] - 1
        if rowval[idx] == col
            map = idx
            break
        end
    end
    if nzval[map] == 0.0
        invdiag[col] = 0
    else
        invdiag[col] = 1.0 / nzval[map]
    end
end

@kernel function _get_diag!(diag, colptr, rowval, nzval)
    col = @index(Global, Linear)
    map = 1
    for idx in colptr[col]:colptr[col+1] - 1
        if rowval[idx] == col
            map = idx
            break
        end
    end
    diag[col] = nzval[map]
end

#unfortunately this cannot run on a GPU so we have to resort to that ugly loop in _update_diag!
@inline map_row_to_diag_element(i, rowval, colptr) =  colptr[i] - 1 + findfirst(rowval[colptr[i]:colptr[i+1]-1] .== i)

@inline function validate_laplacian_direction(N, topo, reduced_dim)
    dim = N > 1 && reduced_dim == false
    if N < 3 && topo == Bounded && dim == true
        throw(ArgumentError("Cannot calculate laplacian in bounded domain with N < 3!"))
    end

    return dim
end

@inline validate_laplacian_size(N, dim) = dim == true ? N : 1

@inline ensure_diagonal_elements_are_present!(A) = fkeep!(A, (i, j, x) -> (i == j || !iszero(x)))


text:
using Javis
using Documenter

makedocs(;
    modules=[Javis],
    authors="Ole Kröger <o.kroeger@opensourc.es> and contributors",
    repo="https://github.com/Wikunia/Javis.jl/blob/{commit}{path}#L{line}",
    sitename="Javis.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Wikunia.github.io/Javis.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
	"Mission" => "mission.md",
        "References" => "references.md",
        "Contributing" => "contributing.md",
	"Tutorials" => [
			"tutorials/tutorial_1.md",
			"tutorials/tutorial_2.md",
			"tutorials/tutorial_3.md",
			],
    ],
)

deploydocs(;
    repo="github.com/Wikunia/Javis.jl",
)


text:
"""
    PosteriorSample

Samples from the poseteriors distribution.
"""
abstract type PosteriorSample end
abstract type SinglePosteriorSample <: PosteriorSample end
abstract type MultiplePosteriorSample <: PosteriorSample end

struct BernoulliPosteriorSample <: SinglePosteriorSample
    θ::Vector{Real}
end

struct ExponentialPosteriorSample <: SinglePosteriorSample
    θ::Vector{Real} # scale, inverse of rate λ in wikipedia
end

struct PoissonPosteriorSample <: SinglePosteriorSample
    λ::Vector{Real}
end

struct NormalPosteriorSample <: MultiplePosteriorSample
    μ::Vector{Real}
    σ²::Vector{Real}
end

struct LogNormalPosteriorSample <: MultiplePosteriorSample
    μ_logx::Vector{Real}
    σ²_logx::Vector{Real}
    μ_x::Vector{Real}
    σ²_x::Vector{Real}
end

function mean(postsample::PosteriorSample, parameter::Symbol)
    return Dict(parameter => mean(getfield(postsample, parameter)))
end
mean(postsample::BernoulliPosteriorSample) = mean(postsample, :θ)
mean(postsample::ExponentialPosteriorSample) = mean(postsample, :θ)

function mean(postsample::MultiplePosteriorSample)
    result = Dict{Symbol,Real}()
    for field in fieldnames(typeof(postsample))
        result[field] = mean(getfield(postsample, field))
    end
    return result
end

"""
    ProbabilisticModel

`ProbabilisticModel` is a a model of the parameters that we are interested in.
The model is defined by its prior distribution and likelihood function.
"""
abstract type ProbabilisticModel end


text:
using Metaheuristics

if VERSION < v"0.7.0"
    using Base.Test
    srand(31415926534)
else
    using Test
    import Random: seed!
    seed!(31415926534)
end

for tests in [
            "common-methods.jl",
            "box-constrained.jl",
            "constrained.jl",
            "multi-objective.jl",
]
    include(tests)
end


text:
using Distributed
@everywhere using SharedArrays, LinearAlgebra

addprocs(4)

## Markov Chain
@everywhere function myrange_burn(q::SharedArray, burn)
    idx = indexpids(q)
    if idx == 0 # This worker is not assigned a piece
       return 1:0, 1:0
    end
    nchunks = length(procs(q))
    splits = [round(Int, s) for s in range(burn, stop=size(q,1), length=nchunks+1)]
    splits[idx]+1:splits[idx+1]
end

@everywhere function fill_chain!(m::SharedArray, p, q, eδ, eq, burn)
    ew=eq-q
    range = myrange_burn(m, burn)
    for r in range
        W, δsim, qsim = jump(p=p, q=q, eδ=eδ, eq=eq)
        logtrydens = norm(p.*W)-norm(p.*ew)

        if logtrydens.>log.(rand())
            ew=W
            eδ=δsim
            eq=qsim
        end

        m[r,:]=gM(p,ew)
    end
    m
end

function shared_fill!(m, p, q, eδ, eq, burn)
    @sync begin
        for w in procs(m)
            @async remotecall_wait(fill_chain!, w, m, p, q, eδ, eq, burn)
        end
    end
    m
end;

function gchain_par(p, q; Par::Par=Param, eδ=eδ, eq=eq)
    @unpack burn, nsim = Par
    T, N = size(q)
    # eδ, eq = warm_start()

    ew=eq-q

    mchain=zeros((nsim-burn),N)
    mchain_par=SharedArray{Float64,2}(mchain)
    shared_fill!(mchain_par, p, q, eδ, eq, burn)

    println("Chain is ready.")
    return mchain_par
end

## Testin Parallel
@fetchfrom 2 myrange_burn(mchain_shared, 1000)
@fetchfrom 3 myrange_burn(mchain_shared, 1000)
mchain_par = gchain_par(p[1:12],q[1:12,1:6])

@time gchain(p[1:12],q[1:12,1:6])
@time gchain_par(p[1:12],q[1:12,1:6])

@fetchfrom 2 fill_chain!(mchain_shared,p[1:12],q[1:12,1:6], $eδ, $eq, 1000 )

## Optmization

@everywhere function myrange(q::SharedArray)
    idx = indexpids(q)
    if idx == 0 # This worker is not assigned a piece
       return 1:0, 1:0
    end
    nchunks = length(procs(q))
    splits = [round(Int, s) for s in range(0, stop=size(q,1), length=nchunks+1)]
    splits[idx]+1:splits[idx+1]
end

@everywhere function hM_chunk!(hM, mchain, gamma)
   # @show (irange, jrange, trange)  # display so we can see what's happening

   # hM=zeros(size(mchain,2))
   T, N, nfast = size(mchain)
   geta=mchain[:,:,1]
   range = myrange(hM)
   for t in range, r=1:nfast
       val=0.0
       gtry=mchain[t,:,r]
       for j=1:N
           val=+gtry[j].*gamma[j]-geta[t,j].*gamma[j]
       end
       geta[t,:]=log.(rand()) < val[1] ? gtry : geta[t,:]
       hM[t,:]=+geta[t,:]./nfast
       # return hM
   end
   hM
end

function shared_hM_fill!(hM, m, g)
    @sync begin
        for w in procs(hM)
            @async remotecall_wait(hM_chunk!, w, hM, m, g)
        end
    end
    hM
end;

@everywhere function shared_hM_chunck!(hM, m, g)
    hM = hM_chunk!(hM, m, g, myrange(hM))
    return hM
end

function shared_dist_hM!(m, gamma) # display so we can see what's happening
    # hM = zeros(size(m,2))
    hM = @distributed (+) for p in workers()
        remotecall_fetch(shared_hM_chunck!, p, m, gamma)
    end
    return hM
end;

function objfun_par(gamma0::Vector...; Par::Par=Param, mchain=mchain_op)
    @unpack nsim, nfast = Par
    @everywhere gamma = $gamma0
    N = size(mchain[1,:])

    mchain_shared = SharedArray{Float64,2}(mchain)
    hM = shared_dist_hM!(mchain_shared, gamma)

    hM=hM./nfast
    mhM=mean(hM)
    var=Statistics.var(hM)

    obj=(mhM^2)/var

    return obj
end

##
@time objfun_par(collect(1:6))
@time objfun(collect(1:6))
@time objfun_bbo(collect(1:6))

@fetchfrom 2 myrange(mchain_shared)

@fetchfrom 4 myrange(mchain_shared)

myrange(mchain_shared)

## parallel obj funs

@everywhere function myrange(q::SharedArray)
    idx = indexpids(q)
    if idx == 0 # This worker is not assigned a piece
       return 1:0, 1:0
    end
    nchunks = length(procs(q))
    splits = [round(Int, s) for s in range(0, stop=size(q,1), length=nchunks+1)]
    splits[idx]+1:splits[idx+1]
end

@everywhere function hM_chunk!(hM, mchain, gamma)
   # @show (irange, jrange, trange)  # display so we can see what's happening

   # hM=zeros(size(mchain,2))
   T, N, nfast = size(mchain)
   geta=mchain[:,:,1]
   range = myrange(hM)
   for t in range, r=1:nfast
       val=0.0
       gtry=mchain[t,:,r]
       for j=1:N
           val=+gtry[j].*gamma[j]-geta[t,j].*gamma[j]
       end
       geta[t,:]=log.(rand()) < val[1] ? gtry : geta[t,:]
       hM[t,:]=+geta[t,:]./nfast
       # return hM
   end
   hM
end

function shared_hM_fill!(hM, m, g)
    @sync begin
        for w in procs(hM)
            @async remotecall_wait(hM_chunk!, w, hM, m, g)
        end
    end
    hM
end;

function objfun_par(gamma0::Vector ...; mchain=mchain_op)
    # @unpack nsim, nfast = Par
    gamma = gamma0[1]
    T, N, nfast = size(mchain)

    hM = SharedArray{Float64,2}(zeros(T,N))
    shared_hM_fill!(hM, mchain, gamma)

    mhM=sum(hM, dims=1)./T
    omega=(hM.-mhM)'*(hM.-mhM)

    obj=mhM*inv(omega)*mhM'

    return obj[1]
end

function objfun_bbo_par(gamma0; mchain=mchain_op)

    gamma = gamma0
    T,N,nfast = size(mchain)

    hM = SharedArray{Float64,2}(zeros(T,N))
    shared_hM_fill!(hM, mchain, gamma)

    mhM=sum(hM, dims=1)./T
    omega=(hM.-mhM)'*(hM.-mhM)

    obj=mhM*inv(omega)*mhM'

    return obj[1]
end

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
module TestMonotonicer

using TSML
using TSML.Utils
using TSML.TSMLTypes
using TSML.TSMLTransformers

using TSML.Monotonicers
using TSML.Monotonicers: ismonotonic
using TSML.DataReaders
using TSML.Statifiers

using DataFrames
using Dates
using Random
using Test

function test_basicmonotonicer()
  Random.seed!(123)
  dat = rand(1:10,5)
  # check monotonic increasing
  @test ismonotonic(dat) == false
  @test ismonotonic(cumsum(dat)) == true

  mdates = DateTime(2017,12,31,1):Dates.Hour(1):DateTime(2017,12,31,10) |> collect
  mvals = rand(length(mdates)) |> cumsum
  df =  DataFrame(Date=mdates ,Value = mvals)

  mono = Monotonicer(Dict())
  fit!(mono,df)
  res=transform!(mono,df)
  @test ismonotonic(res[:Value]) == false

  fname = joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
  csvfilter = DataReader(Dict(:filename=>fname,:dateformat=>"dd/mm/yyyy HH:MM"))
  valgator = DateValgator(Dict(:dateinterval=>Dates.Hour(1)))
  valnner = DateValNNer(Dict(:dateinterval=>Dates.Hour(1)))
  stfier = Statifier(Dict(:processmissing=>true))
  mono = Monotonicer(Dict())

  mpipeline1 = Pipeline(Dict(
      :transformers => [csvfilter,valgator,mono,stfier]
     )
  )
  fit!(mpipeline1)
  respipe1 = transform!(mpipeline1)

  mpipeline2 = Pipeline(Dict(
      :transformers => [csvfilter,valgator,stfier]
     )
  )
  fit!(mpipeline2)
  respipe2 = transform!(mpipeline2)
  @test sum(respipe1[1,:] .== respipe2[1,:]) == ncol(respipe1)

  mpipeline3 = Pipeline(Dict(
      :transformers => [csvfilter,valgator,valnner,mono,stfier]
     )
  )
  fit!(mpipeline3)
  respipe3 = transform!(mpipeline3)
  mpipeline5 = Pipeline(Dict(
      :transformers => [csvfilter,valgator,valnner,stfier]
     )
  )
  fit!(mpipeline5)
  respipe5 = transform!(mpipeline5)
  val1 = respipe3[1,3:end] |> Vector
  val2 = respipe5[1,3:end] |> Vector
  ok1 = (!isnan).(val1)
  ok2 = (!isnan).(val2)
  @test sum( val1[ok1].== val2[ok2]) == length(val1[ok1])

end
@testset "Monotonicer: readcsv |> valgator |> valnner |> mono |> stfier" begin
  test_basicmonotonicer()
end

function test_typesmonotonicer()
  regularfile = joinpath(dirname(pathof(TSML)),"../data/typedetection/regular.csv")
  monofile = joinpath(dirname(pathof(TSML)),"../data/typedetection/monotonic.csv")
  dailymonofile = joinpath(dirname(pathof(TSML)),"../data/typedetection/dailymonotonic.csv")
  regularfilecsv = DataReader(Dict(:filename=>regularfile,:dateformat=>"dd/mm/yyyy HH:MM"))
  monofilecsv = DataReader(Dict(:filename=>monofile,:dateformat=>"dd/mm/yyyy HH:MM"))
  dailymonofilecsv = DataReader(Dict(:filename=>dailymonofile,:dateformat=>"dd/mm/yyyy HH:MM"))

  valgator = DateValgator(Dict(:dateinterval=>Dates.Hour(1)))
  valnner = DateValNNer(Dict(:dateinterval=>Dates.Hour(1)))
  stfier = Statifier(Dict(:processmissing=>true))
  mono = Monotonicer(Dict())

  regpipeline = Pipeline(Dict(
      :transformers => [regularfilecsv,valgator,valnner,mono]
     )
  )
  fit!(regpipeline)
  regulardf=transform!(regpipeline)

  monopipeline = Pipeline(Dict(
      :transformers => [monofilecsv,valgator,valnner,mono]
     )
  )
  fit!(monopipeline)
  monodf=transform!(monopipeline)

  dailymonopipeline = Pipeline(Dict(
      :transformers => [dailymonofilecsv,valgator,valnner,mono]
     )
  )
  fit!(dailymonopipeline)
  dailymonodf=transform!(dailymonopipeline)

  @test round(dailyflips(regulardf),digits=2) == 11.98
  @test round(dailyflips(monodf),digits=2) == 10.09
  @test round(dailyflips(dailymonodf),digits=2) == 9.11
end
@testset "Monotonicer: check type format detection" begin
  test_basicmonotonicer()
end

end


text:
function GetWorkunitForLoop(n::Array{Int64,1},h::Array{Float64,1})
hinv = 1./h;
times = 20;
t_wu = 0.0;
if length(n)==2
	t = ones(n[1]+2,n[2]+2);
	z = zeros(n[1]+2,n[2]+2);
	tic()
	for k=1:times
		for j = 1:n[2]
			for i = 1:n[1]
				i_t = i+1;
				j_t = j+1;
				z[i,j] = ((0.5*hinv[1]*(t[i_t+1,j_t] - t[i_t-1,j_t]))^2 + (0.5*hinv[2]*(t[i_t,j_t+1] - t[i_t,j_t-1]))^2 );
			end
		end
	end
	t_wu = toq();
elseif length(n)==3
	t = ones(n[1]+2,n[2]+2,n[3]+2);
	z = zeros(n[1]+2,n[2]+2,n[3]+2);
	tic()
	for k=1:times
		for k = 1:n[3]
			for j = 1:n[2]
				for i = 1:n[1]
					i_t = i+1;
					j_t = j+1;
					k_t = k+1;
					z[i,j,k] = ((0.5*hinv[1]*(t[i_t+1,j_t,k_t] - t[i_t-1,j_t,k_t]))^2 + (0.5*hinv[2]*(t[i_t,j_t+1,k_t] - t[i_t,j_t-1,k_t]))^2 + (0.5*hinv[3]*(t[i_t,j_t,k_t+1] - t[i_t,j_t,k_t-1]))^2);
				end
			end
		end
	end
	t_wu = toq();
else
	error("unkonwn dimention");
end
WU = t_wu/times;
println("workunit is: ",WU," for n = ",n);
return WU;
end

# function showWorkunit(h::Array{Float64,1},WU::Float64)
# I = [1.6,1.6,0.8];

# n = zeros(Int64,3);
# n[1] = round(Int64,(I[1]/h[1])+1);
# n[2] = round(Int64,(I[2]/h[2])+1);
# n[3] = round(Int64,(I[3]/h[3])+1);

# t = rand(prod(n));
# (D1,D2,D3) = getLongDiffGradOperators3D(h,n);
# times = 30;
# t_wu = 0.0;
# tic()
# for k=1:times
	# y = D1*t;
	# y.*=y;
	# z = D2*t;
	# z.*=z;
	# y+=z;
	# z = D3*t;
	# z.*=z;
	# y+=z;
# end
# t_wu = toq();

# WU = t_wu/times;
# D1 = D2 = D3 = 0.0;
# println("workunit is: ",WU," for n = ",n);
# return;
# end


text:
struct EmulationModelSimulationResults <: OperationModelSimulationResults
    variables::Dict{OptimizationContainerKey, DataFrames.DataFrame}
    duals::Dict{OptimizationContainerKey, DataFrames.DataFrame}
    parameters::Dict{OptimizationContainerKey, DataFrames.DataFrame}
    aux_variables::Dict{OptimizationContainerKey, DataFrames.DataFrame}
    expressions::Dict{OptimizationContainerKey, DataFrames.DataFrame}
    container_key_lookup::Dict{String, OptimizationContainerKey}
end

function SimulationProblemResults(
    ::Type{EmulationModel},
    store::SimulationStore,
    model_name::AbstractString,
    problem_params::ModelStoreParams,
    sim_params::SimulationStoreParams,
    path,
    container_key_lookup;
    kwargs...,
)
    variables = list_emulation_model_keys(store, STORE_CONTAINER_VARIABLES)
    parameters = list_emulation_model_keys(store, STORE_CONTAINER_PARAMETERS)
    duals = list_emulation_model_keys(store, STORE_CONTAINER_DUALS)
    aux_variables = list_emulation_model_keys(store, STORE_CONTAINER_AUX_VARIABLES)
    expressions = list_emulation_model_keys(store, STORE_CONTAINER_EXPRESSIONS)

    return SimulationProblemResults{EmulationModelSimulationResults}(
        store,
        model_name,
        problem_params,
        sim_params,
        path,
        EmulationModelSimulationResults(
            Dict(x => DataFrames.DataFrame() for x in variables),
            Dict(x => DataFrames.DataFrame() for x in duals),
            Dict(x => DataFrames.DataFrame() for x in parameters),
            Dict(x => DataFrames.DataFrame() for x in aux_variables),
            Dict(x => DataFrames.DataFrame() for x in expressions),
            container_key_lookup,
        );
        kwargs...,
    )
end

function Base.empty!(res::SimulationProblemResults{EmulationModelSimulationResults})
    for container in _get_containers(res)
        for df in values(container)
            empty!(df)
        end
    end
end

function Base.isempty(res::SimulationProblemResults{EmulationModelSimulationResults})
    for container in _get_containers(res)
        for df in values(container)
            if !isempty(df)
                return false
            end
        end
    end

    return true
end

function Base.length(res::SimulationProblemResults{EmulationModelSimulationResults})
    count_not_empty = 0
    for container in _get_containers(res)
        for df in values(container)
            if !isempty(df)
                count_not_empty += 1
            end
        end
    end

    return count_not_empty
end

function _get_store_value(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    container_keys::Vector{<:OptimizationContainerKey},
    ::Nothing;
    start_time=nothing,
    len=nothing,
)
    simulation_store_path = joinpath(get_execution_path(res), "data_store")
    return open_store(HdfSimulationStore, simulation_store_path, "r") do store
        _get_store_value(res, container_keys, store, start_time=start_time, len=len)
    end
end

function _get_store_value(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    container_keys::Vector{<:OptimizationContainerKey},
    store::SimulationStore;
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Nothing, Int}=nothing,
)
    base_power = res.base_power
    results = Dict{OptimizationContainerKey, DataFrames.DataFrame}()
    for key in container_keys
        start_time, _len, resolution = _check_offsets(res, key, store, start_time, len)
        start_index = (start_time - first(res.timestamps)) ÷ resolution + 1
        df = read_results(store, key, index=start_index, len=_len)
        if convert_result_to_natural_units(key)
            df .*= base_power
        end
        time_col = range(start_time, length=_len, step=res.resolution)
        DataFrames.insertcols!(df, 1, :DateTime => time_col)
        results[key] = df
    end

    return results
end

function _check_offsets(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    key,
    store,
    start_time,
    len,
)
    dataset_size = get_emulation_model_dataset_size(store, key)
    resolution =
        (last(res.timestamps) - first(res.timestamps) + res.resolution) ÷ dataset_size[1]
    if isnothing(start_time)
        start_time = first(res.timestamps)
    elseif start_time < first(res.timestamps) || start_time > last(res.timestamps)
        throw(
            IS.InvalidValue(
                "start_time = $start_time is not in the results range $(res.timestamps)",
            ),
        )
    elseif (start_time - first(res.timestamps)) % resolution != Dates.Millisecond(0)
        throw(
            IS.InvalidValue(
                "start_time = $start_time is not a multiple of resolution = $resolution",
            ),
        )
    end

    if isnothing(len)
        len = (last(res.timestamps) + resolution - start_time) ÷ resolution
    elseif start_time + resolution * len > last(res.timestamps) + res.resolution
        throw(
            IS.InvalidValue(
                "len = $len resolution = $resolution exceeds the results range $(res.timestamps)",
            ),
        )
    end

    return start_time, len, resolution
end

function read_aux_variables_with_keys(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    aux_variables::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Nothing, Int}=nothing,
)
    return _read_aux_variables(res, aux_variables, nothing, start_time, len)
end

function _read_aux_variables(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    aux_variable_keys,
    store,
    start_time=nothing,
    len=nothing,
)
    isempty(aux_variable_keys) &&
        return Dict{OptimizationContainerKey, DataFrames.DataFrame}()
    if store === nothing && res.store !== nothing
        # In this case we have an InMemorySimulationStore.
        store = res.store
    end

    _validate_keys(keys(get_aux_variables(res)), aux_variable_keys)
    keys_with_values = Set((k for (k, v) in get_aux_variables(res) if !isempty(v)))
    if isempty(setdiff(aux_variable_keys, keys_with_values))
        @debug "reading aux_variables from SimulationsResults"
        vals = filter(p -> (p.first ∈ aux_variable_keys), get_aux_variables(res))
    else
        @debug "reading aux_variables from data store"
        vals =
            _get_store_value(res, aux_variable_keys, store, start_time=start_time, len=len)
    end
    return vals
end

function _read_expressions(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    expression_keys,
    store,
    start_time=nothing,
    len=nothing,
)
    isempty(expression_keys) &&
        return Dict{OptimizationContainerKey, DataFrames.DataFrame}()
    if store === nothing && res.store !== nothing
        # In this case we have an InMemorySimulationStore.
        store = res.store
    end

    _validate_keys(keys(get_expressions(res)), expression_keys)
    keys_with_values = Set((k for (k, v) in get_expressions(res) if !isempty(v)))
    if isempty(setdiff(expression_keys, keys_with_values))
        @debug "reading expressions from SimulationsResults"
        vals = filter(p -> (p.first ∈ expression_keys), get_expressions(res))
    else
        @debug "reading expressions from data store"
        vals = _get_store_value(res, expression_keys, store, start_time=start_time, len=len)
    end
    return vals
end

function read_duals_with_keys(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    duals::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Nothing, Int}=nothing,
)
    return _read_duals(res, duals, nothing, start_time, len)
end

function _read_duals(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    dual_keys,
    store,
    start_time=nothing,
    len=nothing,
)
    isempty(dual_keys) && return Dict{OptimizationContainerKey, DataFrames.DataFrame}()
    if store === nothing && res.store !== nothing
        # In this case we have an InMemorySimulationStore.
        store = res.store
    end

    _validate_keys(keys(get_duals(res)), dual_keys)
    keys_with_values = Set((k for (k, v) in get_duals(res) if !isempty(v)))
    if isempty(setdiff(dual_keys, keys_with_values))
        @debug "reading duals from SimulationsResults"
        vals = filter(p -> (p.first ∈ dual_keys), get_duals(res))
    else
        @debug "reading duals from data store"
        vals = _get_store_value(res, dual_keys, store, start_time=start_time, len=len)
    end
    return vals
end

function read_expressions_with_keys(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    expressions::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Nothing, Int}=nothing,
)
    return _read_expressions(res, expressions, nothing, start_time, len)
end

function read_parameters_with_keys(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    parameters::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Nothing, Int}=nothing,
)
    return _read_parameters(res, parameters, nothing, start_time, len)
end

function _read_parameters(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    parameter_keys,
    store,
    start_time=nothing,
    len=nothing,
)
    isempty(parameter_keys) && return Dict{OptimizationContainerKey, DataFrames.DataFrame}()
    if store === nothing && res.store !== nothing
        # In this case we have an InMemorySimulationStore.
        store = res.store
    end

    _validate_keys(keys(get_parameters(res)), parameter_keys)
    keys_with_values = Set((k for (k, v) in get_parameters(res) if !isempty(v)))
    if isempty(setdiff(parameter_keys, keys_with_values))
        @debug "reading parameters from SimulationsResults"
        vals = filter(p -> (p.first ∈ parameter_keys), get_parameters(res))
    else
        @debug "reading parameters from data store"
        vals = _get_store_value(res, parameter_keys, store, start_time=start_time, len=len)
    end
    return vals
end

function read_variables_with_keys(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    variables::Vector{<:OptimizationContainerKey};
    start_time::Union{Nothing, Dates.DateTime}=nothing,
    len::Union{Nothing, Int}=nothing,
)
    return _read_variables(res, variables, nothing, start_time, len)
end

function _read_variables(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    variable_keys,
    store,
    start_time=nothing,
    len=nothing,
)
    isempty(variable_keys) && return Dict{OptimizationContainerKey, DataFrames.DataFrame}()
    if store === nothing && res.store !== nothing
        # In this case we have an InMemorySimulationStore.
        store = res.store
    end

    _validate_keys(keys(get_variables(res)), variable_keys)
    keys_with_values = Set((k for (k, v) in get_variables(res) if !isempty(v)))
    if isempty(setdiff(variable_keys, keys_with_values))
        @debug "reading variables from SimulationsResults"
        # TODO: This needs to filter by start_time and len.
        vals = filter(p -> (p.first ∈ variable_keys), get_variables(res))
    else
        @debug "reading variables from data store"
        vals = _get_store_value(res, variable_keys, store, start_time=start_time, len=len)
    end
    return vals
end

"""
Load the simulation results into memory for repeated reads. Running this function twice
overwrites the previously loaded results. This is useful when loading results from remote
locations over network connections.

For each variable/parameter/dual, etc., each element must be the name encoded as a string,
like `"ActivePowerVariable__ThermalStandard"`` or a Tuple with its constituent types, like `(ActivePowerVariable, ThermalStandard)`.

# Arguments

  - `aux_variables::Vector{Union{String, Tuple}}`: Optional list of aux variables to load.
  - `duals::Vector{Union{String, Tuple}}`: Optional list of duals to load.
  - `expressions::Vector{Union{String, Tuple}}`: Optional list of expressions to load.
  - `parameters::Vector{Union{String, Tuple}}`: Optional list of parameters to load.
  - `variables::Vector{Union{String, Tuple}}`: Optional list of variables to load.
"""
function load_results!(
    res::SimulationProblemResults{EmulationModelSimulationResults};
    aux_variables=Vector{Tuple}(),
    duals=Vector{Tuple}(),
    expressions=Vector{Tuple}(),
    parameters=Vector{Tuple}(),
    variables=Vector{Tuple}(),
)
    # TODO: consider extending this to support start_time and len
    aux_variable_keys = [_deserialize_key(AuxVarKey, res, x...) for x in aux_variables]
    dual_keys = [_deserialize_key(ConstraintKey, res, x...) for x in duals]
    expression_keys = [_deserialize_key(ExpressionKey, res, x...) for x in expressions]
    parameter_keys = [_deserialize_key(ParameterKey, res, x...) for x in parameters]
    variable_keys = [_deserialize_key(VariableKey, res, x...) for x in variables]
    function merge_results(store)
        merge!(get_aux_variables(res), _read_aux_variables(res, aux_variable_keys, store))
        merge!(get_duals(res), _read_duals(res, dual_keys, store))
        merge!(get_expressions(res), _read_expressions(res, expression_keys, store))
        merge!(get_parameters(res), _read_parameters(res, parameter_keys, store))
        merge!(get_variables(res), _read_variables(res, variable_keys, store))
    end

    if res.store isa InMemorySimulationStore
        merge_results(res.store)
    else
        simulation_store_path = joinpath(res.execution_path, "data_store")
        open_store(HdfSimulationStore, simulation_store_path, "r") do store
            merge_results(store)
        end
    end

    return
end

# TODO: These aren't being written to the store.
function _read_optimizer_stats(
    res::SimulationProblemResults{EmulationModelSimulationResults},
    store::SimulationStore,
)
    return
end


text:
import Base.Unicode: utf8proc_map, normalize, UTF8PROC_STABLE, UTF8PROC_COMPAT, UTF8PROC_DECOMPOSE, UTF8PROC_COMPOSE, UTF8PROC_IGNORE, UTF8PROC_REJECTNA, UTF8PROC_NLF2LS, UTF8PROC_NLF2PS, UTF8PROC_NLF2LF, UTF8PROC_STRIPCC, UTF8PROC_CASEFOLD, UTF8PROC_LUMP, UTF8PROC_STRIPMARK

function _compute_options(
    stable::Bool=false,
    compat::Bool=false,
    compose::Bool=false,
    decompose::Bool=false,
    stripignore::Bool=false,
    rejectna::Bool=false,
    newline2ls::Bool=false,
    newline2ps::Bool=false,
    newline2lf::Bool=false,
    stripcc::Bool=false,
    casefold::Bool=false,
    lump::Bool=false,
    stripmark::Bool=false,
)
    flags = 0
    stable && (flags = flags | UTF8PROC_STABLE)
    compat && (flags = flags | UTF8PROC_COMPAT)
    # TODO: error if compose & decompose?
    if decompose
        flags = flags | UTF8PROC_DECOMPOSE
    elseif compose
        flags = flags | UTF8PROC_COMPOSE
    elseif compat || stripmark
        throw(ArgumentError("compat=true or stripmark=true require compose=true or decompose=true"))
    end
    stripignore && (flags = flags | UTF8PROC_IGNORE)
    rejectna && (flags = flags | UTF8PROC_REJECTNA)
    newline2ls + newline2ps + newline2lf > 1 && throw(ArgumentError("only one newline conversion may be specified"))
    newline2ls && (flags = flags | UTF8PROC_NLF2LS)
    newline2ps && (flags = flags | UTF8PROC_NLF2PS)
    newline2lf && (flags = flags | UTF8PROC_NLF2LF)
    stripcc && (flags = flags | UTF8PROC_STRIPCC)
    casefold && (flags = flags | UTF8PROC_CASEFOLD)
    lump && (flags = flags | UTF8PROC_LUMP)
    stripmark && (flags = flags | UTF8PROC_STRIPMARK)
    flags
end

_compute_options(nf::Symbol) =
    nf == :NFC ? (UTF8PROC_STABLE | UTF8PROC_COMPOSE) :
    nf == :NFD ? (UTF8PROC_STABLE | UTF8PROC_DECOMPOSE) :
    nf == :NFKC ? (UTF8PROC_STABLE | UTF8PROC_COMPOSE | UTF8PROC_COMPAT) :
    nf == :NFKD ? (UTF8PROC_STABLE | UTF8PROC_DECOMPOSE | UTF8PROC_COMPAT) :
    nf == :NFKC_CF ? (UTF8PROC_STABLE | UTF8PROC_COMPOSE | UTF8PROC_COMPAT | UTF8PROC_CASEFOLD) :
    throw(ArgumentError(":$nf is not one of :NFC, :NFD, :NFKC, :NFKD, :NFKC_CF"))

#normalize(s::AbstractString, nf::Symbol) = utf8proc_map(s, _compute_options(nf))
normalize(s::AbstractString, option::Int) = utf8proc_map(s, option)

abstract type AbstractNormalizer end

struct UnNormalizer <: AbstractNormalizer end

normalize(::UnNormalizer, s::AbstractString) = s

struct UtfNormalizer <: AbstractNormalizer
    option::Int
end

UtfNormalizer(nf::Symbol) = UtfNormalizer(_compute_options(nf))
function UtfNormalizer(
    stable::Bool=false,
    compat::Bool=false,
    compose::Bool=false,
    decompose::Bool=false,
    stripignore::Bool=false,
    rejectna::Bool=false,
    newline2ls::Bool=false,
    newline2ps::Bool=false,
    newline2lf::Bool=false,
    stripcc::Bool=false,
    casefold::Bool=false,
    lump::Bool=false,
    stripmark::Bool=false,
)
    return UtfNormalizer(
        _compute_options(
            stable=stable,
            compat=compat,
            compose=compose,
            decompose=decompose,
            stripignore=stripignore,
            rejectna=rejectna,
            newline2ls=newline2ls,
            newline2ps=newline2ps,
            newline2lf=newline2lf,
            stripcc=stripcc,
            casefold=casefold,
            lump=lump,
            stripmark=stripmark,
        )
    )
end

normalize(nr::UtfNormalizer, s::AbstractString) = normalize(s, nr.option)


text:
include("../base.jl")
include("../utils/nodeSetsParser.jl")

function myrandom()
    return rand(["arriba", "abajo", "izquierda", "derecha"])
end

functionset, terminalset = createNodes()

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using FastRefine
using LazySets
import FastRefine: forward_network, forward_affine_map, ishull

nnet = read_nnet("nnet/86442.nnet")

delta = 0.4

solver3 = DimGrid(delta)

in_hyper = Hyperrectangle(fill(1.0, 8), fill(1.0, 8))
out_hyper = Hyperrectangle(fill(0.0, 2), fill(10.0, 2))
problem = Problem(nnet, in_hyper, out_hyper)

(W, b) = (problem.network.layers[1].weights, problem.network.layers[1].bias)
input1 = forward_affine_map(solver, W, b, problem.input)
lower, upper = low(input1), high(input1)
n_hypers_per_dim = BigInt.(max.(ceil.(Int, (upper-lower) / delta), 1))

C = vcat(W, -W)
d = vcat(upper-b, b-lower)

HP = HPolytope(C, d)

inter = intersection(in_hyper, HP)

a = isempty(inter)