text:
using CSV
using DataFrames
using Statistics
using CUDA
using StencilInterpolations

const viz = false

function random_particles(nxcell, x, y, z, dx, dy, dz, nx, ny, nz)
    ncells = (nx-1)*(ny-1)*(nz-1)
    px, py, pz = zeros(nxcell*ncells), zeros(nxcell*ncells), zeros(nxcell*ncells)
    for i in 1:nx-1
        for j in 1:ny-1, k in 1:nz-1
            # lowermost-left corner of the cell
            x0, y0, z0 = x[i], y[j], z[k]
            # cell index
            cell = i + (nx-1)*(j-1) + (nx-1)*(ny-1)*(k-1)
            for l in 1:nxcell
                px[(cell-1)*nxcell + l] = rand()*dx + x0
                py[(cell-1)*nxcell + l] = rand()*dy + y0
                pz[(cell-1)*nxcell + l] = rand()*dz + z0
            end
        end
    end
    return px, py, pz
end

function main(nx, ny, nz, nxcell)

    # model domain
    lx = ly = lz = 1
    dx, dy, dz = lx/(nx-1), ly/(ny-1), lz/(nz-1)
    x = LinRange(0, lx, nx)
    y = LinRange(0, ly, ny)
    z = LinRange(0, lz, nz)

    # random particles
    px, py, pz = random_particles(nxcell, x, y, z, dx, dy, dz, nx, ny, nz)
    particle_coords = (px, py, pz)
    N = length(px)

    # field to interpolate
    F = [-sin(2*zi)*cos(3*π*xi) for xi in x, _ in y, zi in z]
    F0 = deepcopy(F)

    ## CPU -----------------------------------------------------------------------------------------

    # scattering operation (tri-linear interpolation)
    t_scatter_cpu = @elapsed Fp = scattering((x, y, z), F, particle_coords)

    # gathering operation (inverse distance weighting)
    t_gather_cpu = @elapsed gathering!(F, Fp, (x, y, z), particle_coords)

    # Compute error
    sol = [-sin(2*zi)*cos(3*π*xi) for (xi, zi) in zip(px, pz)]
    misfit_scatter = @.(log10(abs(Fp-sol)))
    misfit_gather = @.(log10(abs(F-F0)))

    ## CUDA -----------------------------------------------------------------------------------------

    Fpd = CUDA.zeros(Float64, N)
    Fd = CuArray(F)
    Fd0 = deepcopy(Fd)
    particle_coords_dev = CuArray.(particle_coords)

    # scattering operation (for now just bi-linear interpolation)
    t_scatter_cuda = @elapsed scattering!(Fpd, (x, y, z), Fd, particle_coords_dev)

    # gathering operation (inverse distance weighting)
    fill!(Fd, 0.0)
    t_gather_cuda = @elapsed gathering!(Fd, Fpd, (x, y, z), particle_coords_dev)

    # Compute error
    sol_gpu = CuArray(sol)
    misfit_scatter_cuda = @.(log10(abs(Fpd-sol_gpu)))
    misfit_gather_cuda = @.(log10(abs(Fd-Fd0)))

    println("Finished for Ω ∈ [0,1] × [0,1] × [0,1]; $(nx) × $(ny) × $(nz) nodes; $nxcell particles per cell or $(Float64(N)) particles")

    # Plots CPU
    if viz == true
        xx = vec([x for x in x, y in y, z in z])
        yy = vec([y for x in x, y in y, z in z])
        zz = vec([z for x in x, y in y, z in z])

        f = Figure(resolution=(2400, 1200), fontsize=20)
        a = Axis(f[1,1], aspect=1, title="Analytical")
        scatter!(a, xx, yy, zz, color=vec(F0), colormap=:batlow, markersize=10)
        xlims!(0, lx)
        ylims!(0, ly)
        ylims!(0, lz)
        hidexdecorations!(a)

        a = Axis(f[1, 2], aspect=1, title="Scattered")
        scatter!(a, px, py, pz, color=Fp, colormap=:batlow)
        xlims!(0, lx)
        ylims!(0, ly)
        ylims!(0, lz)
        hidedecorations!(a)

        a = Axis(f[1, 3], aspect=1, title="log10 error")
        s = scatter!(a, px, py, pz, color=misfit_scatter, colormap=:batlow)
        xlims!(0, lx)
        ylims!(0, ly)
        ylims!(0, lz)
        hidedecorations!(a)

        Colorbar(f[1,4], s, height = 400)

        a = Axis(f[1, 5], aspect=1, title="Gathered")
        scatter!(a, xx, yy, zz, color=vec(F), colormap=:batlow, markersize=10)
        xlims!(0, lx)
        ylims!(0, ly)
        ylims!(0, lz)
        hidedecorations!(a)

        a = Axis(f[1, 6], aspect=1, title="log10 error")
        s = scatter!(a, xx, yy, zz, color=misfit_gather[:], colormap=:batlow)
        xlims!(0, lx)
        ylims!(0, ly)
        ylims!(0, lz)
        hidedecorations!(a)

        Colorbar(f[1,7], s, height = 400)

        display(f)

        # Plots GPU
        a = Axis(f[2,1], aspect=1, title="CUDA")
        scatter!(a, xx, yy, zz, color=vec(F), colormap=:batlow, markersize=10)
        xlims!(0, lx)
        ylims!(0, ly)
        ylims!(0, lz)

        a = Axis(f[2, 2], aspect=1, title="CUDA")
        scatter!(a, px, py, color=Array(Fpd), colormap=:batlow)
        xlims!(0, lx)
        ylims!(0, ly)
        ylims!(0, lz)
        hideydecorations!(a)

        a = Axis(f[2, 3], aspect=1, title="CUDA")
        s = scatter!(a, px, py, color=Array(misfit_scatter_cuda), colormap=:batlow)
        xlims!(0, lx)
        ylims!(0, ly)
        ylims!(0, lz)
        hideydecorations!(a)
        Colorbar(f[2,4], s, height = 400)

        a = Axis(f[2, 5], aspect=1, title="CUDA")
        scatter!(a, xx, yy, zz, color=vec(Array(Fd)), colormap=:batlow, markersize=10)
        xlims!(0, lx)
        ylims!(0, ly)
        ylims!(0, lz)
        hideydecorations!(a)

        a = Axis(f[2, 6], aspect=1, title="CUDA")
        s = scatter!(a, xx, yy, zz, color=Array(misfit_gather_cuda[:]), colormap=:batlow)
        xlims!(0, lx)
        ylims!(0, ly)
        ylims!(0, lz)
        hideydecorations!(a)

        Colorbar(f[2,7], s, height = 400)

        display(f)

    end

    e1_cpu, e1_cuda = mean(Fp.-sol), mean(Fpd.-sol_gpu)
    return t_scatter_cpu, t_scatter_cuda, t_gather_cpu, t_gather_cuda,  e1_cpu, e1_cuda
end

function perf_test()

    df = DataFrame(
        threads=Int16[],
        nx=Int64[],
        ny=Int64[],
        nz=Int64[],
        nxcell=Float64[],
        t_scatter_cpu=Float64[],
        t_scatter_cuda=Float64[],
        t_gather_cpu=Float64[],
        t_gather_cuda=Float64[],
        error_cpu=Float64[],
        error_cuda=Float64[]
    )

    nx = ny = nz = 128
    for nxcell in (4, 15, 25, 50)
        out = main(nx, ny, nz, nxcell)
        push!(df, [Threads.nthreads() nx ny nz nxcell out...])
    end

    CSV.write("scatter_perf_$(Threads.nthreads())_nxcell_fma_3D.csv", df)

end

perf_test()


text:
function header(; title = "Julia 语言与数据科学初步", author = "战立侃")
    item = filter(i -> i[2] == split(basename(title), "#")[1], StatisticsWithJulia.NOTEBOOKS)
    title = isempty(item) ? title : item[1][1]
    HTML{String}("""

        $title

        $author

    """)
end

function footer()
    main_page = _linkname(Base.current_project(), "../index.jl", "")
    HTML{String}("""

 返回《 Julia 语言与数据科学初步 》课程主页。

    """)
end

ext: jl


text:
#! /bin/julia
using SequentialFit, Plots
N = 20 # Number of frames

gaussian(x, A, mu, sigma) = A * exp(-sign(sigma) * ((x - mu) / sigma)^2)
p_true = [0.7, -1.5, 0.4, 0.9, 1.2, 0.7]

function slowFunction(x)
    # sleep(2) # slow evaluation, disabled for illustration
    return +(
        gaussian(x, p_true[1:3]...), # Gaussian peak
        gaussian(x, p_true[4:6]...), # Gaussian peak
        0.05sin(15 * x), # noise
    )
end

@. model(x, p) = gaussian(x, p[1:3]...) + gaussian(x, p[4:6]...)
p_0 = [1.0, -1.0, 0.5, 1.0, 1.0, 0.5] # Initial guess of parameters. Has to be decent
#x_0 = [-2.0, 1.0] # Initial sample points

s = SeqFit(slowFunction, model, p_0)
err = zeros(N, length(p_0))

anim = @animate for i in 1:N
    update!(s)
    err[i, :] .= abs.(s.fit.param ./ p_true .- 1)
    plot(s; groundtruth=true, selection=true, xlim=(-5, 5), ylim=(-0.1, 1.0))
end

gif(anim, "doc/example.gif"; fps=1)

savefig(
    plot(err; legend=false, xlabel="Iteration", ylabel="Error", yscale=:log10),
    "doc/example_errors.png",
)


text:
struct RestrictedMeasure{F,M} <: AbstractMeasure
    f::F
    base::M
end

@inline function logdensity(d::RestrictedMeasure, x)
    d.f(x) || return -Inf
    return 0.0
end

function density(d::RestrictedMeasure, x)
    d.f(x) || return 0.0
    return 1.0
end

basemeasure(μ::RestrictedMeasure) = μ.base


text:
#=
Base code for processing the drums dataset.
=#
#######################
# Channeling
#######################
#map pitches to destination channels (INVALID for discarding)
const INVALID = 9
const TOCHANNEL = Dict{UInt8,UInt8}(
    0x24=>1, #"Kick",
    0x26=>2, #"Snare",
    0x28=>2, #"Snare Rimshot",
    0x1a=>3, #"Hihat Rim",
    0x2e=>3, #"Hihat Head",
    0x33=>4, #"Ride Head",
    0x35=>4, #"Ride Bell",
    0x27=>5, #Tom 4 Rimshot,
    0x29=>5, #Tom 4,
    0x2b=>5, #"Tom 3",
    0x2d=>5, #"Tom 2",
    0x2f=>5, #"Tom 2 Rimshot",
    0x30=>5, #"Tom 1",
    0x32=>5, #"Tom 1 Rimshot",
    0x3a=>5, #"Tom 3 Rimshot",
    0x25=>INVALID, #"Snare RimClick",
    0x2c=>INVALID, #"Hihat Foot Close",
    0x31=>INVALID, #"Cymbal 1",
    0x3b=>INVALID, #"Ride Rim",
    0x37=>INVALID, #"Cymbal 1",
    0x39=>INVALID, #"Cymbal 2",
    0x34=>INVALID, #"Cymbal 2",
    0x2a=>INVALID,
    0x16=>INVALID,
)

const NOTANALYZABLE = [x for x in keys(TOCHANNEL) if TOCHANNEL[x] == INVALID]

const CHANNELNAMES = Dict{UInt8,String}(
    1=>"Kick",
    2=>"Snare",
    3=>"Hihat",
    4=>"Ride",
    5=>"Toms",
    6=>"All in one",
    INVALID=>"Not analyzable"
)

function channel_direct!(notes::Notes, map = TOCHANNEL)
    for note in notes
        note.pitch = map[note.pitch]
    end
    return notes
end

"""
    channel_direct(notes::Notes, map = TOCHANNEL)
Channel `notes` into the mapping specified by `map`
(by changing their pitch).
"""
channel_direct(notes, map = TOCHANNEL) = channel_direct!(copy(notes), map)

#######################
# Cleanup
#######################
function drums_cleanup(drums;
    back = 200, forw = 200, cutoff_fc = 30, cutoff_vel = 20)

    Drums.rm_hihatfake!(drums;
    back = back, forw = forw, cutoff_fc = cutoff_fc, cutoff_vel = cutoff_vel,
    cut_pitches=[0x2e])

    remove = NOTANALYZABLE
    drums = removepitches(drums, remove)
    return drums
end

function remove_peaking_velocities!(drums)
    digital = Drums.DIGITAL
    peak(v) = v ∈ digital ? 159 : 127
    remove = Int[]
    for i in 1:length(drums)
        note = drums[i]
        if note.velocity ≥ peak(note.velocity)
            push!(remove, i)
        end
    end
    deleteat!(drums.notes, remove)
    return drums
end

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
function add_initial_sample!(nnpool, labels, oracle, nn_stability, initial_n, query_history, gamma_progress)
    outlier_sample = StatsBase.sample(findall(labels .== :outlier), trunc(Int, initial_n / 2.0), replace=false)
    inlier_sample = StatsBase.sample(findall(labels .== :inlier), trunc(Int, initial_n / 2.0), replace=false)

    for (i,o) in enumerate(inlier_sample ∪ outlier_sample)
        add_to_pools!(nnpool, o, ask_oracle(oracle, o), nn_stability)
        push!(query_history, i, o)
        push!(gamma_progress, i, 1.0)
    end
end

