text:
struct BackwardReach end

import LazySets:HPoly

# This is the main function
function solve(solver::BackwardReach, problem::Problem)
    safe_zone = backward_network(solver, problem.network, problem.output)
    if isempty(safe_zone)
        return BasicResult(:unknown)
    end
    print("safe_zone $(dim(safe_zone))\n")
    return issubset(problem.input, safe_zone) ? BasicResult(:holds) : BasicResult(:violated)
end

function backward_network(solver::BackwardReach, nnet::Network, input::HPolyhedron{Float64})
    reach = input
    for i in length(nnet.layers):-1:1
        reach = backward_layer(solver, nnet.layers[i], nnet.layers[i].activation, reach)
    end
    return reach
end

function backward_layer(solver::BackwardReach, layer::Layer, act::Id, input::HPolyhedron{Float64})
    lc = get_constraints(input)
    (W, b) = (layer.weights, layer.bias)
    lco = backward_affine_map(lc, W, b)
    remove_redundant_constraints!(lco)
    return HPolyhedron(lco)
end

function backward_layer(solver::BackwardReach, layer::Layer, act::ReLU, input::HPolyhedron{Float64})
    n_dim = dim(input)
    lc = get_constraints(input)
    (W, b) = (layer.weights, layer.bias)
    lcs = Vector{LinearConstraint{Float64, Vector{Float64}}}(undef, 0)
    diag = zeros(Float64, (n_dim, n_dim))
    b = zeros(Float64, n_dim)
    act_status = zeros(Int64, n_dim)
    n_part = 2 ^ n_dim
    for i in 1:n_part
        n = i
        lci = Vector{LinearConstraint{Float64, Vector{Float64}}}(undef, 0)
        for j in 1:n_dim
            n, act_status[j] = fldmod1(n, 2)
            diag[j, j] = 1.0
            if act_status[j] == 1
                push!(lci, LinearConstraint(diag[j,:], 0.0))
            end
            push!(lci, LinearConstraint(-diag[j,:], 0.0))
        end #setfield! immutable struct of type HalfSpace cannot be changed
        append!(lci, lc)
        if isempty(HPolyhedron(lci))
            continue
        end
        remove_redundant_constraints!(lci) #LinearConstraint
        lci_to_delete = Vector{Int64}(undef, 0)
        for k in 1:length(lci)
            if !in(lci[k], lc)
                push!(lci_to_delete, k)
            end
        end
        deleteat!(lci, lci_to_delete)
        Wi = diag
        bi = b
        for l in 1:n_dim
            Wi[l, l] = act_status[l] == 1 ? 0.0 : 1.0
            bi[l] = 0.0
        end
        append!(lcs, backward_affine_map(lci, Wi, bi))
    end
    if length(lcs) == 0
        ni = size(W, 2)
        diagi = zeros(Float64, (ni, ni))
        for m in 1:ni
            diagi[m, m] = 1.0
            push!(lcs, LinearConstraint(diagi[m,:], 0.0))
        end
        return HPolyhedron(lcs)
    end
    lco = backward_affine_map(lcs, W, b)
    remove_redundant_constraints!(lco)
    return HPolyhedron(lco)
end

function backward_layer(solver::BackwardReach, layer::Layer, act::PiecewiseLinearActivation, input::HPolyhedron{Float64})
    n_dim = dim(input)
    lc = get_constraints(input)
    (W, b) = (layer.weights, layer.bias)
    inflections = get_inflection(act)
    lcs = Vector{LinearConstraint{Float64, Vector{Float64}}}(undef, 0)
    diag = zeros(Float64, (n_dim, n_dim))
    b = zeros(Float64, n_dim)
    if inflections == []
        W1 = diag
        b1 = b
        for l in 1:n_dim
            W1[l, l] = act.pieces[1].slope
            b1[l] = act.pieces[1].intercept
        end
        append!(lcs, backward_affine_map(lc, W1, b1))
    else
        out_bounds = [-Inf; act.(inflections); Inf]
        n_interval = length(inflections) + 1
        act_status = zeros(Int64, n_dim)
        n_part = n_interval ^ n_dim
        for i in 1:n_part
            n = i
            lci = Vector{LinearConstraint{Float64, Vector{Float64}}}(undef, 0)
            for j in 1:n_dim
                n, act_status[j] = fldmod1(n, n_interval)
                diag[j, j] = 1.0
                act_status[j] == n_interval || push!(lci, LinearConstraint(diag[j,:], out_bounds[act_status[j] + 1])) #right bound
                act_status[j] == 1 || push!(lci, LinearConstraint(-diag[j,:], -out_bounds[act_status[j]])) #left bound
            end #setfield! immutable struct of type HalfSpace cannot be changed
            append!(lci, lc)
            if isempty(HPolyhedron(lci))
                continue
            end
            remove_redundant_constraints!(lci) #LinearConstraint
            lci_to_delete = Vector{Int64}(undef, 0)
            for k in 1:length(lci)
                if !in(lci[k], lc)
                    push!(lci_to_delete, k)
                end
            end
            deleteat!(lci, lci_to_delete)
            Wi = diag
            bi = b
            for l in 1:n_dim
                Wi[l, l] = act.pieces[act_status[l]].slope
                bi[l] = act.pieces[act_status[l]].intercept
            end
            append!(lcs, backward_affine_map(lci, Wi, bi))
        end
    end
    lco = backward_affine_map(lcs, W, b)
    remove_redundant_constraints!(lco)
    return HPolyhedron(lco)
end

#ensure LinearConstraint{Float64, Vector{Float64}}
function get_constraints(input::HPolyhedron{Float64})
    C, d = tosimplehrep(input) #Cannot `convert` an object of type HalfSpace{Float64,LazySets.Arrays.SingleEntryVector{Float64}} to an object of type HalfSpace{Float64,Array{Float64,1}}
    lc = Vector{LinearConstraint{Float64, Vector{Float64}}}(undef, length(d))
    for i in 1:length(d)
        lc[i] = LinearConstraint(C[i,:], d[i])
    end
    return lc
end

