text:
@time @testset "kernels" begin
    @testset begin
        A = [1 2; 3 4]
        B = [5 6; 7 8]
        C1 = zeros(Int, 2, 2)
        C2 = zeros(Int, 2, 2)
        Gaius.add_gemm_kernel!(C1, A, B, Val(1))
        LinearAlgebra.mul!(C2, A, B)
        @test C1 == C2
    end

    @testset begin
        A = [1 2; 3 4]
        v = [5; 6]
        u1 = zeros(Int, 2)
        u2 = zeros(Int, 2)
        Gaius.add_gemm_kernel!(u1, A, v)
        LinearAlgebra.mul!(u2, A, v)
        @test u1 == u2
    end

    @testset begin
        v = adjoint([5; 6])
        A = [1 2; 3 4]
        u1 = adjoint(zeros(Int, 2))
        u2 = adjoint(zeros(Int, 2))
        Gaius.add_gemm_kernel!(u1, v, A)
        LinearAlgebra.mul!(u2, v, A)
        @test u1 == u2
    end

    @testset begin
        A = [1 2; 3 4]
        B = [5 6; 7 8]
        C = zeros(Int, 2, 2)
        Gaius.add_gemm_kernel!(C, A, B, Val(1))
        @test C == A * B
    end

    @testset begin
        A = [1 2; 3 4]
        B = [5 6; 7 8]
        C = zeros(Int, 2, 2)
        Gaius.add_gemm_kernel!(C, A, B, Val(3))
        @test C == 3 * A * B
    end

end


text:
num_part(::UncertainBound{T, N}) where {T, N} = N
ub_eltype(::UncertainBound{T, N}) where {T, N} = T

@testset "UncertainBound (default constructor)" begin
    ub = UncertainBound(1..2, 3..4)
    @test num_part(ub) == NP
    @test ub.left.particles |> length == NP
end

@testset "get_draw()" begin
    ps = Particles([1, 2, 3])
    ub = UncertainBound(ps, ps)
    for i ∈ 1:3
        @test nsi.get_draw(i, ub) == [i, i]
    end
end

@testset "UncertainBound (construct from distributions)" begin
    dleft = Normal(0, 1)
    dright = Uniform(2, 3)
    seed!(1)
    ub = UncertainBound(dleft, dright)
    @test num_part(ub) == 10_000 # default no. of samples from constructor
    @test ub_eltype(ub) == Float64  # Distributions.jl should yield Float64 ny default
    @test ub.left |> std ≈ 1 atol = 1e-4
    @test ub.right |> std ≈ 1/√(12) atol = 1e-4
    ub = UncertainBound(dleft, dright, 1000)
    @test num_part(ub) == 1000
end

@testset "UncertainBound (construct from position(s) and width)" begin

    seed!(1)

    # build an uncertain curve with known position of peaks
    # in each draw
    uc = begin
        x = collect(Float64, 0:0.1:20)
        y = zeros(Float64, 9, length(x))
        # make peaks fall in x-range 10:10.9 and 15:15.9
        # draw 1: peak at 10.0 and 15.0
        # draw 2: peak at 10.1 and 15.1
        # draw 3: peak at 10.2 and 15.2
        # etc.
        for i ∈ 1:9
            y[i, i + 100] = 1.
            y[i, i + 150] = 1.
        end
        UncertainCurve(x, Particles(y))
    end

    dwidth = Uniform(1, 1 + 1e-8) # I want essentially the same width everywhere for testing
    ub = UncertainBound(10.5, dwidth, uc)
    # with the above setup, the values of the `left` and `right` particles should be as follows
    # the deviation comes from the fact that the width distribution (Uniform) hast to have some small width
    @test ub.left.particles[1] ≈ 9.5
    @test ub.left.particles[5] ≈ 9.9
    @test ub.left.particles[9] ≈ 10.3
    @test ub.right.particles[1] ≈ 10.5
    @test ub.right.particles[5] ≈ 10.9
    @test ub.right.particles[9] ≈ 11.3

    # test creation of bounds with correlated width
    dwidth = Uniform(1, 2) # now the width should vary more, to test the correlation of width between the two generated bounds
    ub1, ub2 = UncertainBound([10.5, 15.5], dwidth, uc)
    for i ∈ 1:9
        @test ub1.left.particles[i] - ub1.right.particles[i] ≈ ub2.left.particles[i] - ub2.right.particles[i] atol = 1e-9
    end
end

@testset "UncertainBound creation throws error when distribution has invalid support region" begin
    seed!(1)
    uc = begin # create uncertain curve with one symmetric peak
        x = 0:0.1:10;
        y = @. exp(-(x - 5)^2)
        add_noise(Curve(x, y), GaussianNoiseModel(0.03))
    end;
    @test_throws ArgumentError UncertainBound(5.0, Normal(0, 1), uc)
end

@testset "scale_shift_beta()" begin
    seed!(1)
    β = scale_shift_beta(2, 2, 1, 2)
    xs = rand(β, 1_000_000)
    @test minimum(xs) ≈ 1 atol = 1e-3
    @test maximum(xs) ≈ 2 atol = 1e-3
    @test mean(xs) ≈ 1.5 atol = 1e-3
end

ext: jl


text:
module FH

export finite_horizon, backward_recursion

function backward_recursion(s,n, V, S,A,N,T,R)
    V_a = zeros(Float64, length(A))
    for a in eachindex(A)
        for s_next in eachindex(S)
            if (s, a, s_next, n) in keys(T)
                r = R[s, a, s_next, n]
                p = T[s, a, s_next, n]

                # Backward induction formula
                V_a[a] += p*(r+V[s_next, n+1])
            end
        end
    end
    return maximum(V_a), argmax(V_a)
end

function finite_horizon(S,A,N,T,R; goal_reward = 0.)
    V = Matrix{Float64}(undef, length(S), length(N))
    V[:,N[end]].=goal_reward # Set rewords for reaching goal states
    π = Matrix{Int64}(undef, length(S), length(N))

    # println("Start")
    for n in reverse(N[1:end-1])
        # println(n)
        for s in S
            V[s, n], π[s,n] = backward_recursion(s,n,V, S,A,N,T,R)
        end
    end
    return V, π
end

end


text:
struct Instruction
  register::String
  operation::Function
  condition_register::String
  condition::Function
end

conditions = Dict{String, Function}(
  "==" => ==, ">" => >, "<" => <,
  ">=" => >=, "<=" => <=, "!=" => !=
)

function instruction_from_source(source::String)
  # Split at spaces
  words = split(source, ' ')

  # Instruction parameters
  register = words[1]
  value = parse(Int, words[3])
  words[2] == "inc" && (operation = (x) -> x + value)
  words[2] == "dec" && (operation = (x) -> x - value)

  # Condition
  condition_register = words[5]
  condition_value = parse(Int, words[7])
  condition = (x) -> conditions[words[6]](x, condition_value)

  return Instruction(register, operation, condition_register, condition)
end

function execute(registers, instructions)
  max_value = 0

  for instr in instructions
    # Check if condition is met
    if instr.condition(get(registers, instr.condition_register, 0))
      # Calculate new value
      new_value = instr.operation(get(registers, instr.register, 0))

      # Set to register and compare with old max value
      registers[instr.register] = new_value
      max_value = max(max_value, new_value)
    end
  end

  return max_value
end

# Create an array of all instructions
instructions = map(line -> instruction_from_source(line), eachline(STDIN))

# Execute input
registers = Dict{String, Int}()
max_value = execute(registers, instructions)
println("Max value ever used: $max_value")

# Find max
max_register = reduce((x, y) -> registers[x] >= registers[y] ? x : y, keys(registers))
println("Max value at register $max_register: $(registers[max_register])")


text:
# L2 norm (times a constant)

"""
  NormL2(λ::Real=1.0)

Returns the function `g(x) = λ||x||_2`, for a real parameter `λ ⩾ 0`.
"""

immutable NormL2{R <: Real} <: ProximableConvex
  lambda::R
  function NormL2(lambda::R)
    if lambda < 0
      error("parameter λ must be nonnegative")
    else
      new(lambda)
    end
  end
end

NormL2{R <: Real}(lambda::R=1.0) = NormL2{R}(lambda)

function (f::NormL2)(x::AbstractArray)
  return f.lambda*vecnorm(x)
end

function prox!{T <: RealOrComplex}(f::NormL2, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  vecnormx = vecnorm(x)
  scale = max(0, 1-f.lambda*gamma/vecnormx)
  for i in eachindex(x)
    y[i] = scale*x[i]
  end
  return f.lambda*scale*vecnormx
end

fun_name(f::NormL2) = "Euclidean norm"
fun_dom(f::NormL2) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::NormL2) = "x ↦ λ||x||_2"
fun_params(f::NormL2) = "λ = $(f.lambda)"

function prox_naive{T <: RealOrComplex}(f::NormL2, x::AbstractArray{T}, gamma::Real=1.0)
  vecnormx = vecnorm(x)
  scale = max(0, 1-f.lambda*gamma/vecnormx)
  y = scale*x
  return y, f.lambda*scale*vecnormx
end


text:
abstract type LTISystem <: AbstractSystem end
+(sys1::LTISystem, sys2::LTISystem) = +(promote(sys1, sys2)...)
-(sys1::LTISystem, sys2::LTISystem) = -(promote(sys1, sys2)...)
*(sys1::LTISystem, sys2::LTISystem) = *(promote(sys1, sys2)...)
/(sys1::LTISystem, sys2::LTISystem) = /(promote(sys1, sys2)...)

"""`issiso(sys)`

Returns `true` if `sys` is SISO, else returns `false`."""
function issiso(sys::LTISystem)
    return ninputs(sys) == 1 && noutputs(sys) == 1
end

"""`iscontinuous(sys)`

Returns `true` if `sys` is continuous, else returns `false`."""
function iscontinuous(sys::LTISystem)
    return sys.Ts == 0
end

"""`isstable(sys)`

Returns `true` if `sys` is stable, else returns `false`."""
function isstable(sys::LTISystem)
    if iscontinuous(sys)
        if any(real.(pole(sys)).>=0)
            return false
        end
    else
        if any(abs.(pole(sys)).>=1)
            return false
        end
    end
    return true
end

function _check_consistent_sampling_time(systems::AbstractVector{LTISystem})
    if !all(s.Ts == Ts for s in systems)
       error("Sampling time mismatch")
   end
end
function _check_consistent_sampling_time(sys1::LTISystem, sys2::LTISystem)
    if sys1.Ts != sys2.Ts
       error("Sampling time mismatch")
   end
end


text:
# ---
# title: 1073. Adding Two Negabinary Numbers
# id: problem1073
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/adding-two-negabinary-numbers/description/>
# hidden: true
# ---
#
# Given two numbers `arr1` and `arr2` in base **-2** , return the result of
# adding them together.
#
# Each number is given in _array format_ :  as an array of 0s and 1s, from most
# significant bit to least significant bit.  For example, `arr = [1,1,0,1]`
# represents the number `(-2)^3 + (-2)^2 + (-2)^0 = -3`.  A number `arr` in
# _array, format_ is also guaranteed to have no leading zeros: either `arr ==
# [0]` or `arr[0] == 1`.
#
# Return the result of adding `arr1` and `arr2` in the same format: as an array
# of 0s and 1s with no leading zeros.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr1 = [1,1,1,1,1], arr2 = [1,0,1]
#     Output: [1,0,0,0,0]
#     Explanation: arr1 represents 11, arr2 represents 5, the output represents 16.
#
#
# **Example 2:**
#
#
#
#     Input: arr1 = [0], arr2 = [0]
#     Output: [0]
#
#
# **Example 3:**
#
#
#
#     Input: arr1 = [0], arr2 = [1]
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr1.length, arr2.length <= 1000`
#   * `arr1[i]` and `arr2[i]` are `0` or `1`
#   * `arr1` and `arr2` have no leading zeros
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import Base: <, iszero, ==, <=, >=, sign, signbit, isinteger, isfinite
export <, iszero, ==, <=, >=, sign, signbit, isinteger, isfinite

function isinteger(a::NumberInterval)
    if floor(a) ≺ a ≺ ceil(a)
        return false
    elseif issingleton(a)
        #NOTE given the previous result, singletons must contain single integers
        return true
    end
    return missing_or_exception(a)
end

function isfinite(a::NumberInterval)
    # per IEEE standard, intervals cannot contain infinities
    if !isempty(a) && !isnan(a)
        return true
    end
    return missing_or_exception(a)
end

for (numberf, setf) in ((:<, :strictprecedes), (:<=, :precedes))
    @eval function $numberf(a::NumberInterval, b::NumberInterval)
        if $setf(a, b)
            return true
        elseif $setf(b, a)
            return false
        end
        return missing_or_exception((a, b))
    end
end

>=(a::NumberInterval, b::NumberInterval) = b <= a

function iszero(a::NumberInterval)
    if !contains_zero(a)
        return false
    end
    if a ⊆ zero(typeof(a))
        return true
    end
    return missing_or_exception(a)
end

function ==(a::NumberInterval, b::NumberInterval)
    if isnan(a) || isnan(b)
        return false
    elseif isdisjoint(a, b)
        return false
    elseif issingleton(a) && issingleton(b) && a ⊆ b
        return true
    end
    return missing_or_exception((a, b))
end

function sign(a::NumberInterval)
    z = zero(typeof(a))
    if a ≺ z
        return -1
    elseif z ≺ a
        return +1
    elseif a ⊆ z
        return  0
    end
    return missing_or_exception(a)
end

function signbit(a::NumberInterval)
    z = zero(typeof(a))
    if a ≺ z
        return true
    elseif precedes(z, a)
        return false
    end
    return missing_or_exception(a)
end


text:
module APIUtils

using ..CUDA

using Libdl

# helpers that facilitate working with CUDA APIs
include("call.jl")
include("enum.jl")

end


text:
# src/core.jl
# Core UnitfulCurrency Module

@reexport module _unitfulcurrency

using Unitful

export CURRENCY, Currency, CurrencyFreeUnits, CurrencyUnits, @currencyunit

@dimension CURRENCY "CURRENCY" Currency

"""```
@currencyunit code name
    === Unitful.@refunit code String(code) name CURRENCY false