function search_gamma(data, labels, oracle, nn_stability, stability_window,
                      stability_threshold, max_iteration, initial_n, max_γ, seed, adjusted)

    Random.seed!(seed)
    gamma_history = History(Float64)
    stability_history = History(Float64)
    query_history = History(Int64)

    nnpool = NNPool(data, k=nn_stability)
    add_initial_sample!(nnpool, labels, oracle, nn_stability, initial_n, query_history, gamma_history)

    for i in initial_n+1:max_iteration
        gamma = find_gamma(nnpool, data, nn_stability; max_γ=max_γ, adjusted=adjusted)

        K = kernelmatrix(Val(:col), GaussianKernel(gamma), data)
        current_alignment = alignment(nnpool, data, K, nn_stability; adjusted=adjusted)

        query = next_query(nnpool, data, K, query_history.values, nn_stability; adjusted=adjusted)
        add_to_pools!(nnpool, query, ask_oracle(oracle, query), nn_stability)

        stability = 0.0
        if (i > stability_window + 3)
            # smoothing
            tmp = rollmean(gamma_history.values, stability_window)
            # take average lag-1 difference
            stability = mean(min_max_normalize(abs.(tmp[1:end-1] - tmp[2:end]))[end-2:end])
        end

        push!(query_history, i, query)
        push!(gamma_history, i, gamma)
        push!(stability_history, i, stability)
        @info i, gamma, stability
        if (stability > 0) && (stability < stability_threshold)
            @info "Gamma stable. Stability: $stability"
            break
        end
    end
    return nnpool, gamma_history, stability_history, query_history
end

function calculate_C_grid(data, labels, oracle, gamma, min_C, max_C, steps, solver)
    C_range = range(min_C, max_C, length=steps)
    model_range = VanillaSVDD[]
    for c in C_range
        @debug "Training with C=$c"
        m = VanillaSVDD(data)
        SVDD.initialize!(m, SimpleCombinedStrategy(FixedGammaStrategy(MLKernels.GaussianKernel(gamma)), FixedCStrategy(c)))
        set_adjust_K!(m, true)
        SVDD.fit!(m, solver)
        push!(model_range, m)
    end
    return (collect(C_range), model_range)
end

function calculate_qm(nnpool, data, labels, model_range, qm)
    idx_pools = (nnpool.pools .!= :U)
    class_range = map(m -> SVDD.classify.(SVDD.predict(m, data[:, idx_pools])), model_range)
    cms_pools = map(x -> ConfusionMatrix(x, labels[idx_pools]), class_range)
    return qm.(cms_pools)
end

function find_min_feasible_C(data, gamma, solver, c_max = 0.5, c_min = 0.0, iter=10)
    c = (c_max - c_min) / 2
    c_opt = c

    for i in 1:10
        m = VanillaSVDD(data)
        SVDD.initialize!(m, SimpleCombinedStrategy(FixedGammaStrategy(MLKernels.GaussianKernel(gamma)), FixedCStrategy(c)))
        set_adjust_K!(m, true)
        res = SVDD.fit!(m, solver)
        if res == JuMP.MathOptInterface.OPTIMAL
            c_opt = c
            c_max = c
            c = c_min + (c - c_min) / 2
        else
            c_min = c
            c = c + (c_max - c) / 2
        end
    end
    return c_opt
end

function find_max_feasible_C(data, gamma, solver, c_max = 1.0, c_min = 0.0, iter=10)
    c = (c_max - c_min) / 2
    c_opt = c

    for i in 1:10
        m = VanillaSVDD(data)
        SVDD.initialize!(m, SimpleCombinedStrategy(FixedGammaStrategy(MLKernels.GaussianKernel(gamma)), FixedCStrategy(c)))
        set_adjust_K!(m, true)
        SVDD.fit!(m, solver)

        if any(SVDD.classify.(SVDD.predict(m, data)) .== :outlier)
            c_min = c
            c = c + (c_max - c) / 2
        else
            c_opt = c
            c_max = c
            c = c_min + (c - c_min) / 2
        end
    end
    return c_opt
end

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
module GLRasterization

using ModernGL, GLFW
using ColorVectorSpace, StaticArrays, GeometryTypes
import GLAbstraction, Colors
using GLAbstraction: GLBuffer, Shader, compile_shader, Texture
using Colors: RGBA

import Transpiler
using Transpiler: emit_vertex_shader, emit_geometry_shader, emit_fragment_shader, glsl_gensym
import Transpiler: gli

using FieldTraits
using FieldTraits: @reactivecomposed, @field, @composed, @needs
using FieldTraits: Composable, ComposableLike, Field, UsageError, Fields, Links
using FieldTraits: ReactiveComposable, Partial
import FieldTraits: on, default
using ..Visualize
using Visualize: WindowEvents, Area, Window, Mouse, Keyboard, Visible, Focused
using Visualize: Debugging, Name, AbstractWindow, DroppedFiles, Renderlist, Open
using Visualize: Projection, NativeWindow, ProjectionView, Resolution, Camera, Color
using Visualize: EyePosition, View, Scene, SceneUniforms, NativeWindowEvents, AbstractVertex
using Visualize: LookAt, PerspectiveCamera, Translation, Rotation, Pan, Zoom, poll_actions

import Visualize: add!, IRect, show!, destroy!, swapbuffers!, Drawable, rasterizer, renderloop

include("glutils.jl")
include("glfw.jl")
include("gtk.jl")
include("rasterpipeline.jl")
include("lines.jl")

export UniformBuffer, VertexArray, GLFWWindow, GLRasterizer

end


text:
using RobustAndOptimalControl
using ControlSystems
using Test
using LinearAlgebra
using Random

execute_tests = [true, true, true, true, true, true, true, true, true]

"""
Tests for the public and private methods of the hInfSynthesis function. This
function utilizes the preexisting ControlSystems toolbox, and performs a
H-infinity synthesis using the dual Riccati equation approach. As such,
the synthesis is done in a set of steps.

(1) Re-writing the specifications on an extended state-space form.
(2) Verifying that the resulting extended state-space object satisfies a set of
    assumptions required for proceeding with the synthesis.
(3) A coordinate transform to enable the synthesis.
(4) Synthesis using the γ-iterations, checking if a solution to the H-infinity
    problem exists in each iteration and applying a bisection method.
(5) Re-transforming the system to the original coordinates
(6) Verification that the computed solution is correct.

In addition to these six ponts, the code also enables

(7) A bilinear discretization with an inverse operation to move from continuous
to discrete time, thereby enabling approximate discrete-time synthesis.
(8) Plotting functionality to visualize the H-infinity synthesis.
(9) Three examples which can be used to demonstrate the tool.

Many of the tests are quite intuitive, and all points (1)-(9) are tested
extensively with detailed comments for each test-set, and below is a list
of all the functions tested in this unit test

  hinfsynthesize
  hinfassumptions
  hinfsignals
  hinfpartition
  bilineard2c
  bilinearc2d
  _detectable
  _stabilizable
  _synthesizecontroller
  _assertrealandpsd
  _checkfeasibility
  _solvehamiltonianare
  _solvematrixequations
  _γIterations
  _transformp2pbar
  _scalematrix
  _input2ss
  _coordinateTtansformqr
  _coordinateTtansformsvd

"""

if execute_tests[1]
    @testset "(1) Specifications" begin
        """
        These tests make sure that the specifications are written on a suitable form
        for any combination of none/empty, static gain, SISO, and MIMO
        specificaitons in the weighting functions. Essentially, the tests are made
        to verify that hInf_partition(G, WS, WU, WT) works as expected.
        """

        @testset "Conversion of user input" begin

            @testset "Empty and nothing" begin

                # Check that conversion of nothing is OK
                A, B, C, D = RobustAndOptimalControl._input2ss(nothing)
                @test isa(A, Array{Float64,2})
                @test isa(B, Array{Float64,2})
                @test isa(C, Array{Float64,2})
                @test isa(D, Array{Float64,2})

                @test size(A) == (0, 0)
                @test size(B) == (0, 0)
                @test size(C) == (0, 0)
                @test size(D) == (0, 0)

                # Check that conversion of empty objects are OK
                A, B, C, D = RobustAndOptimalControl._input2ss([])
                @test isa(A, Array{Float64,2})
                @test isa(B, Array{Float64,2})
                @test isa(C, Array{Float64,2})
                @test isa(D, Array{Float64,2})

                @test size(A) == (0, 0)
                @test size(B) == (0, 0)
                @test size(C) == (0, 0)
                @test size(D) == (0, 0)
            end

            @testset "Static gains" begin
                # Fixture
                number = 2.0

                # Check that conversion of numbers are OK
                A, B, C, D = RobustAndOptimalControl._input2ss(number)
                @test isa(A, Array{Float64,2})
                @test isa(B, Array{Float64,2})
                @test isa(C, Array{Float64,2})
                @test isa(D, Array{Float64,2})

                @test D[1, 1] == number

                @test size(A) == (0, 0)
                @test size(B) == (0, 1)
                @test size(C) == (1, 0)
                @test size(D) == (1, 1)
            end

            @testset "LTI models" begin

                # Fixture
                Random.seed!(0)
                M = 3
                N = 2
                SISO_tf = tf([1, 0], [1, 0, 1])
                SISO_ss = ss(SISO_tf)
                MIMO_ss = ss(
                    rand(Float64, (M, M)),
                    rand(Float64, (M, N)),
                    rand(Float64, (N, M)),
                    rand(Float64, (N, N)),
                )
                MIMO_tf = tf(MIMO_ss)

                # check that conversion of SISO tf data is OK
                A, B, C, D = RobustAndOptimalControl._input2ss(SISO_tf)
                @test isa(A, Array{Float64,2})
                @test isa(B, Array{Float64,2})
                @test isa(C, Array{Float64,2})
                @test isa(D, Array{Float64,2})

                @test ss(A, B, C, D) == ss(SISO_tf)

                # check that conversion of SISO tf data is OK
                A, B, C, D = RobustAndOptimalControl._input2ss(SISO_ss)
                @test isa(A, Array{Float64,2})
                @test isa(B, Array{Float64,2})
                @test isa(C, Array{Float64,2})
                @test isa(D, Array{Float64,2})

                @test ss(A, B, C, D) == SISO_ss

                # check that conversion of MIMO tf data is OK
                A, B, C, D = RobustAndOptimalControl._input2ss(MIMO_tf)
                @test isa(A, Array{Float64,2})
                @test isa(B, Array{Float64,2})
                @test isa(C, Array{Float64,2})
                @test isa(D, Array{Float64,2})

                @test ss(A, B, C, D) == ss(MIMO_tf)

                # check that conversion of MIMO tf data is OK
                A, B, C, D = RobustAndOptimalControl._input2ss(MIMO_ss)
                @test isa(A, Array{Float64,2})
                @test isa(B, Array{Float64,2})
                @test isa(C, Array{Float64,2})
                @test isa(D, Array{Float64,2})

                @test ss(A, B, C, D) == MIMO_ss
            end
        end

        #
        # @testset "Dimensionality checks"
        #
        #   # Check sensitivity function weight - must have the same number of inputs
        #   # as outputs, can have an arbitrarily large statespace, and must have the
        #   # same number of inputs as there are outputs in the process which is to be
        #   # controlled.
        #
        #   # Fixture
        #   M = 5; N = 3; L=4;
        #   G = ss(rand(Float64, (M,M)), rand(Float64, (M,L)),
        #          rand(Float64, (N,M)), rand(Float64, (N,L)))
        #
        #   # Create a sensitivity weighting function with ii inputs and ii outputs
        #   for ii = 2:6
        #     for jj = 2:6
        #       for kk = 2:6
        #         # Create some randome weighting functions
        #         WS = ss(rand(Float64, (M,M)), rand(Float64, (M,ii)),
        #                 rand(Float64, (ii,M)), rand(Float64, (ii,ii)))
        #         WU = ss(rand(Float64, (M,M)), rand(Float64, (M,jj)),
        #                 rand(Float64, (jj,M)), rand(Float64, (jj,jj)))
        #         WT = ss(rand(Float64, (M,M)), rand(Float64, (M,kk)),
        #                 rand(Float64, (kk,M)), rand(Float64, (kk,kk)))
        #
        #         println([ii,jj,kk])
        #         # Chech that the specifications can be re-written is possible
        #         if ii == N && jj == L && kk == N
        #           @test isa(hInf_partition(G, WS, WU, WT), ControlSystems.ExtendedStateSpace)
        #         else
        #           @test_throws ErrorException hInf_partition(G, WS, WU, WT)
        #         end
        #       end
        #     end
        #   end
        # end
    end
end