function backward_affine_map(input::Vector{LinearConstraint{Float64, Vector{Float64}}}, M::Matrix{Float64}, v::Vector{Float64})
    output = Vector{LinearConstraint{Float64, Vector{Float64}}}(undef, 0)
    for lc in input
        vp = vprod(lc.a, M)
        if !iszero(vp)
            push!(output, LinearConstraint(vp, lc.b - lc.a' * v))
        end
    end
    return output
end

function vprod(x::Vector{Float64}, y::Matrix{Float64})
    @assert length(x) == size(y, 1) "Vector length $(length(x)) doesn't match Matrix size $(size(y))"
    n = size(y, 2)
    prod = zeros(Float64, n)
    for i in 1:n
        prod[i] = x' * y[:,i]
    end
    return prod
end


text:
__precompile__()

module DiffEqBase

using RecipesBase, SimpleTraits, RecursiveArrayTools, Compat
using Ranges # For plot recipes with units
import Base: length, ndims, size, getindex, setindex!, endof, show, print,
             next, start, done, eltype, eachindex, similar

import Base: resize!, deleteat!

import RecursiveArrayTools: recursivecopy!

# Problems
@compat abstract type DEProblem end
@compat abstract type DEElement end
@compat abstract type DESensitivity end

export DEProblem, DEElement, DESensitivity

@compat abstract type AbstractSteadyStateProblem{uType,isinplace} <: DEProblem end

export AbstractSteadyStateProblem

@compat abstract type AbstractODEProblem{uType,tType,isinplace} <: DEProblem end
@compat abstract type AbstractDiscreteProblem{uType,tType,isinplace} <:
                      AbstractODEProblem{uType,tType,isinplace} end

export AbstractODEProblem, AbstractDiscreteProblem

@compat abstract type AbstractRODEProblem{uType,tType,isinplace,ND} <: DEProblem end

export AbstractRODEProblem

@compat abstract type AbstractSDEProblem{uType,tType,isinplace,ND} <:
                      AbstractRODEProblem{uType,tType,isinplace,ND} end

export AbstractSDEProblem

@compat abstract type AbstractDAEProblem{uType,duType,tType,isinplace} <: DEProblem end

export AbstractDAEProblem

@compat abstract type AbstractDDEProblem{uType,tType,lType,isinplace} <: DEProblem end
@compat abstract type AbstractConstantLagDDEProblem{uType,tType,lType,isinplace} <:
                      AbstractDDEProblem{uType,tType,lType,isinplace} end

export AbstractDDEProblem,
       AbstractConstantLagDDEProblem

@compat abstract type AbstractSecondOrderODEProblem{uType,tType,isinplace} <: AbstractODEProblem{uType,tType,isinplace} end
@compat abstract type AbstractSecondOrderSDEProblem{uType,tType,isinplace} <: AbstractSDEProblem{uType,tType,isinplace} end

export AbstractSecondOrderODEProblem, AbstractSecondOrderSDEProblem

# Algorithms
@compat abstract type DEAlgorithm end
@compat abstract type AbstractSteadyStateAlgorithm <: DEAlgorithm end
@compat abstract type AbstractODEAlgorithm <: DEAlgorithm end
@compat abstract type AbstractSplitODEAlgorithm <: DEAlgorithm end
@compat abstract type AbstractPartitionedODEAlgorithm <: DEAlgorithm end
@compat abstract type AbstractSecondOrderODEAlgorithm <: DEAlgorithm end
@compat abstract type AbstractSplitDAEAlgorithm <: DEAlgorithm end
@compat abstract type AbstractPartitionedDAEAlgorithm <: DEAlgorithm end
@compat abstract type AbstractSplitConstrainedODEAlgorithm <: DEAlgorithm end
@compat abstract type AbstractPartitionedConstrainedODEAlgorithm <: DEAlgorithm end
@compat abstract type AbstractRODEAlgorithm <: DEAlgorithm end
@compat abstract type AbstractSDEAlgorithm <: DEAlgorithm end
@compat abstract type AbstractDAEAlgorithm <: DEAlgorithm end
@compat abstract type AbstractDDEAlgorithm <: DEAlgorithm end

export DEAlgorithm, AbstractSteadyStateAlgorithm, AbstractODEAlgorithm,
       AbstractSplitODEAlgorithm,
       AbstractPartitionedODEAlgorithm, AbstractSecondOrderODEAlgorithm,
       AbstractRODEAlgorithm, AbstractSDEAlgorithm,
       AbstractDAEAlgorithm, AbstractDDEAlgorithm

# Monte Carlo Simulations
@compat abstract type AbstractMonteCarloProblem <: DEProblem end
@compat abstract type AbstractMonteCarloEstimator <: DEProblem end

export AbstractMonteCarloProblem, MonteCarloProblem
export MonteCarloSolution, MonteCarloTestSolution
export AbstractMonteCarloEstimator

@compat abstract type DEOptions end
@compat abstract type DECache end
@compat abstract type DECallback end
@compat abstract type DEDataArray{T} <: AbstractVector{T} end

export DEOptions, DECache, DECallback, DEDataArray

# Integrators
@compat abstract type DEIntegrator end
@compat abstract type AbstractSteadyStateIntegrator <: DEIntegrator end
@compat abstract type AbstractODEIntegrator <: DEIntegrator end
@compat abstract type AbstractSplitODEIntegrator <: DEIntegrator end
@compat abstract type AbstractPartitionedODEIntegrator <: DEIntegrator end
@compat abstract type AbstractSecondOrderODEIntegrator <: DEIntegrator end
@compat abstract type AbstractRODEIntegrator <: DEIntegrator end
@compat abstract type AbstractSDEIntegrator <: DEIntegrator end
@compat abstract type AbstractDDEIntegrator <: DEIntegrator end
@compat abstract type AbstractDAEIntegrator <: DEIntegrator end

export DEIntegrator, AbstractODEIntegrator, AbstractSplitODEIntegrator,
       AbstractPartitionedODEIntegrator, AbstractSecondOrderODEIntegrator,
       AbstractRODEIntegrator, AbstractSDEIntegrator,
       AbstractDDEIntegrator, AbstractDAEIntegrator

# Solutions
@compat abstract type DESolution end
@compat abstract type AbstractSteadyStateSolution <: DESolution end

export DESolution, AbstractSteadyStateSolution

@compat abstract type AbstractODESolution <: DESolution end

export AbstractODESolution

# Needed for plot recipes
@compat abstract type AbstractDDESolution <: AbstractODESolution end
@compat abstract type  AbstractRODESolution <: AbstractODESolution end
@compat abstract type  AbstractDAESolution <: AbstractODESolution end
@compat abstract type AbstractSensitivitySolution end
@compat abstract type AbstractMonteCarloSolution <: DESolution end

export AbstractRODESolution,
       AbstractDAESolution,
       AbstractDDESolution,
       AbstractSensitivitySolution, AbstractMonteCarloSolution

# Misc
@compat abstract type Tableau end
@compat abstract type ODERKTableau <: Tableau end
@compat abstract type DECostFunction end

export Tableau, ODERKTableau, DECostFunction

@compat abstract type AbstractParameterizedFunction{isinplace} <: Function end
@compat abstract type ConstructedParameterizedFunction{isinplace} <: AbstractParameterizedFunction{isinplace} end
export AbstractParameterizedFunction, ConstructedParameterizedFunction

@compat abstract type AbstractNoiseProcess{isinplace} end
export AbstractNoiseProcess

include("utils.jl")
include("extended_functions.jl")
include("solutions/steady_state_solutions.jl")
include("solutions/ode_solutions.jl")
include("solutions/rode_solutions.jl")
include("solutions/dae_solutions.jl")
include("solutions/monte_solutions.jl")
include("solutions/solution_interface.jl")
include("tableaus.jl")
include("problems/discrete_problems.jl")
include("problems/steady_state_problems.jl")
include("problems/ode_problems.jl")
include("problems/constrained_ode_problems.jl")
include("problems/rode_problems.jl")
include("problems/sde_problems.jl")
include("problems/dae_problems.jl")
include("problems/dde_problems.jl")
include("problems/monte_problems.jl")
include("problems/problem_traits.jl")
include("callbacks.jl")
include("integrator_interface.jl")
include("parameters_interface.jl")
include("constructed_parameterized_functions.jl")
include("linear_nonlinear.jl")
include("data_array.jl")

type ConvergenceSetup{P,C}
    probs::P
    convergence_axis::C
end

function solve end
function solve! end
function init end
function step!(d::DEIntegrator) error("Integrator stepping is not implemented") end

export DEParameters, Mesh, ExplicitRKTableau, ImplicitRKTableau

export isinplace, noise_class

export solve, solve!, init, step!

export tuples, intervals

export resize!,deleteat!,addat!,full_cache,user_cache,u_cache,du_cache,
       resize_non_user_cache!,deleteat_non_user_cache!,addat_non_user_cache!,
       terminate!,
       add_tstop!,add_saveat!,set_abstol!,
       set_reltol!,get_du,get_dt,get_proposed_dt,set_proposed_dt!,u_modified!,
       savevalues!

export numargs, @def, @muladd

export construct_correlated_noisefunc

export HasJac, HastGrad, HasParamFuncs, HasParamDeriv, HasParamJac,
       HasInvJac,HasInvW, HasInvW_t, HasHes, HasInvHes, HasSyms

export has_jac, has_invjac, has_invW, has_invW_t, has_hes, has_invhes,
       has_tgrad, has_paramfuncs, has_paramderiv, has_paramjac,
       has_syms, has_analytic

export DiscreteProblem

export SteadyStateProblem, SteadyStateSolution
export ODEProblem, ODESolution
export SplitODEProblem, PartitionedODEProblem
export RODEProblem, RODESolution, SDEProblem
export SplitSDEProblem, PartitionedSDEProblem, SecondOrderSDEProblem
export DAEProblem, DAESolution
export SplitDAEProblem, PartitionedDAEProblem, ConstrainedODEProblem,
       SplitConstrainedODEProblem, PartitionedConstrainedODEProblem
export ConstantLagDDEProblem, DDEProblem

export ParameterizedFunction, DAEParameterizedFunction, DDEParameterizedFunction,
       ProbParameterizedFunction, OutputParameterizedFunction

export calculate_monte_errors

export build_solution, calculate_solution_errors!

export ConvergenceSetup

export ContinuousCallback, DiscreteCallback, CallbackSet

export initialize!

export copy_non_array_fields

export param_values, num_params, problem_new_parameters

export white_noise_func_wrapper, white_noise_func_wrapper!

export is_diagonal_noise, is_sparse_noise

export LinSolveFactorize, DEFAULT_LINSOLVE

export isautodifferentiable

end # module

ext: jl


text:
module BCC

using ..CLE
using JuLIP
using JuLIPMaterials: Vec3, Mat3, Ten43, cluster

function lattice_constants_111(sym::Symbol; calc = nothing)
   # See cell_111 for description of outputs
   if calc == nothing
      a0 = rnn(sym)
   else
      at0 = bulk(sym, cubic=true, pbc=true)
      variablecell!(at0)
      set_calculator!(at0, calc)
      minimise!(at0, verbose=0)
      a0 = sqrt(3)/2 * cell(at0)[1,1]
   end
   b0 = sqrt(8/9) * a0   # b0^2 + (a0/3)^2 = a0^2
   c0 = sqrt(3/4) * b0   # (b0/2)^2 + c0^2 = b0^2
   return a0, b0, c0
end

function cell_111(sym::Symbol; calc=nothing)
   # a0 is z-direction cell height
   # b0 is x-spacing between atoms in cell
   # c0 is y-spacing between atoms in cell
   # Full cell is size = 3b0 × 2c0 × a0
   a0, b0, c0 = lattice_constants_111(sym, calc=calc)
   X = [ [0.0, 0.0, 0.0] [b0, 0.0, a0/3] [2*b0, 0.0, 2*a0/3] [b0/2, c0, 2*a0/3] [3*b0/2, c0, 0.0] [5*b0/2, c0, a0/3] ] |> vecs
   F = diagm([3*b0, 2*c0, a0])
   at = Atoms(sym, X)
   set_defm!(at, F)
   set_pbc!(at, true)
   return at
end

"""
`screw_111(sym::Symbol, R::Float64; x0=:center, layers=1)`

Construct a circular cluster with a screw dislocation in the centre.
* `sym`: chemical symbol
* `R`: radius of cluster
* `x0`: position of screw dislocation core, relative to a lattice site
"""
function screw_111(sym::Symbol, R::Float64;
            x0 = :center, layers=1, soln = :anisotropic, calc = nothing,
            bsign = 1, Nquad=30)::AbstractAtoms
   a0, b0, c0 = lattice_constants_111(sym; calc=calc)
   x00 = JVecF( ([b0, 0, a0/3] + [b0/2, c0, 2*a0/3]) / 3 )
   if (x0 == :center) || (x0 == :centre)
      # center of mass of a triangle. (the z coordinate is irrelevant!)
      x0 = x00
   end
   # create a cluster
   atu = cell_111(sym, calc=calc)
   at = cluster(atu, R, dims=[1,2], x0=x0)
   at = at * (1,1,layers)
   # get positions  to manipulate them
   X = positions(at) |> mat
   # reference positions
   X0 = copy(X)
   X0[1,:] -= X0[1,1]+x0[1]
   X0[2,:] -= X0[2,1]+x0[2]
   set_data!(at, :X0, copy(X))
   set_data!(at, :core, X[:,1]+x0)
   set_positions!(at, X0)

   # get the screw displacement (Burgers vector = (0, 0, a0))
   if soln == :isotropic
      disl = IsoScrewDislocation3D( bsign*a0, remove_singularity = true )
      # apply to `at` and return
      X0 = vecs(X0)
      X = vecs(X)
      X .+= disl.(X0)
   elseif soln == :anisotropic
      set_pbc!(atu, true)
      set_calculator!(atu, calc)
      t = Vec3([0.0, 0.0, 1.0])
      b = Vec3([0.0, 0.0, a0*bsign])
      m0 = Vec3([1.0, 0.0, 0.0])
      ℂ = CLE.elastic_moduli(atu)
      ℂ = Ten43{Float64}(ℂ[:])
      disl = CLE.Dislocation(b, t, ℂ, Nquad = Nquad, cut = m0, remove_singularity = true)
      X0 = vecs(X0)
      X = vecs(X)
      X .+= disl.(X0)
   else
      error("unknown `soln`")
   end

   set_positions!(at, X)
   return at
end

end


text:
using Test

@test true

#@testset
#end # testset

max_stars_repo_name: Entze/jlArrowSolver

max_issues_repo_name: Entze/jlArrowSolver

max_forks_repo_name: Entze/jlArrowSolver


text:
## KL Divergence between the GP Prior and the variational distribution
function GaussianKL(model::AbstractGPModel, state)
    return mapreduce(GaussianKL, +, model.f, Zviews(model), state.kernel_matrices)
end

function GaussianKL(gp::AbstractLatent, X::AbstractVector, k_mat)
    return GaussianKL(mean(gp), pr_mean(gp, X), cov(gp), k_mat.K)
end

## See https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence#Multivariate_normal_distributions ##
function GaussianKL(
    μ::AbstractVector,
    μ₀::AbstractVector,
    Σ::Symmetric{T,Matrix{T}},
    K::Cholesky{T,Matrix{T}},
) where {T<:Real}
    return 0.5 * (logdet(K) - logdet(Σ) + tr(K \ Σ) + invquad(K, μ - μ₀) - length(μ))
end

function GaussianKL(
    μ::AbstractVector{T},
    μ₀::AbstractVector,
    Σ::Symmetric{T,Matrix{T}},
    K::AbstractMatrix{T},
) where {T<:Real}
    K
    return 0.5 * (logdet(K) - logdet(Σ) + tr(K \ Σ) + dot(μ - μ₀, K \ (μ - μ₀)) - length(μ))
end

extraKL(::AbstractGPModel{T}, ::Any) where {T} = zero(T)

"""
    extraKL(model::OnlineSVGP)

Extra KL term containing the divergence with the GP at time t and t+1
"""
function extraKL(model::OnlineSVGP{T}, state) where {T}
    return mapreduce(
        +, model.f, state.opt_state, state.kernel_matrices
    ) do gp, opt_state, kernel_mat
        prev_gp = opt_state.previous_gp
        κₐμ = kernel_mat.κₐ * mean(gp)
        KLₐ = prev_gp.prev𝓛ₐ
        KLₐ +=
            -0.5 * sum(
                trace_ABt.(
                    Ref(prev_gp.invDₐ),
                    [kernel_mat.K̃ₐ, kernel_mat.κₐ * cov(gp) * transpose(kernel_mat.κₐ)],
                ),
            )
        KLₐ += dot(prev_gp.prevη₁, κₐμ) - 0.5 * dot(κₐμ, prev_gp.invDₐ * κₐμ)
        return KLₐ
    end
end

InverseGammaKL(α, β, αₚ, βₚ) = GammaKL(α, β, αₚ, βₚ)
"""
    GammaKL(α, β, αₚ, βₚ)

KL(q(ω)||p(ω)), where q(ω) = Ga(α,β) and p(ω) = Ga(αₚ,βₚ)
"""
function GammaKL(α, β, αₚ, βₚ)
    return sum(
        (α - αₚ) .* digamma(α) .- log.(gamma.(α)) .+ log.(gamma.(αₚ)) .+
        αₚ .* (log.(β) .- log.(βₚ)) .+ α .* (βₚ .- β) ./ β,
    )
end

"""
    PoissonKL(λ, λ₀)

KL(q(ω)||p(ω)), where q(ω) = Po(ω|λ) and p(ω) = Po(ω|λ₀)
"""
function PoissonKL(λ::AbstractVector{T}, λ₀::Real) where {T}
    return λ₀ * length(λ) - (one(T) + log(λ₀)) * sum(λ) + sum(xlogx, λ)
end

"""
    PoissonKL(λ, λ₀, ψ)

KL(q(ω)||p(ω)), where q(ω) = Po(ω|λ) and p(ω) = Po(ω|λ₀) with ψ = E[log(λ₀)]
"""
function PoissonKL(
    λ::AbstractVector{<:Real}, λ₀::AbstractVector{<:Real}, ψ::AbstractVector{<:Real}
)
    # sum(λ₀) - sum(λ) + sum(xlogx, λ) - dot(λ, ψ)
    # sum(λ₀) - sum(λ) + mapreduce(xlogx, +, λ) - dot(λ, ψ)
    return sum(λ₀) - sum(λ) + sum(xlogx.(λ)) - dot(λ, ψ)
end

"""
    PolyaGammaKL(b, c, θ)

KL(q(ω)||p(ω)), where q(ω) = PG(b,c) and p(ω) = PG(b,0). θ = 𝑬[ω]
"""
function PolyaGammaKL(b, c, θ)
    return dot(b, logcosh.(0.5 * c)) - 0.5 * dot(abs2.(c), θ)
end

"""
    GIGEntropy(a, b, p)

Entropy of GIG variables with parameters a,b and p and omitting the derivative d/dpK_p cf <https://en.wikipedia.org/wiki/Generalized_inverse_Gaussian_distribution#Entropy>
"""
function GIGEntropy(a, b, p)
    sqrt_ab = sqrt.(a .* b)
    return 0.5 * (sum(log, a) - sum(log, b)) +
           mapreduce((p, s) -> log(2 * besselk(p, s)), +, p, sqrt_ab) +
           sum(
               0.5 * sqrt_ab ./ besselk.(p, sqrt_ab) .*
               (besselk.(p + 1, sqrt_ab) + besselk.(p - 1, sqrt_ab)),
           )
end


text:
@testset "Summary" begin
    gfa_result = read_from_gfa("data/gfa_sample_2.gfa")
    res = find_graph_component(gfa_result)
    result = @capture_out get_summary(res)
    s_test = "No of simple graphs: 2\nNo of lone cycles: 1\nNo of lone nodes: 1\n"
    @test result == s_test
    result_2 = @capture_out get_summary(res.graph[1])
    s_test_2 = "No of vertices: 5\nNo of edges: 4\nNo of source_node: 2\nNo of sink_node: 1\nNo of path: 0\n"
    @test result_2 == s_test_2
    result_3 = @capture_out get_summary(res.lone_cycle[1])
    s_test_3 = "No of vertices: 3\nNo of edges: 3\nNo of source_node: 1\nNo of sink_node: 0\nNo of path: 0\n"
    @test result_3 == s_test_3
    result_4 = @capture_out get_summary(res.lone_node[1])
    s_test_4 = "No of vertices: 1\nNo of edges: 0\nNo of source_node: 1\nNo of sink_node: 1\nNo of path: 0\n"
    @test result_4 == s_test_4
end


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

using .Blake3Hash

data = rand(UInt8, 2 * 8192)

sa     = Blake3Hash.SAImpl.Blake3Ctx()
sa_out = Vector{UInt8}(undef, 64)

Blake3Hash.SAImpl.update!(sa, data)
Blake3Hash.SAImpl.digest(sa, sa_out)

t     = Blake3Hash.Simd128.Blake3Ctx()
t_out = Vector{UInt8}(undef, 64)

Blake3Hash.Simd128.update!(t, data)
Blake3Hash.Simd128.digest(t, t_out)

println("E: $sa_out")
println("R: $t_out")


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

function gnuplot_init()
    global gnuplot_state

    pr = (0,0,0,0) # stdin, stdout, stderr, pid
    try
		pr = popen3(`gnuplot`)
    catch
        error("There was a problem starting up gnuplot.")
    end
    # It's possible that `popen3` runs successfully, but gnuplot exits
    # immediately. Double-check that gnuplot is running at this point.
    if Base.process_running(pr[4])
	    gnuplot_state.running = true
	    gnuplot_state.fid = pr
		# Start tasks to read and write gnuplot's pipes
		yield()  # get async tasks started (code blocks without this line)
		notify(StartPipes)
	else
        error("There was a problem starting up gnuplot.")
	end
end

# Async tasks to read/write to gnuplot's pipes.
const StartPipes = Condition()  # signal to start reading pipes

# This task reads all characters available from gnuplot's stdout.
@async while true
	wait(StartPipes)
	pout = gnuplot_state.fid[2]
	while true
		if !isopen(pout)
			break
		end
		gnuplot_state.gp_stdout = ascii(readavailable(pout))
	end
end

# This task reads all characters available from gnuplot's stderr.
@async while true
	wait(StartPipes)
	perr = gnuplot_state.fid[3]
	while true
		if !isopen(perr)
			break
		end
		gnuplot_state.gp_stderr = ascii(readavailable(perr))
	end
end

# close gnuplot pipe
function gnuplot_exit(x...)
    global gnuplot_state

    if gnuplot_state.running
        # close pipe
        close(gnuplot_state.fid[1])
        close(gnuplot_state.fid[2])
        close(gnuplot_state.fid[3])
    end
    # reset gnuplot_state
    gnuplot_state.running = false
    gnuplot_state.current = 0
    gnuplot_state.figs = Any[]
    return 0
end

# Return index to figure with handle 'c'. If no such figure exists, returns 0.
function findfigure(c)
    global gnuplot_state
    i = 0
    for j = 1:length(gnuplot_state.figs)
        if gnuplot_state.figs[j].handle == c
            i = j
            break
        end
    end
    return i
end

# convert marker string description to gnuplot's expected number
function pointtype(x::AbstractString)
    if x == "+"
        return 1
    elseif x == "x"
        return 2
    elseif x == "*"
        return 3
    elseif x == "esquare"
        return 4
    elseif x == "fsquare"
        return 5
    elseif x == "ecircle"
        return 6
    elseif x == "fcircle"
        return 7
    elseif x == "etrianup"
        return 8
    elseif x == "ftrianup"
        return 9
    elseif x == "etriandn"
        return 10
    elseif x == "ftriandn"
        return 11
    elseif x == "edmd"
        return 12
    elseif x == "fdmd"
        return 13
    end
    return 1
end

# return configuration string for a single plot
function linestr_single(conf::CurveConf)
    s = ""
    if conf.legend != ""
        s = string(s, " title '", conf.legend, "' ")
    else
        s = string(s, "notitle ")
    end
    s = string(s, " with ", conf.plotstyle, " ")
    if conf.color != ""
        s = string(s, "linecolor rgb '", conf.color, "' ")
    end
    s = string(s, "lw ", string(conf.linewidth), " ")
    # some plotstyles don't allow point specifiers
    cp = conf.plotstyle
    if cp != "lines" && cp != "impulses" && cp != "pm3d" && cp != "image" &&
        cp != "rgbimage" && cp != "boxes" && cp != "dots" && cp != "steps" &&
        cp != "fsteps" && cp != "fillsteps" && cp != "financebars"
        if conf.marker != ""
            s = string(s, "pt ", string(pointtype(conf.marker)), " ")
        end
        s = string(s, "ps ", string(conf.pointsize), " ")
    end
    return s
end

# build a string with plot commands according to configuration
function linestr(curves::Vector{CurveData}, cmd::AbstractString, file::AbstractString)
    # We have to insert "," between plot commands. One easy way to do this
    # is create the first plot command, then the rest
    # We also need to keep track of the current index (starts at zero)
    index = 0
    s = string(cmd, " '", file, "' ", " i 0 ", linestr_single(curves[1].conf))
    if length(curves) > 1
        for i in curves[2:end]
            index += 1
            s = string(s, ", '", file, "' ", " i ", string(index)," ",
                linestr_single(i.conf))
        end
    end
    return s
end

# create a Z-coordinate matrix from x, y coordinates and a function
function meshgrid(x,y,f)
    Z = zeros(length(x),length(y))
    for k = 1:length(x)
        Z[k,:] = [ f(i,j) for i=x[k], j=y ]
    end
    return Z
end

# dereference CurveConf, by adding a method to copy()
function copy(conf::CurveConf)
    new = CurveConf()
    new.legend = conf.legend
    new.plotstyle = conf.plotstyle
    new.color = conf.color
    new.marker = conf.marker
    new.linewidth = conf.linewidth
    new.pointsize = conf.pointsize
    return new
end

# dereference AxesConf
function copy(conf::AxesConf)
    new = AxesConf()
    new.title = conf.title
    new.xlabel = conf.xlabel
    new.ylabel = conf.ylabel
    new.zlabel = conf.zlabel
    new.fill = conf.fill
    new.grid = conf.grid
    new.box = conf.box
    new.axis = conf.axis
    new.xrange = conf.xrange
    new.yrange = conf.yrange
    new.zrange = conf.zrange
    return new
end

# Build a "set term" string appropriate for the terminal type
function termstring(term::AbstractString)
    global gnuplot_state
    global gaston_config

    gc = gaston_config

    if is_term_screen(term)
        ts = "set term $term $(gnuplot_state.current)"
    else
        if term == "pdf"
            s = "set term pdfcairo $(gc.print_color) "
            s = "$s font \"$(gc.print_fontface),$(gc.print_fontsize)\" "
            s = "$s fontscale $(gc.print_fontscale) "
            s = "$s linewidth $(gc.print_linewidth) "
            s = "$s size $(gc.print_size)"
        elseif term == "eps"
            s = "set term epscairo $(gc.print_color) "
            s = "$s font \"$(gc.print_fontface),$(gc.print_fontsize)\" "
            s = "$s fontscale $(gc.print_fontscale) "
            s = "$s linewidth $(gc.print_linewidth) "
            s = "$s size $(gc.print_size)"
        elseif term == "png"
            s = "set term pngcairo $(gc.print_color) "
            s = "$s font \"$(gc.print_fontface),$(gc.print_fontsize)\" "
            s = "$s fontscale $(gc.print_fontscale) "
            s = "$s linewidth $(gc.print_linewidth) "
            s = "$s size $(gc.print_size)"
        elseif term == "gif"
            s = "set term gif "
            s = "$s font $(gc.print_fontface) $(gc.print_fontsize) "
            s = "$s fontscale $(gc.print_fontscale) "
            s = "$s linewidth $(gc.print_linewidth) "
            s = "$s size $(gc.print_size)"
        elseif term == "svg"
            s = "set term svg "
            s = "$s font \"$(gc.print_fontface),$(gc.print_fontsize)\" "
            s = "$s linewidth $(gc.print_linewidth) "
            s = "$s size $(gc.print_size)"
        end
        ts = "$s \nset output \"$(gc.outputfile)\""
    end
    return ts
end

# send gnuplot the current figure's configuration
function gnuplot_send_fig_config(config)
	# fill style
	if config.fill != ""
		gnuplot_send(string("set style fill ",config.fill))
	end
	# grid
	if config.grid != ""
		if config.grid == "on"
			gnuplot_send(string("set grid"))
		else
			gnuplot_send(string("unset grid"))
		end
	end
    # legend box
    if config.box != ""
        gnuplot_send(string("set key ",config.box))
    end
    # plot title
    if config.title != ""
        gnuplot_send(string("set title '",config.title,"' "))
    end
    # xlabel
    if config.xlabel != ""
        gnuplot_send(string("set xlabel '",config.xlabel,"' "))
    end
    # ylabel
    if config.ylabel != ""
        gnuplot_send(string("set ylabel '",config.ylabel,"' "))
    end
    # zlabel
    if config.zlabel != ""
        gnuplot_send(string("set zlabel '",config.zlabel,"' "))
    end
    # axis log scale
    if config.axis != "" || config.axis != "normal"
        if config.axis == "semilogx"
            gnuplot_send("set logscale x")
        end
        if config.axis == "semilogy"
            gnuplot_send("set logscale y")
        end
        if config.axis == "loglog"
            gnuplot_send("set logscale xy")
        end
    end
    # ranges
    gnuplot_send("set autoscale")
    if config.xrange != ""
        gnuplot_send(string("set xrange ",config.xrange))
    end
    if config.yrange != ""
        gnuplot_send(string("set yrange ",config.yrange))
    end
    if config.zrange != ""
        gnuplot_send(string("set zrange ",config.zrange))
    end
end

# Validation functions.
# These functions validate that configuration parameters are valid and
# supported. They return true iff the argument validates.

# Validate terminal type.
function validate_terminal(s::AbstractString)
    supp_terms = ["qt", "wxt", "x11", "svg", "gif", "png", "pdf", "aqua", "eps"]
    if s == "aqua" && OS_NAME != :Darwin
        return false
    end
    if in(s, supp_terms)
        return true
    end
    return false
end

# Identify terminal by type: file or screen
function is_term_screen(s::AbstractString)
    screenterms = ["qt", "wxt", "x11", "aqua"]
    if in(s, screenterms)
        return true
    end
    return false
end

function is_term_file(s::AbstractString)
    screenterms = ["svg", "gif", "png", "pdf", "eps"]
    if in(s, screenterms)
        return true
    end
    return false
end

# Valid plotstyles supported by gnuplot's plot
function validate_2d_plotstyle(s::AbstractString)
    valid = ["lines", "linespoints", "points", "impulses", "boxes",
        "errorlines", "errorbars", "dots", "steps", "fsteps", "fillsteps",
        "financebars"]
    if in(s, valid)
        return true
    end
    return false
end

# Valid plotstyles supported by gnuplot's splot
function validate_3d_plotstyle(s::AbstractString)
    valid = ["lines", "linespoints", "points", "impulses", "pm3d",
            "image", "rgbimage", "dots"]
    if in(s, valid)
        return true
    end
    return false
end

# Valid axis types
function validate_axis(s::AbstractString)
    valid = ["", "normal", "semilogx", "semilogy", "loglog"]
    if in(s, valid)
        return true
    end
    return false
end

# Valid markers supported by Gaston
function validate_marker(s::AbstractString)
    valid = ["", "+", "x", "*", "esquare", "fsquare", "ecircle", "fcircle",
    "etrianup", "ftrianup", "etriandn", "ftriandn", "edmd", "fdmd"]
    if in(s, valid)
        return true
    end
    return false
end

# Validate fill style
function validate_fillstyle(s::AbstractString)
	valid = ["","empty","solid","pattern"]
	if in(s,valid)
		return true
	end
	return false
end

# Validate grid styles
function validate_grid(s::AbstractString)
	valid = ["", "on", "off"]
	if in(s,valid)
		return true
	end
	return false
end

function validate_range(s::AbstractString)

    # floating point, starting with a dot
    f1 = "[-+]?\\.\\d+([eE][-+]?\\d+)?"
    # floating point, starting with a digit
    f2 = "[-+]?\\d+(\\.\\d*)?([eE][-+]?\\d+)?"
    # floating point
    f = "($f1|$f2)"
    # autoscale directive (i.e. `*` surrounded by
    # optional bounds lb < * < ub)
    as = "(($f\\s*<\\s*)?\\*(\\s*<\\s*$f)?)"
    # full range item: a floating point, or an
    # autoscale directive, or nothing
    it = "(\\s*($as|$f)?\\s*)"

    # empty range
    er = "\\[\\s*\\]"
    # full range: two colon-separated items
    fr = "\\[$it:$it\\]"

    # range regex
    rx = Regex("^\\s*($er|$fr)\\s*\$")

    if isempty(s) || ismatch(rx, s)
        return true
    end
    return false
end

function writemime(io::IO, ::MIME"image/png", x::Figure)
	# The plot is written to /tmp/gaston-ijula.png. Read the file and
	# write it to io.
	data = open(readbytes, "$(gnuplot_state.tmpdir)gaston-ijulia.png","r")
	write(io,data)
end

# Execute command `cmd`, and return a tuple `(in, out, err, r)`, where
# `in`, `out`, `err` are pipes to the process' STDIN, STDOUT, and STDERR, and
# `r` is a process descriptor.
function popen3(cmd::Cmd)
    pin = Base.Pipe()
    out = Base.Pipe()
    err = Base.Pipe()
    r = spawn(cmd, (pin, out, err))
    Base.close_pipe_sync(out.in)
    Base.close_pipe_sync(err.in)
    Base.close_pipe_sync(pin.out)
    Base.start_reading(out.out)
    Base.start_reading(err.out)
    return (pin.in, out.out, err.out, r)
end


text:
import Pkg

atreplinit() do _
    for pkg in (:BenchmarkTools, :Revise)
        @info "Loading $pkg"
        pkgstr = String(pkg)
        if Base.identify_package(pkgstr) === nothing
            Pkg.add(pkgstr)
        end
        :(using $pkg) |> eval
    end
end

max_stars_repo_name: AnAverageHuman/dotfiles

max_issues_repo_name: AnAverageHuman/dotfiles

max_forks_repo_name: AnAverageHuman/dotfiles


text:
tests = ["testini.jl",
         "testhttp.jl",
         "testsimple.jl",
         "testmerge.jl"]

for test in tests
    include(test)
end

outfile = joinpath(@__DIR__, "confs", "out.conf")
if isfile(outfile)
    rm(outfile)
end


text:
### A Pluto.jl notebook ###
# v0.11.8

using Markdown
using InteractiveUtils

# ╔═╡ f9030e18-e484-11ea-24e0-3132afabdf83
begin
	using Pkg
	Pkg.activate(".")
	Pkg.add(PackageSpec(name="MixedModels", rev="master"))
	using MixedModels
end

# ╔═╡ 1357b670-e484-11ea-2275-cba956358da4
md"# Analysis of the sleepstudy data"

# ╔═╡ 2a6a0106-e484-11ea-2d1c-65ea0604a1e3
md"""
The `sleepstudy` data are from a study on the effects of sleep deprivation on
response time.  Eighteen subjects were allowed only 3 hours of time to sleep each night for 9 successive nights.  Their reaction time was measured each day, starting the day before the first night of sleep deprivation, when the subjects were on their regular sleep schedule.
"""

# ╔═╡ bf6eab08-e484-11ea-21b3-79b2acabd22f
md"## Loading the data"

# ╔═╡ 44d5f548-e4ae-11ea-13c2-ddb1bbed274e
md"""
First attach the MixedModels package.  In Pluto this must be done within a `begin ... end` block in a cell.
"""

# ╔═╡ 786e3abe-e4ae-11ea-0965-6bbd681d0f72
md"""
The `sleepstudy` data are one of the datasets available with recent versions of the `MixedModels` package.
"""

# ╔═╡ 66a6f3bc-e485-11ea-2b47-13a48b5d7e52
sleepstudy = MixedModels.dataset("sleepstudy")

# ╔═╡ f3ed171c-e4ab-11ea-1eac-8b0265a6c6ea
md"## Fitting an intial model"

# ╔═╡ 8a902898-e485-11ea-2d50-b3c44d666f9c
m1 = fit(MixedModel, @formula(reaction ~ 1+days+(1+days|subj)), sleepstudy)

# ╔═╡ 06560c54-e4ac-11ea-3b36-e57ee4bc3134
md"""
This model includes fixed effects for the intercept, representing the typical reaction time at the beginning of the experiment with zero days of sleep deprivation, and the slope w.r.t. days of sleep deprivation.  In this case about 250 ms. typical reaction time without deprivation and a typical increase of 10.5 ms. per day of sleep deprivation.

The random effects represent shifts from the typical behavior for each subject.The shift in the intercept has a standard deviation of about 24 ms. which would indicate a range of about 200 ms. to 300 ms. in the intercepts.  Similarly within-subject slopes would be expected to have a range of about 0 ms./day up to 20 ms./day.

The within-subject correlation of the random effects for intercept and slope is small, +0.8, indicating that a simpler model with uncorrelated random effects may be sufficient.
"""

# ╔═╡ c01a9062-e4ac-11ea-383c-471462a1e649
md"## A model with uncorrelated random effects"

# ╔═╡ 3e3b227c-e4ad-11ea-3c24-fb23d26e1171
md"""
The `zerocorr` function applied to a random-effects term creates uncorrelated vector-valued per-subject random effects.
"""

# ╔═╡ d68c3c14-e485-11ea-0b41-6d234ab4b676
m2 = fit(MixedModel, @formula(reaction ~ 1+days+zerocorr(1+days|subj)),sleepstudy)

# ╔═╡ 8e066b36-e4ad-11ea-2655-47a9a8d3c031
md"""
This model has a slightly lower log-likelihood than does `m1` but one less parameter than `m1`.   A likelihood-ratio test can be used to compare these nested models.
"""

# ╔═╡ ec436546-e485-11ea-19b8-a15f2f8247df
MixedModels.likelihoodratiotest(m2, m1)

# ╔═╡ d0856ae0-e4ab-11ea-0016-2f170c65768f
md"""
Alternatively, the AIC or BIC values can be compared.  They are on a scale where "smaller is better".  Both these model-fit statistics prefer `m2`, the simpler model.
"""

# ╔═╡ Cell order:
# ╟─1357b670-e484-11ea-2275-cba956358da4
# ╟─2a6a0106-e484-11ea-2d1c-65ea0604a1e3
# ╟─bf6eab08-e484-11ea-21b3-79b2acabd22f
# ╟─44d5f548-e4ae-11ea-13c2-ddb1bbed274e
# ╠═f9030e18-e484-11ea-24e0-3132afabdf83
# ╟─786e3abe-e4ae-11ea-0965-6bbd681d0f72
# ╠═66a6f3bc-e485-11ea-2b47-13a48b5d7e52
# ╟─f3ed171c-e4ab-11ea-1eac-8b0265a6c6ea
# ╠═8a902898-e485-11ea-2d50-b3c44d666f9c
# ╟─06560c54-e4ac-11ea-3b36-e57ee4bc3134
# ╟─c01a9062-e4ac-11ea-383c-471462a1e649
# ╟─3e3b227c-e4ad-11ea-3c24-fb23d26e1171
# ╠═d68c3c14-e485-11ea-0b41-6d234ab4b676
# ╟─8e066b36-e4ad-11ea-2655-47a9a8d3c031
# ╠═ec436546-e485-11ea-19b8-a15f2f8247df
# ╟─d0856ae0-e4ab-11ea-0016-2f170c65768f


text:
export Convection

"""
    struct Convection{N} <: Application{N}
Defines first-order linear convection equation u_t + c⋅∇u = 0
with optional Dirichlet boundary condition: u(boundary) = bc_d
"""
struct Convection{N} <: Application{N}
    c    # Convection constant
    bc_d # Dirichlet BC
    Convection{N}(c,bc_d=0.0,bc_n=0.0) where {N} = new{N}(c,bc_d,bc_n)
    function Convection{N}(param_dict::Dict{String,Any}) where {N}
        c = param_dict["c"]
        bc_d = get!(param_dict, "bc_d", 0.0)
        new{N}(c, bc_d, bc_n)
    end
end

"""
    flux_c!(flux, app::Convection, u)
Computes the convective flux for convection equation: f_c(u) = cu
"""
flux_c!(flux, app::Convection, u) = (flux .= u*app.c')

function boundary_flux_c!(flux, app::Convection, u, bc, normal_k)
    flux_k = similar(u)
    if bc == -1 # Dirichlet BC: u = bc_d
        boundary_u = fill(app.bc_d, size(u))
        flux_c!(flux_k, app, boundary_u)
    elseif bc == -2 # Neumann BC: ∇u⋅n = 0
        flux_c!(flux_k, app, u)
    end
    flux .= flux_k.*normal_k
end

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
using StatsModels
using GLM
using DataFrames
using CSV
using Lathe.preprocess: TrainTestSplit
using Plots
using Statistics

# 1. Priprema i provera podataka
# Ucitavanje podataka
data = DataFrame(CSV.File("tacke10.csv"))
# Napomena, testirati i za fajlove tacke40.csv, tacke40a.csv,
# tacke40b.csv, tacke40c.csv, tacke40d.csv i tacke40e.csv,
c =  cor(data.x, data.y)
println("Koeficijent korelacije je: $c")
if c>0.9
   println("Postoji veoma jaka veza izmedju podataka")
elseif c>0.7
   println("Postoji jaka veza izmedju podataka ")
elseif c>0.5
   println("Postoji umerena veza izmedju podataka ")
else
   println("Veza izmedju podataka nije dovoljno jaka")
end

# Podela na skup za obuku i testiranja
dataTrain, dataTest = TrainTestSplit(data,.80)
# Prikaz ulaznih podataka
# Prikaz opisa podataka
# Za prikaz strukture podataka, odkomentarisati sledeci red
# describe(data)
# Prikaz podataka u obliku grafikona
myInputPlot = scatter(dataTrain.x,dataTrain.y, title = "Tacke, pre linearne regresije", ylabel = "Y values", xlabel = "X value")
scatter!(myInputPlot, dataTest.x, dataTest.y)

# 2. Linearna regresija
# Priprema za linearnu regresiju sa skupom za obuku
fm = @formula(y ~ x)
linearRegressor = lm(fm, dataTrain)
# Testranje podataka linearnom regresijom
dataPredictedTrain = predict(linearRegressor, dataTrain)
dataPredictedTest = predict(linearRegressor, dataTest)

# Ispis podataka sa testiranja
for i in 1:length(dataPredictedTest)
   println("(X$i, Y$i) = ($(dataTest.x[i]), $(dataTest.y[i])) Predicted: Y$i $(dataPredictedTest[i])")
end
# Grafikon podataka sa testiranja
myOutputPlot = scatter(dataTrain.x, dataTrain.y, title = "Tacke, tacke posle linearne regresije", ylabel = "Y values", xlabel = "X value", legend=:bottomright)
scatter!(myOutputPlot, dataTest.x, dataTest.y)
scatter!(myOutputPlot, dataTest.x, dataPredictedTest)
scatter(myOutputPlot)
# 3. Analiza rezultata
# racunanje vrednosti r2 (R squared)
println()
rSquared = r2(linearRegressor)
println("Vrednost r squared iznosi $rSquared")
if (rSquared>0.9)
   println("Ovaj model je jako dobar za predvidjanje")
elseif (rSquared>0.7)
   println("Ovaj model je veoma dobar za predvidjanje")
elseif (rSquared>0.5)
   println("Ovaj model je relativno dobar za predvidjanje")
else
   println("Ovaj model nije dobar za predvidjanje")
end

# Racunanje gresaka za skup za obuku
# Racunanje gresaka za svaku vrednost posebno
errorsTrain = dataTrain.y-dataPredictedTrain
println()
println("Spisak svih gresaka pri obuci je: $(round.(errorsTrain; digits = 3))")
# Racunanje proseka greske
absMeanErrorTrain = mean(abs.(errorsTrain))
# Racunanje prosecne relativne greske u procentima MAPE
mapeTrain = mean(abs.(errorsTrain./dataTrain.y))
# Racunanje kvadarata gresaka za svaku vrednost posebno
errorTrainSquared = errorsTrain.*errorsTrain
# Racunanje MSE i RMSE (Root Mean Square Error)
mseTrain = mean(errorTrainSquared)
rmseTrain = sqrt(mean(errorTrainSquared))

println("Prosecna absolutna greska pri obuci je: $absMeanErrorTrain")
println("Prosecna relativna greska pri obuci je: $mapeTrain")
println("Prosek kvadarata greske pri obuci je: $mseTrain")
println("Koren proseka kvadarata greske pri obuci je: $rmseTrain")

# Racunanje gresaka za test skup
# Racunanje gresaka za svaku vrednost posebno
errorsTest = dataTest.y-dataPredictedTest
println()
println("Spisak svih gresaka pri testiranju je: $(round.(errorsTest; digits = 3))")
# Racunanje proseka greske
absMeanErrorTest = mean(abs.(errorsTest))
# Racunanje prosecne relativne greske u procentima MAPE
mapeTest = mean(abs.(errorsTest./dataTest.y))
# Racunanje kvadarata gresaka za svaku vrednost posebno
errorTestSquared = errorsTest.*errorsTest
# Racunanje MSE i RMSE (Root Mean Square Error)
mseTest = mean(errorTestSquared)
rmseTest = sqrt(mean(errorTestSquared))

println("Prosecna absolutna greska pri testiranju je: $absMeanErrorTest")
println("Prosecna relativna greska pri testiranju je: $mapeTest")
println("Prosek kvadarata greske pri testiranju je: $mseTest")
println("Koren proseka kvadarata greske pri testiranju je: $rmseTest")

# Provera da li je sistem pogodan za testiranje
if (rmseTrain<rmseTest)
   println("Sistem je dobro istreniran")
else
   println("Sistem nije dobro istreniran")
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
using DelimitedFiles

"""
This file contains a selection of tools to interact with this module to make optimisation
easier.

These include export functions for pulses an plotting functions
"""
function export_pulse(filename, pulse)
    open(filename, "w") do io
        writedlm(io, [pulse.T pulse.Ω pulse.ϕ])
    end
end

function export_pulse(filename, T, Ω, ϕ)
    open(filename, "w") do io
        writedlm(io, [T Ω ϕ])
    end
end

"""
Should really offer in place versions too but need to see how that works!
"""
function plot_pulse(pulse)
    # noice
end

function plot_pulse(T, Ω, ϕ; res = 1000)

end


text:
"""
Queue type on which a computation can be carried.
"""
abstract type QueueType end

struct QueueCompute <: QueueType end
struct QueueGraphics <: QueueType end
struct QueueTransfer <: QueueType end
struct QueueSparseBinding <: QueueType end
struct QueueVideoDecode <: QueueType end
struct QueueVideoEncode <: QueueType end

"""
Render pass execution specification for commands.
"""
abstract type RenderPassRequirement end

"""
The command can be executed inside a render pass.
"""
struct RenderPassInside <: RenderPassRequirement end

"""
The command can be executed outside a render pass.
"""
struct RenderPassOutside <: RenderPassRequirement end

"""
Function type classification.

Types:
- `FTYPE_CREATE`: constructor (functions that begin with `vkCreate`).
- `FTYPE_DESTROY`: destructor (functions that begin with `vkDestroy`).
- `FTYPE_ALLOCATE`: allocator (functions that begin with `vkAllocate`).
- `FTYPE_FREE`: deallocator (functions that begin with `vkFree`).
- `FTYPE_COMMAND`: Vulkan command (functions that begin with `vkCmd`).
- `FTYPE_QUERY`: used to query parameters, returned directly or indirectly through pointer mutation (typically, functions that begin with `vkEnumerate` and `vkGet`, but not all of them and possibly others).
- `FTYPE_OTHER`: no identified type.
"""
@enum FunctionType FTYPE_CREATE = 1 FTYPE_DESTROY FTYPE_ALLOCATE FTYPE_FREE FTYPE_COMMAND FTYPE_QUERY FTYPE_OTHER

"""
Structure type classification.

Types:
- `STYPE_CREATE_INFO`: holds constructor parameters (structures that end with `CreateInfo`).
- `STYPE_ALLOCATE_INFO`: holds allocator parameters (structures that end with `AllocateInfo`).
- `STYPE_GENERIC_INFO`: holds parameters for another function or structure (structures that end with `Info`, excluding those falling into the previous types).
- `STYPE_DATA`: usually represents user or Vulkan data.
- `STYPE_PROPERTY`: is a property returned by Vulkan in a `returnedonly` structure, usually done through `FTYPE_QUERY` type functions.
"""
@enum StructType STYPE_CREATE_INFO = 1 STYPE_ALLOCATE_INFO STYPE_GENERIC_INFO STYPE_DATA STYPE_PROPERTY

"""
Parameter requirement. Applies both to struct members and function parameters.

Requirement types:
- `OPTIONAL`: may have its default zero (or nullptr) value, acting as a sentinel value (similar to `Nothing` in Julia).
- `REQUIRED`: must be provided, no sentinel value is allowed.
- `POINTER_OPTIONAL`: is a pointer which may be null, but must have valid elements if provided.
- `POINTER_REQUIRED`: must be a valid pointer, but its elements are optional (e.g. are allowed to be sentinel values).
"""
@enum PARAM_REQUIREMENT OPTIONAL = 1 REQUIRED POINTER_OPTIONAL POINTER_REQUIRED

PARAM_REQUIREMENT(node::Node) =
    !haskey(node, "optional") || node["optional"] == "false" ? REQUIRED :
    PARAM_REQUIREMENT(findfirst(node["optional"] .== ["true", "false", "true,false", "false,true"]))

"""
Everything that a Vulkan specification can apply to: data structures, functions, parameters...
"""
abstract type Spec end

Base.broadcastable(spec::Spec) = Ref(spec)

"""
Specification for a function parameter.
"""
struct SpecFuncParam <: Spec
    "Name of the parent function."
    func::Symbol
    "Identifier."
    name::Symbol
    "Expression of its idiomatic Julia type."
    type::ExprLike
    "If constant, cannot be mutated by Vulkan functions."
    is_constant::Bool
    "Whether it must be externally synchronized before calling the function."
    is_externsync::Bool
    "[`PARAM_REQUIREMENT`](@ref) classification."
    requirement::PARAM_REQUIREMENT
    "Name of the parameter (of the same function) which represents its length. `Nothing` for non-vector types."
    len::Optional{Symbol}
    "Name of the parameters (of the same function) it is a length of."
    arglen::Vector{Symbol}
    "Whether automatic validity documentation is enabled. If false, this means that the parameter may be an exception to at least one Vulkan convention."
    autovalidity::Bool
end

"""
Specification for a function.
"""
struct SpecFunc <: Spec
    "Name of the function."
    name::Symbol
    "[`FunctionType`](@ref) classification."
    type::FunctionType
    "Return type (void if `Nothing`)."
    return_type::Optional{ExprLike}
    "Whether the function can be executed inside a render pass, outside, or both. Empty if not specified, in which case it is equivalent to both inside and outside."
    render_pass_compatibility::Vector{RenderPassRequirement}
    "Type of queues on which the function can be executed. Empty if not specified, in which case it is equivalent to being executable on all queues."
    queue_compatibility::Vector{QueueType}
    "Function parameters."
    params::StructVector{SpecFuncParam}
    success_codes::Vector{Symbol}
    error_codes::Vector{Symbol}
end

"""
Specification for a structure parameter.
"""
struct SpecStructMember <: Spec
    "Name of the parent structure."
    parent::Symbol
    "Identifier."
    name::Symbol
    "Expression of its idiomatic Julia type."
    type::ExprLike
    "If constant, cannot be mutated by Vulkan functions."
    is_constant::Bool
    "Whether it must be externally synchronized before calling any function which uses the parent structure."
    is_externsync::Bool
    "[`PARAM_REQUIREMENT`](@ref) classification."
    requirement::PARAM_REQUIREMENT
    "Name of the member (of the same structure) which represents its length. `Nothing` for non-vector types."
    len::Optional{ExprLike}
    "Name of the members (of the same structure) it is a length of."
    arglen::Vector{ExprLike}
    "Whether automatic validity documentation is enabled. If false, this means that the member may be an exception to at least one Vulkan convention."
    autovalidity::Bool
end

"""
Specification for a structure.
"""
struct SpecStruct <: Spec
    "Name of the structure."
    name::Symbol
    "[`StructType`](@ref) classification."
    type::StructType
    "Whether the structure is only returned by Vulkan functions (and never requested as input)."
    is_returnedonly::Bool
    "Name of the structures it extends, usually done through the original structures' `pNext` argument."
    extends::Vector{Symbol}
    "Structure members."
    members::StructVector{SpecStructMember}
end

"""
Specification for a union type.
"""
struct SpecUnion <: Spec
    "Name of the union type."
    name::Symbol
    "Possible types for the union."
    types::Vector{ExprLike}
    "Fields which cast the struct into the union types"
    fields::Vector{Symbol}
    "Selector values, if any, to determine the type of the union in a given context (function call for example)."
    selectors::Vector{Symbol}
    "Whether the type is only returned by Vulkan functions (and never requested as input)."
    is_returnedonly::Bool
end

"""
Specification for a constant.
"""
struct SpecConstant <: Spec
    "Name of the constant."
    name::Symbol
    "Value of the constant."
    value::Any
end

"""
Specification for an enumeration type.
"""
struct SpecEnum <: Spec
    "Name of the enumeration type."
    name::Symbol
    "Vector of possible enumeration values."
    enums::StructVector{SpecConstant}
end

"""
Specification for a bit used in a bitmask.
"""
struct SpecBit <: Spec
    "Name of the bit."
    name::Symbol
    "Position of the bit."
    position::Int
end

"""
Specification for a bitmask type that must be formed through a combination of `bits`.

Is usually an alias for a `UInt32` type which carries meaning through its bits.
"""
struct SpecBitmask <: Spec
    "Name of the bitmask type."
    name::Symbol
    "Valid bits that can be combined to form the final bitmask value."
    bits::StructVector{SpecBit}
    width::Integer
end

"""
Specification for a flag type `name` that is a type alias of `typealias`. Can be associated with a bitmask structure, in which case the `bitmask` number is set to the corresponding `SpecBitmask`.
"""
struct SpecFlag <: Spec
    "Name of the flag type."
    name::Symbol
    "The type it aliases."
    typealias::Symbol
    "Bitmask, if applicable."
    bitmask::Optional{SpecBitmask}
end

"""
Specification for handle types.

A handle may possess a parent. In this case, the handle can only be valid if its parent is valid.

Some handles are dispatchable, which means that they are represented as opaque pointers.
Non-dispatchable handles are 64-bit integer types, and may encode information directly into their value.
"""
struct SpecHandle <: Spec
    "Name of the handle type."
    name::Symbol
    "Name of the parent handle, if any."
    parent::Optional{Symbol}
    "Whether the handle is dispatchable or not."
    is_dispatchable::Bool
end

"""
Specification for an alias of the form `const <name> = <alias>`.
"""
struct SpecAlias{S<:Spec} <: Spec
    "Name of the new alias."
    name::Symbol
    "Aliased specification element."
    alias::S
end

"""
Function `func` that creates a `handle` from a create info structure `create_info_struct` passed as the value of the parameter `create_info_param`.

If `batch` is true, then `func` expects a list of multiple create info structures and will create multiple handles at once.
"""
struct CreateFunc <: Spec
    func::SpecFunc
    handle::SpecHandle
    create_info_struct::Optional{SpecStruct}
    create_info_param::Optional{SpecFuncParam}
    batch::Bool
end

"""
Function `func` that destroys a `handle` passed as the value of the parameter `destroyed_param`.

If `batch` is true, then `func` expects a list of multiple handles and will destroy all of them at once.
"""
struct DestroyFunc <: Spec
    func::SpecFunc
    handle::SpecHandle
    destroyed_param::SpecFuncParam
    batch::Bool
end

@enum ExtensionType begin
    EXTENSION_TYPE_INSTANCE
    EXTENSION_TYPE_DEVICE
    EXTENSION_TYPE_ANY
end

@enum PlatformType::UInt32 begin
    PLATFORM_NONE        = 0
    PLATFORM_XCB         = 1
    PLATFORM_XLIB        = 2
    PLATFORM_XLIB_XRANDR = 4
    PLATFORM_WAYLAND     = 8
    PLATFORM_METAL       = 16
    PLATFORM_MACOS       = 32
    PLATFORM_IOS         = 64
    PLATFORM_WIN32       = 128
    PLATFORM_ANDROID     = 256
    PLATFORM_GGP         = 512
    PLATFORM_VI          = 1024
    PLATFORM_FUCHSIA     = 2048
    PLATFORM_DIRECTFB    = 4096
    PLATFORM_SCREEN      = 8192
    PLATFORM_PROVISIONAL = 16384
end

struct SpecPlatform <: Spec
    type::PlatformType
    description::String
end

struct SpecExtension <: Spec
    name::Symbol
    type::ExtensionType
    requirements::Vector{Symbol}
    is_disabled::Bool
    author::Optional{String}
    symbols::Vector{Symbol}
    platform::PlatformType
    is_provisional::Bool
end

function Base.show(io::IO, ::MIME"text/plain", ext::SpecExtension)
    @match ext.type begin
        &EXTENSION_TYPE_INSTANCE => print(io, "Instance extension ")
        &EXTENSION_TYPE_DEVICE => print(io, "Device extension ")
        &EXTENSION_TYPE_ANY => print(io, "Extension ")
    end
    print(io, ext.name)
    inline_infos = String[]
    ext.is_provisional && push!(inline_infos, "provisional")
    ext.is_disabled && push!(inline_infos, "disabled")
    !isempty(inline_infos) && print(io, " (", join(inline_infos, ", "), ')')
    println(io)
    ext.platform ∉ (PLATFORM_NONE, PLATFORM_PROVISIONAL) && println(io, "• Platform: ", replace(string(ext.platform), "PLATFORM_" => ""), " (", first(spec_by_field(spec_platforms, :type, ext.platform)).description, ')')
    !isempty(ext.requirements) && println(io, "• Depends on: ", join(ext.requirements, ", "))
    n = length(ext.symbols)
    if n > 0
        println(io, "• $n symbols: ")
        limit = 8
        foreach(enumerate(sort(ext.symbols))) do (i, symbol)
            (i ≤ limit / 2 || i > n - limit ÷ 2) && (i == n ? print : println)(io, "  ∘ ", symbol)
            i == limit ÷ 2 && println(io, "  ⋮")
        end
    end
    if !isnothing(ext.author)
        print(io, "\n• From: ", ext.author)
        authors = spec_by_field(author_tags, :tag, ext.author)
        if !isempty(authors)
            @assert length(authors) == 1
            print(io, " (", first(authors).author, ')')
        end
    end
end

struct AuthorTag
    tag::String
    author::String
end


text:
module Utils

function pause()
    ccall(:jl_cpu_pause, Cvoid, ())
    GC.safepoint()
end

"""
    foreachthread(f, xs)

Like `foreach(f, xs)` but run `f(x)` in its own thread and try to make different
invocations of `f` as "concurrent" as possible for testing race-related bugs.
"""
function foreachthread(f, xs)
    if length(xs) > Threads.nthreads()
        error(
            "length (=",
            length(xs),
            ") of the input vector must be smaller than or equal to nthread (=",
            Threads.nthreads(),
            ")",
        )
    end
    if Threads.nthreads() == 1
        foreach(f, xs)
        return
    end

    nactives = Threads.Atomic{Int}(length(xs))
    go = Threads.Atomic{Bool}(false)
    Threads.@threads for x in xs
        if Threads.atomic_sub!(nactives, 1) > 1
            go[] = true
        else
            for _ in 1:1000_000
                go[] && break
                pause()
            end
            if !go[]
                @error "Failed to wait for go signal (threadid: $(Threads.threadid()))"
            end
        end

        f(x)
    end
end

end  # module


text:
using Oceananigans.Fields: AbstractField
using Oceananigans.Grids:
    interior_indices,
    left_halo_indices, right_halo_indices,
    underlying_left_halo_indices, underlying_right_halo_indices

# TODO: Move to Grids/grid_utils.jl

#####
##### Viewing halos
#####

west_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, left_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx), :, :) :
                      view(f.data, left_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx),
                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))