```
Build a `Currency` type using the `Unitful.@refunit` with given `code` and `name`.
"""
macro currencyunit(code, name)
    abbreviation = String(code)
    esc(quote
        Unitful.@refunit $code $abbreviation $name CURRENCY false
    end)
    # FIXME: allow base10 prefixes
end

end  # module _unitfulcurrency


text:
@plottest begin
    # SRM0 params
    η₀ = 5.0
    τᵣ = 1.0

    # Input spike train params
    rate = 0.01
    T = 20
    ∂t = 0.01
    n = convert(Int, ceil(T / ∂t))

    srm = Neuron(QueuedSynapse(Synapse.Delta()), SRM0(η₀, τᵣ), Threshold.Poisson(5.0, 0.5, 0.1))
    input = ConstantRate(rate)
    spikes = excite!(srm, input, n)

    # callback to record voltages
    voltages = Float64[]
    record = function ()
        push!(voltages, getvoltage(srm))
    end

    # simulate
    output = simulate!(srm, n; dt = ∂t, cb = record, dense = true)

    # plot raster plot
    raster_plot = rasterplot(∂t .* spikes, ∂t .* output, label = ["Input", "Output"], xlabel = "Time (sec)",
                    title = "Raster Plot (δ response)")
    xlims!(0, T)

    # plot dense voltage recording
    plot(∂t .* collect(1:n), voltages,
        title = "SRM Membrane Potential with Varying Presynaptic Responses", xlabel = "Time (sec)", ylabel = "Potential (V)", label = "δ response")

    # resimulate using presynaptic response
    srm = Neuron(QueuedSynapse(Synapse.Alpha()), SRM0(η₀, τᵣ), Threshold.Poisson(5.0, 0.5, 0.1))
    voltages = Float64[]
    excite!(srm, spikes)
    simulate!(srm, n; dt = ∂t, cb = record, dense = true)

    # plot voltages with response function
    voltage_plot = plot!(∂t .* collect(1:n), voltages, label = "α response")
    xlims!(0, T)

    plot(raster_plot, voltage_plot, layout = grid(2, 1), xticks = 0:T)
end joinpath(datadir, "Srm0-Threshold-Test.png") !isci


text:
procs = addprocs(12)

using OPCSPs
using MCTS

N = 1000

problems = [gen_problem(noise=10.0,
                        p=0.2,
                        n_nodes=10,
                        rng=MersenneTwister(i+1000))
            for i in 1:N]

iss = initial_states(problems, rng_offset=100)

naive = evaluate_performance(problems, iss, GurobiExactSolver(multithreaded=false), rng_offset=1000)
@show mean(naive)

cheating = evaluate_performance(problems, iss, Cheater(), rng_offset=1000)
@show mean(cheating)

mean_feedback = evaluate_performance(problems, iss, FeedbackSolver(GurobiExactSolver(multithreaded=false)), rng_offset=1000)
@show mean(mean_feedback)

heur_feedback = evaluate_performance(problems, iss, FeedbackSolver(HeuristicSolver()), rng_offset=1000)
@show mean(heur_feedback)

s = GurobiExactSolver(multithreaded=false)

srng = MersenneTwister(1947)
solver = AgUCTSolver(
    aggregator=VoxelOPCSPAg(2.0),
    rollout_solver=FeedbackSolver(s),
    exploration_constant=30.0,
    n_iterations=50000,
    rng=srng
)

@time ag = evaluate_performance(problems, iss, solver, rng_offset=1000, parallel=true)
@show mean(ag)

rmprocs(procs)


text:
# Custom plugin for Documentation workflow
struct DocumenterCI <: PkgTemplates.FilePlugin
    documentation_yml::String

    # Default constructor; don't use this
    DocumenterCI(documentation_yml::AbstractString) = new(documentation_yml)
end

# Construtor with kwargs; use this. Don't use the default value
DocumenterCI(; documentation_yml::AbstractString = "") = DocumenterCI(documentation_yml)

# Specifying source template and destination
PkgTemplates.source(p::DocumenterCI) = p.documentation_yml
PkgTemplates.destination(::DocumenterCI) = "./.github/workflows/documentation.yml"

# view definition
PkgTemplates.view(::DocumenterCI, ::Template, ::AbstractString) = Dict(
    "secrets.GITHUB_TOKEN" => "{{ secrets.GITHUB_TOKEN }}",
    "secrets.DOCUMENTER_KEY" => "{{ secrets.DOCUMENTER_KEY }}"
)


text:
include("lattices.jl")

mutable struct ElectronTinder
    """ Database to match lonely up spins with lonely down spins by site"""
    lonely_up::Array{Int}
    lonely_down::Array{Int}
end

function addUser!(t::ElectronTinder,site::Int,spin::Int)
    spin==1 ? append!(t.lonely_up,site) : append!(t.lonely_down,site)
end

function SwipeRight(database::ElectronTinder)
    """ Creates list of matches (moves) generated by
    electron tinder's database """
    @assert length(database.lonely_up) == length(database.lonely_down)
    lu = database.lonely_up
    ld = database.lonely_down
    nmatches = length(lu)
    matches = Array{ExchangeMove}(undef,nmatches)
    for i in 1:nmatches
        matches[i] = ExchangeMove((lu[i],ld[i]))
    end
    return matches
end

using Statistics: mean
using LightGraphs: neighborhood
function Swap(chain::DoubleGutzwillerChain,regionsize::Int)
    """ regionsize ≈ r, the radius of the region selected """

    states = get_State(chain) # (state1, state2)
    lg = get_Model(chain)["lattice"].graph
    ns = prod(get_Model(chain)["dims"])

    type = typeof(states[1].det_A_up)
    swap_per_site = zeros(type,ns)

    # redundancy for mixed BC lattice
    lg0 = SimpleGraph(adjacency_matrix(lg))

    for v = 1:ns
        sites = unique(neighborhood(lg0,v,regionsize))
        swap_per_site[v] = compute_swapregion(chain,sites)
    end

    return mean(swap_per_site)
end

function compute_swapregion(chain::DoubleGutzwillerChain,
                            sites::Array)
    """ This function is somewhat nontrivial, the following is the
    layout of how it works. First, a list of exchange moves is created
    for each replica that maps the current configuration to the new
    one after swapping. This happens in a two stages: first the spin
    spin configurations inside the regions are compared to see where they
    disagree. Every site where they disagree, an entry is added the
    appropriate spin list in Electron Tinder. Finally, the completed
    ElectronTinder database is fed into SwipeRight to generate a list of moves.

    After this, the list of moves are submitted to preexisting functions
    to compute the determinant ratio produced by that sequence of moves.

    Finally what is returned is the expectation value of the swap operator
    acting on this region which is,
    det(swapped copy 1)*det(swapped copy 2) / (det(og copy 1) * det(og copy 2))

    Later, I may modify this to be more general/support more than 2 copies,
    but for now it seems like thats more trouble than its worth.
    """
    s0 = get_State(chain)
    (spinc1, spinc2) = (s0[1].spin_config.sc,
                        s0[2].spin_config.sc)

    if sum(spinc1[sites])!=sum(spinc2[sites]) # regions must have same quantum #
        return 0
    end
    if spinc1[sites] == spinc2[sites]
        return 1
    end

    Tinder = (ElectronTinder([],[]),ElectronTinder([],[]))
    for site in sites
        if spinc1[site] != spinc2[site]
            addUser!(Tinder[1],site,spinc1[site])
            addUser!(Tinder[2],site,spinc2[site])
        end
    end

    (MoveList1, MoveList2) = SwipeRight.(Tinder)
    @assert size(MoveList1) == size(MoveList2)

    # now because update_chain! modifies the state, make copies of
    # the two states to update
    free_chains = deepcopy(get_Replicas(chain))
    free_states = get_State.(free_chains)

    @assert length(MoveList1) == length(MoveList2) "compute swapregion: move lists
    for two copies are of different length"

    # try to do "fast" determinant update with all moves simultaneously
    n_moves = length(MoveList1)
    num_particles = length(free_states[1].r_up) # not very pretty

    type = typeof(s0[1].det_A_up) # complexf64 or f64

    U_up_copy1 = zeros(type,(num_particles,n_moves))
    U_down_copy1 = zeros(type,(num_particles,n_moves))
    V_up_copy1 = zeros(type,(num_particles,n_moves))
    V_down_copy1 = zeros(type,(num_particles,n_moves))

    U_up_copy2 = zeros(type,(num_particles,n_moves))
    U_down_copy2 = zeros(type,(num_particles,n_moves))
    V_up_copy2 = zeros(type,(num_particles,n_moves))
    V_down_copy2 = zeros(type,(num_particles,n_moves))

    for i in 1:n_moves
        # update_Rs!(free_states[1],MoveList1[i])
        # update_Rs!(free_states[2],MoveList2[i])

        U_up_copy1[:,i], V_up_copy1[:,i], U_down_copy1[:,i], V_down_copy1[:,i] =
        get_update_vectors(free_states[1],MoveList1[i])

        U_up_copy2[:,i], V_up_copy2[:,i], U_down_copy2[:,i], V_down_copy2[:,i] =
        get_update_vectors(free_states[2],MoveList2[i])

    end

    # its important to do the determinant ratio updates all at once
    # brute force reevaluation of determinants code left in for debugging

    #
    # filled_states = get_Wavefunctions(free_chains[1])
    #
    # det1up = get_Determinant(free_states[1].r_up,filled_states)
    # det1down = get_Determinant(free_states[1].r_down,filled_states)
    # det2up = get_Determinant(free_states[2].r_up,filled_states)
    # det2down = get_Determinant(free_states[2].r_down,filled_states)

    det_ratio_1u = det_ratio_factor(free_states[1].A_inv_up,U_up_copy1,V_up_copy1)
    det_ratio_1d = det_ratio_factor(free_states[1].A_inv_down,U_down_copy1,V_down_copy1)
    det_ratio_2u = det_ratio_factor(free_states[2].A_inv_up,U_up_copy2,V_up_copy2)
    det_ratio_2d = det_ratio_factor(free_states[2].A_inv_down,U_down_copy2,V_down_copy2)

    fast_ratio = det_ratio_1u * det_ratio_1d * det_ratio_2u * det_ratio_2d

    # ss1 = get_test_state(free_chains[1],free_states[1].r_up, free_states[1].r_down)
    # ss2 = get_test_state(free_chains[2],free_states[2].r_up, free_states[2].r_down)

    # sdp = det1up * det1down * det2up * det2down
    # odp = s0[1].det_A_up * s0[1].det_A_down * s0[2].det_A_up * s0[2].det_A_down

    # print("sdp: ",sdp," odp: ",odp,"\n")

    return fast_ratio
end


text:
import PyCall
# check if python has the requird Packages installed

required_python_packages = ["numpy","pandas","xlsxwriter"] # for testing this functionality you can add some package like "nose" "SymPy" "nltk" "Scapy" ....
try
    for x in required_python_packages
    	PyCall.pywrap(PyCall.pyimport(x))
    end
catch
    	@info "One of the required python packages was not found on your system:"
    	println(required_python_packages)
    for x in required_python_packages
        println(x)
    end
    	@info "Trying to install those packages..."
    		# Use eventual proxy info
    		proxy_arg = Array{AbstractString}(undef, 0)
    		if haskey(ENV, "http_proxy")
    			push!(proxy_arg, "--proxy")
    			push!(proxy_arg, ENV["http_proxy"])
    		end

    		# Import pip
    		try
    			@pyimport pip
    		catch
    			# If it is not found, install it
    			println("Pip not found on your sytstem. Downloading it.")
    			get_pip = joinpath(dirname(@__FILE__), "get-pip.py")
    			download("https://bootstrap.pypa.io/get-pip.py", get_pip)
    			run(`$(PyCall.python) $(proxy_arg) $get_pip --user`)
    		end

    		println("Installing required python packages using pip")
    		run(`$(PyCall.python) $(proxy_arg) -m pip install --user --upgrade pip setuptools`)
    		run(`$(PyCall.python) $(proxy_arg) -m pip install --user --upgrade pip`)
    		run(`$(PyCall.python) $(proxy_arg) -m pip install --user $(required_python_packages)`)

    	@info "installation of python packages finished."
end


text:
using BinaryBuilder

name = "kubectl"
version = v"1.20.0"

sources = [
    FileSource("https://storage.googleapis.com/kubernetes-release/release/v1.20.0/bin/linux/amd64/kubectl", "a5895007f331f08d2e082eb12458764949559f30bcc5beae26c38f3e2724262c"; filename="x86_64-linux-gnu-kubectl"),
    FileSource("https://storage.googleapis.com/kubernetes-release/release/v1.20.0/bin/darwin/amd64/kubectl", "82046a4abb056005edec097a42cc3bb55d1edd562d6f6f38c07318603fcd9fca"; filename="x86_64-apple-darwin14-kubectl"),
    FileSource("https://storage.googleapis.com/kubernetes-release/release/v1.20.0/bin/windows/amd64/kubectl.exe", "ee7be8e93349fb0fd1db7f5cdb5985f5698cef69b7b7be012fc0e6bed06b254d"; filename="x86_64-w64-mingw32-kubectl"),
    FileSource("https://www.apache.org/licenses/LICENSE-2.0.txt", "cfc7749b96f63bd31c3c42b5c471bf756814053e847c10f3eb003417bc523d30"; filename="LICENSE.txt")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/
mkdir -p ${bindir}

install_license LICENSE.txt

mv "${target}-kubectl" "${bindir}/kubectl${exeext}"
chmod 755 "${bindir}/kubectl${exeext}"

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("kubectl", :kubectl),
]

dependencies = Dependency[]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the acoustic perturbation equations

equations = AcousticPerturbationEquations2D(v_mean_global=(0.5, 0.0), c_mean_global=1.0,
                                            rho_mean_global=1.0)

# Create DG solver with polynomial degree = 5 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=5, surface_flux=flux_lax_friedrichs)

coordinates_min = (-100.0, 0.0) # minimum coordinates (min(x), min(y))
coordinates_max = (100.0, 200.0) # maximum coordinates (max(x), max(y))

# Create a uniformly refined mesh with periodic boundaries
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4,
                n_cells_max=100_000,
                periodicity=false)

"""
    initial_condition_gauss_wall(x, t, equations::AcousticPerturbationEquations2D)

A Gaussian pulse, used in the `gauss_wall` example elixir in combination with
[`boundary_condition_wall`](@ref). Uses the global mean values from `equations`.
"""
function initial_condition_gauss_wall(x, t, equations::AcousticPerturbationEquations2D)
  v1_prime = 0.0
  v2_prime = 0.0
  p_prime = exp(-log(2) * (x[1]^2 + (x[2] - 25)^2) / 25)

  prim = SVector(v1_prime, v2_prime, p_prime, global_mean_vars(equations)...)

  return prim2cons(prim, equations)
end
initial_condition = initial_condition_gauss_wall

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_condition_wall)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 30.0
tspan = (0.0, 30.0)
ode = semidiscretize(semi, tspan)

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval=100)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval=100, solution_variables=cons2state)

# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl=0.7)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)

###############################################################################
# run the simulation

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks)

# Print the timer summary
summary_callback()


text:
# adopted from https://github.com/JuliaLang/julia/blob/release-0.6/contrib/julia-config.jl

isdebugbuild() = ccall(:jl_is_debugbuild, Cint, ()) != 0

function shell_escape(str)
    str = replace(str, "'" => "'\''")
    return "'$str'"
end

function julia_libdir()
    libname = isdebugbuild() ? "libjulia-debug" :
                               "libjulia"
    return dirname(abspath(Libdl.dlpath(libname)))
end

function julia_private_libdir()
    if Base.DARWIN_FRAMEWORK # taken from Libdl tests
        if ccall(:jl_is_debugbuild, Cint, ()) != 0
            dirname(abspath(Libdl.dlpath(Base.DARWIN_FRAMEWORK_NAME * "_debug")))
        else
            joinpath(dirname(abspath(Libdl.dlpath(Base.DARWIN_FRAMEWORK_NAME))),"Frameworks")
        end
    elseif ccall(:jl_is_debugbuild, Cint, ()) != 0
        dirname(abspath(Libdl.dlpath("libjulia-internal-debug")))
    else
        dirname(abspath(Libdl.dlpath("libjulia-internal")))
    end
end

julia_includedir() = abspath(Sys.BINDIR, Base.INCLUDEDIR, "julia")

function ldflags()
    fl = "-L$(shell_escape(julia_libdir())) -L$(shell_escape(julia_private_libdir()))"
    if Sys.iswindows()
        fl = fl * " -Wl,--stack,8388608"
        fl = fl * " -Wl,--export-all-symbols"
    elseif Sys.islinux()
        fl = fl * " -Wl,--export-dynamic"
    end
    return fl
end

function ldlibs()
    libnames = isdebugbuild() ? "-ljulia-debug -ljulia-internal-debug" :
                                "-ljulia       -ljulia-internal"
    if Sys.islinux()
        return "-Wl,-rpath-link,$(shell_escape(julia_libdir())) -Wl,-rpath-link,$(shell_escape(julia_private_libdir())) $libnames"
    elseif Sys.iswindows()
        return "$libnames -lopenlibm"
    else
        return "$libnames"
    end
end

function cflags()
    flags = IOBuffer()
    print(flags, "-std=gnu99")
    include = shell_escape(julia_includedir())
    print(flags, " -I", include)
    if Sys.isunix()
        print(flags, " -fPIC")
    end
    return String(take!(flags))
end

function rpath_executable()
    Sys.iswindows() ? `` :
    Sys.isapple()   ? `-Wl,-rpath,'@executable_path/../lib' -Wl,-rpath,'@executable_path/../lib/julia'` :
                      `-Wl,-rpath,\$ORIGIN/../lib:\$ORIGIN/../lib/julia`
end

function rpath_sysimage()
    Sys.iswindows() ? `` :
    Sys.isapple()   ? `-Wl,-rpath,'@executable_path' -Wl,-rpath,'@executable_path/julia'` :
                      `-Wl,-rpath,\$ORIGIN:\$ORIGIN/julia`
end


text:
"""
Internal module used for testing.
"""
module _Testing

using Test

struct XFailTestSet <: Test.AbstractTestSet
    ts::Test.AbstractTestSet

    XFailTestSet(ts::Test.AbstractTestSet) = new(ts)
end

XFailTestSet(args...; kwargs...) =
    XFailTestSet(typeof(Test.get_testset())(args...; kwargs...))

Test.finish(ts::XFailTestSet) = Test.finish(ts.ts)
Test.record(ts::XFailTestSet, result::Test.Pass) = Test.record(
    ts.ts,
    Test.Fail(
        Symbol(:xfail_, result.test_type),
        result.orig_expr,
        result.data,
        result.value,
        LineNumberNode(0, Symbol("<unknown>")),
    ),
)
Test.record(ts::XFailTestSet, result::Test.Fail) = Test.record(
    ts.ts,
    Test.Pass(
        Symbol(:xfail_, result.test_type),
        result.orig_expr,
        result.data,
        result.value,
    ),
)
Test.record(ts::XFailTestSet, result::Test.Error) = Test.record(ts.ts, result)

end  # module


text:
#=
Decomposition of Return
    rₜ=√hₜ zₜ, zₜ~N(0,1) so rₜ|Fₜ₋₁~N(0,hₜ). hₜ=Var(rₜ|Fₜ₋₁)

GARCH(p,q)(Bollerslev 1986)

    hₜ₊₁=ω+Σβⱼhₜ₋ⱼ+Σαᵢ(rₜ₋ᵢ²-hₜ₋ᵢ)

    αᵢ,βⱼ,ω≧0

EGARCH(p,q)(Nelson 1991)

    loghₜ₊₁=ω+Σαᵢg(zₜ)+Σβⱼloghₜ₋ⱼ

    where g(z)=θ₁z+θ₂(|z|-E|z|) (Leverage function)

REGARCH(1,1,K)(Hansen and Zhuo 2015)

    GARCH equation

        loghₜ₊₁=ω + τ(zₜ) + βloghₜ + γ′uₜ

        where τ(z)=τ₁z+τ₂(z²-1)

    Measurement Equation

        log xₖ,ₜ = ξₖ + ϕₖ log hₜ + δ₍ₖ₎(zₜ) + uₖ,ₜ for k=1,...,K (we assume ϕₖ=1)

        where δ₍ₖ₎(z)=δₖ,₁z + δₖ,₂(z²-1)

=#

module GARCHModels

    using Optim
    using Optim: TwiceDifferentiable, optimize, minimizer
    using FiniteDiff, ForwardDiff
    using LinearAlgebra, Statistics, Distributions
    using DataFrames

    export REGARCH, GARCH, EGARCH

    mutable struct Volatilityestimation
        parameter
        robustcov
        conditionalvariance
        LL
    end

    subscript(i::Integer) = i<0 ? error("$i is negative") : join('₀'+d for d in reverse(digits(i)))

    function GARCH_coefnames(p,q)
        names = Array{String, 1}()
        push!(names,"ω")
        for k in 1:p
            push!(names, "β"*subscript(k))
        end
        for k in 1:q
            push!(names, "α"*subscript(k))
        end
        push!(names,"h0")
        return names
    end
    function GARCH_startingvals(p,q,rt)
        x0 = zeros(p+q+2)
        x0[1]=0.01
        x0[2:p+1] = 0.7/p
        x0[p+2:p+1+q] .= 0.2/q
        x0[end] = (var(rt))
        return x0
    end
    function GARCH_constraints(p,q)
        lower = -ones(p+q+2)
        upper = ones(p+q+2)
        lower[end] = 0
        return lower, upper
    end
    function GARCH_loglik(params::AbstractVector{T},p,q,rt,output) where T
    # Return data must be a Vector, not a Matrix.
        lower, upper = GARCH_constraints(p,q)
        all(lower.<params.<upper) || return -Inf
        n = size(rt)[1]
        presample = max(p,q)
        names = Array{String, 1}()

        ω = params[1]
        β = params[2:p+1]
        α = params[p+2:p+q+1]
        h0 = params[end]
        ht = zeros(T,n)
        zt = zeros(T,n)
        ht[1:presample] .= h0
        zt[1:presample] = rt[1:presample]./sqrt.(ht[1:presample])
        for t in presample+1:n
            mht = ω
            for i in 1:p
                mht += β[i]*ht[t-i]#log(xt[t,k]) - ξ[k] - log(ht[t]) - δ₁[k]*zt[t] - δ₂[k]*(zt[t]^2 - 1)
            end
            for j in 1:q
                mht += α[j]*(rt[t-j]^2 - ht[t-j])
            end
            ht[t] = mht
            zt[t] = rt[t]./sqrt.(ht[t])
        end
        LL = -log.(ht) .- zt.^2
        if output=="LL"
            return sum(LL)
        elseif output=="ht"
            return ht
        elseif output=="LLs"
            return LL
        end
    end
    function GARCH_scores(params,p,q,rt)
        f = x -> GARCH_loglik(x,p,q,rt,"LLs")
        S = ForwardDiff.jacobian(f, params)
    end
    function GARCH_InformationMatrix(params,p,q,rt)
        g = x -> GARCH_loglik(x,p,q,rt,"LL")
        H = ForwardDiff.hessian(g, params)
        J = -H
    end

    function GARCH_fit(p,q,rt,xt)
        obj = params -> -GARCH_loglik(params,p,q,rt,"LL")
        params = GARCH_startingvals(p,q,rt)
        res = optimize(obj, params, autodiff=:forward,iterations =10000,show_trace=true)
        params .= Optim.minimizer(res)
    end

    function GARCH(p,q,rt)
        params = GARCH_fit(p,q,rt)
        I = GARCH_InformationMatrix(params,p,q,rt)
        S = GARCH_scores(params,p,q,rt)
        ht = GARCH_loglik(params,p,q,rt,"ht")
        J = S'*S
        robustcov = I\J/I
        LL = GARCH_loglik(params,p,q,rt,"LL")
        coefnames = GARCH_coefnames(K)
        z = params./sqrt.(diag(robustcov))
        pvalue = (1 .- cdf.(Normal(),abs.(z)))*2
        parameters = DataFrame([coefnames,params,sqrt.(diag(robustcov)),z,pvalue])
        DataFrames.rename!(parameters,["Parameter","Estimate","Std","z","P(>|z|)"])
        res = Volatilityestimation(parameters,robustcov,ht,LL)
    end

    function EGARCH_coefnames(p,q)
        names = Array{String, 1}()
        push!(names,"ω")
        for k in 1:p
            push!(names, "β"*subscript(k))
        end
        for k in 1:q
            push!(names, "α"*subscript(k))
        end
        for i in 1:2
            push!(θ, "θ"*subscript(k))
        end
        push!(names,"logh0")
        return names
    end

    function EGARCH_startingvals(p,q,rt)
        x0 = zeros(p+q+4)
        x0[1]=0
        x0[2:p+1] = 0.7/p
        x0[p+2:p+1+q] .= 0.2/q
        x0[p+q+2:p+q+3] .= 0
        x0[end] = (var(rt))
        return x0
    end

    function EGARCH_constraints(p,q)
        lower = -ones(p+q+4)
        upper = ones(p+q+4)
        lower[end] = -10
        return lower, upper
    end

    function EGARCH_loglik(params::AbstractVector{T},p,q,rt,output) where T
    # Return data must be a Vector, not a Matrix.
        lower, upper = EGARCH_constraints(p,q)
        all(lower.<params.<upper) || return -Inf

        n = size(rt)[1]
        presample = max(p,q)
        names = Array{String, 1}()

        ω = params[1]
        β = params[2:p+1]
        α = params[p+2:p+q+1]
        θ = params[p+q+2:p+q+3]
        logh0 = params[end]
        ht = zeros(T,n)
        zt = zeros(T,n)
        ht[1:presample] .= exp(logh0)
        zt[1:presample] = rt[1:presample]./sqrt.(ht[1:presample])
        for t in presample+1:n
            mht = ω
            for i in 1:p
                mht += β[i]*log(ht[t-i])#log(xt[t,k]) - ξ[k] - log(ht[t]) - δ₁[k]*zt[t] - δ₂[k]*(zt[t]^2 - 1)
            end
            for j in 1:q
                mht += α[j]*(abs(zt[t-j]) - sqrt(2/pi))
            end
            ht[t] = mht
            zt[t] = rt[t]./sqrt.(ht[t])
        end
        LL = -log.(ht) .- zt.^2
        if output=="LL"
            return sum(LL)
        elseif output=="ht"
            return ht
        elseif output=="LLs"
            return LL
        end
    end
    function EGARCH_scores(params,p,q,rt)
        f = x -> EGARCH_loglik(x,p,q,rt,"LLs")
        S = ForwardDiff.jacobian(f, params)
    end
    function EGARCH_InformationMatrix(params,p,q,rt)
        g = x -> EGARCH_loglik(x,p,q,rt,"LL")
        H = ForwardDiff.hessian(g, params)
        J = -H
    end

    function EGARCH_fit(p,q,rt,xt)
        obj = params -> -EGARCH_loglik(params,p,q,rt,"LL")
        params = EGARCH_startingvals(p,q,rt)
        res = optimize(obj, params, autodiff=:forward,iterations =10000,show_trace=true)
        params .= Optim.minimizer(res)
    end

    function EGARCH(p,q,rt)
        params = EGARCH_fit(p,q,rt)
        I = EGARCH_InformationMatrix(params,p,q,rt)
        S = EGARCH_scores(params,p,q,rt)
        ht = EGARCH_loglik(params,p,q,rt,"ht")
        J = S'*S
        robustcov = I\J/I
        LL = EGARCH_loglik(params,p,q,rt,"LL")
        coefnames = EGARCH_coefnames(K)
        z = params./sqrt.(diag(robustcov))
        pvalue = (1 .- cdf.(Normal(),abs.(z)))*2
        parameters = DataFrame([coefnames,params,sqrt.(diag(robustcov)),z,pvalue])
        DataFrames.rename!(parameters,["Parameter","Estimate","Std","z","P(>|z|)"])
        res = Volatilityestimation(parameters,robustcov,ht,LL)
    end

    function REGARCH_coefnames(p,q)
        names = Array{String, 1}()
        push!(names,"ω")
        push!(names, "β")
        push!(names, "τ₁")
        push!(names, "τ₂")
        for k in 1:K
            push!(names, "γ"*subscript(k))
        end
        for k in 1:K
            push!(names, "ξ"*subscript(k))
        end
        for k in 1:K
            push!(names,"δ₁"*subscript(k))
        end
        for k in 1:K
            push!(names,"δ₂"*subscript(k))
        end
        push!(names,"logh0")
        return names
    end
    function REGARCH_startingvals(K,rt)
        x0 = zeros(5+4*K)
        x0[1]=0.01
        x0[2] = 0.8
        x0[5:4+K] .= 0.3/K
        x0[5+K:4+2*K] .= -0.3
        x0[end] = log(var(rt))
        return x0
    end
    function REGARCH_constraints(K)
        lower = -ones(4*K+5)
        upper = ones(4*K+5)
        upper[5+2*K:4+3*K] .= 2
        lower[end] = -50
        return lower, upper
    end
    function REGARCH_loglik(params::AbstractVector{T},K,rt,xt,output) where T
    # Return data must be a Vector, not a Matrix.
        size(rt)[1]==size(xt)[1] || return error("Check the length of the data")
        lower, upper = REGARCH_constraints(K)
        all(lower.<params.<upper) || return -Inf
        n = size(rt)[1]
        presample = 1
        nparam = 5+5*K
        ω = params[1]
        β = params[2]
        τ₁ = params[3]
        τ₂ = params[4]
        γ = params[5:4+K]
        ξ = params[5+K:4+2*K]
        δ₁ = params[5+2*K:4+3*K]
        δ₂ = params[5+3*K:4+4*K]
        logh0 = params[end]
        ht = zeros(T,n)
        ut = zeros(T,n,K)
        zt = zeros(T,n)
        ht[1:presample] .= exp(logh0)
        ut[1:presample] .= 0
        zt[1:presample] = rt[1:presample]./sqrt.(ht[1:presample])
        for t in presample+1:n
            mlht = ω + β*log(ht[t-1]) .+ τ₁*zt[t-1] .+ τ₂*(zt[t-1]^2-1) .+ sum(γ.*(ut[t-1,:]))
            ht[t] = exp(mlht)
            zt[t] = rt[t]./sqrt.(ht[t])
            for k in 1:K
                ut[t,k] = log(xt[t,k]) - ξ[k] - log(ht[t]) - δ₁[k]*zt[t] - δ₂[k]*(zt[t]^2 - 1)
            end
        end
        Σ = ut'*ut/n
        LL = -log.(ht) .- zt.^2 .- log(det(Σ))
        if output=="LL"
            return sum(LL)
        elseif output=="ht"
            return ht
        elseif output=="LLs"
            return LL
        end
    end
    function REGARCH_scores(params,K,rt,xt)
        f = x -> REGARCH_loglik(x,K,rt,xt,"LLs")
        S = ForwardDiff.jacobian(f, params)
    end
    function REGARCH_InformationMatrix(params,K,rt,xt)
        g = x -> REGARCH_loglik(x,K,rt,xt,"LL")
        H = ForwardDiff.hessian(g, params)
        J = -H
    end

    function REGARCH_fit(K,rt,xt)
        obj = params -> -REGARCH_loglik(params,K,rt,xt,"LL")
        params = REGARCH_startingvals(K,rt)
        res = optimize(obj, params, autodiff=:forward,iterations =10000,show_trace=true)
        params .= Optim.minimizer(res)
    end

    function REGARCH(K,rt,xt)
        params = REGARCH_fit(K,rt,xt)
        I = REGARCH_InformationMatrix(params,K,rt,xt)
        S = REGARCH_scores(params,K,rt,xt)
        ht = REGARCH_loglik(params,K,rt,xt,"ht")
        J = S'*S
        robustcov = I\J/I
        LL = REGARCH_loglik(params,K,rt,xt,"LL")
        coefnames = REGARCH_coefnames(K)
        z = params./sqrt.(diag(robustcov))
        pvalue = (1 .- cdf.(Normal(),abs.(z)))*2
        parameters = DataFrame([coefnames,params,sqrt.(diag(robustcov)),z,pvalue])
        DataFrames.rename!(parameters,["Parameter","Estimate","Std","z","P(>|z|)"])
        res = Volatilityestimation(parameters,robustcov,ht,LL)
    end

end


text:
using PreprocessData
using PreprocessData: rename_cols, place_target, find_in, getpath, getdata, splits
using PreprocessData: Iris, df_to_array, df_or_array, getfilename, getpath
using Test
using DataFrames
using CSV
using DataDeps
using LinearAlgebra
using Statistics

# rename_cols
dfOrig = DataFrame(A=[1,2,3],B=["a","b","c"],C=[5,6,7])
df_renamed = rename(dfOrig, [Symbol("Column 1"), Symbol("Column 2"), :C])
@test rename_cols(dfOrig) == df_renamed

# place_target
dfOrig = DataFrame(A=[1,2,3],B=["a","b","c"],C=[5,6,7])
dfNew = DataFrame(A=[1,2,3],C=[5,6,7],B=["a","b","c"])
rename!(dfNew, [Symbol("Column 1"), Symbol("Column 2"), Symbol("Target")])
@test place_target(2,dfOrig,"") == dfNew

# place_target
dfOrig = DataFrame(A=[1,2,3],B=["a","b","c"],C=[5,6,7])
dfNew = DataFrame(A=[1,2,3],B=["a","b","c"], C=[5,6,7], Target=["first","second","third"])
typeSplit = "train"
labels = DataFrame(NO=["first","second","third"])
file = joinpath(pwd(), string("labels-", typeSplit))
CSV.write(file, labels, header=false) => file
@test place_target(PreprocessData.Labels,dfOrig,typeSplit) == dfNew

# find_in
word1 = "isthereword.valid?"
@test find_in(word1) == "valid"
word2 = "there_is_nothing.really"
@test find_in(word1) == "valid"

# extract
# The Zip archive file.csv.zip must be in the folder in order to test
# pathNotExtracted = joinpath(pwd(), string("file.csv.zip"))
# pathExtracted = joinpath(pwd(), string("file.csv"))
# @test PreprocessData.extract(pathNotExtracted) == pathExtracted

irispath = getpath(Iris())

# getdata
iris = CSV.File(joinpath(irispath, "data-train.csv"), header=true) |> DataFrame
@test getdata(Iris(), PreprocessData.Train) == iris

# load
@test PreprocessData.load(Iris()) == iris

# splits
@test splits(Iris(), iris, [1,3], [2,4]) == (iris[[1,3],:], iris[[2,4],:])

# loadheader
header = ["sepal length", "sepal width", "petal length", "petal width", "class"]
@test PreprocessData.loadheader(joinpath(irispath, "header.csv"),Iris()) == header

# getheader
@test PreprocessData.getheader(Iris()) == header

#getfilename
@test PreprocessData.getfilename(joinpath(irispath, "data-train.csv")) == "data-train"

# df_to_array
dfOrig = DataFrame(A=[1,2,3],B=["a","b","c"],C=[5,6,7], D=[1,1,1])
array1 = [1 "a";2 "b";3 "c"]
array2 = [5 1;6 1;7 1]
@test df_to_array(dfOrig, cols=2) == (array1, array2)

# df_or_array
@test PreprocessData.df_or_array(dfOrig, true, cols=2) == (array1, array2)

# labels
@test PreprocessData.labels(Iris()) == iris[:,end]

# classes
@test classes(Iris()) == ["Iris-setosa", "Iris-versicolor", "Iris-virginica"]

# binarize
dfOrig = DataFrame(A=[1,2,3],B=["a","b","c"],Target=[5,6,7])
@test binarize(dfOrig[:,:Target], [5]) == [1; 0; 0]
@test binarize(dfOrig[:,:Target], [5,6]) == [true; true; false]

# makerow
@test PreprocessData.makerow([1;2;3]) == [1 2 3]

# selectnumeric (and ignore Target)
dfOrig = DataFrame(A=[1,2,3],B=["a","b","c"],Target=[5,6,7])
@test PreprocessData.selectnumeric(dfOrig) == [1]

# selectnumeric
array = [1 "a";2 "b";3 "c"]
@test PreprocessData.selectnumeric(array) == (:,[1])

# l2norm
dfOrig = DataFrame(A=[1,2,3],B=["a","b","c"],C=[5.4,6.3,7.2],Target=[5,6,7])
array = [1 "a" 5.4;2 "b" 6.3;3 "c" 7.2]
Norm = [norm([1,2,3]); norm([5.4,6.3,7.2])]
@test l2norm(dfOrig) == Norm
@test l2norm(array) == PreprocessData.makerow(Norm)

# minmax
min = [minimum([1,2,3]); minimum([5.4,6.3,7.2])]
max = [maximum([1,2,3]); maximum([5.4,6.3,7.2])]
@test PreprocessData.minmax(dfOrig) == (min, max)
@test PreprocessData.minmax(array) == (PreprocessData.makerow(min),PreprocessData.makerow(max))

# _select
dfOrig = DataFrame(A=[1,2,3],B=["a","b","c"],C=[5.4,6.3,7.2],Target=[5,6,7])
dfNew = DataFrame(A=[1,2,3],B=["a","b","c"])
@test PreprocessData._select(dfOrig, 1:2) == dfNew
array = [1 "a" 5.4;2 "b" 6.3;3 "c" 7.2]
arrayNew = [1 "a";2 "b";3 "c"]
@test PreprocessData._select(array, (:,1:2)) == arrayNew

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
"""
    BGKandStream!(fout, feq, ftemp, Fx, Fy, τ)

Performs a BGK collision operation with a WFM forcecorrection and a subsequent streaming of the resulting populations.

# Arguments

- `fout :: Array{<:Number,3}`: Streamed distribution after the collision processes
- `feq :: Array{<:Number,3}`: Equilibrium distribution, computed with equilibrium!
- `ftemp :: Array{<:Number,3}`: Temporary distribution from the time step before, only useful if `` \\tau \\neq 1``
- `Fx :: Array{<:Number,2}`: Sum of forces acting on the fluid in x-direction
- `Fy :: Array{<:Number,2}`: Sum of forces acting on the fluid in y-direction
- `τ <: Number`: Relaxtion time, if not supplied `` \\tau = 1 `` assumed

# Mathematics

The lattice Boltzmann equation in its discretized format is relatively simple to write down

`` f_{\\alpha}(\\mathbf{x}+\\mathbf{e}_{\\alpha}\\Delta t, t+\\Delta t) - f_{\\alpha}(\\mathbf{x}, t) = -\\frac{1}{\\tau}(f_{\\alpha}(\\mathbf{x}, t) - f^{\\text{eq}}_{\\alpha}(\\mathbf{x}, t)) + \\Delta t \\mathcal{S}_{\\alpha}, ``

where the collision kernel is approximated with a BKG single relaxation time (SRT)

`` \\Omega_{\\alpha} = \\frac{1}{\\tau}(f_{\\alpha}(\\mathbf{x}, t) - f^{\\text{eq}}_{\\alpha}(\\mathbf{x}, t), ``

and a source term `` \\mathcal{S} `` which is given by

`` \\mathcal{S}_{\\alpha} = \\frac{3 w_{\\alpha}}{e_{\\alpha x}e_{\\alpha x}+e_{\\alpha y}e_{\\alpha y}}\\mathbf{e_{\\alpha}}\\cdot\\mathbf{F}_{\\alpha} . ``

The term `` e_{\\alpha x}e_{\\alpha x}+e_{\\alpha y}e_{\\alpha y} `` is either zero for the zeroth population, 1 for the first four populations or 2 for the remaining ones.

# Examples

```jldoctest
julia> using Swalbe