if execute_tests[2]
    @testset "(2) Assumptions" begin
        """
        Tests the methods used to check that the assumptions are satisfied,
        incorporating two separate algorithms to determine detectability and
        observability, and a rectangular pseudoinverse to check the rank conditions.
        """

        @testset "Stabilizability check" begin
            """
            Test the check for stabilizability using the Hautus Lemma
            """
            ### Fixture
            Random.seed!(0)
            N = 10
            M = 5
            R = rand(Float64, (N, N))
            Q = eigvecs(R + R')
            L = rand(Float64, N) .- 0.5
            A = Q * Diagonal(L) * Q'
            E = eigvals(A)

            # Here we should return false if for any positive eigenvalue λ of A, we get
            # that rank(A-λI, B) < N. We define B as conlums which are linearly
            # dependent with columns in A-λI, and check that the system returns false
            # in every case where λ > 0 and never when λ≦0.
            for ii = 1:N
                for jj = 1:(N-M)
                    Ahat = A - E[ii] * Matrix{Float64}(I, N, N)
                    B = Ahat[:, jj:(jj+M)]
                    if E[ii] >= 0
                        @test !RobustAndOptimalControl._stabilizable(A, B)
                    else
                        @test RobustAndOptimalControl._stabilizable(A, B)
                    end
                end
            end

            # Check common input types and ensure that a method error is thrown when
            # not using two abstract matrices, but some other common type
            N = 10
            M = 5
            A = rand(Float64, (N, N))
            B = rand(Float64, (N, M))
            @test_throws MethodError RobustAndOptimalControl._stabilizable(A, nothing)
            @test_throws MethodError RobustAndOptimalControl._stabilizable(nothing, B)
            @test_throws MethodError RobustAndOptimalControl._stabilizable(A, [])
            @test_throws MethodError RobustAndOptimalControl._stabilizable([], B)
            @test_throws MethodError RobustAndOptimalControl._stabilizable(A, ss(1))
            @test_throws MethodError RobustAndOptimalControl._stabilizable(ss(1), B)
            @test_throws MethodError RobustAndOptimalControl._stabilizable(A, tf(1))
            @test_throws MethodError RobustAndOptimalControl._stabilizable(tf(1), B)
        end

        @testset "Detectability check" begin
            """
            Test the check for detectability using the Hautus Lemma
            """
            ### Fixture
            Random.seed!(0)
            N = 10
            M = 5
            R = rand(Float64, (N, N))
            Q = eigvecs(R + R')
            L = rand(Float64, N) .- 0.5
            A = Q * Diagonal(L) * Q'
            E = eigvals(A)

            # Here we should return false if for any positive eigenvalue λ of A, we get
            # that rank(A-λI; C) < N. We define C as rows which are linearly
            # dependent with rows in A-λI, and check that the system returns false
            # in every case where λ > 0 and never when λ≦0.
            for ii = 1:N
                for jj = 1:(N-M)
                    Ahat = A - E[ii] * Matrix{Float64}(I, N, N)
                    C = Ahat[jj:(jj+M), :]
                    if E[ii] >= 0
                        @test !RobustAndOptimalControl._detectable(A, C)
                    else
                        @test RobustAndOptimalControl._detectable(A, C)
                    end
                end
            end

            # Check common input types and ensure that a method error is thrown when
            # not using two abstract matrices, but some other common type
            N = 10
            M = 5
            A = rand(Float64, (M, M))
            C = rand(Float64, (N, M))
            @test_throws MethodError RobustAndOptimalControl._detectable(A, nothing)
            @test_throws MethodError RobustAndOptimalControl._detectable(nothing, C)
            @test_throws MethodError RobustAndOptimalControl._detectable(A, [])
            @test_throws MethodError RobustAndOptimalControl._detectable([], C)
            @test_throws MethodError RobustAndOptimalControl._detectable(A, ss(1))
            @test_throws MethodError RobustAndOptimalControl._detectable(ss(1), C)
            @test_throws MethodError RobustAndOptimalControl._detectable(A, tf(1))
            @test_throws MethodError RobustAndOptimalControl._detectable(tf(1), C)
        end

        # TODO: write tests using the above submethods directly in hInf_assumptions
    end
end

if execute_tests[3]
    @testset "(3) Coordinate transform" begin
        """
        Computes the coordinate transfomration to write the system on the form of
        assumption A4 in [1], i.e. take an extended statespace model, P, satisfying
        assumption A2, and transform is to a system P_bar, where the matrices
        D12_bar = [0,I] and D21_bar = [0;I].
        """

        @testset "Find transformation for full rank matrix with QR" begin
            """
            Test computation of the transformation in square and rectangular cases.
            """
            # Fixture
            tolerance = 1e-10

            for N = 1:5
                for M = 1:5
                    A = rand(Float64, (M, N))
                    I_mat = Matrix{Float64}(I, min(M, N), min(M, N))

                    if M == N
                        # Square case
                        A_bar_true = I_mat
                    elseif M > N
                        # Rectangular case (more rows than columns)
                        Z_mat = zeros(Float64, (max(M, N) - min(M, N), min(M, N)))
                        A_bar_true = [Z_mat; I_mat]
                    else
                        # Rectangular case (more columns than rows)
                        Z_mat = zeros(Float64, (min(M, N), max(M, N) - min(M, N)))
                        A_bar_true = [Z_mat I_mat]
                    end
                    Tl, Tr = RobustAndOptimalControl._coordinatetransformqr(A)
                    @test opnorm(Tl * A * Tr - A_bar_true) < tolerance
                end
            end
        end

        @testset "Find transformation for full rank matrix with SVD" begin
            """
            Test computation of the transformation in square and rectangular cases.
            """
            # Fixture
            tolerance = 1e-10

            for N = 1:5
                for M = 1:5
                    A = rand(Float64, (M, N))
                    I_mat = Matrix{Float64}(I, min(M, N), min(M, N))

                    if M == N
                        # Square case
                        A_bar_true = I_mat
                    elseif M > N
                        # Rectangular case (more rows than columns)
                        Z_mat = zeros(Float64, (max(M, N) - min(M, N), min(M, N)))
                        A_bar_true = [Z_mat; I_mat]
                    else
                        # Rectangular case (more columns than rows)
                        Z_mat = zeros(Float64, (min(M, N), max(M, N) - min(M, N)))
                        A_bar_true = [Z_mat I_mat]
                    end
                    Tl, Tr = RobustAndOptimalControl._coordinatetransformsvd(A)
                    @test opnorm(Tl * A * Tr - A_bar_true) < tolerance
                end
            end
        end

        @testset "Find transformation for arbitrary matrix inputs" begin
            """
            Test computation of the transformation in square and rectangular cases.
            """
            # Fixture
            test_matrices_full = [
                rand(Float64, (1, 1)),
                rand(Float64, (5, 5)),
                rand(Float64, (2, 5)),
                rand(Float64, (5, 2)),
            ]

            # Compute rank deficient matrices
            test_matrices_rank_deficient = test_matrices_full
            for ii = 1:length(test_matrices_full)
                U, S, V = svd(test_matrices_full[ii])
                S[1] = 0
                test_matrices_rank_deficient[ii] = U * Diagonal(S) * V'
            end

            # Make sure that an exception is raised if assumption A2 is violated,
            # any of the matrices which are to be decomposed are rank deficient
            for A in test_matrices_rank_deficient
                println(A)
                @test_throws ErrorException RobustAndOptimalControl._scalematrix(A)
            end

            # Various bad inputs
            @test_throws MethodError RobustAndOptimalControl._scalematrix(tf(1))
            @test_throws MethodError RobustAndOptimalControl._scalematrix(ss(1))
            @test_throws MethodError RobustAndOptimalControl._scalematrix([])
            @test_throws MethodError RobustAndOptimalControl._scalematrix(nothing)
            @test_throws MethodError RobustAndOptimalControl._scalematrix(1)

            # Check that errors are thrown if the matric has zise zero
            @test_throws ErrorException RobustAndOptimalControl._scalematrix(zeros(Float64, (0, 1)))

            # Various bad methods
            A = test_matrices_full[1]
            @test_throws ErrorException RobustAndOptimalControl._scalematrix(
                A;
                method = "bad method",
            )
            @test_throws ErrorException RobustAndOptimalControl._scalematrix(A; method = 1)
        end

        @testset "Test application of the transformation to an ESS object" begin
            """
            The wa
            """
        end
    end
end

if execute_tests[4]
    @testset "(4) Gamma iterations" begin
        """
        Tests the core methods of the γ-iteration bisection method, including
        the ARE hamiltonian solver, the eigenvalue solver, the feasibility checks.
        """

        @testset "Solution feasibility check" begin
            """
            Check that a solution to the dual riccati equations is correctly reported
            as being feasible if satisfying the conditions of positive definiteness on
            the X and Y solutions are met, and the spectral radius condition of X*Y is
            also met
            """
            # Fixture
            tolerance = 1e-10
            iteration = 1
            Random.seed!(0)
            N = 10
            M = 5
            R = rand(Float64, N, N)
            Q = eigvecs(R + R')
            ρX = 0.2
            ρY = 0.3
            LX = rand(Float64, N)
            LX = ρX * sort(LX / maximum(LX))
            LY = rand(Float64, N)
            LY = ρY * sort(LY / maximum(LY))
            Xinf = Q * Diagonal(LX) * Q'
            Yinf = Q * Diagonal(LY) * Q'
            γ = 1

            # Test that the fesibility is true if ρ(Xinf*Yinf) < γ^2, that is the
            # check should be true for any
            #
            #    γ = sqrt(ρX*ρY) + ϵ
            #
            # with any epsilon greater than or equal to zero.
            @test !RobustAndOptimalControl._checkfeasibility(
                Xinf,
                Yinf,
                sqrt(ρX * ρY) - tolerance,
                tolerance,
                iteration;
                verbose = false,
            )
            @test RobustAndOptimalControl._checkfeasibility(
                Xinf,
                Yinf,
                sqrt(ρX * ρY) + tolerance,
                tolerance,
                iteration;
                verbose = false,
            )

            # Test that errors are thrown if the matrix Xinf and Yinf are not PSD down
            # to the numerical tolerance.
            L = LX
            L[1] += -L[1] + 2 * tolerance
            Xpos = Q * Diagonal(L) * Q' # slightly positive eigenvalue
            L[1] += -L[1]
            Xzero = Q * Diagonal(LX) * Q'               # exactly one zero eigenvalue
            L[1] += -L[1] - 2 * tolerance
            Xneg = Q * Diagonal(L) * Q'   # slightly negative eigenvalue
            @test RobustAndOptimalControl._checkfeasibility(
                Xpos,
                Yinf,
                sqrt(ρX * ρY) + tolerance,
                tolerance,
                iteration;
                verbose = false,
            )
            @test RobustAndOptimalControl._checkfeasibility(
                Xzero,
                Yinf,
                sqrt(ρX * ρY) + tolerance,
                tolerance,
                iteration;
                verbose = false,
            )
            @test !RobustAndOptimalControl._checkfeasibility(
                Xneg,
                Yinf,
                sqrt(ρX * ρY) + tolerance,
                tolerance,
                iteration;
                verbose = false,
            )

            L = LY
            L[1] += -L[1] + 2 * tolerance
            Ypos = Q * Diagonal(L) * Q' # slightly positive eigenvalue
            L[1] += -L[1]
            Yzero = Q * Diagonal(L) * Q'               # exactly one zero eigenvalue
            L[1] += -L[1] - 2 * tolerance
            Yneg = Q * Diagonal(L) * Q'   # slightly negative eigenvalue
            @test RobustAndOptimalControl._checkfeasibility(
                Xinf,
                Ypos,
                sqrt(ρX * ρY) + tolerance,
                tolerance,
                iteration;
                verbose = false,
            )
            @test RobustAndOptimalControl._checkfeasibility(
                Xinf,
                Yzero,
                sqrt(ρX * ρY) + tolerance,
                tolerance,
                iteration;
                verbose = false,
            )
            @test !RobustAndOptimalControl._checkfeasibility(
                Xinf,
                Yneg,
                sqrt(ρX * ρY) + tolerance,
                tolerance,
                iteration;
                verbose = false,
            )
        end

        # TODO: Include a check to verify that the bisection works as intended.
        # TODO: Check to verify that the Hamiltonian Shur-solver is working
        # TODO: Check to verify that the Hamiltonian Eigenvalue-solver is working
    end
end

if execute_tests[7]
    @testset "(7) Bilinear discretization" begin
        """
        This tests the bilinear method of discretizing a continuous time StateSpace
        or ExtendedStateSpace object, moving from the Laplace-domain to the Z-domain.
        However, importantly, the method also enables approximating discrete-time
        systems with continuous-time state-space domain, moving from the Z-domain back
        to the laplace-domain. Since the L∞-norm is invariant over the bilinear
        discretization and "continuization", this effectively allows us to approximate
        discrete-time systems with a continuous-time equivalent, design an H∞-optimal
        controller for the continuous-time plant, and then re-discretizing this
        controller.

        The structure of these tests is to simply discretize a set of plants with
        various dimensions, transport them back into continuous-time domain, and then
        to discrete-time again. By then comparing the resulting two pairs of discrete-
        time and continuous-time systems in the frequency domain, we can test if the
        bilinear discretization operates as expected.
        """

        @testset "SS data" begin
            # Fixture for the SS-type data
            Random.seed!(0)
            N = [1, 5, 10]
            M = [1, 5, 10]
            H = [0.1, 0.01, 0.001]
            tolerance = 1e-7

            # Tests for the SS-type data
            for (ii, m) in enumerate(N)
                for (jj, n) in enumerate(M)
                    for (kk, h) in enumerate(H)
                        freq = [10^i for i in range(-6, stop = log10(pi / h), length = 101)]

                        AcTrue, BcTrue, CcTrue, DcTrue =
                            rand(m, m), rand(m, n), rand(n, m), rand(n, n)

                        ev = eigvals(AcTrue)
                        if !isempty(ev[real(ev).<0])
                            AcTrue =
                                AcTrue -
                                one(AcTrue) * 2 * maximum(abs.(real(ev[real(ev).<0])))
                        end

                        Gtrue = ss(AcTrue, BcTrue, CcTrue, DcTrue)
                        valA, fA = sigma(Gtrue, freq)
                        sysB = bilinearc2d(Gtrue, h)
                        valB, fB = sigma(sysB, freq)
                        sysC = bilineard2c(bilinearc2d(Gtrue, h))
                        valC, fC = sigma(sysC, freq)
                        sysD = bilinearc2d(bilineard2c(bilinearc2d(Gtrue, h)), h)
                        valD, fD = sigma(sysD, freq)

                        @test abs(maximum(svd(sysB.B).S) - maximum(svd(sysB.C).S)) <
                              tolerance
                        @test abs(maximum(svd(sysC.B).S) - maximum(svd(sysC.C).S)) <
                              tolerance
                        @test abs(maximum(svd(sysD.B).S) - maximum(svd(sysD.C).S)) <
                              tolerance

                        # Test that the C->D->C and D->C->D results in the same
                        @test norm(valA - valC, Inf) < tolerance
                        @test norm(valB - valD, Inf) < tolerance
                    end
                end
            end
        end

        @testset "Extended SS data" begin
            # Fixture for the extended SS-type data
            Random.seed!(0)
            N = [1, 4]
            P1 = [1, 4]
            P2 = [1, 4]
            M1 = [1, 4]
            M2 = [1, 4]
            H = [0.1, 0.01, 0.001]
            tolerance = 1e-5

            for (in1, n) in enumerate(N)
                for (im1, m1) in enumerate(M1)
                    for (im2, m2) in enumerate(M2)
                        for (ip1, p1) in enumerate(P1)
                            for (ip2, p2) in enumerate(P2)
                                for (ih, h) in enumerate(H)
                                    freq = [
                                        10^i for i in
                                        range(-4, stop = log10(pi / h), length = 101)
                                    ]

                                    A, B1, B2 = 0.01 * rand(n, n), rand(n, m1), rand(n, m2)
                                    C1, D11, D12 = rand(p1, n), rand(p1, m1), rand(p1, m2)
                                    C2, D21, D22 = rand(p2, n), rand(p2, m1), rand(p2, m2)

                                    ev = eigvals(A)
                                    if !isempty(ev[real(ev).>0])
                                        A =
                                            A -
                                            one(A) *
                                            2 *
                                            maximum(abs.(real(ev[real(ev).>0])))
                                    end

                                    EsysA = ss(A, B1, B2, C1, C2, D11, D12, D21, D22)

                                    valA, fA = sigma(
                                        ss(
                                            EsysA.A,
                                            [EsysA.B1 EsysA.B2],
                                            [EsysA.C1; EsysA.C2],
                                            [EsysA.D11 EsysA.D12; EsysA.D21 EsysA.D22],
                                        ),
                                        freq,
                                    )

                                    # To discrete time
                                    EsysB = bilinearc2d(EsysA, h)
                                    valB, fB = sigma(
                                        ss(
                                            EsysB.A,
                                            [EsysB.B1 EsysB.B2],
                                            [EsysB.C1; EsysB.C2],
                                            [EsysB.D11 EsysB.D12; EsysB.D21 EsysB.D22],
                                            h,
                                        ),
                                        freq,
                                    )

                                    # To continuous time
                                    EsysC = bilineard2c(bilinearc2d(EsysA, h))
                                    valC, fC = sigma(
                                        ss(
                                            EsysC.A,
                                            [EsysC.B1 EsysC.B2],
                                            [EsysC.C1; EsysC.C2],
                                            [EsysC.D11 EsysC.D12; EsysC.D21 EsysC.D22],
                                        ),
                                        freq,
                                    )

                                    # To discrete time
                                    EsysD =
                                        bilinearc2d(bilineard2c(bilinearc2d(EsysA, h)), h)
                                    valD, fD = sigma(
                                        ss(
                                            EsysD.A,
                                            [EsysD.B1 EsysD.B2],
                                            [EsysD.C1; EsysD.C2],
                                            [EsysD.D11 EsysD.D12; EsysD.D21 EsysD.D22],
                                            h,
                                        ),
                                        freq,
                                    )

                                    # Test that the C->D->C and D->C->D results in the same
                                    @test norm(valA - valC, Inf) < tolerance
                                    @test norm(valB - valD, Inf) < tolerance
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

if execute_tests[9]
    @testset "(9) Synthesis Examples" begin

        @testset "DC motor example" begin
            # Fixture
            tolerance = 1e-4

            # Make sure that the code runs
            @test_nowarn include("../examples/hinf_example_DC.jl")
            Ω = [10^i for i in range(-3, stop = 3, length = 201)]

            # Check that the optimal gain is correct
            @test abs(γ - 4.463211059) < tolerance

            # Check that the closed loop satisfies ||F_l(P(jω), C(jω)||_∞ < γ,  ∀ω ∈ Ω
            valPcl = sigma(Pcl, Ω)[1]
            @test all(valPcl .< (γ + tolerance))

            # Check that ||S(jω)/WS(jω)||_∞ < γ,  ∀ω ∈ Ω
            if isa(WS, LTISystem) || isa(WS, Number)
                valSWS = sigma(S * WS, Ω)[1]
                @test all(valSWS .< (γ + tolerance))
            end

            # Check that ||C(jω)S(jω)/W_U(jω)||_∞ < γ,  ∀ω ∈ Ω
            if isa(WU, LTISystem) || isa(WU, Number)
                valKSWU = sigma(CS * WU, Ω)[1]
                @test all(valKSWU .< (γ + tolerance))
            end

            # Check that ||T(jω)/W_T(jω)||_∞ < γ,  ∀ω ∈ Ω
            if isa(WT, LTISystem) || isa(WT, Number)
                valTWT = sigma(T * WT, Ω)[1]
                @test all(valTWT .< (γ + tolerance))
            end

            # C_test = let
            #     _A = [-1.99970880438311e-7 -3.28465977959578e-12 -1.46200072004469e-10; 5213.86459190004 -599.125731164477 -26175.7449027424; 0.000132853900209503 0.249985014202232 -0.000667020485484034]
            #     _B = [3.9999999999998; -5.48401518999795e-14; -9.03638143773611e-7]
            #     _C = [651.733073987504 -74.8757163955595 -3271.9681128428]
            #     _D = [0.0]
            #     ss(_A, _B, _C, _D)
            # end
            # @test C ≈ C_test

            # Pcl_test = let
            #     _A = [-2.0e-7 0.0 -22.4 0.0 0.0 0.0; 0.0 -0.12 0.0 5213.86459190003 -599.005731164476 -26175.7449027424; 0.0 0.25 0.0 0.0 0.0 0.0; 0.0 0.0 -22.3999999999989 -1.99970880438311e-7 -3.28465977959578e-12 -1.46200072004469e-10; 0.0 0.0 3.07104850639885e-13 5213.86459190004 -599.125731164477 -26175.7449027424; 0.0 0.0 5.06037360513222e-6 0.000132853900209503 0.249985014202232 -0.000667020485484034]
            #     _B = [4.0; 0.0; 0.0; 3.9999999999998; -5.48401518999795e-14; -9.03638143773611e-7]
            #     _C = [4.99999996666667 0.0 -3.73333333333333 0.0 0.0 0.0; 0.0 0.0 0.0 651.733073987504 -74.8757163955595 -3271.9681128428]
            #     _D = [0.666666666666667; 0.0]
            #     ss(_A, _B, _C, _D)
            # end
            # @test Pcl ≈ Pcl_test

            @test_nowarn include("../examples/hinf_example_DC_discrete.jl")
        end

        @testset "MIT open courseware example" begin

            # Fixture
            tolerance = 1e-4

            # Make sure that the code runs
            @test_nowarn include("../examples/hinf_example_MIT.jl")
            Ω = [10^i for i in range(-7, stop = 7, length = 201)]

            # Check that the optimal gain is correct
            @test abs(γ - 0.923430124918) < tolerance

            # Check that the closed loop satisfies ||F_l(P(jω), C(jω)||_∞ < γ,  ∀ω ∈ Ω
            valPcl = sigma(Pcl, Ω)[1]
            @test all(valPcl .< (γ + tolerance))

            # Check that ||S(jω)/WS(jω)||_∞ < γ,  ∀ω ∈ Ω
            if isa(WS, LTISystem) || isa(WS, Number)
                valSWS = sigma(S * WS, Ω)[1]
                @test all(valSWS .< (γ + tolerance))
            end

            # Check that ||C(jω)S(jω)/W_U(jω)||_∞ < γ,  ∀ω ∈ Ω
            if isa(WU, LTISystem) || isa(WU, Number)
                valKSWU = sigma(CS * WU, Ω)[1]
                @test all(valKSWU .< (γ + tolerance))
            end

            # Check that ||T(jω)/W_T(jω)||_∞ < γ,  ∀ω ∈ Ω
            if isa(WT, LTISystem) || isa(WT, Number)
                valTWT = sigma(T * WT, Ω)[1]
                @test all(valTWT .< (γ + tolerance))
            end
            @test_nowarn include("../examples/hinf_example_MIT_discrete.jl")
        end

        @testset "Quad tank example" begin

            # Fixture
            tolerance = 1e-2

            # Make sure that the code runs
            @test_nowarn include("../examples/hinf_example_tank.jl")
            Ω = [10^i for i in range(-7, stop = 7, length = 201)]

            # Check that the optimal gain is correct
            @test abs(γ - 0.9534467) < tolerance

            # Check that the closed loop satisfies ||F_l(P(jω), C(jω)||_∞ < γ,  ∀ω ∈ Ω
            valPcl = sigma(Pcl, Ω)[1]
            @test all(valPcl[:, 1] .< (γ + tolerance))

            # Check that ||S(jω)/WS(jω)||_∞ < γ,  ∀ω ∈ Ω
            if isa(WS, LTISystem) || isa(WS, Number)
                valSWS = sigma(S * WS, Ω)[1]
                @test all(valSWS[:, 1] .< (γ + tolerance))
            end

            # Check that ||C(jω)S(jω)/W_U(jω)||_∞ < γ,  ∀ω ∈ Ω
            if isa(WU, LTISystem) || isa(WU, Number)
                valKSWU = sigma(CS * WU, Ω)[1]
                @test all(valKSWU[:, 1] .< (γ + tolerance))
            end

            # Check that ||T(jω)/W_T(jω)||_∞ < γ,  ∀ω ∈ Ω
            if isa(WT, LTISystem) || isa(WT, Number)
                valTWT = sigma(T * WT, Ω)[1]
                @test all(valTWT[:, 1] .< (γ + tolerance))
            end
            @test_nowarn include("../examples/hinf_example_tank_discrete.jl")
        end
    end
end


text:
module BigWaterfall

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Big Waterfall (FG)" => Ahorn.EntityPlacement(
        Maple.BigWaterfall,
        "rectangle",
        Dict{String, Any}(
            "layer" => "FG"
        )
    ),
    "Big Waterfall (BG)" => Ahorn.EntityPlacement(
        Maple.BigWaterfall,
        "rectangle",
        Dict{String, Any}(
            "layer" => "BG"
        )
    ),
)

const fillColor = Ahorn.XNAColors.LightBlue .* 0.3
const surfaceColor = Ahorn.XNAColors.LightBlue .* 0.8

const waterSegmentLeftMatrix = [
    1 1 1 0 1 0;
    1 1 1 0 1 0;
    1 1 1 0 1 0;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 1 0 1;
    1 1 1 0 1 0;
    1 1 1 0 1 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
    1 1 0 1 0 0;
]

const waterSegmentLeft = Ahorn.matrixToSurface(
    waterSegmentLeftMatrix,
    [
        fillColor,
        surfaceColor
    ]
)

const waterSegmentRightMatrix = [
    0 1 0 1 1 1;
    0 1 0 1 1 1;
    0 1 0 1 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 0 1 0 1 1;
    0 1 0 1 1 1;
    0 1 0 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
    1 0 1 1 1 1;
]

const waterSegmentRight = Ahorn.matrixToSurface(
    waterSegmentRightMatrix,
    [
        fillColor,
        surfaceColor
    ]
)

Ahorn.minimumSize(entity::Maple.BigWaterfall) = 8, 8
Ahorn.resizable(entity::Maple.BigWaterfall) = true, true

Ahorn.selection(entity::Maple.BigWaterfall) = Ahorn.getEntityRectangle(entity)

Ahorn.editingOptions(entity::Maple.BigWaterfall) = Dict{String, Any}(
    "layer" => String["FG", "BG"]
)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.BigWaterfall, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 16))
    height = Int(get(entity.data, "height", 64))

    segmentHeightLeft, segmentWidthLeft = size(waterSegmentLeftMatrix)
    segmentHeightRight, segmentWidthRight = size(waterSegmentRightMatrix)

    Ahorn.Cairo.save(ctx)

    Ahorn.rectangle(ctx, 0, 0, width, height)
    Ahorn.clip(ctx)

    for i in 0:segmentHeightLeft:ceil(Int, height / segmentHeightLeft) * segmentHeightLeft
        Ahorn.drawImage(ctx, waterSegmentLeft, 0, i)
        Ahorn.drawImage(ctx, waterSegmentRight, width - segmentWidthRight, i)
    end

    # Drawing a rectangle normally doesn't guarantee that its the same color as above
    if height >= 0 && width >= segmentWidthLeft + segmentWidthRight
        fillRectangle = Ahorn.matrixToSurface(fill(0, (height, width - segmentWidthLeft - segmentWidthRight)), [fillColor])
        Ahorn.drawImage(ctx, fillRectangle, segmentWidthLeft, 0)
    end

    Ahorn.restore(ctx)