east_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, right_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx), :, :) :
                      view(f.data, right_halo_indices(LX, topology(f, 1), f.grid.Nx, f.grid.Hx),
                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))

south_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, :, left_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy), :) :
                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
                                   left_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy),
                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))

north_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, :, right_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy), :) :
                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
                                   right_halo_indices(LY, topology(f, 2), f.grid.Ny, f.grid.Hy),
                                   interior_indices(LZ, topology(f, 3), f.grid.Nz))

bottom_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, :, :, left_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz)) :
                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
                                   left_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz))

top_halo(f::AbstractField{LX, LY, LZ}; include_corners=true) where {LX, LY, LZ} =
    include_corners ? view(f.data, :, :, right_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz)) :
                      view(f.data, interior_indices(LX, topology(f, 1), f.grid.Nx),
                                   interior_indices(LY, topology(f, 2), f.grid.Ny),
                                   right_halo_indices(LZ, topology(f, 3), f.grid.Nz, f.grid.Hz))

underlying_west_halo(f, grid, location) =
    view(f.parent, underlying_left_halo_indices(location, topology(grid, 1), grid.Nx, grid.Hx), :, :)

underlying_east_halo(f, grid, location) =
    view(f.parent, underlying_right_halo_indices(location, topology(grid, 1), grid.Nx, grid.Hx), :, :)