julia> feq = ones(5,5,9); ftemp = zeros(5,5,9); fout = zeros(5,5,9);

julia> feq[1,1,:] .= 2.0 # To check the streaming process
9-element view(::Array{Float64, 3}, 1, 1, :) with eltype Float64:
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0
 2.0

julia> Swalbe.BGKandStream!(fout, feq, ftemp, zeros(5,5), zeros(5,5))

julia> fout[:,:,6] # The value 2 should have moved one down and one to the right!
5×5 Matrix{Float64}:
 1.0  1.0  1.0  1.0  1.0
 1.0  2.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0

```

# References

- [Salmon](https://www.ingentaconnect.com/contentone/jmr/jmr/1999/00000057/00000003/art00005#)
- [Dellar](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.65.036309)
- [Peng et al.](https://onlinelibrary.wiley.com/doi/full/10.1002/fld.4726)

See also: [`Swalbe.equilibrium`](@ref)
"""
function BGKandStream!(fout, feq, ftemp, Fx, Fy, τ)
    # All distribution functions
    fe0, fe1, fe2, fe3, fe4, fe5, fe6, fe7, fe8 = viewdists(feq)
    ft0, ft1, ft2, ft3, ft4, ft5, ft6, ft7, ft8 = viewdists(ftemp)
    fo0, fo1, fo2, fo3, fo4, fo5, fo6, fo7, fo8 = viewdists(fout)

    omeg = 1-1/τ
    # Collision for the nine populations
    # Zeroth, no forcing!
    fo0 .= omeg .* ft0 .+ 1/τ .* fe0
    # Straight ones, force correction is 1/3, 3*1/9
    fo1 .= omeg .* ft1 .+ 1/τ .* fe1 .+ 1/3 .* Fx
    fo2 .= omeg .* ft2 .+ 1/τ .* fe2 .+ 1/3 .* Fy
    fo3 .= omeg .* ft3 .+ 1/τ .* fe3 .- 1/3 .* Fx
    fo4 .= omeg .* ft4 .+ 1/τ .* fe4 .- 1/3 .* Fy
    # Diagonal ones, force correction 1/24 -> 3/2*1/36
    fo5 .= omeg .* ft5 .+ 1/τ .* fe5 .+ 1/24 .* (Fx .+ Fy)
    fo6 .= omeg .* ft6 .+ 1/τ .* fe6 .+ 1/24 .* (Fy .- Fx)
    fo7 .= omeg .* ft7 .+ 1/τ .* fe7 .- 1/24 .* (Fx .+ Fy)
    fo8 .= omeg .* ft8 .+ 1/τ .* fe8 .+ 1/24 .* (Fx .- Fy)

    # This is the streaming step with implicite periodic boundarys
    circshift!(ft0, fo0, (0, 0))
    circshift!(ft1, fo1, (1, 0))
    circshift!(ft2, fo2, (0, 1))
    circshift!(ft3, fo3, (-1, 0))
    circshift!(ft4, fo4, (0, -1))
    circshift!(ft5, fo5, (1, 1))
    circshift!(ft6, fo6, (-1, 1))
    circshift!(ft7, fo7, (-1, -1))
    circshift!(ft8, fo8, (1, -1))

    # Overwrite fout with ftemp
    fout .= ftemp
    return nothing
end
# With State struct
function BGKandStream!(state::LBM_state_2D, sys::SysConst)
    # All distribution functions
    fe0, fe1, fe2, fe3, fe4, fe5, fe6, fe7, fe8 = viewdists(state.feq)
    ft0, ft1, ft2, ft3, ft4, ft5, ft6, ft7, ft8 = viewdists(state.ftemp)
    fo0, fo1, fo2, fo3, fo4, fo5, fo6, fo7, fo8 = viewdists(state.fout)

    omeg = 1-1/sys.τ

    # Collision for the nine populations
    # Zeroth, no forcing!
    fo0 .= omeg .* ft0 .+ 1/sys.τ .* fe0
    # Straight ones, force correction is 1/3, 3*1/9
    fo1 .= omeg .* ft1 .+ 1/sys.τ .* fe1 .+ 1/3 .* state.Fx
    fo2 .= omeg .* ft2 .+ 1/sys.τ .* fe2 .+ 1/3 .* state.Fy
    fo3 .= omeg .* ft3 .+ 1/sys.τ .* fe3 .- 1/3 .* state.Fx
    fo4 .= omeg .* ft4 .+ 1/sys.τ .* fe4 .- 1/3 .* state.Fy
    # Diagonal ones, force correction 1/24 -> 3/2*1/36
    fo5 .= omeg .* ft5 .+ 1/sys.τ .* fe5 .+ 1/24 .* (state.Fx .+ state.Fy)
    fo6 .= omeg .* ft6 .+ 1/sys.τ .* fe6 .+ 1/24 .* (state.Fy .- state.Fx)
    fo7 .= omeg .* ft7 .+ 1/sys.τ .* fe7 .- 1/24 .* (state.Fx .+ state.Fy)
    fo8 .= omeg .* ft8 .+ 1/sys.τ .* fe8 .+ 1/24 .* (state.Fx .- state.Fy)

    # This is the streaming step with implicite periodic boundarys
    circshift!(ft0, fo0, (0, 0))
    circshift!(ft1, fo1, (1, 0))
    circshift!(ft2, fo2, (0, 1))
    circshift!(ft3, fo3, (-1, 0))
    circshift!(ft4, fo4, (0, -1))
    circshift!(ft5, fo5, (1, 1))
    circshift!(ft6, fo6, (-1, 1))
    circshift!(ft7, fo7, (-1, -1))
    circshift!(ft8, fo8, (1, -1))

    # Overwrite fout with ftemp
    state.fout .= state.ftemp
    return nothing
end

# Explicit version for the case τ = 1, which simplifies the computation quite a lot.
function BGKandStream!(fout, feq, ftemp, Fx, Fy)
    # All distribution functions
    fe0, fe1, fe2, fe3, fe4, fe5, fe6, fe7, fe8 = viewdists(feq)
    ft0, ft1, ft2, ft3, ft4, ft5, ft6, ft7, ft8 = viewdists(ftemp)
    fo0, fo1, fo2, fo3, fo4, fo5, fo6, fo7, fo8 = viewdists(fout)

    # Collision for the nine populations
    # Zeroth, no forcing!
    fo0 .= fe0
    # Straight ones, force correction is 1/3, 3*1/9
    fo1 .= fe1 .+ 1/3 .* Fx
    fo2 .= fe2 .+ 1/3 .* Fy
    fo3 .= fe3 .- 1/3 .* Fx
    fo4 .= fe4 .- 1/3 .* Fy
    # Diagonal ones, force correction 1/24 -> 3/2*1/36
    fo5 .= fe5 .+ 1/24 .* (Fx .+ Fy)
    fo6 .= fe6 .+ 1/24 .* (Fy .- Fx)
    fo7 .= fe7 .- 1/24 .* (Fx .+ Fy)
    fo8 .= fe8 .+ 1/24 .* (Fx .- Fy)

    # This is the streaming step with implicite periodic boundarys
    circshift!(ft0, fo0, (0, 0))
    circshift!(ft1, fo1, (1, 0))
    circshift!(ft2, fo2, (0, 1))
    circshift!(ft3, fo3, (-1, 0))
    circshift!(ft4, fo4, (0, -1))
    circshift!(ft5, fo5, (1, 1))
    circshift!(ft6, fo6, (-1, 1))
    circshift!(ft7, fo7, (-1, -1))
    circshift!(ft8, fo8, (1, -1))

    # Overwrite fout with ftemp
    fout .= ftemp
    return nothing
end
# Explicit version for the case τ = 1, with state struct
function BGKandStream!(state::LBM_state_2D)
    # All distribution functions
    fe0, fe1, fe2, fe3, fe4, fe5, fe6, fe7, fe8 = viewdists(state.feq)
    ft0, ft1, ft2, ft3, ft4, ft5, ft6, ft7, ft8 = viewdists(state.ftemp)
    fo0, fo1, fo2, fo3, fo4, fo5, fo6, fo7, fo8 = viewdists(state.fout)

    # Collision for the nine populations
    # Zeroth, no forcing!
    fo0 .= fe0
    # Straight ones, force correction is 1/3, 3*1/9
    fo1 .= fe1 .+ 1/3 .* state.Fx
    fo2 .= fe2 .+ 1/3 .* state.Fy
    fo3 .= fe3 .- 1/3 .* state.Fx
    fo4 .= fe4 .- 1/3 .* state.Fy
    # Diagonal ones, force correction 1/24 -> 3/2*1/36
    fo5 .= fe5 .+ 1/24 .* (state.Fx .+ state.Fy)
    fo6 .= fe6 .+ 1/24 .* (state.Fy .- state.Fx)
    fo7 .= fe7 .- 1/24 .* (state.Fx .+ state.Fy)
    fo8 .= fe8 .+ 1/24 .* (state.Fx .- state.Fy)

    # This is the streaming step with implicite periodic boundarys
    circshift!(ft0, fo0, (0, 0))
    circshift!(ft1, fo1, (1, 0))
    circshift!(ft2, fo2, (0, 1))
    circshift!(ft3, fo3, (-1, 0))
    circshift!(ft4, fo4, (0, -1))
    circshift!(ft5, fo5, (1, 1))
    circshift!(ft6, fo6, (-1, 1))
    circshift!(ft7, fo7, (-1, -1))
    circshift!(ft8, fo8, (1, -1))

    # Overwrite fout with ftemp
    state.fout .= state.ftemp
    return nothing
end

function BGKandStream!(fout, feq, ftemp, F::Vector, τ)
    # All distribution functions
    fe0, fe1, fe2 = viewdists_1D(feq)
    ft0, ft1, ft2 = viewdists_1D(ftemp)
    fo0, fo1, fo2 = viewdists_1D(fout)

    omeg = 1-1/τ
    # Collision for the nine populations
    # Zeroth, no forcing!
    fo0 .= omeg .* ft0 .+ 1/τ .* fe0
    # Straight ones, force correction is 1/3, 3*1/9
    fo1 .= omeg .* ft1 .+ 1/τ .* fe1 .+ 1/2 .* F
    fo2 .= omeg .* ft2 .+ 1/τ .* fe2 .- 1/2 .* F

    # This is the streaming step with implicite periodic boundarys
    circshift!(ft0, fo0, 0)
    circshift!(ft1, fo1, 1)
    circshift!(ft2, fo2, -1)

    # Overwrite fout with ftemp
    fout .= ftemp
    return nothing
end

# Case τ=1
function BGKandStream!(fout, feq, ftemp, F::Vector)
    # All distribution functions
    fe0, fe1, fe2 = viewdists_1D(feq)
    ft0, ft1, ft2 = viewdists_1D(ftemp)
    fo0, fo1, fo2 = viewdists_1D(fout)

    # Collision for the three populations
    # Zeroth, no forcing!
    fo0 .= fe0
    # Straight ones, force correction is 1/3, 3*1/9
    fo1 .= fe1 .+ 1/2 .* F
    fo2 .= fe2 .- 1/2 .* F

    # This is the streaming step with implicite periodic boundarys
    circshift!(ft0, fo0, 0)
    circshift!(ft1, fo1, 1)
    circshift!(ft2, fo2, -1)

    # Overwrite fout with ftemp
    fout .= ftemp
    return nothing
end
# with state struct
function BGKandStream!(state::State_1D, sys::SysConst_1D)
    # All distribution functions
    fe0, fe1, fe2 = viewdists_1D(state.feq)
    ft0, ft1, ft2 = viewdists_1D(state.ftemp)
    fo0, fo1, fo2 = viewdists_1D(state.fout)

    omeg = 1-1/sys.τ
    # Collision for the nine populations
    # Zeroth, no forcing!
    fo0 .= omeg .* ft0 .+ 1/sys.τ .* fe0
    # Straight ones, force correction is 1/3, 3*1/9
    fo1 .= omeg .* ft1 .+ 1/sys.τ .* fe1 .+ 1/2 .* state.F
    fo2 .= omeg .* ft2 .+ 1/sys.τ .* fe2 .- 1/2 .* state.F

    # This is the streaming step with implicite periodic boundarys
    circshift!(ft0, fo0, 0)
    circshift!(ft1, fo1, 1)
    circshift!(ft2, fo2, -1)

    # Overwrite fout with ftemp
    state.fout .= state.ftemp
    return nothing
end

function BGKandStream!(state::State_1D)
    # All distribution functions
    fe0, fe1, fe2 = viewdists_1D(state.feq)
    ft0, ft1, ft2 = viewdists_1D(state.ftemp)
    fo0, fo1, fo2 = viewdists_1D(state.fout)

    # Collision for the three populations
    # Zeroth, no forcing!
    fo0 .= fe0
    # Straight ones, force correction is 1/3, 3*1/9
    fo1 .= fe1 .+ 1/2 .* state.F
    fo2 .= fe2 .- 1/2 .* state.F

    # This is the streaming step with implicite periodic boundarys
    circshift!(ft0, fo0, 0)
    circshift!(ft1, fo1, 1)
    circshift!(ft2, fo2, -1)

    # Overwrite fout with ftemp
    state.fout .= state.ftemp
    return nothing
end

"""
    viewdists(f)

Generates a view for all nine populations of a **D2Q9** distribution function.

# Examples
```jldoctest
julia> using Swalbe, Test