end

end


text:
module LazyArraysTests

using Test
using Gridap.Arrays
using FillArrays
using Gridap.Helpers

a = rand(3,2,4)
test_array(a,a)

a = rand(3,2,4)
a = CartesianIndices(a)
test_array(a,a)

a = rand(16,32)
a = CartesianIndices(a)
c = lazy_map(-,a)
test_array(c,-a)

a = rand(12)
c = lazy_map(-,a)
test_array(c,-a)

a = rand(12)
b = rand(12)
c = lazy_map(-,a,b)
test_array(c,a.-b)

c = lazy_map(-,Float64,a,b)
test_array(c,a.-b)

a = rand(0)
b = rand(0)
c = lazy_map(-,a,b)
test_array(c,a.-b)

a = fill(rand(Int,2,3),12)
b = fill(rand(Int,1,3),12)
c = map(array_cache,(a,b))
i = 1
v = map((ci,ai) -> getindex!(ci,ai,i),c,(a,b))
# v = getitems!(c,(a,b),i)
@test c == (nothing,nothing)
@test v == (a[i],b[i])

a = fill(rand(Int,2,3),12)
b = fill(rand(Int,1,3),12)
ai = testitem(a)
@test ai == a[1]
ai, bi = map(testitem,(a,b))
@test ai == a[1]
@test bi == b[1]