underlying_south_halo(f, grid, location) =
    view(f.parent, :, underlying_left_halo_indices(location, topology(grid, 2), grid.Ny, grid.Hy), :)

underlying_north_halo(f, grid, location) =
    view(f.parent, :, underlying_right_halo_indices(location, topology(grid, 2), grid.Ny, grid.Hy), :)

underlying_bottom_halo(f, grid, location) =
    view(f.parent, :, :, underlying_left_halo_indices(location, topology(grid, 3), grid.Nz, grid.Hz))

underlying_top_halo(f, grid, location) =
    view(f.parent, :, :, underlying_right_halo_indices(location, topology(grid, 3), grid.Nz, grid.Hz))

#####
##### Viewing boundary grid points (used to fill other halos)
#####

left_boundary_indices(loc, topo, N, H) = 1:H
left_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty

right_boundary_indices(loc, topo, N, H) = N-H+1:N
right_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty

underlying_left_boundary_indices(loc, topo, N, H) = 1+H:2H
underlying_left_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty

underlying_right_boundary_indices(loc, topo, N, H) = N+1:N+H
underlying_right_boundary_indices(::Type{Nothing}, topo, N, H) = 1:0 # empty

underlying_west_boundary(f, grid, location) =
    view(f.parent, underlying_left_boundary_indices(location, topology(grid, 1), grid.Nx, grid.Hx), :, :)