julia> ftest = reshape(collect(1.0:225.0),5,5,9);

julia> f0, f1, f2, f3, f4, f5, f6, f7, f8 = Swalbe.viewdists(ftest);

julia> @test all(f3 .== ftest[:,:,4])
Test Passed
  Expression: all(f3 .== ftest[:, :, 4])

```

See also: [`Swalbe.BGKandStream`](@ref)
"""
function viewdists(f)
    f0 = view(f, :, :, 1)
    f1 = view(f, :, :, 2)
    f2 = view(f, :, :, 3)
    f3 = view(f, :, :, 4)
    f4 = view(f, :, :, 5)
    f5 = view(f, :, :, 6)
    f6 = view(f, :, :, 7)
    f7 = view(f, :, :, 8)
    f8 = view(f, :, :, 9)

    return f0, f1, f2, f3, f4, f5, f6, f7, f8
end

"""
    viewdists_1D(f)

Generates a view for all three populations of a **D1Q3** distribution function.

# Examples
```jldoctest
julia> using Swalbe, Test

julia> ftest = reshape(collect(1.0:30.0),10,3);

julia> f0, f1, f2 = Swalbe.viewdists_1D(ftest);

julia> @test all(f1 .== ftest[:,2])
Test Passed
  Expression: all(f1 .== ftest[:, 2])