a = fill(rand(Int,2,3),0)
b = fill(1,0)
ai = testitem(a)
@test ai == Array{Int,2}(undef,0,0)
ai, bi = map(testitem,(a,b))
@test ai == Array{Int,2}(undef,0,0)
@test bi == zero(Int)

a = fill(+,10)
x = rand(10)
y = rand(10)
v = lazy_map(+,x,y)
r = [(xi+yi) for (xi,yi) in zip(x,y)]
test_array(v,r)
v = lazy_map(+,Float64,x,y)
test_array(v,r)

p = 4
p1 = [1,2]
p2 = [2,1]
p3 = [4,3]
p4 = [6,1]

x = [p1,p2,p3,p4]

fa(x) = 2*x
fb(x) = sqrt.(x)

aa = Fill(fa,4)
r = lazy_map(fa,x)
@test all([ r[i] ≈ 2*x[i] for i in 1:4])

bb = Fill(fb,4)
r = lazy_map(fb,x)
@test all([ r[i] ≈ sqrt.(x[i]) for i in 1:4])

aaop = lazy_map(Operation,aa)
cm = lazy_map(evaluate,aaop,bb)
r = lazy_map(evaluate,cm,x)
@test all([ r[i] ≈ 2*(sqrt.(x[i])) for i in 1:4])

kk = cm[1]
ckk = return_cache(kk,p)
evaluate!(ckk,kk,p)

aop = lazy_map(Operation(+),aa,bb)
lazy_map(evaluate,aa,x)+lazy_map(evaluate,bb,x)
lazy_map(evaluate,aop,x)
@test lazy_map(evaluate,aop,x) == lazy_map(evaluate,aa,x)+lazy_map(evaluate,bb,x)

# Broadcasting

a = fill(rand(2,3),12)
b = rand(12)
c = lazy_map(Broadcasting(-),a,b)
test_array(c,[ai.-bi for (ai,bi) in zip(a,b)])

a = fill(rand(2,3),0)
b = rand(0)
c = lazy_map(Broadcasting(-),a,b)
test_array(c,[ai.-bi for (ai,bi) in zip(a,b)])

a = fill(rand(2,3),12)
b = rand(12)
c = lazy_map(Broadcasting(-),a,b)
d = lazy_map(Broadcasting(+),a,c)
e = lazy_map(Broadcasting(*),d,c)
test_array(e,[((ai.-bi).+ai).*(ai.-bi) for (ai,bi) in zip(a,b)])

a = Fill(Broadcasting(+),10)
x = [rand(2,3) for i in 1:10]
y = [rand(1,3) for i in 1:10]
v = lazy_map(evaluate,a,x,y)
r = [(xi.+yi) for (xi,yi) in zip(x,y)]
test_array(v,r)

a = Fill(Broadcasting(+),10)
x = [rand(mod(i-1,3)+1,3) for i in 1:10]
y = [rand(1,3) for i in 1:10]
v = lazy_map(evaluate,a,x,y)
r = [(xi.+yi) for (xi,yi) in zip(x,y)]
test_array(v,r)

# Operations

x = fill(4,3,3)

ax = Fill(x,4)
aa = Fill(fa,4)
bb = Fill(fb,4)

aop = lazy_map(Operation(Broadcasting(+)),aa,bb)
aax = lazy_map(evaluate,aa,ax)
bbx = lazy_map(evaluate,bb,ax)
aopx = lazy_map(evaluate,aop,ax)
@test aopx == aax+bbx

aop = lazy_map(Operation(Broadcasting(*)),aa,bb)
aopx = lazy_map(evaluate,aop,ax)
@test aopx[1] == aax[1].*bbx[1]

# Allocations

x = fill(4,3,3)

ax = Fill(x,4)
aa = Fill(Operation(fa),4)
bb = Fill(fb,4)
cm = lazy_map(evaluate,aa,bb)
r = lazy_map(evaluate,cm,ax)
@test all([ r[i] ≈ 2*(sqrt.(ax[i])) for i in 1:4])

nn = 2
an = Fill(nn,4)
ap = Fill(Broadcasting(*),4)
cm = lazy_map(evaluate,ap,ax,an)
@test all([cm[i] == nn*ax[i] for i in 1:4])

c_cm = array_cache(cm)
@allocated getindex!(c_cm,cm,1)
nalloc = 0
for i in length(cm)
  global nalloc
  nalloc += @allocated getindex!(c_cm,cm,i)
end
# @test nalloc == 0

as = Fill(Broadcasting(sqrt),4)
cs = lazy_map(evaluate,as,ax)
@test all([cs[i] == sqrt.(ax[i]) for i in 1:4])

c_cs = array_cache(cs)
@allocated getindex!(c_cs,cs,1)
nalloc = 0
for i in length(cs)
  global nalloc
  nalloc += @allocated getindex!(c_cs,cs,i)
end
# @test nalloc == 0

asm = lazy_map(Operation,as)
ah = lazy_map(evaluate,asm,ap)
ch = lazy_map(evaluate,ah,ax,an)
@test all([ ch[i] ≈ sqrt.(nn*ax[i]) for i in 1:4])
c_ch = array_cache(ch)
@allocated getindex!(c_ch,ch,1)
nalloc = 0
for i in length(ch)
  global nalloc
  nalloc += @allocated getindex!(c_ch,ch,i)
end
# @test nalloc == 0

# Empty arrays

g(x) = 1.0*x
a = lazy_map(g,Int[])
@test test_array(a,Float64[])

myf(x) = sqrt(x-1)
Arrays.testargs(::typeof(myf),x) = one(x)
@test test_array(lazy_map(myf,Int[]),Float64[])

# Fill optimizations

l = 10
v = 3
a = lazy_map(g,Fill(v,l))
test_array(a,Fill(g(v),l))
@test isa(a,Fill)

a = lazy_map(g,Float64,Fill(v,l))
test_array(a,Fill(g(v),l))
@test isa(a,Fill)

# array_cache

a = [rand(3) for i in 1:l]
b = [rand(3) for i in 1:l]
c = lazy_map(Broadcasting(+),a,b)
d = lazy_map(Broadcasting(*),a,c)

cache = array_cache(d)
@test getindex!(cache,d,1) === getindex!(cache,d,2)

i = 1
cache = array_cache(d)
@test getindex!(cache,d,i) === getindex!(cache,d,i)

ci = CartesianIndex((1))
cache = array_cache(d)
@test getindex!(cache,d,ci) === getindex!(cache,d,ci)

# Shapes and indices

s = (3,4)
a = rand(s...)
b = rand(s...)
c = lazy_map(+,a,b)
@test size(c) == s
d = a.+b
test_array(c,d)
@test c[1] == d[1]
@test c[2,3] == d[2,3]
@test c[CartesianIndex(2,3)] == d[CartesianIndex(2,3)]
@test c[2,3,1,1,1] == d[2,3,1,1,1]
@test c[CartesianIndex(2,3),1,1] == d[CartesianIndex(2,3),1,1]

s = (3,4)
a = rand(s...)
b = rand(prod(s))
c = lazy_map(+,a,b)
d = map(+,a,b)
@test size(c) == size(d)
@test  c == d

a = fill(x->2*x,3,3)
b = rand(9)
c = lazy_map(evaluate,a,b)
@test size(c) == (9,)
@test ndims(c) == 1
test_array(c,2*b)

# Test the intermediate results caching mechanism

a = Arrays.ArrayWithCounter(fill(rand(2,3),12))
b = Arrays.ArrayWithCounter(rand(12))
c = lazy_map(Broadcasting(-),a,b)
d = lazy_map(Broadcasting(+),a,c)
e = lazy_map(Broadcasting(*),d,c)
r = [ (ai.-bi).*(ai.+(ai.-bi)) for (ai,bi) in zip(a,b)]
test_array(e,r)
cache = array_cache(e)
Arrays.resetcounter!(a)
Arrays.resetcounter!(b)
for i in 1:length(e)
  ei = getindex!(cache,e,i)
  ei = getindex!(cache,e,i)
  ei = getindex!(cache,e,i)
end
@test all(a.counter .== 2)
@test all(b.counter .== 1)

# using lazy_map and test_map at the same time

n = 4
k  = Fill(+,n)
a = fill(2,n)
b = fill(4,n)
c = a + b
@test all(lazy_map(test_map,c,k,a,b))
@test all(lazy_map((x...)->test_map(x...,cmp=≈),c,k,a,b))

end # module


text:
include("turbines.jl")

function get_turbine()
    coord = Coord(0.0,0.0,0.0)
    return nrel_5mw_turbine = Turbine(coord, 126.4, 90.0, 1/3, 2.5, 10.0, 0.8)
end

turbine = get_turbine()
println(turbine)

# NREL 5MW turbine)
# rotor diameter: 126.4 (m)
# Hub Height: 90.0 (m)
# cut_in_speed: 3.0  ( m/s)
# rated_power: 5000.0  (kW)
# generator_efficiency: 0.944


text:
function ex1()

	ex1=JuMP.Model()
	JuMP.@variable(ex1,x[1:2],start=1.0)
	JuMP.@NLobjective(ex1,Min,x[1]-x[2])
	JuMP.@constraint(ex1,1-x[2]>=0)
	ex1=MathOptNLPModel(ex1)#MathProgNLPModel(ex1)
	G=JuMP.Model()
	JuMP.@variable(G,x[1:2],start=1.0)
	JuMP.@constraint(G,x[1]>=0)
	JuMP.@NLobjective(G,Min,0.0)
	G=MathOptNLPModel(G)#MathProgNLPModel(G)
	H=JuMP.Model()
	JuMP.@variable(H,x[1:2],start=1.0)
	JuMP.@constraint(H,x[2]>=0)
	JuMP.@NLobjective(H,Min,0.0)
	H=MathOptNLPModel(H)#MathProgNLPModel(H)

 return MPCCNLPs(ex1,G = G,H = H)
end

function ex1bd()

	ex1bd=JuMP.Model()
	ux(i)=[Inf;1][i]
	JuMP.@variable(ex1bd,x[i=1:2], upper_bound=ux(i),start=1.0)
	JuMP.@NLobjective(ex1bd,Min,x[1]-x[2])
	ex1bd=MathOptNLPModel(ex1bd)#MathProgNLPModel(ex1bd)
	G=JuMP.Model()
	JuMP.@variable(G,x[1:2],start=1.0)
	JuMP.@constraint(G,x[1]>=0)
	JuMP.@NLobjective(G,Min,0.0)
	G=MathOptNLPModel(G) #MathProgNLPModel(G)
	H=JuMP.Model()
	JuMP.@variable(H,x[1:2],start=1.0)
	JuMP.@constraint(H,x[2]>=0)
	JuMP.@NLobjective(H,Min,0.0)
	H=MathOptNLPModel(H) #MathProgNLPModel(H)

 return MPCCNLPs(ex1bd,G = G,H = H)
end

function ex2()

	ex2=JuMP.Model()
	JuMP.@variable(ex2,x[1:2],start=-1.0)
	JuMP.@NLobjective(ex2,Min,0.5*((x[1]-1)^2+(x[2]-1)^2))
	ex2=MathOptNLPModel(ex2) #MathProgNLPModel(ex2)
	G=JuMP.Model()
	JuMP.@variable(G,x[1:2],start=1.0)
	JuMP.@constraint(G,x[1]>=0)
	JuMP.@NLobjective(G,Min,0.0)
	G=MathOptNLPModel(G) #MathProgNLPModel(G)
	H=JuMP.Model()
	JuMP.@variable(H,x[1:2],start=1.0)
	JuMP.@constraint(H,x[2]>=0)
	JuMP.@NLobjective(H,Min,0.0)
	H=MathOptNLPModel(H) #MathProgNLPModel(H)

 return MPCCNLPs(ex2,G = G,H = H)
end