underlying_east_boundary(f, grid, location) =
    view(f.parent, underlying_right_boundary_indices(location, topology(grid, 1), grid.Nx, grid.Hx), :, :)

underlying_south_boundary(f, grid, location) =
    view(f.parent, :, underlying_left_boundary_indices(location, topology(grid, 2), grid.Ny, grid.Hy), :)

underlying_north_boundary(f, grid, location) =
    view(f.parent, :, underlying_right_boundary_indices(location, topology(grid, 2), grid.Ny, grid.Hy), :)

underlying_bottom_boundary(f, grid, location) =
    view(f.parent, :, :, underlying_left_boundary_indices(location, topology(grid, 3), grid.Nz, grid.Hz))

underlying_top_boundary(f, grid, location) =
    view(f.parent, :, :, underlying_right_boundary_indices(location, topology(grid, 3), grid.Nz, grid.Hz))


text:
function save_jld2(
    file_path::AbstractString,
    group_path::AbstractString,
    groups...,
)
    path, file = splitdir(file_path)
    mkpath(path)
    jldopen(file_path, "a+") do f
        for (n, v) in Iterators.partition(groups, 2)
            try
                f[group_path * "/" * n] = v
            catch
                @warn "Data name: $(n) already exists. Data not saved."
            end
        end
    end