```

See also: [`Swalbe.BGKandStream`](@ref)
"""
function viewdists_1D(f)
    f0 = view(f, :, 1)
    f1 = view(f, :, 2)
    f2 = view(f, :, 3)

    return f0, f1, f2
end


text:
module ModifiedBoosting

using Distributions,SpecialFunctions
using LinearAlgebra, MultivariateStats
using StatsBase

mutable struct Node

    is_leaf     :: Bool
    is_root     :: Bool
    left        :: Node
    right       :: Node
    parent      :: Node

    # branch parameters
    var         :: Int64
    val         :: Float64

    lower       :: Float64
    upper       :: Float64

    # leaf parameters
    μ           :: Float64

    ## Parameters for computing weights efficiently
    current_weight  :: Float64

    # Node(is_leaf, is_root, left, right, parent, var, val, lower, upper, μ, current_weight)    = new(is_leaf, is_root, left, right, parent, var, val, lower, upper, μ, current_weight)
    Node() = new()

end

function Root()
    x = Node()

    x.is_leaf = true
    x.is_root = true
    x.left = x
    x.right = x
    x.parent = x

    x.var = 0
    x.val = 0.0
    x.lower = 0.0
    x.upper = 1.0
    x.μ = 0.0
    x.current_weight = 1.0

    return x
end

function InitRoot(X::Array{Float64,2})
    x = Root()
    return x
end

# binary split
function AddLeaves!(x::Node)
    x.left = Node()
    x.right = Node()
    x.is_leaf = false

    ## Initialize the leafs when splitting
    x.left.is_leaf = true
    x.left.parent = x
    x.left.var = 0
    x.left.val = 0.0
    x.left.is_root = false
    x.left.lower = 0.0
    x.left.upper = 1.0
    x.left.μ = 0.0
    x.left.current_weight = 0.0
    x.right.is_leaf = true
    x.right.parent = x
    x.right.var = 0
    x.right.val = 0.0
    x.right.is_root = false
    x.right.lower = 0.0
    x.right.upper = 1.0
    x.right.μ = 0.0
    x.right.current_weight = 0.0

    return nothing
end

mutable struct Hypers
    α       :: Float64
    β        :: Float64
    γ       :: Float64
    σ       :: Float64
    σᵤ    :: Float64
    shape       :: Float64
    width       :: Float64
    num_tree    :: Int64
    s           :: Vector{Float64}

    ## Double hypers
    ̂σ   :: Float64
    ̂σᵤ :: Float64
    α_scale :: Float64
    α_shape_1 :: Float64
    α_shape_2 :: Float64

    Hypers()    = new()
end

# hyper-parameters involved
function Hypers(X::Array{Float64, 2},
            Y::Vector{Float64};
            α = 0.95,
            β = 2.0,
            γ = 0.95,
            k = 2.0,
            width = 0.1,
            shape = 1.0,
            num_tree = 50,
            α_scale = 0.0,
            α_shape_1 = 0.5,
            α_shape_2 = 1.0)

            out = Hypers()
            out.α = α
            out.β = β
            out.γ = γ
            out.σ = GetSigma(X,Y)
            out.σᵤ = 0.5 / (k * sqrt(num_tree))
            out.shape = shape
            out.width = width
            out.num_tree = num_tree
            out.s = ones(size(X,2)) / size(X,2)

            out.̂σ = out.σ
            out.̂σᵤ = out.σᵤ

            out.α_scale = α_scale == 0.0 ? size(X,2) : α_scale
            out.α_shape_1 = α_shape_1
            out.α_shape_2 = α_shape_2

            return out;

end

function GetLimits(x::Node)
    y = x
    lower = 0.0
    upper = 1.0
    my_bool = y.is_root ? false : true # my_bool is not defined, so it is initially false
    while my_bool
        is_left = (y.parent.left == y)
        y = y.parent
        my_bool = y.is_root ? false : true # when y is not root
        if y.var == x.var
            my_bool = false
            if is_left
                upper = y.val
                lower = y.lower
            else
                upper = y.upper
                lower = y.val
            end
        end
    end
    return lower, upper
end

function BirthLeaves!(x::Node, hypers::Hypers)
    if !x.is_leaf
        error("Node must be a leaf")
    end
    AddLeaves!(x)
    wv = StatsBase.Weights(hypers.s)
    x.var = sample(wv)
    # x.var = rand(1:P) ## TODO change this to allow non-iid selection
    lower, upper = GetLimits(x)
    x.lower = lower
    x.upper = upper
    x.val = (upper - lower) * rand() + lower
end

function is_left(n::Node)
    if n == n.parent.left
        return true
    else
        return false
    end
end

function DeleteLeaves!(x::Node)
    x.left = x
    x.right = x
    x.is_leaf = true
    return nothing
end

function depth(x::Node)
    x.is_root ? 0 : 1 + depth(x.parent)
end

function leaves!(x::Node, leaves::Vector{Node})
    if x.is_leaf
        push!(leaves, x);
    else
        leaves!(x.left, leaves);
        leaves!(x.right, leaves);
    end
end

function leaves(x::Node)
    leaves = Node[];
    leaves!(x, leaves);
    return leaves;
end

function fried(X::Array{Float64, 2})
    n = size(X,1)
    p = size(X,2)

    Y = zeros(n)
    for i in 1:n
        Y[i] = 10 * sin(2 * pi * X[i,1] * X[i,2]) +
            20 * (X[i,3] - 0.5)^2 + 10 * X[i,4] + 5 * X[i,5]
    end
    return Y
end

function test_data(n::Int64; sigma = 1.0, p = 10)
    X = rand(n,p)
    Y = fried(X) + sigma * randn(n)
    return X, Y
end

function logit(x::Float64)
    return log(x) - log(1 - x)
end

function normalize(y::Vector{Float64})
    a = minimum(y)
    b = maximum(y)
    return (y .- a) ./ (b - a) .- 0.5
end

function expit(x::Float64)
    return 1.0 / (1 + exp(-x))
end

function activation(x::Float64, c::Float64, hypers::Hypers)
    return 1- expit((x - c) / hypers.width)
end

function get_var_counts(forest::Vector{Node}, hypers::Hypers)
    counts = zeros(Int64, size(hypers.s,1))
    num_tree = length(forest)
    for t in 1:num_tree
        get_var_counts!(forest[t], counts)
    end
    return counts
end

function get_var_counts!(node::Node, counts::Vector{Int64})
    if !node.is_leaf
        counts[node.var] += 1
        get_var_counts!(node.left, counts)
        get_var_counts!(node.right, counts)
    end
end

end  # module

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
using Plots
using FLoops
using BenchmarkTools
using Octavian

struct Spins
	N::Int32
	M::Int32
	A::Array{Float32}
end

function Init(N = 100::Int,M = 100::Int)
	S = Spins(N, M, rand(Float32, N+2 ,M+2)*2π .- π)
	S.A[1,:] .= 0.0
	S.A[N + 2, :] .= 0.0
	S.A[:,1] .= 0.0
	S.A[:,M + 2] .= 0.0

	boolarr1 = S.A[2:N + 1,2:M + 1] .== 0.0
	temparr = rand(Float32,N,M)
	boolarr2 = temparr .> 0.7
	boolarr = boolarr1 .& boolarr2
	S.A[2:N + 1,2:M + 1] -= boolarr
	boolarr = boolarr1 .& .~boolarr2
	S.A[2:N + 1,2:M + 1] += boolarr

	return S
end

function EnergyParallel(S::Spins)
	E = zeros(Float32,S.N+2,S.M+2)
	p = [(-1,0),(1,0),(0,-1),(0,1)]
	@floop ThreadedEx(basesize = 1, nestlevel=Val(:inf), simd = true) for i in p
			Epart = matmul(S.A[2:S.N+1,2:S.M+1],
					S.A[2+i[1]:S.N+1+i[1],
						2+i[2]:S.M+1+i[2]])
			@reduce(E += Epart)
	end
	return sum(E)
end

function Magnetisation(S::Spins)
	mx = sum(cos.(S.A))/ *(S.N,S.M)
	my = sum(sin.(S.A))/ *(S.N,S.M)
	return √(mx^2 + my^2)
end

function Advance(S::Spins,max_steps=10000::Int,T=5.0::Float32)
	eps = 0.001

	E = zeros(Float32,max_steps + 1)
	m = zeros(Float32,max_steps + 1)

	Enow = EnergyParallel(S)
	E[1] = Enow
	m[1] = Magnetisation(S)

	step = 1

	#while *(step<max_steps, m[step] > eps, 1 - m[step] > eps)
	for _ in 1:max_steps

		step += 1
		i = rand(2:S.N+1)
		j = rand(2:S.M+1)

		InitVal = S.A[i,j]

		S.A[i,j]= (S.A[i,j] + 2π*(rand()-0.5))%π

		Enext = EnergyParallel(S)
		dE = Enext - Enow

        if dE < 0
			Enow = Enext
		elseif	rand() < exp(-dE/T)
            Enow = Enext
        else
            S.A[i,j] = InitVal
		end

		E[step] = Enow
		#m[step] = Magnetisation(S)
	end
	return E,m,S.A
end

function Run(N=100::Int,M=100::Int,max_steps=10000::Int,T=Float32(5.0)::Float32)
	S = Init(N,M)
	h1 = heatmap(S.A)
	E,m,S=Advance(S,max_steps,T)

	"""l = @layout [a{0.6h};b{0.5h}]
	p1 = plot(m, linecolor = :blue, label="m")
	p2 = plot(E, lc = :red, label="E")
	plot(p1,p2,layout = l)"""

	h2 = heatmap(S)
	savefig(plot(h1,h2), "./Out")
end

@time Run(100,100,10)
println(Threads.nthreads())

max_stars_repo_name: anirbrhm/curious-eavesdropper

max_issues_repo_name: anirbrhm/curious-eavesdropper

max_forks_repo_name: anirbrhm/curious-eavesdropper


text:
# SCIP defines for character, string, and numeric constants
{% for def_name, def_ref in parser.defines.items() %}const {{ def_name }} = {{ def_ref }}
{% endfor %}
# SCIP defines for known C types and pointers to SCIP structs.
{% for def_name, def_ref in parser.typealiases.items() %}typealias {{ def_name }} {{ def_ref }}
{% endfor %}


text:
using JuliaEuler
using Dates
using Test

include("helper.jl")
include("ex1_24.jl")
include("ex25_50.jl")

max_stars_repo_name: binnisb/JuliaEuler

max_issues_repo_name: binnisb/JuliaEuler

max_forks_repo_name: binnisb/JuliaEuler


text:
abstract type MissionNode ; end
abstract type MissionNodeIterator <: MissionNode; end

mutable struct Agent <: MissionNode
    fn                  ::Union{ Thunk, Function }
    machines            ::Vector{String}
    priority            ::Int
    min_workers         ::Int
    max_workers         ::Int
    dispatched_workers  ::Vector{Int}
end
Agent( fn::Union{Thunk,Function}, machines::Vector{String} ) = Agent( fn, machines, 1, 1, 1, [] )
Agent( fn::Union{Thunk,Function}, machines::Vector{String}, maxworkers::Int ) = Agent( fn, machines, 1, 1, maxworkers, [] )

mutable struct Stash <: MissionNode
    src                 ::String
    fn                  ::Union{Nothing, Thunk, Function}
    machines            ::Vector{ String }
    priority            ::Int
    min_workers         ::Int
    max_workers         ::Int
    dispatched_workers  ::Vector{ Int }
end

Stash( src::String, machines::Vector{String} ) = Stash( src, nothing, machines, 1, 1, 1, [] )
Stash( src::String, machines::Vector{String}, maxworkers::Int ) = Stash( src, nothing, machines, 1, 1, maxworkers,[] )
Stash( src::String, fn::Union{Thunk,Function}, machines::Vector{String} ) = Stash( src, fn, machines, 1, 1, 1, [] )
Stash( src::String, fn::Union{Thunk,Function}, machines::Vector{String}, maxworkers::Int ) = Stash( src, fn, machines, 1, 1, maxworkers,[] )

mutable struct StashIterator <: MissionNodeIterator
    iterator            ::Union{FileIterator, Vector{ String } }
    iteratorstate       ::Any
    #below matches exactly with Stash constructor...
    src                 ::String
    fn                  ::Union{Nothing, Thunk, Function}
    machines            ::Vector{ String }
    priority            ::Int
    min_worker_constant ::Int
    min_workers         ::Int
    max_workers         ::Int
    dispatched_workers  ::Vector{ Int }
end

function StashIterator( src_iterator::Union{FileIterator,Vector{String}},
                        fn::Union{Thunk,Function},
                        machines::Vector{String},
                        priority::Int = 1,
                        min_workers_per_iter::Int = 1,
                        max_workers_per_iter::Int = 1 )
    result, state   = iterate( src_iterator )
    return StashIterator(   src_iterator, state, result, fn, machines,
                            priority,
                            copy(min_workers_per_iter),
                            min_workers_per_iter, max_workers_per_iter,
                            Vector{Int}(undef, 0) )
end

"""
    next!( si::MissionNodeIterator )

Updates a MissionNodeIterator's source parameter by updating the internal iterator

"""
function next!( si::MissionNodeIterator )
    result = iterate( si.iterator, si.iteratorstate )
    if isnothing( result )
        si.src = nothing
    else
        si.src, si.iteratorstate = result
    end
end

max_stars_repo_name: caseykneale/MrPhelps

max_issues_repo_name: caseykneale/MrPhelps

max_forks_repo_name: caseykneale/MrPhelps


text:
using LinearAlgebra: det
using PyCall: pyimport, PyObject

"""
    identify_seam_end_inds(points::AbstractVector{T}) where T

Find all indices that are at the second half of a seam.
"""
function identify_seam_end_inds(points::AbstractVector{T}) where T
    point_set = Set{T}()
    seam_end_inds = Int[]
    for (i, p) in enumerate(points)
        if p in point_set
            push!(seam_end_inds, i)
        else
            push!(point_set, p)
        end
    end
    return seam_end_inds
end

"""
    find_latest_ending_seam(points::AbstractVector{T},
        seam_end_inds::AbstractVector{Int}) where T

Find the latest ending seam.
"""
function find_latest_ending_seam(points::AbstractVector{T},
    seam_end_inds::AbstractVector{Int}) where T
    begin_inds, end_inds = Int[], Int[]
    for ind in reverse(1:length(seam_end_inds))
        j = seam_end_inds[ind]
        if (length(end_inds) != 0) && (j != end_inds[end]-1)
            return begin_inds, end_inds, ind
        end
        i = findfirst(isequal(points[j]), points)
        if (length(begin_inds) != 0) && (i != begin_inds[end]+1)
            return begin_inds, end_inds, ind
        end
        push!(begin_inds, i); push!(end_inds, j)
    end
    return begin_inds, end_inds, 0
end

"""
    remove_seams(points::AbstractVector{T}) where T
    remove_seams(points::AbstractVector{T}, seam_end_inds::AbstractVector{Int}) where T

Find the sublists of `points` that result from removing all seams.

# Examples
```jldoctest
DocTestSetup = quote
    using DiscretePDEs