# Exemple 3 :
# minimize x1+x2-x3
# s.t.     0<=x1 _|_ x2>=0
#          -4x1+x3<=0
#          -4x2+x3<=0
function ex3()

	ex3=JuMP.Model()
	JuMP.@variable(ex3,x[1:3],start=1.0)
	JuMP.@NLobjective(ex3,Min,x[1]+x[2]-x[3])
	JuMP.@NLconstraint(ex3,c1,-4*x[1]+x[3]<=0)
	JuMP.@NLconstraint(ex3,c2,-4*x[2]+x[3]<=0)
	ex3=MathOptNLPModel(ex3) #MathProgNLPModel(ex3)
	G=JuMP.Model()
	JuMP.@variable(G,x[1:3],start=1.0)
	JuMP.@constraint(G,x[1]>=0)
	JuMP.@NLobjective(G,Min,0.0)
	G=MathOptNLPModel(G) #MathProgNLPModel(G)
	H=JuMP.Model()
	JuMP.@variable(H,x[1:3],start=1.0)
	JuMP.@constraint(H,x[2]>=0)
	JuMP.@NLobjective(H,Min,0.0)
	H=MathOptNLPModel(H) #MathProgNLPModel(H)

 return MPCCNLPs(ex3,G = G,H = H)
end

#Bard1 (MacMPEC)
function bard1()

 	ex3=JuMP.Model()
	JuMP.@variable(ex3,x[1:5],start=1.0)
	JuMP.@NLobjective(ex3,Min,(x[1]-5)^2+(2*x[2]+1)^2)
	JuMP.@NLconstraint(ex3,2*(x[2]-1)-1.5*x[1]+x[3]-0.5*x[4]+x[5]==0)
	ex3=MathOptNLPModel(ex3) #MathProgNLPModel(ex3)
	G=JuMP.Model()
	JuMP.@variable(G,x[1:5],start=1.0)
	JuMP.@constraint(G,3*x[1]-x[2]-3>=0)
	JuMP.@constraint(G,-x[1]+0.5*x[2]+4>=0)
	JuMP.@constraint(G,-x[1]-x[2]+7>=0)
	JuMP.@NLobjective(G,Min,0.0)
	G=MathOptNLPModel(G) #MathProgNLPModel(G)
	H=JuMP.Model()
	JuMP.@variable(H,x[1:5],start=1.0)
	JuMP.@constraint(H,x[3]>=0)
	JuMP.@constraint(H,x[4]>=0)
	JuMP.@constraint(H,x[5]>=0)
	JuMP.@NLobjective(H,Min,0.0)
	H=MathOptNLPModel(H) #MathProgNLPModel(H)

 return MPCCNLPs(ex3,G = G,H = H)
end


text:
isplanar(tri::Triangle3D) = true
function isplanar(quad::Quadrilateral3D)
    # If the surface normals of two triangles, composed of vertices (1,2,3) and (1,3,4),
    # face the same direction, the quad is planar
    𝗻₁ = (quad[2] - quad[1]) × (quad[3] - quad[1])
    𝗻₂ = (quad[4] - quad[1]) × (quad[3] - quad[1])
    return norm(𝗻₁ ⋅ 𝗻₂) ≈ norm(𝗻₁)*norm(𝗻₂)
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
export iauUtctai
"""
Time scale transformation:  Coordinated Universal Time, UTC, to
International Atomic Time, TAI.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical.

Given:
   utc1,utc2  double   UTC as a 2-part quasi Julian Date (Notes 1-4)

Returned:
   tai1,tai2  double   TAI as a 2-part Julian Date (Note 5)

Returned (function value):
               int      status: +1 = dubious year (Note 3)
                                 0 = OK
                              -1 = unacceptable date

Notes:

   1. utc1+utc2 is quasi Julian Date (see Note 2), apportioned in any
   convenient way between the two arguments, for example where utc1
   is the Julian Day Number and utc2 is the fraction of a day.

   2. JD cannot unambiguously represent UTC during a leap second unless
   special measures are taken.  The convention in the present
   function is that the JD day represents UTC days whether the
   length is 86399, 86400 or 86401 SI seconds.  In the 1960-1972 era
   there were smaller jumps (in either direction) each time the
   linear UTC(TAI) expression was changed, and these "mini-leaps"
   are also included in the SOFA convention.

   3. The warning status "dubious year" flags UTCs that predate the
   introduction of the time scale or that are too far in the future
   to be trusted.  See iauDat for further details.

   4. The function iauDtf2d converts from calendar date and time of day
   into 2-part Julian Date, and in the case of UTC implements the
   leap-second-ambiguity convention described above.

   5. The returned TAI1,TAI2 are such that their sum is the TAI Julian
   Date.

Called:
   iauJd2cal    JD to Gregorian calendar
   iauDat       delta(AT) = TAI-UTC
   iauCal2jd    Gregorian calendar to JD

References:

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

   Explanatory Supplement to the Astronomical Almanac,
   P. Kenneth Seidelmann (ed), University Science Books (1992)

This revision:  2013 July 26

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauUtctai(double utc1, double utc2, double *tai1, double *tai2)

function iauUtctai(utc1::Real, utc2::Real)
   ref_tai1 = Ref{Float64}(0.0)
   ref_tai2 = Ref{Float64}(0.0)

   status = ccall((:iauUtctai, libsofa_c), Cint,
                   (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                   convert(Float64, utc1), convert(Float64, utc2),
                   ref_tai1, ref_tai2)

   return status, ref_tai1[], ref_tai2[]
end


text:
model = :CapitalAdjustModel
dipar = Dict(:nK => 40, :nz => 5, :γ => 0.5)
p_1_All = createmodel(model; dipar..., intdim = :All)
p_1_Separable = createmodel(model; dipar..., intdim = :Separable)
p_1_Separable_States = createmodel(model; dipar..., intdim = :Separable_States)
p_1_Separable_ExogStates = createmodel(model; dipar..., intdim = :Separable_ExogStates)

diopt = Dict(:concavity => false, :monotonicity => false, :rewardcall=>:pre)
sol_1_All = solve(p_1_All; diopt...)
sol_1_Separable = solve(p_1_Separable; diopt...)
sol_1_Separable_States = solve(p_1_Separable_States; diopt...)
sol_1_Separable_ExogStates = solve(p_1_Separable_ExogStates; diopt...)

model = :CapitalAdjustModel2
dipar = Dict(:nK => 10, :nN => 7, :nz => 3)
p_2_All = createmodel(model; dipar..., intdim = :All)
p_2_Sep = createmodel(model; dipar..., intdim = :Separable)
p_2_Sep_States = createmodel(model; dipar..., intdim = :Separable_States)
p_2_Sep_ExogStates = createmodel(model; dipar..., intdim = :Separable_ExogStates)

optdict = Dict(
	:monotonicity=>[false,false], :concavity=>[false,false],
	:rewardcall=>:jit)
sol_2_All = solve(p_2_All; optdict...)
sol_2_Separable = solve(p_2_Sep; optdict...)
sol_2_Separable_States = solve(p_2_Sep_States; optdict...)
sol_2_Separable_ExogStates = solve(p_2_Sep_ExogStates; optdict...)

@testset "intdims" begin

	@testset "compare solution methods" begin

        @test isapprox(sol_1_All, sol_1_Separable, rtol=mytol)
        @test isapprox(sol_1_Separable, sol_1_Separable_States, rtol=mytol)
        @test isapprox(sol_1_Separable_States, sol_1_Separable_ExogStates, rtol=mytol)

		@test isapprox(sol_2_All, sol_2_Separable, rtol=mytol)
		@test isapprox(sol_2_Separable, sol_2_Separable_States, rtol=mytol)
        @test isapprox(sol_2_Separable_States, sol_2_Separable_ExogStates,rtol=mytol)

    end

end # intdims


text:
###
### Bipartite graph utilities
###

"""
    find_augmenting_path(g::BipartiteGraph, eq, var_eq_matching, varwhitelist, vcolor=falses(ndsts(g)), ecolor=falses(nsrcs(g))) -> path_found::Bool

Try to find augmenting paths.
"""
function find_augmenting_path(g, eq, var_eq_matching, varwhitelist, vcolor=falses(ndsts(g)), ecolor=falses(nsrcs(g)))
    ecolor[eq] = true

    # if a `var` is unassigned and the edge `eq <=> var` exists
    for var in 𝑠neighbors(g, eq)
        if (varwhitelist === nothing || varwhitelist[var]) && var_eq_matching[var] === unassigned
            var_eq_matching[var] = eq
            return true
        end
    end

    # for every `var` such that edge `eq <=> var` exists and `var` is uncolored
    for var in 𝑠neighbors(g, eq)
        ((varwhitelist === nothing || varwhitelist[var]) && !vcolor[var]) || continue
        vcolor[var] = true
        if find_augmenting_path(g, var_eq_matching[var], var_eq_matching, varwhitelist, vcolor, ecolor)
            var_eq_matching[var] = eq
            return true
        end
    end
    return false
end

"""
    matching(s::Union{SystemStructure,BipartiteGraph}, varwhitelist=nothing, eqwhitelist=nothing) -> assign

Find equation-variable bipartite matching. `s.graph` is a bipartite graph.
"""
matching(s::SystemStructure, varwhitelist=nothing, eqwhitelist=nothing) = matching(s.graph, varwhitelist, eqwhitelist)
function matching(g::BipartiteGraph, varwhitelist=nothing, eqwhitelist=nothing)
    var_eq_matching = Matching(ndsts(g))
    for eq in 𝑠vertices(g)
        if eqwhitelist !== nothing
            eqwhitelist[eq] || continue
        end
        find_augmenting_path(g, eq, var_eq_matching, varwhitelist)
    end
    return var_eq_matching
end

function error_reporting(sys, bad_idxs, n_highest_vars, iseqs)
    io = IOBuffer()
    if iseqs
        error_title = "More equations than variables, here are the potential extra equation(s):\n"
        out_arr = equations(sys)[bad_idxs]
    else
        error_title = "More variables than equations, here are the potential extra variable(s):\n"
        out_arr = structure(sys).fullvars[bad_idxs]
    end

    Base.print_array(io, out_arr)
    msg = String(take!(io))
    neqs = length(equations(sys))
    if iseqs
        throw(ExtraEquationsSystemException(
            "The system is unbalanced. "
            * "There are $n_highest_vars highest order derivative variables "
            * "and $neqs equations.\n"
            * error_title
            * msg
        ))
    else
        throw(ExtraVariablesSystemException(
            "The system is unbalanced. "
            * "There are $n_highest_vars highest order derivative variables "
            * "and $neqs equations.\n"
            * error_title
            * msg
        ))
    end
end

###
### Structural check
###
function check_consistency(sys::AbstractSystem)
    s = structure(sys)
    @unpack graph, var_to_diff, fullvars = s
    n_highest_vars = count(v->length(outneighbors(s.var_to_diff, v)) == 0, vertices(s.var_to_diff))
    neqs = nsrcs(graph)
    is_balanced = n_highest_vars == neqs

    if neqs > 0 && !is_balanced
        varwhitelist = var_to_diff .== nothing
        var_eq_matching = matching(graph, varwhitelist) # not assigned
        # Just use `error_reporting` to do conditional
        iseqs = n_highest_vars < neqs
        if iseqs
            eq_var_matching = invview(complete(var_eq_matching)) # extra equations
            bad_idxs = findall(isnothing, @view eq_var_matching[1:nsrcs(graph)])
        else
            bad_idxs = findall(isequal(unassigned), var_eq_matching)
        end
        error_reporting(sys, bad_idxs, n_highest_vars, iseqs)
    end

    # This is defined to check if Pantelides algorithm terminates. For more
    # details, check the equation (15) of the original paper.
    extended_graph = (@set graph.fadjlist = Vector{Int}[graph.fadjlist; map(collect, edges(var_to_diff))])
    extended_var_eq_matching = matching(extended_graph)

    unassigned_var = []
    for (vj, eq) in enumerate(extended_var_eq_matching)
        if eq === unassigned
            push!(unassigned_var, fullvars[vj])
        end
    end

    if !isempty(unassigned_var) || !is_balanced
        io = IOBuffer()
        Base.print_array(io, unassigned_var)
        unassigned_var_str = String(take!(io))
        errmsg = "The system is structurally singular! " *
                 "Here are the problematic variables: \n" *
                 unassigned_var_str
        throw(InvalidSystemException(errmsg))
    end

    return nothing
end

###
### BLT ordering
###

"""
    find_scc(g::BipartiteGraph, assign=nothing)