end

function load_jld2(
    file_path::AbstractString,
    group_path::AbstractString,
    groups...,
)
    names = [group_path * "/" * i for i in groups]
    FileIO.load(file_path, names...)
end

function check_jld2(
    file_path::AbstractString,
    group_path::AbstractString,
    group_name,
)
    group = group_path * "/" * group_name
    res = jldopen(file_path, "r") do f
        haskey(f, group)
    end
end


text:
using DataStructures

abstract gameState

abstract function draw(g::gameState,dt::Real)
abstract function update(g::gameState,dt::Real)
abstract function handleInput(g::gameState)

type game
    stack::Stack
    window
end

abstract function pushState(g::game,state::gameState)
abstract function popState(g::game)
abstract function changeState(g:game,state::gameState)
abstract function peekState(g:game)

abstract function gameLoop(g:game)

max_stars_repo_name: JKSano/netgameproj

max_issues_repo_name: JKSano/netgameproj

max_forks_repo_name: JKSano/netgameproj


text:
using LinearAlgebra
using SparseArrays
using Plots
using ForwardDiff

"This code applies the method of manufactured solutions to the 2D Laplace's equation
with nonzero Dirichlet BCs using a finite difference method."

m = 100 # number of points

uexact(x,y) = log(2+sin(pi*x)*sin(pi*y))
dudx_exact(x,y) = ForwardDiff.derivative(x->uexact(x,y),x)
dudy_exact(x,y) = ForwardDiff.derivative(y->uexact(x,y),y)
du2dx2_exact(x,y) = ForwardDiff.derivative(x->dudx_exact(x,y),x)
du2dy2_exact(x,y) = ForwardDiff.derivative(y->dudy_exact(x,y),y)
f(x,y) = -(du2dx2_exact(x,y) + du2dy2_exact(x,y))

u_left(y)   = uexact(-1,y)
u_right(y)  = uexact(1,y)
u_bottom(x) = uexact(x,-1)
u_top(x)    = uexact(x,1)
uBC = (u_left,u_right,u_bottom,u_top)