end
points = [0,1,2,3,4,5,6,5,7,3,2,8,9]
DiscretePDEs.remove_seams(points)

# output

3-element Array{Array{Int64,1},1}:
 [0, 1, 2, 8, 9]
 [3, 4, 5, 7]
 [5, 6]
```
"""
remove_seams(points::AbstractVector{T}) where T = remove_seams(points,
    identify_seam_end_inds(points))

# this method is a helper that performs the recursion
function remove_seams(points::AbstractVector{T}, seam_end_inds::AbstractVector{Int}) where T
    if length(seam_end_inds) == 0
        return [points]
    else
        begin_inds, end_inds, ind = find_latest_ending_seam(points, seam_end_inds)
        outside = [points[1:begin_inds[1]]; points[(end_inds[1]+1):end]]
        inside = points[begin_inds[end]:(end_inds[end]-1)]
        outside_seam_end_inds = [filter(x -> x <= begin_inds[1], seam_end_inds[1:ind]);
                                 filter(x -> end_inds[1] < x, seam_end_inds[1:ind]) .-
                                    (end_inds[1] - begin_inds[1])]
        inside_seam_end_inds = filter(x -> begin_inds[end] <= x < end_inds[end],
            seam_end_inds[1:ind]) .- (begin_inds[end]-1)
        return [remove_seams(outside, outside_seam_end_inds);
            remove_seams(inside, inside_seam_end_inds)]
    end
end

"""
    demarcate_self_intersections(points_array::AbstractVector{<:Tuple{<:Real, <:Real}})

In a gds file, if a polygon intersects itself at a vertex, that vertex is not necessarily
repeated. This causes means that `remove_seams` will not correctly identify the seams.
Return an array of points with the missing points added.
"""
function demarcate_self_intersections(points_array::AbstractVector{<:Tuple{<:Real, <:Real}})
    output = typeof(points_array)()
    for i1 in 1:length(points_array)
        i2 = mod1(i1+1, length(points_array))
        p1, p2 = points_array[i1], points_array[i2]
        # find all points on the line from p1 to p2 using determinant
        collinears = filter(p -> det([[p1..., 1] [p2..., 1] [p..., 1]]) == 0, points_array)
        # p = t * p1 + (1-t) * p2 implies t = (p - p2) ⋅ (p1 - p2)/|p1 - p2|^2
        p_diff = collect(p1) - collect(p2)
        p_norm_square = transpose(p_diff) * p_diff
        fractional_dists = [transpose(p_diff) * (collect(p) - collect(p2))/p_norm_square for p in collinears]
        # remove anything on the line but not strictly between p1 and p2
        between_inds = 0 .< fractional_dists .< 1
        collinears = collinears[between_inds]
        fractional_dists = fractional_dists[between_inds]
        # add all the points in order from p1 to p2, not including p2
        append!(output, [p1; collinears[sortperm(fractional_dists)]])
    end
    return output
end

"""
    identify_outer_polygon!(
        points_array::AbstractVector{<:AbstractVector{<:Tuple{<:Real, <:Real}}})

Find the outermost polygon, which is the one containing an extremal point. Move
this outermost polygon to the end of the list.
"""
function identify_outer_polygon!(
    points_array::AbstractVector{<:AbstractVector{<:Tuple{<:Real, <:Real}}})
    all_points = vcat(points_array...)
    extremal_point = all_points[argmax([p[1] for p in all_points])]
    ind = findfirst(arr -> extremal_point in arr, points_array)
    extremal_arr = points_array[ind]
    return push!(deleteat!(points_array, ind), extremal_arr)
end

"""
    GeoCode(cell::PyObject;
        layer_translation::Dict{Int, <:Real}=Dict{Int, Float64}(),
        layer_extrusion::Dict{Int, <:Real}=Dict{Int, Float64}(),
        layer_volume_difference::AbstractVector{Tuple{Int, Int, Bool}}=
            Tuple{Int, Int, Bool}[],
        layer_surface_groups::Dict{String, <:AbstractVector{Int}}=
            Dict{String, Vector{Int}}(),
        layer_volume_groups::Dict{String, <:AbstractVector{Int}}=
            Dict{String, Vector{Int}}(),
        layer_mesh_size::Dict{Int, <:Real}=Dict{Int, Float64}())

Construct a GeoCode object from a gds cell.

# Arguments
    - `cell::PyObject`: the gds cell.
    - `layer_translation::Dict{Int, <:Real}`: a mapping from layer to a translation in z.
    - `layer_extrusion::Dict{Int, <:Real}`: a mapping from layer to an extrusion in z,
        after the translation in z.
    - `layer_volume_difference::AbstractVector{Tuple{Int, Int, Bool}}`: a vector of tuples
        of the form `(layer_1, layer_2, keep_layer_2)` indicating that the extruded volumes
        on layer `layer_2` should be subtracted from the extruded volumes on layer
        `layer_1`. The extruded volumes on layer `layer_2` will be kept if `keep_layer_2`.
    - `layer_surface_groups::Dict{String, <:AbstractVector{Int}}`: a mapping from group name
        to a vector of layers indicating that the surfaces on those layers will be collected
        into a physical group.
    - `layer_volume_groups::Dict{String, <:AbstractVector{Int}}`: a mapping from group name
        to a vector of layers indicating that the volumes extruded from those layers will be
        collected into a physical group. The layers must be a sublist of the keys in
        `layer_extrusion`.
    - `layer_mesh_size::Dict{Int, <:Real}`: a mapping from layer to mesh size.
"""
function GeoCode(cell::PyObject;
    layer_translation::Dict{Int, <:Real}=Dict{Int, Float64}(),
    layer_extrusion::Dict{Int, <:Real}=Dict{Int, Float64}(),
    layer_volume_difference::AbstractVector{Tuple{Int, Int, Bool}}=Tuple{Int, Int, Bool}[],
    layer_surface_groups::Dict{String, <:AbstractVector{Int}}=Dict{String, Vector{Int}}(),
    layer_volume_groups::Dict{String, <:AbstractVector{Int}}=Dict{String, Vector{Int}}(),
    layer_mesh_size::Dict{Int, <:Real}=Dict{Int, Float64}())
    # flattening ensures all are PolygonSet objects
    elements = cell.copy("", exclude_from_current=true).flatten().elements
    geo_code = GeoCode()
    surface_group_layers = unique(vcat(values(layer_surface_groups)...))
    volume_group_layers = unique(vcat(values(layer_volume_groups)...))
    @assert issubset(volume_group_layers, keys(layer_extrusion))
    layer_to_planes = Dict(l => GeoPlane[] for l in surface_group_layers)
    layer_to_extrusions = Dict(l => GeoSurfaceExtrusion[] for l in keys(layer_extrusion))
    for polygon_set in elements
        for (layer, points_array) in zip(polygon_set.layers, polygon_set.polygons)
            layer = convert(Int, layer) # convert from PyObject
            point_tuples = [(points_array[i,1], points_array[i,2]) for i in 1:size(points_array,1)]
            # find seams and remove them, then create inner polygons and outer polygon
            seamless_polygons = filter(p -> length(p) > 0, remove_seams(demarcate_self_intersections(point_tuples)))
            identify_outer_polygon!(seamless_polygons)
            planes = GeoPlane[] # will need to subtract all planes 1:end-1 from the last
            z = get(layer_translation, layer, 0.0)
            for points in seamless_polygons
                mesh_size = get(layer_mesh_size, layer, 0.0)
                g = GeoCode([(1.0 * p[1], 1.0 * p[2], 1.0 * z, 1.0 * mesh_size)
                    for p in points]) # p is just an x and y coordinate
                append!(geo_code, g)
                push!(planes, g.planes[1])
            end
            # subtract inner polygons from outer polygon
            for p in planes[1:end-1]
                push!(geo_code.surface_differences,
                    GeoSurfaceDifference(planes[end], p, false))
            end
            # add last plane to a physical group
            if layer in keys(layer_to_planes)
                push!(layer_to_planes[layer], planes[end])
            end
            # put extrusion code later since it creates points
            if layer in keys(layer_extrusion)
                extrusion = GeoSurfaceExtrusion(0, 0, layer_extrusion[layer], planes[end])
                push!(geo_code.extrusions, extrusion)
                push!(layer_to_extrusions[layer], extrusion)
            end
        end
    end
    # add any volume differences
    for (l1, l2, keep_tool) in layer_volume_difference
        for v1 in layer_to_extrusions[l1]
            for v2 in layer_to_extrusions[l2]
                push!(geo_code.volume_differences,
                    GeoVolumeDifference(v1, v2, keep_tool))
            end
        end
    end
    # add any physical groups
    for name in keys(layer_surface_groups)
        planes = vcat([layer_to_planes[l] for l in layer_surface_groups[name]]...)
        push!(geo_code.surface_groups, GeoSurfaceGroup(name, planes))
    end
    for name in keys(layer_volume_groups)
        extrusions = vcat([layer_to_extrusions[l] for l in layer_volume_groups[name]]...)
        push!(geo_code.volume_groups, GeoVolumeGroup(name, extrusions))
    end
    return geo_code
end

export gdspy_open
"""
    gdspy_open(file_name::AbstractString, cell_name::AbstractString)

Use gdspy to open a gds file. Return the cell corresponding to `cell_name`.
"""
function gdspy_open(file_name::AbstractString, cell_name::AbstractString)
    gds = pyimport("gdspy").GdsLibrary().read_gds(file_name)
    return gds.cell_dict[cell_name]
end


text:
using ModelingToolkit, DiffEqBase, LinearAlgebra

# Define some variables
@parameters t x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2
eq = Dt(u(t, x)) ~ Dxx(u(t, x))
bcs = [u(0, x) ~ -x * (x - 1) * sin(x),
    u(t, 0) ~ 0, u(t, 1) ~ 0]

domains = [t ∈ (0.0, 1.0),
    x ∈ (0.0, 1.0)]

@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u])
@show pdesys

@test all(isequal.(independent_variables(pdesys), [t, x]))


text:
""" `sys = tf(num, den[, Ts]), sys = tf(gain[, Ts])`

Create as a fraction of polynomials:

`sys::TransferFunction{SisoRational{T,TR}} = numerator/denominator`
where T is the type of the coefficients in the polynomial.

`num`: the coefficients of the numerator polynomial. Either scalar or vector to create SISO systems
or an array of vectors to create MIMO system.

`den`: the coefficients of the denominator polynomial. Either vector to create SISO systems
or an array of vectors to create MIMO system.

`Ts`: Sample time if discrete time system.

Other uses:
`tf(sys)`: Convert `sys` to `tf` form.
`tf("s")`, `tf("z")`: Create the continuous transferfunction `s`.

See also: `zpk`, `ss`
"""
function tf(num::AbstractVecOrMat{<:AbstractVector{T1}}, den::AbstractVecOrMat{<:AbstractVector{T2}}, Ts::TE) where {TE<:TimeEvolution, T1<:Number, T2<:Number}
    # Validate input and output dimensions match
    ny, nu = size(num, 1), size(num, 2)
    if (ny, nu) != (size(den, 1), size(den, 2))
        error("num and den dimensions must match")
    end

    T = promote_type(T1,T2)
    matrix = Matrix{SisoRational{T}}(undef, ny, nu)
    for o=1:ny
        for i=1:nu
            matrix[o, i] = SisoRational{T}(Vector{T}(num[o, i]), Vector{T}(den[o, i]))
        end
    end
    return TransferFunction{TE,SisoRational{T}}(matrix, Ts)
end
tf(num::AbstractVecOrMat{<:AbstractVector{T1}}, den::AbstractVecOrMat{<:AbstractVector{T2}}, Ts::Number) where {T1,T2} =
    tf(num, den, Discrete(Ts))
tf(num::AbstractVecOrMat{<:AbstractVector{T1}}, den::AbstractVecOrMat{<:AbstractVector{T2}}) where {T1,T2} =
    tf(num, den, Continuous())

function tf(num::AbstractVector{T1}, den::AbstractVector{T2}, timeevol::TE) where {TE<:TimeEvolution,T1<:Number, T2<:Number}
    T = promote_type(T1, T2)
    return TransferFunction{TE,SisoRational{T}}(fill(SisoRational{T}(num, den), 1, 1), timeevol)
end
tf(num::AbstractVector{T1}, den::AbstractVector{T2}, Ts::Number) where {T1<:Number, T2<:Number} =
    tf(num, den, Discrete(Ts))
tf(num::AbstractVector{T1}, den::AbstractVector{T2}) where {T1<:Number, T2<:Number} =
    tf(num, den, Continuous())

tf(num::Number, den::Vector, args...) = tf([num], den, args...)

# Cases for just static gain
function tf(D::AbstractArray{T}, timeevol::TE) where {TE<:TimeEvolution, T<:Number}
    ny, nu = size(D, 1), size(D, 2)

    matrix = Matrix{SisoRational{T}}(undef, ny, nu)
    for i in eachindex(D)
        matrix[i] = SisoRational{T}([D[i]], [one(T)])
    end
    return TransferFunction{TE,SisoRational{T}}(matrix, timeevol)
end
tf(D::AbstractArray{T}, Ts::Number) where T = tf(D, Discrete(Ts))
tf(D::AbstractArray{T}) where T = tf(D, Continuous())

tf(n::Number, args...; kwargs...) = tf([n], args...; kwargs...)

tf(sys::AbstractStateSpace) = convert(TransferFunction, sys) # NOTE: Would perhaps like to write TransferFunction{SisoRational}, but couldn't get this to work...

function tf(G::TransferFunction{TE,<:SisoTf{T}}) where {TE<:TimeEvolution,T<:Number}
    convert(TransferFunction{TE,SisoRational{T}}, G)
end

# Function for creation of 's' or 'z' var
function tf(var::AbstractString)
    var != "s" && error("var must be 's' for continuous time tf.")
    return tf([1, 0], [1], Continuous())
end
function tf(var::AbstractString, Ts::Real)
    var != "z" && error("var must be 'z' for discrete time tf.")
    Ts == 0 && error("Ts must not be 0 for discrete time tf.")
    return tf([1, 0], [1], Discrete(Ts))
end

## Constructors for polynomial inputs
function tf(num::AbstractArray{PT}, den::AbstractArray{PT},  timeevol::TE) where {TE<:TimeEvolution,T<:Number, PT <: Polynomials.Polynomial{T}}
    ny, nu = size(num, 1), size(num, 2)
    if (ny, nu) != (size(den, 1), size(den, 2))
        error("num and den dimensions must match")
    end

    matrix = Matrix{SisoRational{T}}(undef, ny, nu)
    for o=1:ny
        for i=1:nu
            matrix[o, i] = SisoRational{T}(num[o, i], den[o, i])
        end
    end
    return TransferFunction{TE,SisoRational{T}}(matrix, timeevol)
end
tf(num::AbstractArray{PT}, den::AbstractArray{PT}, Ts::Number) where {T<:Number, PT <: Polynomials.Polynomial{T}} =
    tf(num, den, Discrete(Ts))
tf(num::AbstractArray{PT}, den::AbstractArray{PT}) where {T<:Number, PT <: Polynomials.Polynomial{T}} =
    tf(num, den, Continuous())

function tf(num::PT, den::PT, timeevol::TE) where {TE<:TimeEvolution, T<:Number, PT <: Polynomials.Polynomial{T}}
    tf(fill(num,1,1), fill(den,1,1), timeevol)
end
tf(num::PT, den::PT, Ts::Number) where {T<:Number, PT <: Polynomials.Polynomial{T}} =
    tf(num, den, Discrete(Ts))
tf(num::PT, den::PT) where {T<:Number, PT <: Polynomials.Polynomial{T}} =
    tf(num, den, Continuous())


text:
"""
    Collocation_points

	Type with vector of collocation points xi and respective quadrature
	integration weights wi.
"""
struct Collocation_points
	xi         #  collocation points
	wi         #  quadrature integration weights
end
function show(io ::IO, object :: Collocation_points)
	println(io, "Collocation points and quadrature weights")
	println(io, "...", length(object.xi), " collocation points")
	println(io, ".xi"); println(io, ".wi")
end

"""
    Discrete_domain

	Define the flow and effort collocation points.
"""
struct Discrete_domain
	flow     :: Collocation_points # zi
	effort   :: Collocation_points # xi
end
function show(io ::IO, object :: Discrete_domain)
	println(io, "Discretization domain points and quadrature weights")
	println(io, "...", length(object.flow.xi), " collocation points for the flow variables")
	println(io, "...", length(object.effort.xi), " collocation points for the effort variables")
	println(io, ".flow"); println(io, ".effort")
end
"""
    Phs

	Port-Hamiltonian object. Describe the interconnection matrices,
	Hamiltonian, constraints, etc.

"""
mutable struct Phs
	J :: Array; # interconnection matrix
	B :: Array; # control/output matrix
	Bd :: Array; # distributed control matrix (optional)
	D :: Array; # direct matrix
	Q  # Q matrix, optional
	R :: Array; # damping matrix (optional)
	G  :: Array # constraint matrix (optional)
	G_D :: Array # direct term constraint matrix (optional)
	Hamiltonian :: Function; # Hamiltonian
	GradHam :: Function; # Hamiltonian gradient
	hessian :: Function # Hamiltonian hessian
	disc_data ::  Discrete_domain # discretization data
	TransfMatrix :: Array # transformation matrix after reduction change of variables
							# always start as identity
	StatesNames :: Dict
	InputsNames :: Array
	auxvars :: Dict # auxiliary variables dictionary

	"""
	    Phs object constructor: interconnection matrices and Hamiltonian
		                        function as inputs.
	"""
	function Phs(J :: Array, B :: Array, D :: Array, Ham :: Function)
		this = new()
		this.J = J
		this.B = B
		this.D = D
		this.R = zeros(size(J))
		this.Hamiltonian = Ham
		this.TransfMatrix = Matrix(1I, size(J,1), size(J,1))
		this.StatesNames = Dict()
		this.auxvars = Dict()
		this.InputsNames = repeat([""],size(B,2))
		this
	end
	"The Hamiltonian matrix Q can be used instead of the Hamiltonian function"
	function Phs(J :: Array, B :: Array, D :: Array, Q)
		this = new()
		this.J = J;	this.B = B;	this.D = D;	this.Q = Q;
		this.R = zeros(size(J))
		Ham(x :: Array) = (0.5* transpose(x)*Q*x)[1]
		this.Hamiltonian = Ham
		GHam(x :: Array) = Q*x
		this.GradHam = GHam
		this.hessian = x->Q
		this.TransfMatrix = Matrix(1I, size(J,1), size(J,1))
		this.StatesNames = Dict()
		this.auxvars = Dict()
		this.InputsNames = repeat([""],size(B,2))
		this
	end
end

"""
    set_constraint!(ph :: Phs, G, G_D)

	Set a constraint G to the Phs. If a constraint already exists, it is
	overwritten. The feedthrough matrix G_D is optional.

"""
function set_constraint!(ph :: Phs, G, G_D)
	if isdefined(ph, :G) warn("Constraint was already defined! overwriting!"); end
	ph.G = G
	ph.G_D = G_D
end
function set_constraint!(ph ::Phs, G)
	set_constraint!(ph, G, zeros(size(G,2),size(G,2)))
end
function show(io ::IO, object :: Phs)
	println(io, "Port Hamiltonian system")
	println(io, "...", size(object.J,2), " energy variables")
	println(io, "...", size(object.B,2), " inputs/outputs")
	println(io, ".J"); println(io, ".B");
	println(io, ".D");
	if isdefined(object, :R) println(io, ".R"); end
	if isdefined(object, :Q) println(io, ".Q"); end
	if isdefined(object, :G) println(io, ".G (constrained system)"); end
	if isdefined(object, :Hamiltonian) println(io, ".Hamiltonian"); end
	if isdefined(object, :GradHam) println(io, ".GradHam"); end
	if isdefined(object, :disc_data) println(io, ".disc_data"); end
	if isdefined(object, :Bd) println(io, ".Bd (distributed ports)"); end
end

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
module pool
export Pool, mapWorker
using JLD
using eqjl

    mutable struct Pool
        tmp_dir::String
        rank_type::String
        clean_tmp::Bool
        step::Int64
        function Pool(tmp_dir, rank_type)
            if !(rank_type in ["workers", "leaders"])
                error("rank_type must be one of 'workers' or 'leaders'")
            end

            if !(ispath(tmp_dir))
                mkdir(tmp_dir)
            end

            new(tmp_dir, rank_type, false, 1)
        end
    end

    function __read_result(pool, fname)
        r = jldopen(fname, "r") do f_in
            read(f_in, "data")
        end
        return r
    end

    function mapWorker(pool, func, args)

        func_path = string("func_", pool.step, ".jld")
        func_f = joinpath(pool.tmp_dir, func_path)
        jldopen(func_f, "w") do f_out
            write(f_out, "data", func)
        end

        args_path = string("args_", pool.step, ".jld")
        args_f = joinpath(pool.tmp_dir, args_path)
        jldopen(args_f, "w") do f_out
            write(f_out, "data", args)
        end

        cmd = string("jlmap|", pool.step, "|", length(args), "|", pool.tmp_dir, "|", pool.rank_type)
        println(cmd)
        OUT_put(cmd)

        result = IN_get()
        pool.step += 1
        results = split(result, ";")
        lists = [__read_result(pool, x) for x in results]
        return [y for x in lists for y in x]
    end

end

max_stars_repo_name: emews/EQ-JL

max_issues_repo_name: emews/EQ-JL

max_forks_repo_name: emews/EQ-JL


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

function SCIPincludeConshdlrOr(scip)
    ccall((:SCIPincludeConshdlrOr, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsOr(scip, cons, name, resvar, nvars, vars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsOr, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Cint, Ptr{Ptr{SCIP_VAR}}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, resvar, nvars, vars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicOr(scip, cons, name, resvar, nvars, vars)
    ccall((:SCIPcreateConsBasicOr, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Cint, Ptr{Ptr{SCIP_VAR}}), scip, cons, name, resvar, nvars, vars)
end

function SCIPgetNVarsOr(scip, cons)
    ccall((:SCIPgetNVarsOr, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVarsOr(scip, cons)
    ccall((:SCIPgetVarsOr, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetResultantOr(scip, cons)
    ccall((:SCIPgetResultantOr, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end


text:
_basic_psa_opts(zoom,ps_dict) = Dict{Symbol,Any}(
    #:levels=>expandlevels(zoom.levels),
    :recompute_levels=>zoom.autolev,
    :proj_lev=>zoom.proj_lev,
    :scale_equal=>zoom.scale_equal,
    :real_matrix=>ps_dict[:Aisreal],
    :verbosity=>0)


text:
using BinDeps

# Libdl has been moved to a standard library module
using Compat.Libdl

@BinDeps.setup

# Add library dependency for direct method
osqp = library_dependency("osqp", aliases=["libosqp"])

# Current version
version = "0.3.0"

# Get current operating system
osqp_platform =
if Sys.is_linux()
    "linux"
elseif Sys.is_apple()
    "mac"
elseif Sys.is_windows()
    "windows"
else
    error("Platform not supported!")
end

# Provide binaries for each operating system
archive_name="osqp-$version-$osqp_platform$(Sys.WORD_SIZE)"

provides(Binaries, URI("https://dl.bintray.com/bstellato/generic/OSQP/$version/$archive_name.tar.gz"), [osqp], unpacked_dir="$archive_name/lib", os=:Darwin)
provides(Binaries, URI("https://dl.bintray.com/bstellato/generic/OSQP/$version/$archive_name.tar.gz"), [osqp], unpacked_dir="$archive_name/lib", os=:Windows)

# Build from sources on Linux
provides(Sources, URI("https://github.com/oxfordcontrol/osqp/archive/v$version.tar.gz"),
    [osqp], unpacked_dir="osqp-$version", os = :Linux)

# Define directories locations deps/usr and deps/src
prefix = joinpath(BinDeps.depsdir(osqp),"usr")
srcdir = joinpath(BinDeps.depsdir(osqp),"src","osqp-$version")
blddir = joinpath(srcdir, "build")

# Define library name
libname = "libosqp.$(dlext)"

provides(SimpleBuild,
    (@build_steps begin
    GetSources(osqp)
    CreateDirectory(joinpath(prefix, "lib"))
    @build_steps begin
        # ChangeDirectory(srcdir)
        CreateDirectory(blddir)
        @build_steps begin
        ChangeDirectory(blddir)
        FileRule(joinpath(prefix, "lib", libname), @build_steps begin
            `cmake -G "Unix Makefiles" -DUNITTESTS=OFF ..`
            `make osqp`
            `mv out/$libname $prefix/lib`
        end)
        end
    end
    end),
[osqp], os = :Linux)

@BinDeps.install Dict(:osqp => :osqp)


text:
using MarkdownLiteral
using Test

import HypertextLiteral
import MarkdownLiteral: @markdown

@testset "Not much" begin
    @test Symbol("@markdown") ∉ names(MarkdownLiteral)
    @test Symbol("@htl") ∈ names(HypertextLiteral)

    code_snippet = """
    xs = [1:10..., 20]
    map(xs) do x
    	f(x^2)
    end
    """

    plot(x, y) = HypertextLiteral.@htl("""

    """)

    result = @markdown("""
    # Hello!
    This is *Markdown* but **supercharged**!!
    <marquee style=$((color="purple", font_family="cursive"))>Inline HTML supported!</marquee>

    Here is a list, created using simple string interpolation:
    $((
    	"- item $i\n" for i in 1:3
    ))

    Another list, interpolated as HTML:

    $((
    	HypertextLiteral.@htl("
item $i
") for i in 1:3
    ))

    ![](https://media.giphy.com/media/JmUfwENE6i4Jxig27n/giphy.gif)

    Hello ``world``

    ```math
    \\sqrt{1+1}
    ```

    ## Intepolating a plotly plot
    It works!
    $(plot(1:10, rand(10)))

    ## Code block
    ```julia
    function f(x::Int64)
    	"hello \$(x)"
    end

    # we can interpolate into code blocks!
    $(code_snippet)
    ```
    """)

    htmle = repr(MIME"text/html"(), result)

    @test htmle isa String
    @test occursin("Hello", htmle)
end