Find strongly connected components of the equations defined by `g`. `assign`
gives the undirected bipartite graph a direction. When `assign === nothing`, we
assume that the ``i``-th variable is assigned to the ``i``-th equation.
"""
function find_scc(g::BipartiteGraph, assign=nothing)
    cmog = DiCMOBiGraph{false}(g, Matching(assign === nothing ? Base.OneTo(nsrcs(g)) : assign))
    sccs = Graphs.strongly_connected_components(cmog)
    foreach(sort!, sccs)
    return sccs
end

function sorted_incidence_matrix(sys, val=true; only_algeqs=false, only_algvars=false)
    sys = algebraic_equations_scc(sys)
    s = structure(sys)
    @unpack var_eq_matching, fullvars, scc, graph = s
    g = graph
    varsmap = zeros(Int, ndsts(graph))
    eqsmap = zeros(Int, nsrcs(graph))
    varidx = 0
    eqidx = 0
    for c in scc, eq in c
        var = invview(var_eq_matching)[eq]
        if var != 0
            eqsmap[eq] = (eqidx += 1)
            varsmap[var] = (varidx += 1)
        end
    end
    for i in diffvars_range(s)
        varsmap[i] = (varidx += 1)
    end
    for i in dervars_range(s)
        varsmap[i] = (varidx += 1)
    end
    for i in 1:nsrcs(graph)
        if eqsmap[i] == 0
            eqsmap[i] = (eqidx += 1)
        end
    end

    I = Int[]
    J = Int[]
    for eq in 𝑠vertices(g)
        only_algeqs && (isalgeq(s, eq) || continue)
        for var in 𝑠neighbors(g, eq)
            only_algvars && (isalgvar(s, var) || continue)
            i = eqsmap[eq]
            j = varsmap[var]
            (iszero(i) || iszero(j)) && continue
            push!(I, i)
            push!(J, j)
        end
    end
    #sparse(I, J, val, nsrcs(g), ndsts(g))
    sparse(I, J, val)
end

###
### Structural and symbolic utilities
###

function find_solvables!(sys)
    s = structure(sys)
    @unpack fullvars, graph, solvable_graph = s
    eqs = equations(sys)
    empty!(solvable_graph)
    for (i, eq) in enumerate(eqs)
        isdiffeq(eq) && continue
        term = value(eq.rhs - eq.lhs)
        for j in 𝑠neighbors(graph, i)
            isalgvar(s, j) || continue
            var = fullvars[j]
            isinput(var) && continue
            a, b, islinear = linear_expansion(term, var)
            a = unwrap(a)
            if islinear && (!(a isa Symbolic) && a isa Number && a != 0)
                add_edge!(solvable_graph, i, j)
            end
        end
    end
    s
end

# debugging use
function reordered_matrix(sys, partitions=structure(sys).partitions)
    s = structure(sys)
    @unpack graph = s
    eqs = equations(sys)
    nvars = ndsts(graph)
    I, J = Int[], Int[]
    ii = 0
    M = Int[]
    for partition in partitions
        append!(M, partition.v_solved)
        append!(M, partition.v_residual)
    end
    M = invperm(vcat(M, setdiff(1:nvars, M)))
    for partition in partitions
        for es in partition.e_solved
            isdiffeq(eqs[es]) && continue
            ii += 1
            js = [M[x] for x in 𝑠neighbors(graph, es) if isalgvar(s, x)]
            append!(I, fill(ii, length(js)))
            append!(J, js)
        end

        for er in partition.e_residual
            isdiffeq(eqs[er]) && continue
            ii += 1
            js = [M[x] for x in 𝑠neighbors(graph, er) if isalgvar(s, x)]
            append!(I, fill(ii, length(js)))
            append!(J, js)
        end
    end
    # only plot algebraic variables and equations
    sparse(I, J, true)
end

###
### Nonlinear equation(s) solving
###

@noinline nlsolve_failure(rc) = error("The nonlinear solver failed with the return code $rc.")

function numerical_nlsolve(f, u0, p)
    prob = NonlinearProblem{false}(f, u0, p)
    sol = solve(prob, NewtonRaphson())
    rc = sol.retcode
    rc === :DEFAULT || nlsolve_failure(rc)
    # TODO: robust initial guess, better debugging info, and residual check
    sol.u
end


text:
module TestExamples1DEulerGravity

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1d_dgsem")

@testset "Compressible Euler with self-gravity" begin
  @trixi_testset "elixir_eulergravity_convergence.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"),
      l2   = [0.0002170799126638106, 0.0002913792848717502, 0.0006112320856262327],
      linf = [0.0004977401033188222, 0.0013594223337776157, 0.002041891084400227])
  end
end

end # module


text:
## LICENSE
##   Copyright (c) 2018 GEOPHYSTECH LLC
##
##   Licensed under the Apache License, Version 2.0 (the "License");
##   you may not use this file except in compliance with the License.
##   You may obtain a copy of the License at
##
##       http://www.apache.org/licenses/LICENSE-2.0
##
##   Unless required by applicable law or agreed to in writing, software
##   distributed under the License is distributed on an "AS IS" BASIS,
##   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##   See the License for the specific language governing permissions and
##   limitations under the License.
##
## initial release by Andrey Stepnov, email: a.stepnov@geophsytech.ru
module GroundMotion
  # This source contains a module settings. Please find types and functions in nested source files (see include).
  # The distance() function is used, see https://github.com/JuliaGeo/Geodesy.jl for further read
  using
    DelimitedFiles,
    Geodesy

  export
    # types
    Point_vs30,
    Point_pga_out,
    Earthquake,
    Params_as2008,
    Params_simidorikawa1999,
    Params_mf2013,
    # functions
    read_vs30_file,
    convert_to_point_vs30,
    read_vs30_dl_file,
    convert_to_float_array,
    gmpe_as2008,
    gmpe_simidorikawa1999,
    gmpe_mf2013

  # Types
  include("types.jl")
  # AS2008 GMPE modeling
  include("as2008.jl")
  # Si-Midorikawa 1999 GMPE modeling
  include("simidorikawa1999.jl")
  # Morikawa-Fujiwara 2013 GMPE modeling
  include("morikawafujiwara2013.jl")
  # read/write functions
  include("auxiliary.jl")
end


text:
module Examples

  using ..AlphaZero

  include("../games/tictactoe/main.jl")
  export Tictactoe

  include("../games/connect-four/main.jl")
  export ConnectFour

  include("../games/grid-world/main.jl")
  export GridWorld

  include("../games/fanorona/main.jl")
  export Fanorona

  const games = Dict(
    "grid-world" => GridWorld.GameSpec(),
    "tictactoe" => Tictactoe.GameSpec(),
    "fanorona" => Fanorona.GameSpec(),
    "connect-four" => ConnectFour.GameSpec())

  const experiments = Dict(
    "grid-world" => GridWorld.Training.experiment,
    "tictactoe" => Tictactoe.Training.experiment,
    "fanorona" => Fanorona.Training.experiment,
    "connect-four" => ConnectFour.Training.experiment)

end


text:
struct Simulation{DT <: Number, FT <: Function}
    Δt::DT
    f::FT
    x::OffsetArray{DT,3,Array{DT,3}}

    function Simulation(f::FT, x₀::AbstractArray{DT,2}, Δt::DT, nt::Int) where {DT, FT}
        x = OffsetArray(zeros(DT, size(x₀)..., nt+1), axes(x₀)..., 0:nt)
        x[:,:,0] .= x₀
        new{DT,FT}(Δt, f, x)
    end
end

ntimesteps(sim::Simulation) = lastindex(sim.x,3)
nsamples(sim::Simulation) = length(axes(sim.x,2))
ndims(sim::Simulation) = length(axes(sim.x,1))

eachtimestep(sim::Simulation) = axes(sim.x,3)[1:end]
eachsample(sim::Simulation) = axes(sim.x,2)

function run!(sim::Simulation{DT}) where {DT}
    ẋ = zeros(DT, ndims(sim))
    for n in eachtimestep(sim)
        for i in eachsample(sim)
            sim.f(ẋ, sim.x[:,i,n-1])
            sim.x[:,i,n] .= sim.x[:,i,n-1] .+ sim.Δt .* ẋ
        end
    end
    return sim.x
end

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
using URIParser
using Test

urls = ["hdfs://user:password@hdfshost:9000/root/folder/file.csv#frag",
    "https://user:password@httphost:9000/path1/path2;paramstring?q=a&p=r#frag",
    "https://user:password@httphost:9000/path1/path2?q=a&p=r#frag",
    "https://user:password@httphost:9000/path1/path2;paramstring#frag",
    "https://user:password@httphost:9000/path1/path2#frag",
    "file:///path/to/file/with%3fshould%3dwork%23fine",
    "ftp://ftp.is.co.za/rfc/rfc1808.txt", "http://www.ietf.org/rfc/rfc2396.txt",
    "ldap://[2001:db8::7]/c=GB?objectClass?one", "mailto:John.Doe@example.com",
    "news:comp.infosystems.www.servers.unix", "tel:+1-816-555-1212", "telnet://192.0.2.16:80/",
    "urn:oasis:names:specification:docbook:dtd:xml:4.1.2"]

failed = 0
for url in urls
    global failed
    u = URI(url)
    if !(string(u) == url) || !isvalid(u)
        failed += 1
        println("Test failed for ",url)
    end
end
if failed != 0
    exit(failed)
end

@test URI("hdfs://user:password@hdfshost:9000/root/folder/file.csv") == URI("hdfs","hdfshost",9000,"/root/folder/file.csv","","","user:password")
@test URI("google.com","/some/path") == URI("http://google.com:80/some/path")
g = URI("google.com","/some/path")
@test URI(g,port=160) == URI("http://google.com:160/some/path")

@test escape("abcdef αβ 1234-=~!@#\$()_+{}|[]a;") == "abcdef%20%CE%B1%CE%B2%201234-%3D~%21%40%23%24%28%29_%2B%7B%7D%7C%5B%5Da%3B"
@test unescape(escape("abcdef 1234-=~!@#\$()_+{}|[]a;")) == "abcdef 1234-=~!@#\$()_+{}|[]a;"
@test unescape(escape("👽")) == "👽"

@test escape_form("abcdef 1234-=~!@#\$()_+{}|[]a;") == "abcdef+1234-%3D~%21%40%23%24%28%29_%2B%7B%7D%7C%5B%5Da%3B"
@test unescape_form(escape_form("abcdef 1234-=~!@#\$()_+{}|[]a;")) == "abcdef 1234-=~!@#\$()_+{}|[]a;"

@test ("user", "password") == userinfo(URI("https://user:password@httphost:9000/path1/path2;paramstring?q=a&p=r#frag"))
@test URI("https://user:password@httphost:9000/path1/path2;paramstring?q=a&p=r") == defrag(URI("https://user:password@httphost:9000/path1/path2;paramstring?q=a&p=r#frag"))

@test ["dc","example","dc","com"] == path_params(URI("ldap://ldap.example.com/dc=example,dc=com"))[1]
@test ["servlet","jsessionid","OI24B9ASD7BSSD"] == path_params(URI("http://www.mysite.com/servlet;jsessionid=OI24B9ASD7BSSD"))[1]

@test Dict("q"=>"a","p"=>"r") == query_params(URI("https://httphost/path1/path2;paramstring?q=a&p=r#frag"))
@test Dict("q"=>"a","malformed"=>"") == query_params(URI("https://foo.net/?q=a&malformed"))

@test false == isvalid(URI("file:///path/to/file/with?should=work#fine"))
@test true == isvalid(URI("file:///path/to/file/with%3fshould%3dwork%23fine"))

@test URI("s3://bucket/key") == URI("s3","bucket",0,"/key")

@test sprint(show, URI("http://google.com")) == "URI(http://google.com/)"

# Error paths
# Non-ASCII characters
@test_throws ErrorException URI("http://🍕.com")
# Unexpected start of URL
@test_throws ErrorException URI(".google.com")
# Unexpected character after scheme
@test_throws ErrorException URI("ht!tp://google.com")

#  Issue #27
@test URIParser.escape("t est\n") == "t%20est%0A"

# Issue #2
@test sprint(show, MIME("text/html"), URI("http://google.com")) ==
    """http://google.com/"""


text:
using DashBootstrapComponents

inline_radioitems = html_div([
    dbc_label("Choose one"),
    dbc_radioitems(
        options = [
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
        ],
        value = 1,
        id = "radioitems-inline-input",
        inline = true,
    ),
]);

inline_checklist = html_div([
    dbc_label("Choose a bunch"),
    dbc_checklist(
        options = [
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
        ],
        value = [],
        id = "checklist-inline-input",
        inline = true,
    ),
]);

inline_switches = html_div([
    dbc_label("Toggle a bunch"),
    dbc_checklist(
        options = [
            Dict("label" => "Option 1", "value" => 1),
            Dict("label" => "Option 2", "value" => 2),
        ],
        value = [],
        id = "switches-inline-input",
        inline = true,
        switch = true,
    ),
]);

inline_inputs = dbc_form([inline_radioitems, inline_checklist, inline_switches]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
using Documenter, PyCMA

makedocs(
    modules = [PyCMA],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Johanni Brea",
    sitename = "PyCMA.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/jbrea/PyCMA.jl.git",
    push_preview = true
)


text:
# Autogenerated wrapper script for MPSolve_jll for i686-w64-mingw32-cxx11
export libmps

using GMP_jll
JLLWrappers.@generate_wrapper_header("MPSolve")
JLLWrappers.@declare_library_product(libmps, "libmps-3.dll")
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libmps,
        "bin\\libmps-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using GLWindow, ModernGL, GLAbstraction, Romeo, Color, Images, ImmutableArrays
function genrenderbuffer(format, dimensions, attachment)
    renderbuffer = GLuint[0]
    glGenRenderbuffers(1, renderbuffer)
    glBindRenderbuffer(GL_RENDERBUFFER, renderbuffer[1])
    glRenderbufferStorage(GL_RENDERBUFFER, format, dimensions...)
    glFramebufferRenderbuffer(GL_FRAMEBUFFER, attachment, GL_RENDERBUFFER, renderbuffer[1])
    renderbuffer[1]
end
begin
local const BOUNDINGBOX_FRAMEBUFFER = glGenFramebuffers()
glBindFramebuffer(GL_FRAMEBUFFER, BOUNDINGBOX_FRAMEBUFFER)
framebuffsize = [1,1]

minbuffer = Texture(Vec4, framebuffsize)
maxbuffer = Texture(Vec4, framebuffsize)
glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, minbuffer.id, 0)
glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT1, GL_TEXTURE_2D, maxbuffer.id, 0)
println(GLENUM(glCheckFramebufferStatus(GL_FRAMEBUFFER)).name)
glBindFramebuffer(GL_FRAMEBUFFER, Romeo.RENDER_FRAMEBUFFER)