function solve(m,f,uBC)
    u_left,u_right,u_bottom,u_top = uBC
    x = LinRange(-1,1,m+2)
    y = LinRange(-1,1,m+2)
    xint = x[2:end-1]
    yint = y[2:end-1]
    h = x[2]-x[1]

    id(i,j,m) = i + (j-1)*m # return global index from i,j

    function FD_matrix_2D(h,m)
        A = spzeros(m*m,m*m)
        for i = 1:m, j = 1:m # loop thru x,y indices
            A[id(i,j,m),id(i,j,m)] = 4/h^2
            if i > 1 # avoids leftmost line of nodes
                A[id(i,j,m),id(i-1,j,m)] = -1/h^2 # x-derivative
            end
            if i < m # avoids rightmost line of nodes
                A[id(i,j,m),id(i+1,j,m)] = -1/h^2
            end
            if j > 1
                A[id(i,j,m),id(i,j-1,m)] = -1/h^2 # y-derivative
            end
            if j < m
                A[id(i,j,m),id(i,j+1,m)] = -1/h^2
            end
        end
        return A
    end

    A = FD_matrix_2D(h,m)

    b = vec(f.(xint,yint'))

    # impose non-zero Dirichlet BCs
    for i = 1:m, j = 1:m
        if i==1
            b[id(i,j,m)] += u_left(yint[j])/h^2
        end
        if i==m
            b[id(i,j,m)] += u_right(yint[j])/h^2
        end
        if j==1
            b[id(i,j,m)] += u_bottom(xint[i])/h^2
        end
        if j==m
            b[id(i,j,m)] += u_top(xint[i])/h^2
        end
    end

    # solve Au = f(xi,yj)
    u = A\b
    return u,h,xint,yint
end

errvec = []
hvec = []
for m in [4 8 16 32 64 128]
    u,h,xint,yint = solve(m,f,uBC)
    err = maximum(abs.(u .- vec(uexact.(xint,yint'))))
    append!(errvec,err)
    append!(hvec,h)
end
plot(hvec,errvec,mark=:dot,xaxis=:log,yaxis=:log)
plot!(hvec,hvec.^2,ls=:dash,xaxis=:log,yaxis=:log)
# contourf(xint,yint,))


text:
include("common.jl")

@testset "Thermodynamics - Performance pθq constructor" begin

    trial = BenchmarkTools.@benchmark $thermo_state_ρeq_newton()
    show(stdout, MIME("text/plain"), trial)
    trial = BenchmarkTools.@benchmark $thermo_state_ρeq_regula_falsi()
    show(stdout, MIME("text/plain"), trial)
    trial = BenchmarkTools.@benchmark $thermo_state_ρeq_secant()
    show(stdout, MIME("text/plain"), trial)
    trial = BenchmarkTools.@benchmark $thermo_state_pθq()
    show(stdout, MIME("text/plain"), trial)
    trial = BenchmarkTools.@benchmark $thermo_state_pTq()
    show(stdout, MIME("text/plain"), trial)

end


text:
@testset "points" begin
    @testset "order" begin
        p1 = Point2D(1.0, 2.0)
        p2 = Point2D(1.1, 2.1)
        p3 = Point2D(1.1, 2.2)

        @test p1 < p2
    end

    # @testset "conversion" begin
    #     point = [1.0 ; 2.0]
    #     p = convert(Point2D, point)
    #     p1 = Point2D(1.0, 2.0)
    #     p2 = Point2D(1.1, 2.1)
    #     @test p == p1
    #     @test p != p2
    # end
end

max_stars_repo_name: helgith/GeometricStructures

max_issues_repo_name: helgith/GeometricStructures

max_forks_repo_name: helgith/GeometricStructures


text:
using HorizonSideRobots

function chess_marking(r::Robot)
    back_track = go_to_bottom_left_corner_return_path(r)
    mark_current = !(isodd(length(back_track))) # Если длина back_track чётная, красим левую нижнюю клетку
    side = Ost
    while !(isborder(r, Nord) && isborder(r, side))
        if mark_current
            putmarker!(r)
        end

        if (isborder(r, Ost) || isborder(r, West)) && !(isborder(r, Sud) && isborder(r, West)) && !(isborder(r, Nord))
            move!(r, Nord)
            mark_current = !mark_current
            if mark_current
                putmarker!(r)
            end
            side = HorizonSide(mod(Int(side) + 2, 4))
        end

        if !(isborder(r, side))
            move!(r, side)
            mark_current = !mark_current
        end
    end
    if mark_current
        putmarker!(r)
    end

    go_to_bottom_left_corner_return_path(r)
    back_by_path_o_sides(r, back_track)
end

function go_to_bottom_left_corner_return_path(r::Robot)::Array
    path = []
    while (isborder(r, Sud) == false) || (isborder(r, West) == false)
        if isborder(r, Sud) == false
            move!(r, Sud)
            push!(path, Nord)
        end
        if isborder(r, West) == false
            move!(r, West)
            push!(path, Ost)
        end
    end
    return path
end

function back_by_path_o_sides(r, path::Array)
    n = length(path)
    while n > 0
        move!(r, path[n])
        n = n - 1
    end
end


text:
struct SRIterativeCache{Tc,T,G,E,TT}
    S::T
    S_c::Tc
    F::G
    ∇out::E

    time_evolving::TT
end

function _sr_iterative_cache(algo::SR, prob, net, par_cache)
    T = eltype(trainable_first(net))
    g = grad_cache(T, net)
    gv = vec_data(g)
    if algo.use_fullmat
        S = tuple([similar(gv, T, length(gv), length(gv))
                        for gv in vec_data(g)]...)
        if T <: Real
            Sc = tuple([similar(gv, Complex{T}, length(gv), length(gv))
                            for gv in vec_data(g)]...)
        else
            Sc = S
        end
    else # use O
        S = tuple([SrMatrix(T, similar(gv, T, length(gv), 1), 0, par_cache)
                    for gv in vec_data(g)]...)
        Sc = tuple([nothing for gv in vec_data(g)]...)
    end
    return SRIterativeCache(S, Sc, g, grad_cache(T, net),
                            algo.time_evo)
end

function setup_algorithm!(g::SRIterativeCache, ∇C, Ô, par_cache)
    α = g.time_evolving isa Val{true} ? true * im : true

    O = Ô
    for (Sc, S, F) in zip(g.S_c, g.S, vec_data(g.F))
        if S isa SrMat
            init!(S, O)
            if eltype(F) <: Real
                F .= real.(α .* ∇C)
            else
                F .= α .* ∇C
            end
        else
            T = eltype(S)
            N = size(Ô, 2)

            # TODO MPI Sync Sc
            mul!(Sc, O, O')

            if T <: Real
                S .= real.(Sc ./ N)
                F .= real.(α  .* ∇C)
            else
                # We take the conjugate because O is actually Oconj matrix respect
                # to the standard SR implementation..
                S .= conj.(Sc) ./ N
                F .= α .* ∇C
            end

            workers_mean!(S, par_cache)
        end
    end

end

function precondition!(data::SRIterativeCache, params::SR, iter_n)
    return _precondition!(data, params, iter_n)
end

function _precondition!(data::SRIterativeCache, params::SR, iter_n)
    ϵ = params.sr_diag_shift
    success = true

    for (Δw, S, F) in zip(vec_data(data.∇out), data.S, vec_data(data.F))
        # new matrix
        if params.precondition_type == sr_none
            Sprecond = S
        elseif params.precondition_type == sr_shift
            Sprecond = S + convert(eltype(S), ϵ) * I
        elseif params.precondition_type == sr_multiplicative
            λ0 = params.λ0
            b = params.b
            λmin = params.λmin
            #vv = eigvals(S)
            #println("regol $(max(λ0*b^iter_n, λmin)) --> $(minimum(vv)), a $(maximum(vv))")
            λ = convert(eltype(S), max(λ0 * b^iter_n, λmin))
            Sprecond = S + λ * Diagonal(diag(S))

        end

        if params.algorithm == sr_qlp
            x, hist = minresqlp(
                Sprecond,
                F,
                maxiter = size(S, 2) * 10,
                log = true,
                verbose = false,
                tol = params.sr_precision,
            )
        elseif params.algorithm == sr_minres
            x, hist = minres(
                Sprecond,
                F,
                maxiter = size(S, 2) * 10,
                log = true,
                verbose = false,
                tol = params.sr_precision)
        elseif params.algorithm == sr_lsq
            x, hist = lsqr(
                Sprecond,
                F,
                maxiter = size(S, 2) * 10,
                log = true,
                verbose = false,
            )
        elseif params.algorithm == sr_cg
            x, hist = cg(
                Sprecond,
                F,
                maxiter = size(S, 2) * 10,
                log = true,
                verbose = false,
                tol = params.sr_precision,
            )
        end
            #x, hist = cg(S.+ ϵ*I, F, maxiter=size(S,2)*10, log=true, verbose=true, tol=10e-10)
        if eltype(Δw) <: Real
            Δw .= real.(x)
        else
            Δw .= x
        end
        add_iters = 1
        while !hist.isconverged
            println("minresqlp not conerged. Additional $(size(S,2)*10) iters for the $add_iters time.")
            x, hist = minresqlp(
                Δw,
                Sprecond,
                F,
                maxiter = size(S, 2) * 10,
                log = true,
            )
            Δw .= x
            add_iters += 1
            add_iters > 5 && break
        end
        if add_iters > 5 && !hist.isconverged
            success = false
            Δw .= 0.0
        end
    end
    return data.∇out
end


text:
# definition of samplers and random generation

# allows to call `Sampler` only when the the arg isn't a Sampler itself
sampler(::Type{RNG}, X,          n::Repetition=Val(Inf)) where {RNG<:AbstractRNG} =
    Sampler(RNG, X, n)

sampler(::Type{RNG}, ::Type{X},  n::Repetition=Val(Inf)) where {RNG<:AbstractRNG,X} =
    Sampler(RNG, X, n)

sampler(::Type{RNG}, X::Sampler, n::Repetition=Val(Inf)) where {RNG<:AbstractRNG} = X

sampler(rng::AbstractRNG, X, n::Repetition=Val(Inf)) = sampler(typeof(rng), X, n)

## defaults

### 0-arg

make() = make(Float64)

### type: handles e.g. rand(make(Int))

# rand(rng, ::SamplerType{Make0{X}}) should not be overloaded, as make(T)
# has this special pass-thru Sampler defined below

Sampler(::Type{RNG}, ::Make0{X}, n::Repetition) where {RNG<:AbstractRNG,X} =
    Sampler(RNG, X, n)

### object: handles e.g. rand(make(1:3))

# make(x) where x isn't a type should create a distribution d such that rand(x) is
# equivalent to rand(d); so we need to overload make(x) to give a different type
# than Make1{gentype(x)}, as we can't simply define a generic pass-thru Sampler for
# that Make1 type (because users expect the default to be a SamplerTrivial{<:Make1{...}},
# and might want to define only a rand method on this SamplerTrivial)

# like Make1
struct MakeWrap{T,X} <: Distribution{T}
    x::X
end

# make(::Type) is intercepted in distribution.jl
make(x) = MakeWrap{gentype(x),typeof(x)}(x)

Sampler(::Type{RNG}, x::MakeWrap, n::Repetition) where {RNG<:AbstractRNG} =
    Sampler(RNG, x.x, n)

## Uniform

Sampler(::Type{RNG}, d::Union{UniformWrap,UniformType}, n::Repetition
        ) where {RNG<:AbstractRNG} =
            Sampler(RNG, d[], n)

## floats

### fall-back on Random definitions

for CO in (:CloseOpen01, :CloseOpen12)
    @eval begin
        Sampler(::Type{RNG}, ::$CO{T}, n::Repetition) where {RNG<:AbstractRNG,T} =
            Sampler(RNG, Random.$CO{T}(), n)

        Sampler(::Type{<:AbstractRNG}, ::$CO{BigFloat}, ::Repetition) =
            Random.SamplerBigFloat{Random.$CO{BigFloat}}(precision(BigFloat))
    end
end

### new intervals 01

# TODO: optimize for BigFloat

for CO = (:OpenClose01, :OpenOpen01, :CloseClose01)
    @eval Sampler(::Type{RNG}, I::$CO{T}, n::Repetition) where {RNG<:AbstractRNG,T} =
              SamplerSimple(I, CloseOpen01(T))
end

rand(r::AbstractRNG, sp::SamplerSimple{OpenClose01{T}}) where {T} =
    one(T) - rand(r, sp.data)

rand(r::AbstractRNG, sp::SamplerSimple{OpenOpen01{T}}) where {T} =
    while true
        x = rand(r, sp.data)
        x != zero(T) && return x
    end

# optimizations (TODO: optimize for BigFloat too)

rand(r::AbstractRNG, sp::SamplerSimple{OpenOpen01{Float64}}) =
    reinterpret(Float64, reinterpret(UInt64, rand(r, sp.data)) | 0x0000000000000001)

rand(r::AbstractRNG, sp::SamplerSimple{OpenOpen01{Float32}}) =
    reinterpret(Float32, reinterpret(UInt32, rand(r, sp.data)) | 0x00000001)

rand(r::AbstractRNG, sp::SamplerSimple{OpenOpen01{Float16}}) =
    reinterpret(Float16, reinterpret(UInt16, rand(r, sp.data)) | 0x0001)

# prevfloat(T(2)) - 1 for IEEEFloat
upper01(::Type{Float64}) = 0.9999999999999998
upper01(::Type{Float32}) = 0.9999999f0
upper01(::Type{Float16}) = Float16(0.999)
upper01(::Type{BigFloat}) = prevfloat(one(BigFloat))

rand(r::AbstractRNG, sp::SamplerSimple{CloseClose01{T}}) where {T} =
    rand(r, sp.data) / upper01(T)

### CloseOpenAB

for (CO, CO01) = (CloseOpenAB => CloseOpen01,
                  OpenCloseAB => OpenClose01,
                  CloseCloseAB => CloseClose01,
                  OpenOpenAB => OpenOpen01)

    @eval Sampler(::Type{RNG}, d::$CO{T}, n::Repetition) where {RNG<:AbstractRNG,T} =
        SamplerTag{$CO{T}}((a=d.a, d=d.b - d.a, sp=Sampler(RNG, $CO01{T}(), n)))

    @eval rand(rng::AbstractRNG, sp::SamplerTag{$CO{T}}) where {T} =
        sp.data.a + sp.data.d  * rand(rng, sp.data.sp)
end

## Normal & Exponential

rand(rng::AbstractRNG, ::SamplerTrivial{Normal01{T}}) where {T<:NormalTypes} =
    randn(rng, T)

Sampler(::Type{RNG}, d::Normalμσ{T}, n::Repetition) where {RNG<:AbstractRNG,T} =
    SamplerSimple(d, Sampler(RNG, Normal(T), n))

rand(rng::AbstractRNG, sp::SamplerSimple{Normalμσ{T},<:Sampler}) where {T} =
    sp[].μ + sp[].σ  * rand(rng, sp.data)

rand(rng::AbstractRNG, ::SamplerTrivial{Exponential1{T}}) where {T<:AbstractFloat} =
    randexp(rng, T)

Sampler(::Type{RNG}, d::Exponentialθ{T}, n::Repetition) where {RNG<:AbstractRNG,T} =
    SamplerSimple(d, Sampler(RNG, Exponential(T), n))

rand(rng::AbstractRNG, sp::SamplerSimple{Exponentialθ{T},<:Sampler}) where {T} =
    sp[].θ * rand(rng, sp.data)

## Bernoulli

Sampler(::Type{RNG}, b::Bernoulli, n::Repetition) where {RNG<:AbstractRNG} =
    SamplerTag{typeof(b)}(b.p+1.0)

rand(rng::AbstractRNG, sp::SamplerTag{Bernoulli{T}}) where {T} =
    ifelse(rand(rng, CloseOpen12()) < sp.data, one(T), zero(T))

## Categorical

Sampler(::Type{RNG}, c::Categorical, n::Repetition) where {RNG<:AbstractRNG} =
    SamplerSimple(c, Sampler(RNG, CloseOpen(), n))

# unfortunately requires @inline to avoid allocating
@inline rand(rng::AbstractRNG, sp::SamplerSimple{Categorical{T}}) where {T} =
    let c = rand(rng, sp.data)
        T(findfirst(x -> x >= c, sp[].cdf))
    end

# NOTE:
# if length(cdf) is somewhere between 150 and 200, the following gets faster:
#   T(searchsortedfirst(sp[].cdf, rand(rng, sp.data)))

## random elements from pairs

#= disabled in favor of a special meaning for pairs

Sampler(RNG::Type{<:AbstractRNG}, t::Pair, n::Repetition) =
    SamplerSimple(t, Sampler(RNG, Bool, n))

rand(rng::AbstractRNG, sp::SamplerSimple{<:Pair}) =
    @inbounds return sp[][1 + rand(rng, sp.data)]
=#

## composite types

### sampler for pairs and complex numbers

maketype(::Type{Pair},              x, y)             = Pair{val_gentype(x), val_gentype(y)}
maketype(::Type{Pair{X}},           _, y) where {X}   = Pair{X, val_gentype(y)}
maketype(::Type{Pair{X,Y} where X}, x, _) where {Y}   = Pair{val_gentype(x), Y}
maketype(::Type{Pair{X,Y}},         _, _) where {X,Y} = Pair{X,Y}

maketype(::Type{Complex},    x) = Complex{val_gentype(x)}
maketype(::Type{T}, _) where {T<:Complex} = T

maketype(::Type{Complex},    x, y) = Complex{promote_type(val_gentype(x), val_gentype(y))}
maketype(::Type{T}, _, _) where {T<:Complex} = T

function Sampler(::Type{RNG}, u::Make2{T}, n::Repetition
                 ) where RNG<:AbstractRNG where T <: Union{Pair,Complex}
    sp1 = sampler(RNG, u[1], n)
    sp2 = u[1] == u[2] ? sp1 : sampler(RNG, u[2], n)
    SamplerTag{Cont{T}}((sp1, sp2))
end

rand(rng::AbstractRNG, sp::SamplerTag{Cont{T}}) where {T<:Union{Pair,Complex}} =
    T(rand(rng, sp.data[1]), rand(rng, sp.data[2]))

#### additional convenience methods

# rand(Pair{A,B}) => rand(make(Pair{A,B}, A, B))
Sampler(::Type{RNG}, ::Type{Pair{A,B}}, n::Repetition) where {RNG<:AbstractRNG,A,B} =
    Sampler(RNG, make(Pair{A,B}, A, B), n)

# rand(make(Complex, x)) => rand(make(Complex, x, x))
Sampler(::Type{RNG}, u::Make1{T}, n::Repetition) where {RNG<:AbstractRNG,T<:Complex} =
    Sampler(RNG, make(T, u[1], u[1]), n)

# rand(Complex{T}) => rand(make(Complex{T}, T, T)) (redundant with implem in Random)
Sampler(::Type{RNG}, ::Type{Complex{T}}, n::Repetition) where {RNG<:AbstractRNG,T<:Real} =
    Sampler(RNG, make(Complex{T}, T, T), n)

### sampler for tuples

#### "simple scalar" (non-make) version

Sampler(::Type{RNG}, ::Type{T}, n::Repetition
        ) where {RNG<:AbstractRNG,T<:Union{Tuple,NamedTuple}} =
            Sampler(RNG, make(T), n)

#### make

# implement make(Tuple, S1, S2...), e.g. for rand(make(Tuple, Int, 1:3)),
# and       make(NTuple{N}, S)

_maketype(::Type{T}) where {T<:Tuple} =
    T === Tuple ?
        Tuple{} :
    T === NTuple ?
        Tuple{} :
    T isa UnionAll && Type{T} <: Type{NTuple{N}} where N ?
        T{default_gentype(Tuple)} :
        T

function _maketype(::Type{T}, args...) where T <: Tuple
    types = [t <: Type ? t.parameters[1] : gentype(t) for t in args]
    TT = T === Tuple ?
        Tuple{types...} :
    _isNTuple(T, args...) ?
        (T isa UnionAll ? Tuple{fill(types[1], fieldcount(T))...} : T ) :
        T
    TT
end

_isNTuple(::Type{T}, args...) where {T<:Tuple} =
    length(args) == 1 && T !== Tuple && (
        T <: NTuple || !isa(T, UnionAll)) # !isa(Tuple, UnionAll) !!

@generated function _make(::Type{T}, args...) where T <: Tuple
    isempty(args) && return :(Make0{$(_maketype(T))}())
    TT = _maketype(T, args...)
    samples = [t <: Type ? :(UniformType{$(t.parameters[1])}()) :
               :(args[$i]) for (i, t) in enumerate(args)]
    if _isNTuple(T, args...)
        :(Make1{$TT}($(samples[1])))
    else
        quote
            if T !== Tuple && fieldcount(T) != length(args)
                throw(ArgumentError("wrong number of provided argument with $T (should be $(fieldcount(T)))"))
            else
                Make1{$TT}(tuple($(samples...)))
            end
        end
    end
end

make(::Type{T}, args...) where {T<:Tuple} = _make(T, args...)

# make(Tuple, X, n::Integer)

default_sampling(::Type{Tuple}) = Uniform(Float64)

make(::Type{Tuple}, X,         n::Integer)           = make(NTuple{Int(n)}, X)
make(::Type{Tuple}, ::Type{X}, n::Integer) where {X} = make(NTuple{Int(n)}, X)

make(::Type{Tuple}, n::Integer) = make(Tuple, default_sampling(Tuple), Int(n))

# NTuple{N,T} where N
make(::Type{NTuple{N,T} where N},            n::Integer) where {T}   = make(NTuple{Int(n),T})
make(::Type{NTuple{N,T} where N}, X,         n::Integer) where {T}   = make(NTuple{Int(n),T}, X)
make(::Type{NTuple{N,T} where N}, ::Type{X}, n::Integer) where {T,X} = make(NTuple{Int(n),T}, X)

# disambiguate

for Tupl = (Tuple, NamedTuple)
    @eval begin
        make(::Type{T}, X)         where {T<:$Tupl}   = _make(T, X)
        make(::Type{T}, ::Type{X}) where {T<:$Tupl,X} = _make(T, X)

        make(::Type{T}, X,         Y)         where {T<:$Tupl}     = _make(T, X, Y)
        make(::Type{T}, ::Type{X}, Y)         where {T<:$Tupl,X}   = _make(T, X, Y)
        make(::Type{T}, X,         ::Type{Y}) where {T<:$Tupl,Y}   = _make(T, X, Y)
        make(::Type{T}, ::Type{X}, ::Type{Y}) where {T<:$Tupl,X,Y} = _make(T, X, Y)

        make(::Type{T}, X,         Y,         Z)         where {T<:$Tupl}       = _make(T, X, Y, Z)
        make(::Type{T}, ::Type{X}, Y,         Z)         where {T<:$Tupl,X}     = _make(T, X, Y, Z)
        make(::Type{T}, X,         ::Type{Y}, Z)         where {T<:$Tupl,Y}     = _make(T, X, Y, Z)
        make(::Type{T}, ::Type{X}, ::Type{Y}, Z)         where {T<:$Tupl,X,Y}   = _make(T, X, Y, Z)
        make(::Type{T}, X,         Y,         ::Type{Z}) where {T<:$Tupl,Z}     = _make(T, X, Y, Z)
        make(::Type{T}, ::Type{X}, Y,         ::Type{Z}) where {T<:$Tupl,X,Z}   = _make(T, X, Y, Z)
        make(::Type{T}, X,         ::Type{Y}, ::Type{Z}) where {T<:$Tupl,Y,Z}   = _make(T, X, Y, Z)
        make(::Type{T}, ::Type{X}, ::Type{Y}, ::Type{Z}) where {T<:$Tupl,X,Y,Z} = _make(T, X, Y, Z)
    end
end

#### Sampler for general tuples

@generated function Sampler(::Type{RNG}, c::Make1{T,X}, n::Repetition
                            ) where {RNG<:AbstractRNG,T<:Tuple,X<:Tuple}
    @assert fieldcount(T) == fieldcount(X)
    sps = [:(sampler(RNG, c[1][$i], n)) for i in 1:length(T.parameters)]
    :(SamplerTag{Cont{T}}(tuple($(sps...))))
end

@generated function Sampler(::Type{RNG}, ::Make0{T}, n::Repetition
                            ) where {RNG<:AbstractRNG,T<:Tuple}
    d = Dict{DataType,Int}()
    sps = []
    for t in T.parameters
        i = get(d, t, nothing)
        if i === nothing
            push!(sps, :(Sampler(RNG, $t, n)))
            d[t] = length(sps)
        else
            push!(sps, Val(i))
        end
    end
    :(SamplerTag{Cont{T}}(tuple($(sps...))))
end

@generated function rand(rng::AbstractRNG, sp::SamplerTag{Cont{T},S}) where {T<:Tuple,S<:Tuple}
    @assert fieldcount(T) == fieldcount(S)
    rands = []
    for i = 1:fieldcount(T)
        j = fieldtype(S, i) <: Val ?
              fieldtype(S, i).parameters[1] :
              i
        push!(rands, :(convert($(fieldtype(T, i)),
                               rand(rng, sp.data[$j]))))
    end
    :(tuple($(rands...)))
end

#### for "NTuple-like"

# should catch Tuple{Integer,Integer} which is not NTuple, or even Tuple{Int,UInt}, when only one sampler was passed

Sampler(::Type{RNG}, c::Make1{T,X}, n::Repetition) where {RNG<:AbstractRNG,T<:Tuple,X} =
    SamplerTag{Cont{T}}(sampler(RNG, c[1], n))

@generated function rand(rng::AbstractRNG, sp::SamplerTag{Cont{T},S}) where {T<:Tuple,S<:Sampler}
    rands = [:(convert($(T.parameters[i]), rand(rng, sp.data))) for i in 1:fieldcount(T)]
    :(tuple($(rands...)))
end

### named tuples

make(::Type{T}, args...) where {T<:NamedTuple} = _make(T, args...)

_make(::Type{NamedTuple{}}) = Make0{NamedTuple{}}()

@generated function _make(::Type{NamedTuple{K}}, X...) where {K}
    if length(X) <= 1
        NT = NamedTuple{K,_maketype(NTuple{length(K)}, X...)}
        :(Make1{$NT}(make(NTuple{length(K)}, X...)))
    else
        NT = NamedTuple{K,_maketype(Tuple, X...)}
        :(Make1{$NT}(make(Tuple, X...)))
    end
end

function _make(::Type{NamedTuple{K,V}}, X...) where {K,V}
    Make1{NamedTuple{K,V}}(make(V, X...))
end

# necessary to avoid circular defintions
Sampler(::Type{RNG}, m::Make0{NamedTuple}, n::Repetition) where {RNG<:AbstractRNG} =
    SamplerType{NamedTuple}()

Sampler(::Type{RNG}, m::Make1{T}, n::Repetition) where {RNG<:AbstractRNG, T <: NamedTuple} =
    SamplerTag{Cont{T}}(Sampler(RNG, m[1] , n))

rand(rng::AbstractRNG, sp::SamplerType{NamedTuple{}}) = NamedTuple()

rand(rng::AbstractRNG, sp::SamplerTag{Cont{T}}) where T <: NamedTuple =
    T(rand(rng, sp.data))

## collections

### sets/dicts

const SetDict = Union{AbstractSet,AbstractDict}

make(::Type{T}, X,         n::Integer) where {T<:SetDict}   = Make2{maketype(T, X, n)}(X , Int(n))
make(::Type{T}, ::Type{X}, n::Integer) where {T<:SetDict,X} = Make2{maketype(T, X, n)}(X , Int(n))
make(::Type{T},            n::Integer) where {T<:SetDict}   = make(T, default_sampling(T), Int(n))

Sampler(::Type{RNG}, c::Make2{T}, n::Repetition) where {RNG<:AbstractRNG,T<:SetDict} =
    SamplerTag{Cont{T}}((sp = sampler(RNG, c[1], n),
                         len = c[2]))

function rand(rng::AbstractRNG, sp::SamplerTag{Cont{S}}) where {S<:SetDict}
    # assuming S() creates an empty set/dict
    s = sizehint!(S(), sp.data[2])
    _rand!(rng, s, sp.data.len, sp.data.sp)
end

### sets

default_sampling(::Type{<:AbstractSet}) = Uniform(Float64)
default_sampling(::Type{<:AbstractSet{T}}) where {T} = Uniform(T)

#### Set

maketype(::Type{Set},    X, _)           = Set{val_gentype(X)}
maketype(::Type{Set{T}}, _, _) where {T} = Set{T}

### BitSet

default_sampling(::Type{BitSet}) = Uniform(Int8) # almost arbitrary, may change

maketype(::Type{BitSet}, _, _) = BitSet

### dicts

# K,V parameters are necessary here
maketype(::Type{D}, _,      ::Integer) where {K,V,D<:AbstractDict{K,V}} = D
maketype(::Type{D}, ::Type, ::Integer) where {K,V,D<:AbstractDict{K,V}} = D

#### Dict/ImmutableDict

for D in (Dict, Base.ImmutableDict)
    @eval begin
        # again same inference bug
        # TODO: extend to AbstractDict ? (needs to work-around the inderence bug)
        default_sampling(::Type{$D{K,V}}) where {K,V} = Uniform(Pair{K,V})
        default_sampling(D::Type{<:$D})               = throw(ArgumentError("under-specified scalar type for $D"))

        maketype(::Type{$D{K}},           X, ::Integer) where {K} = $D{K,fieldtype(val_gentype(X), 2)}
        maketype(::Type{$D{K,V} where K}, X, ::Integer) where {V} = $D{fieldtype(val_gentype(X), 1),V}
        maketype(::Type{$D},              X, ::Integer)           = $D{fieldtype(val_gentype(X), 1),fieldtype(val_gentype(X), 2)}
    end
end

rand(rng::AbstractRNG, sp::SamplerTag{Cont{S}}) where {S<:Base.ImmutableDict} =
    foldl((d, _) -> Base.ImmutableDict(d, rand(rng, sp.data.sp)),
          1:sp.data.len,
          init=S())

### AbstractArray

default_sampling(::Type{<:AbstractArray{T}}) where {T} = Uniform(T)
default_sampling(::Type{<:AbstractArray})              = Uniform(Float64)

make(::Type{A}, X,         d1::Integer, dims::Integer...) where {A<:AbstractArray} =
    make(A, X, Dims((d1, dims...)))

make(::Type{A}, ::Type{X}, d1::Integer, dims::Integer...) where {A<:AbstractArray,X} =
    make(A, X, Dims((d1, dims...)))

make(::Type{A}, dims::Dims)                    where {A<:AbstractArray} =
    make(A, default_sampling(A), dims)

make(::Type{A}, d1::Integer, dims::Integer...) where {A<:AbstractArray} =
    make(A, default_sampling(A), Dims((d1, dims...)))

if VERSION < v"1.1.0"
     # to resolve ambiguity
    make(A::Type{<:AbstractArray}, X, d1::Integer)              = make(A, X, Dims((d1,)))
    make(A::Type{<:AbstractArray}, X, d1::Integer, d2::Integer) = make(A, X, Dims((d1, d2)))
end

Sampler(::Type{RNG}, c::Make2{A}, n::Repetition) where {RNG<:AbstractRNG,A<:AbstractArray} =
    SamplerTag{A}((sampler(RNG, c[1], n), c[2]))

rand(rng::AbstractRNG, sp::SamplerTag{A}) where {A<:AbstractArray} =
    rand!(rng, A(undef, sp.data[2]), sp.data[1])

#### Array

# cf. inference bug https://github.com/JuliaLang/julia/issues/28762
# we have to write out all combinations for getting proper inference
maketype(::Type{Array{T}},           _, ::Dims{N}) where {T, N} = Array{T, N}
maketype(::Type{Array{T,N}},         _, ::Dims{N}) where {T, N} = Array{T, N}
maketype(::Type{Array{T,N} where T}, X, ::Dims{N}) where {N}    = Array{val_gentype(X), N}
maketype(::Type{Array},              X, ::Dims{N}) where {N}    = Array{val_gentype(X), N}

#### BitArray

default_sampling(::Type{<:BitArray}) = Uniform(Bool)

maketype(::Type{BitArray{N}}, _, ::Dims{N}) where {N} = BitArray{N}
maketype(::Type{BitArray},    _, ::Dims{N}) where {N} = BitArray{N}

#### sparse vectors & matrices

maketype(::Type{SparseVector},    X, p::AbstractFloat, dims::Dims{1}) = SparseVector{   val_gentype(X), Int}
maketype(::Type{SparseMatrixCSC}, X, p::AbstractFloat, dims::Dims{2}) = SparseMatrixCSC{val_gentype(X), Int}

maketype(::Type{SparseVector{X}},    _, p::AbstractFloat, dims::Dims{1}) where {X} = SparseVector{   X, Int}
maketype(::Type{SparseMatrixCSC{X}}, _, p::AbstractFloat, dims::Dims{2}) where {X} = SparseMatrixCSC{X, Int}

# need to be explicit and split these defs in 2 (or 4) to avoid ambiguities
# TODO: check that using T instead of SparseVector in the RHS doesn't have perfs issues
make(T::Type{SparseVector},    X,         p::AbstractFloat, d1::Integer)                        = make(T, X, p, Dims((d1,)))
make(T::Type{SparseVector},    ::Type{X}, p::AbstractFloat, d1::Integer)              where {X} = make(T, X, p, Dims((d1,)))
make(T::Type{SparseMatrixCSC}, X,         p::AbstractFloat, d1::Integer, d2::Integer)           = make(T, X, p, Dims((d1, d2)))
make(T::Type{SparseMatrixCSC}, ::Type{X}, p::AbstractFloat, d1::Integer, d2::Integer) where {X} = make(T, X, p, Dims((d1, d2)))

make(T::Type{SparseVector},    p::AbstractFloat, d1::Integer)              = make(T, default_sampling(T), p, Dims((d1,)))
make(T::Type{SparseMatrixCSC}, p::AbstractFloat, d1::Integer, d2::Integer) = make(T, default_sampling(T), p, Dims((d1, d2)))

make(T::Type{SparseVector},    p::AbstractFloat, dims::Dims{1}) = make(T, default_sampling(T), p, dims)
make(T::Type{SparseMatrixCSC}, p::AbstractFloat, dims::Dims{2}) = make(T, default_sampling(T), p, dims)

Sampler(::Type{RNG}, c::Make3{A}, n::Repetition) where {RNG<:AbstractRNG,A<:AbstractSparseArray} =
    SamplerTag{Cont{A}}((sp = sampler(RNG, c[1], n),
                         p = c[2],
                         dims = c[3]))

rand(rng::AbstractRNG, sp::SamplerTag{Cont{A}}) where {A<:SparseVector} =
    sprand(rng, sp.data.dims[1], sp.data.p, (r, n)->rand(r, sp.data.sp, n))

rand(rng::AbstractRNG, sp::SamplerTag{Cont{A}}) where {A<:SparseMatrixCSC} =
    sprand(rng, sp.data.dims[1], sp.data.dims[2], sp.data.p, (r, n)->rand(r, sp.data.sp, n), gentype(sp.data.sp))

#### StaticArrays

function random_staticarrays()
    @eval using StaticArrays: tuple_length, tuple_prod, SArray, MArray
    for Arr = (:SArray, :MArray)
        @eval begin
            maketype(::Type{<:$Arr{S}}  , X) where {S<:Tuple}   = $Arr{S,val_gentype(X),tuple_length(S),tuple_prod(S)}
            maketype(::Type{<:$Arr{S,T}}, _) where {S<:Tuple,T} = $Arr{S,T,tuple_length(S),tuple_prod(S)}

            Sampler(::Type{RNG}, c::Make1{A}, n::Repetition) where {RNG<:AbstractRNG,A<:$Arr} =
                SamplerTag{Cont{A}}(Sampler(RNG, c[1], n))

            rand(rng::AbstractRNG, sp::SamplerTag{Cont{$Arr{S,T,N,L}}}) where {S,T,N,L} =
                $Arr{S,T,N,L}(rand(rng, make(NTuple{L}, sp.data)))

            @make_container(T::Type{<:$Arr})
        end
    end
end

### String as a scalar

let b = UInt8['0':'9';'A':'Z';'a':'z'],
    s = Sampler(MersenneTwister, b, Val(Inf)) # cache for the likely most common case

    global Sampler, rand, make

    make(::Type{String})                                   = Make2{String}(8, b)
    make(::Type{String}, chars)                            = Make2{String}(8, chars)
    make(::Type{String}, ::Type{C}) where C                = Make2{String}(8, C)
    make(::Type{String}, n::Integer)                       = Make2{String}(Int(n), b)
    make(::Type{String}, chars,      n::Integer)           = Make2{String}(Int(n), chars)
    make(::Type{String}, ::Type{C},  n::Integer) where {C} = Make2{String}(Int(n), C)
    make(::Type{String}, n::Integer, chars)                = Make2{String}(Int(n), chars)
    make(::Type{String}, n::Integer, ::Type{C}) where {C}  = Make2{String}(Int(n), C)

    Sampler(::Type{RNG}, ::Type{String}, n::Repetition) where {RNG<:AbstractRNG} =
        SamplerTag{Cont{String}}((RNG === MersenneTwister ? s : Sampler(RNG, b, n)) => 8)

    function Sampler(::Type{RNG}, c::Make2{String}, n::Repetition) where {RNG<:AbstractRNG}
        sp = RNG === MersenneTwister && c[2] === b ?
            s : sampler(RNG, c[2], n)
        SamplerTag{Cont{String}}(sp => c[1])
    end

    rand(rng::AbstractRNG, sp::SamplerTag{Cont{String}}) = String(rand(rng, sp.data.first, sp.data.second))
end

## X => a / X => (a, as...) syntax as an alternative to make(X, a) / make(X, a, as...)

# this is experimental

pair_to_make((a, b)::Pair) =
    b isa Tuple ?
        make(a, map(pair_to_make, b)...) :
        make(a, pair_to_make(b))

pair_to_make(x) = x

@inline Sampler(::Type{RNG}, p::Pair, r::Repetition) where {RNG<:AbstractRNG} =
    Sampler(RNG, pair_to_make(p), r)

# nothing can be inferred when only the pair type is available
@inline gentype(::Type{<:Pair}) = Any

@inline gentype(p::Pair) = gentype(pair_to_make(p))