const FLOAT32MIN = -999999999999999999999999f0

function boundingbox(renderobject)
    renderobject.boundingboxfunc(renderobject)
end

function boundingbox(renderobject)
    minbuffer[1:end,1:end] = Vec4[Vec4(FLOAT32MIN) for i=1:framebuffsize[1], j=1:framebuffsize[2]]
    maxbuffer[1:end,1:end] = Vec4[Vec4(FLOAT32MIN) for i=1:framebuffsize[1], j=1:framebuffsize[2]]
    glBindFramebuffer(GL_FRAMEBUFFER, BOUNDINGBOX_FRAMEBUFFER)
    glViewport(0,0, framebuffsize...)
    glClampColor(GL_CLAMP_VERTEX_COLOR, GL_FALSE)
    glClampColor(GL_CLAMP_READ_COLOR, GL_FALSE)
    glClampColor(GL_CLAMP_FRAGMENT_COLOR, GL_FALSE)
    glDrawBuffers(2, [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1])
    glClearColor(FLOAT32MIN, FLOAT32MIN, FLOAT32MIN, 0)
    glClear(GL_COLOR_BUFFER_BIT | GL_COLOR_BUFFER_BIT)

    glDisable(GL_DEPTH_TEST)
    glDisable(GL_CULL_FACE)
    glDisable(GL_ALPHA_TEST)
    glEnable(GL_BLEND)
    glBlendFunc(GL_ONE, GL_ONE)
    glBlendEquation(GL_MIN)

    program = renderobject.vertexarraybb.program
    glUseProgram(program.id)
    glBindVertexArray(renderobject.vertexarraybb.id)
    glDrawElements(GL_POINTS, renderobject.vertexarraybb.indexlength, GL_UNSIGNED_INT, GL_NONE)

    println(data(minbuffer))
    println(-data(maxbuffer))
    glBindFramebuffer(GL_FRAMEBUFFER, Romeo.RENDER_FRAMEBUFFER)
    glClampColor(GL_CLAMP_VERTEX_COLOR, GL_TRUE)
    glClampColor(GL_CLAMP_READ_COLOR, GL_TRUE)
    glClampColor(GL_CLAMP_FRAGMENT_COLOR, GL_TRUE)
end
end

shaderdir = Pkg.dir("Romeo", "src", "shader")
homedir = Pkg.dir("Romeo", "src", "boundingbox")

bbprogram = TemplateProgram(joinpath(homedir, "shader.vert"), joinpath(shaderdir,"boundingbox", "boundingbox.frag"),fragdatalocation=[(0, "minbuffer"),(1, "maxbuffer")])

indexes = indexbuffer(GLuint[i=0:49])
v = Vec3[Vec3(rand(-100f0:100f0), rand(-100f0:100f0), rand(-100f0:100f0))for i=1:50]
println(v)
function Base.max{T, NDIM}(x::Array{Vector3{T},NDIM})
    reduce(x) do v0, v1
        Vector3(v0[1] > v1[1] ? v0[1] : v1[1],
            v0[2] > v1[2] ? v0[2] : v1[2],
            v0[3] > v1[3] ? v0[3] : v1[3])
    end
end
function Base.min{T, NDIM}(x::Array{Vector3{T},NDIM})
    reduce(x) do v0, v1
        Vector3(v0[1] < v1[1] ? v0[1] : v1[1],
            v0[2] < v1[2] ? v0[2] : v1[2],
            v0[3] < v1[3] ? v0[3] : v1[3])
    end
end
const triangle = RenderObject(Dict(
        :vertex => GLBuffer(v),
        :index => indexes
    ),
    bbprogram, bbprogram)

postrender!(triangle, render, triangle.vertexarray)

boundingbox(triangle)
println(min(v))
println(max(v))
#=
glClearColor(0,0,0,0)
while Romeo.ROOT_SCREEN.inputs[:open].value
    glBindFramebuffer(GL_FRAMEBUFFER, 0)
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0,0, 1920, 1280)
    render(triangle)
    GLFW.SwapBuffers(Romeo.ROOT_SCREEN.nativewindow)
    GLFW.PollEvents()
    sleep(0.01)
end
=#
GLFW.Terminate()


text:
@testset "Expression Manipulation" begin
    a, b, c = LogicalSymbol.([:a, :b, :c])
    x, y, z = LogicalSymbol.([:x, :y, :z])

    @testset "Evaluation" begin
        st1 = a ∧ b
        @test evaluate(st1, Dict(a => true, b => true))
        @test !evaluate(st1, Dict(a => true, b => false))
        @test !evaluate(st1, Dict(a => false, b => false))

        st2 = (a ∧ b) → c
        @test evaluate(st2, Dict(a => true, b => true, c => true))
        @test evaluate(st2, Dict(a => false, b => true, c => true))
        @test evaluate(st2, Dict(a => false, b => false, c => true))
        @test !evaluate(st2, Dict(a => true, b => true, c => false))
    end

    @testset "Matching" begin
        st1 = a ∧ ¬b
        @test matches(st1, x)
        @test matches(st1, x ∧ y)
        @test !matches(st1, x ∨ y)
        @test !matches(st1, ¬x ∧ y)
        @test matches(st1, x ∧ ¬y)

        st2 = a → (b ∨ ¬c)
        @test matches(st2, x)
        @test matches(st2, x → y)
        @test matches(st2, x → (y ∨ z))
        @test matches(st2, x → (y ∨ ¬z))
        @test !matches(st2, ¬x)
        @test !matches(st2, ¬x → y)
        @test !matches(st2, x → ¬y)
        @test !matches(st2, x → (y ∧ z))

        st3 = (a ∧ b) ∨ (a ∧ b)
        @test matches(st3, x)
        @test matches(st3, x ∨ x)
        @test matches(st3, x ∨ y)
        @test matches(st3, (x ∧ y) ∨ (x ∧ y))
        @test !matches(st3, (x ∧ y) ∨ (y ∧ x))
        @test !matches(st3, ¬x)
    end

    # TODO: Rename this from replace to something more relevant to rule application
    @testset "Replacement" begin
        @test isequal(replace(a, x => ¬x), ¬a)
        @test isequal(replace(a ∧ b, x => ¬x), ¬(a ∧ b))
        @test isequal(replace(a ∧ a, (x ∧ x) => x), a)
        @test isequal(replace(a → b, (x → y) => (¬x ∨ y)), ¬a ∨ b)
        @test isequal(replace(a → ¬b, (x → y) => (¬x ∨ y)), ¬a ∨ ¬b)
        @test isequal(replace(¬(a ∨ b) ∨ ¬(a ∧ b), (¬x ∨ ¬y) => ¬(x ∧ y)), ¬((a ∨ b) ∧ (a ∧ b)))
    end
end


text:
##############################################################################
##
## We use SubDataFrame's to maintain a reference to a subset of a DataFrame
## without making copies.
##
##############################################################################

# a SubDataFrame is a lightweight wrapper around a DataFrame used most
# frequently in split/apply sorts of operations.

immutable SubDataFrame{T <: AbstractVector{Int}} <: AbstractDataFrame
    parent::DataFrame
    rows::T # maps from subdf row indexes to parent row indexes

    function SubDataFrame(parent::DataFrame, rows::T)
        if length(rows) > 0
            rmin, rmax = extrema(rows)
            if rmin < 1 || rmax > size(parent, 1)
                throw(BoundsError())
            end
        end
        new(parent, rows)
    end
end

function SubDataFrame{T <: AbstractVector{Int}}(parent::DataFrame, rows::T)
    return SubDataFrame{T}(parent, rows)
end

function SubDataFrame(parent::DataFrame, row::Integer)
    return SubDataFrame(parent, [row])
end

function SubDataFrame{S <: Integer}(parent::DataFrame,
                                    rows::AbstractVector{S})
    return sub(parent, int(rows))
end

function Base.getindex(df::SubDataFrame, colinds::Any)
    return df.parent[df.rows, colinds]
end

function Base.getindex(df::SubDataFrame, rowinds::Any, colinds::Any)
    return df.parent[df.rows[rowinds], colinds]
end

function Base.setindex!(df::SubDataFrame, val::Any, colinds::Any)
    df.parent[df.rows, colinds] = val
    return df
end

function Base.setindex!(df::SubDataFrame, val::Any, rowinds::Any, colinds::Any)
    df.parent[df.rows[rowinds], colinds] = val
    return df
end

function Base.size(df::SubDataFrame)
    return length(df.rows), size(df.parent, 2)
end

function Base.size(df::SubDataFrame, i::Integer)
    if i == 1
        return length(df.rows)
    elseif i == 2
        return size(df.parent, 2)
    else
        throw(ArgumentError("Invalid size index"))
    end
end

Base.names(df::SubDataFrame) = names(df.parent)

# TODO: Remove this?
index(df::SubDataFrame) = index(df.parent)

function Base.sub{S <: Real}(df::DataFrame, rowinds::AbstractVector{S})
    return SubDataFrame(df, rowinds)
end

function Base.sub{S <: Real}(sdf::SubDataFrame, rowinds::AbstractVector{S})
    return SubDataFrame(sdf.parent, sdf.rows[rowinds])
end

function Base.sub(df::DataFrame, rowinds::AbstractVector{Bool})
    return sub(df, getindex(SimpleIndex(size(df, 1)), rowinds))
end

function Base.sub(sdf::SubDataFrame, rowinds::AbstractVector{Bool})
    return sub(sdf, getindex(SimpleIndex(size(sdf, 1)), rowinds))
end

function Base.sub(adf::AbstractDataFrame, rowinds::Integer)
    return SubDataFrame(adf, Int[rowinds])
end

function Base.sub(adf::AbstractDataFrame, rowinds::Any)
    return sub(adf, getindex(SimpleIndex(size(adf, 1)), rowinds))
end

function Base.sub(adf::AbstractDataFrame, rowinds::Any, colinds::Any)
    return sub(adf[[colinds]], rowinds)
end

function Base.delete!(df::SubDataFrame, c::Any)
    return SubDataFrame(delete!(df.parent, c), df.rows)
end

without(df::SubDataFrame, c::Any) = SubDataFrame(without(df.parent, c), df.rows)

Base.similar(df::SubDataFrame, dims) =
    DataFrame([similar(df[x], dims) for x in names(df)], names(df))

nas(df::SubDataFrame, dims) =
    DataFrame([nas(df[x], dims) for x in names(df)], names(df))


text:
function get_max_growth(obj::MultistepFieldIntegrator)
    return jcall(obj, "getMaxGrowth", jdouble, ())
end

function get_max_step(obj::AdaptiveStepsizeFieldIntegrator)
    return jcall(obj, "getMaxStep", RealFieldElement, ())
end

function get_min_reduction(obj::MultistepFieldIntegrator)
    return jcall(obj, "getMinReduction", jdouble, ())
end

function get_min_step(obj::AdaptiveStepsizeFieldIntegrator)
    return jcall(obj, "getMinStep", RealFieldElement, ())
end

function get_n_steps(obj::MultistepFieldIntegrator)
    return jcall(obj, "getNSteps", jint, ())
end

function get_safety(obj::MultistepFieldIntegrator)
    return jcall(obj, "getSafety", jdouble, ())
end

function get_starter_integrator(obj::MultistepFieldIntegrator)
    return jcall(obj, "getStarterIntegrator", FieldODEIntegrator, ())
end

function initialize_step(obj::AdaptiveStepsizeFieldIntegrator, arg0::jboolean, arg1::jint, arg2::Vector{RealFieldElement}, arg3::FieldODEStateAndDerivative, arg4::FieldEquationsMapper)
    return jcall(obj, "initializeStep", RealFieldElement, (jboolean, jint, Vector{RealFieldElement}, FieldODEStateAndDerivative, FieldEquationsMapper), arg0, arg1, arg2, arg3, arg4)
end

function set_initial_step_size(obj::AdaptiveStepsizeFieldIntegrator, arg0::RealFieldElement)
    return jcall(obj, "setInitialStepSize", void, (RealFieldElement,), arg0)
end

function set_max_growth(obj::MultistepFieldIntegrator, arg0::jdouble)
    return jcall(obj, "setMaxGrowth", void, (jdouble,), arg0)
end

function set_min_reduction(obj::MultistepFieldIntegrator, arg0::jdouble)
    return jcall(obj, "setMinReduction", void, (jdouble,), arg0)
end

function set_safety(obj::MultistepFieldIntegrator, arg0::jdouble)
    return jcall(obj, "setSafety", void, (jdouble,), arg0)
end

function set_starter_integrator(obj::MultistepFieldIntegrator, arg0::FieldODEIntegrator)
    return jcall(obj, "setStarterIntegrator", void, (FieldODEIntegrator,), arg0)
end

function set_step_size_control(obj::AdaptiveStepsizeFieldIntegrator, arg0::jdouble, arg1::jdouble, arg2::Vector{jdouble}, arg3::Vector{jdouble})
    return jcall(obj, "setStepSizeControl", void, (jdouble, jdouble, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2, arg3)
end

function set_step_size_control(obj::AdaptiveStepsizeFieldIntegrator, arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return jcall(obj, "setStepSizeControl", void, (jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end