text:
using Test: @test, @testset
using Random: randn!
import Random
import Zygote
import Flux
import RestrictedBoltzmannMachines as RBMs

@testset "Binary-Binary RBM partition function (brute force)" begin
    rbm = RBMs.RBM(RBMs.Binary(3), RBMs.Binary(2), randn(3,2))
    randn!(rbm.visible.θ)
    randn!(rbm.hidden.θ)

    Z = 0.0
    for v1 in (0,1), v2 in (0,1), v3 in (0,1)
        for h1 in (0,1), h2 in (0,1)
            v = [v1;v2;v3;;]
            h = [h1;h2;;]
            Z += exp(-only(RBMs.energy(rbm, v, h)))
        end
    end
    @test log(Z) ≈ RBMs.log_partition(rbm)

    rbm.w .= 0
    gs = Zygote.gradient(rbm) do rbm
        RBMs.log_partition(rbm)
    end
    @test only(gs).visible.θ ≈ Flux.sigmoid.(rbm.visible.θ)
    @test only(gs).hidden.θ  ≈ Flux.sigmoid.(rbm.hidden.θ)
end

@testset "Spin-Binary RBM partition function (brute force)" begin
    rbm = RBMs.RBM(RBMs.Spin(3), RBMs.Binary(2), randn(3,2))
    randn!(rbm.visible.θ)
    randn!(rbm.hidden.θ)

    Z = 0.0
    for v1 in (-1,1), v2 in (-1,1), v3 in (-1,1)
        for h1 in (0,1), h2 in (0,1)
            v = [v1;v2;v3;;]
            h = [h1;h2;;]
            Z += exp(-only(RBMs.energy(rbm, v, h)))
        end
    end
    @test log(Z) ≈ RBMs.log_partition(rbm)

    rbm.w .= 0
    gs = Zygote.gradient(rbm) do rbm
        RBMs.log_partition(rbm)
    end
    @test only(gs).visible.θ ≈ tanh.(rbm.visible.θ)
    @test only(gs).hidden.θ ≈ Flux.sigmoid.(rbm.hidden.θ)
end


text:
"""
    gmtinfo(cmd0::String="", arg1=[]; kwargs...)

Reads files and finds the extreme values in each of the columns.

Full option list at [`gmtinfo`](http://gmt.soest.hawaii.edu/doc/latest/gmtinfo.html)

Parameters
----------

- **A** : -- Str --

    Specify how the range should be reported.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/gmtinfo.html#a)
- **C** : **per_column** : -- Bool or [] --

    Report the min/max values per column in separate columns [Default uses <min/max> format].
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/gmtinfo.html#c)
- **D** : **center** : -- Bool or [] --

    Modifies results obtained by -I by shifting the region to better align with the center of the data.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/gmtinfo.html#d)
- **E** : **get_record** : -- Str or [] --

    Returns the record whose column col contains the minimum (l) or maximum (h) value.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/gmtinfo.html#e)
- **F** : **counts** : -- Str or [] --

    Returns the counts of various records depending on the appended mode.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/gmtinfo.html#f)
- **I** : **report_region** : -- Number or Str --

    Report the min/max of the first n columns to the nearest multiple of the provided increments
    and output results in the form -Rw/e/s/n
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/gmtinfo.html#i)
- **L** : **common_limits** : -- Bool or [] --

    Determines common limits across tables or segments.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/gmtinfo.html#l)
- **S** : **for_error_bars** : -- Str or [] --

    Add extra space for error bars. Useful together with I option and when later plotting with psxy E.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/gmtinfo.html#s)
- **T** : **nearest_multiple** : -- Number or Str --

    Report the min/max of the first (0’th) column to the nearest multiple of dz and output this as
    the string -Tzmin/zmax/dz.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/gmtinfo.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function gmtinfo(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && !isa(arg1, GMTdataset) && return monolitic("gmtinfo", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	cmd = parse_V_params("", d)
	cmd, opt_bi = parse_bi(cmd, d)
	cmd, opt_di = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, opt_i = parse_i(cmd, d)
	cmd, = parse_o(cmd, d)
	cmd, = parse_r(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'A', d, [:A])
	cmd = add_opt(cmd, 'C', d, [:C :per_column])
	cmd = add_opt(cmd, 'D', d, [:D :center])
	cmd = add_opt(cmd, 'E', d, [:E :get_record])
	cmd = add_opt(cmd, 'F', d, [:F :counts])
	cmd = add_opt(cmd, 'I', d, [:I :report_region])
	cmd = add_opt(cmd, 'L', d, [:L :common_limits])
	cmd = add_opt(cmd, 'S', d, [:S :for_error_bars])
	cmd = add_opt(cmd, 'T', d, [:T :nearest_multiple])

	# If file name sent in, read it.
	cmd, arg1, = read_data(d, cmd0, cmd, arg1, " ", opt_i, opt_bi, opt_di)

	(haskey(d, :Vd)) && println(@sprintf("\tgmtinfo %s", cmd))

	if (!isempty_(arg1))  return gmt("gmtinfo " * cmd, arg1)
	else                  return gmt("gmtinfo " * cmd)
	end
end

# ---------------------------------------------------------------------------------------------------
gmtinfo(arg1, cmd0::String=""; kw...) = gmtinfo(cmd0, arg1; kw...)

ext: jl


text:
function fac(x::BigInt)
    if x < 2
        return 1
    end
    return x * fac(x - 1)
end

print(fac(BigInt(100999)))

max_stars_repo_name: vagoff/shootout

max_issues_repo_name: vagoff/shootout

max_forks_repo_name: vagoff/shootout


text:
#-------------------------------------------------------------------------------
# File: vns_pcc.jl
# Description: This files contains all function that are used in VNS_PCC.
# Date: November 03, 2019
# Author: Jonathan Fontaine, Killian Fretaud, Rémi Garcia,
#         Boualem Lamraoui, Benoît Le Badezet, Benoit Loger
#-------------------------------------------------------------------------------

"""
    perturbation_VNS_PCC_exchange(solution::Solution, k::Int, instance::Instance)

Returns a new solution obtain by `k` random exchanges between cars. Exchanged
cars must have the same HPRC constraints.
"""
function perturbation_VNS_PCC_exchange!(solution::Solution, k::Int, instance::Instance)
    b0 = instance.nb_late_prec_day+1
    # Array that contain only key defining more than one car
    valid_key_HPRC = Array{Int,1}(undef, 0)
    for group in instance.same_HPRC
        if length(group.second) >= 2
            push!(valid_key_HPRC, group.first)
        end
    end

    for _ in 1:k
        same_HPRC_array = instance.same_HPRC[rand(valid_key_HPRC)]
        car_pos_a = rand(same_HPRC_array)
        car_pos_b = rand(same_HPRC_array)
        # Cannot be the same
        while car_pos_a == car_pos_b
            car_pos_b = rand(same_HPRC_array)
        end
        solution.sequence[car_pos_a], solution.sequence[car_pos_b] = solution.sequence[car_pos_b], solution.sequence[car_pos_a]
        if is_sequence_valid(solution.sequence, instance.nb_cars, instance)
            solution.sequence[car_pos_a], solution.sequence[car_pos_b] = solution.sequence[car_pos_b], solution.sequence[car_pos_a]
            move_exchange!(solution, car_pos_a, car_pos_b, instance)
        else
            solution.sequence[car_pos_a], solution.sequence[car_pos_b] = solution.sequence[car_pos_b], solution.sequence[car_pos_a]
        end
    end

    return solution
end

"""
    perturbation_VNS_PCC_insertion(solution::Solution, k::Int,
                                   critical_cars::Array{Int, 1},
                                   instance::Instance)

Delete k vehicles from the sequence and add them back in the sequence according
to a greedy criterion.
"""
function perturbation_VNS_PCC_insertion!(solution::Solution, k::Int,
                                         critical_cars::Array{Int, 1},
                                         instance::Instance)
    remove!(solution, instance, k, critical_cars)
    for i in 1:k
        greedy_add!(solution, instance, instance.nb_cars, 3, true)
    end
    @assert instance.nb_cars == solution.length
    return solution
end

"""
    perturbation_VNS_PCC(solution_init::Solution, p::Int, k::Int, instance::Instance)

Calls both perturbations and return a new solution. This function does not
modify the solution given in parameters.
"""
function perturbation_VNS_PCC(solution_init::Solution, p::Int, k::Int, instance::Instance)
    solution = deepcopy(solution_init)
    if p == 1
        perturbation_VNS_PCC_exchange!(solution, k, instance)
    else
        critical_cars = find_critical_cars(solution, instance, 2)
        k = min(k, length(critical_cars))
        perturbation_VNS_PCC_insertion!(solution, k, critical_cars, instance)
    end

    return solution
end

"""
    local_search_intensification_VNS_PCC_exchange!(solution::Solution, instance::Instance)

Optimizes the weighted sum of three objectives using `move_exchange!`.
"""
function local_search_intensification_VNS_PCC_exchange!(solution::Solution, instance::Instance, start_time::UInt)
    # useful variable
    b0 = instance.nb_late_prec_day+1
    n = instance.nb_cars

    improved = true
    sequence = copy(solution.sequence)
    while improved && TIME_LIMIT > (time_ns() - start_time) / 1.0e9
        improved = false
        critical_cars = find_critical_cars(solution, instance, 2)
        for index_car_a in critical_cars
            best_delta = 0
            best_positions = Array{Int, 1}()
            for index_car_b in b0:instance.nb_cars
                if index_car_a != index_car_b
                    sequence[index_car_a], sequence[index_car_b] = sequence[index_car_b], sequence[index_car_a]
                    if is_sequence_valid(sequence, instance.nb_cars, instance)
                        delta = weighted_sum(cost_move_exchange(solution, index_car_a, index_car_b, instance, 3))
                        if delta < best_delta
                            best_positions = Array{Int, 1}([index_car_b])
                            best_delta = delta
                        elseif delta == best_delta
                            push!(best_positions, index_car_b)
                        end
                    end
                    sequence[index_car_a], sequence[index_car_b] = sequence[index_car_b], sequence[index_car_a]
                end
            end
            if !isempty(best_positions)
                index_car_b = rand(best_positions)
                move_exchange!(solution, index_car_a, index_car_b, instance)
                sequence[index_car_a], sequence[index_car_b] = sequence[index_car_b], sequence[index_car_a]
                if best_delta < 0
                    improved = true
                end
            end
        end
    end

    return solution
end

"""
    local_search_intensification_VNS_PCC_insertion!(solution::Solution, instance::Instance)

Optimizes the weighted sum of three objectives using `move_insertion!`.
"""
function local_search_intensification_VNS_PCC_insertion!(solution::Solution, instance::Instance, start_time::UInt)
    # useful variable
    b0 = instance.nb_late_prec_day+1
    improved = true
    while improved && TIME_LIMIT > (time_ns() - start_time) / 1.0e9
        improved = false
        for index_car in b0:instance.nb_cars
            best_delta = 0
            best_positions = Array{Int, 1}()
            matrix_deltas = cost_move_insertion(solution, index_car, instance, 3)
            penalize_costs!(matrix_deltas, index_car, solution, instance)
            for position in b0:instance.nb_cars
                if position != index_car
                    delta = weighted_sum(matrix_deltas[position, :])
                    if delta < best_delta
                        best_positions = Array{Int, 1}([position])
                        best_delta = delta
                    elseif delta == best_delta
                        push!(best_positions, position)
                    end
                end
            end

            if !isempty(best_positions)
                index_insert = rand(best_positions)
                move_insertion!(solution, index_car, index_insert, instance)
                if best_delta < 0
                    improved = true
                end
            end
        end
    end

    return solution
end

"""
    intensification_VNS_PCC!(solution::Solution, instance::Instance)

Calls both intensification.
"""
function intensification_VNS_PCC!(solution::Solution, instance::Instance, start_time::UInt)
    local_search_intensification_VNS_PCC_insertion!(solution, instance, start_time)
    local_search_intensification_VNS_PCC_exchange!(solution, instance, start_time)
    return solution
end

"""

"""
function local_search_VNS_PCC!(solution::Solution, instance::Instance, start_time::UInt)
    b0 = instance.nb_late_prec_day+1

    improved = true
    sequence = copy(solution.sequence)
    while improved && TIME_LIMIT > (time_ns() - start_time) / 1.0e9
        improved = false
        for index_car_a in b0:instance.nb_cars
            best_delta = 0
            best_positions = Array{Int, 1}()
            hprc_value = instance.HPRC_keys[solution.sequence[index_car_a]]
            for index_car_b in instance.same_HPRC[hprc_value]
                if index_car_a != index_car_b
                    sequence[index_car_a], sequence[index_car_b] = sequence[index_car_b], sequence[index_car_a]
                    if is_sequence_valid(sequence, instance.nb_cars, instance)
                        delta = weighted_sum(cost_move_exchange(solution, index_car_a, index_car_b, instance, 3))
                        if delta < best_delta
                            best_positions = Array{Int, 1}([index_car_b])
                            best_delta = delta
                        elseif delta == best_delta
                            push!(best_positions, index_car_b)
                        end
                    end
                    sequence[index_car_a], sequence[index_car_b] = sequence[index_car_b], sequence[index_car_a]
                end
            end
            if !isempty(best_positions)
                index_car_b = rand(best_positions)
                move_exchange!(solution, index_car_a, index_car_b, instance)
                sequence[index_car_a], sequence[index_car_b] = sequence[index_car_b], sequence[index_car_a]
                if best_delta < 0
                    improved = true
                end
            end
        end
    end

    return solution
end

"""

"""
function VNS_PCC(solution_init::Solution, instance::Instance, start_time::UInt)
    solution = deepcopy(solution_init)
    p = 1
    k = VNS_PCC_MINMAX[p+1][1]
    costs_solution = cost(solution, instance, 3)
    cost_solution = weighted_sum(costs_solution)
    cost_HPRC_solution = costs_solution[1]
    while TIME_LIMIT > (time_ns() - start_time) / 1.0e9
        while (k <= VNS_PCC_MINMAX[p+1][2]) && (TIME_LIMIT > (time_ns() - start_time) / 1.0e9)
            solution_perturbation = perturbation_VNS_PCC(solution, p, k, instance)
            if cost_HPRC_solution < cost(solution_perturbation, instance, 1)[1]
                solution_perturbation = deepcopy(solution)
            end
            local_search_VNS_PCC!(solution_perturbation, instance, start_time)
            cost_solution_perturbation = weighted_sum(solution_perturbation, instance)
            if cost_solution_perturbation < cost_solution
                k = VNS_PCC_MINMAX[p+1][1]
            else
                k += 1
            end
            if cost_solution_perturbation <= cost_solution
                solution = deepcopy(solution_perturbation)
                costs_solution = cost(solution, instance, 3)
                cost_solution = copy(cost_solution_perturbation)
                cost_HPRC_solution = costs_solution[1]
            end
        end
        intensification_VNS_PCC!(solution, instance, start_time)
        costs_solution = cost(solution, instance, 3)
        cost_solution = weighted_sum(costs_solution)
        cost_HPRC_solution = costs_solution[1]
        p = 1 - p
        k = VNS_PCC_MINMAX[p+1][1]
    end

    return solution
end


text:
Optional{T} = Union{Nothing,T}

mutable struct Address
    street::String
    house_nr::String
    zip_code::String
    town::String
    iid::DbId
end

function Address()
    return Address("","","","",DbId())
end

mutable struct Employee
    name::String
    contact_person_first::String
    contact_person_second::String
    job_title::String
    salutation::String
    call_text::String
    done::Bool
    address_rel::Optional{Address}
    callback_date::Optional{DateTime}
    iid::DbId
end

max_issues_repo_name: Humans-of-Julia/TypeDBORM

max_forks_repo_name: Humans-of-Julia/TypeDBORM


text:
function RCCD{T}(ref::Fermi.HartreeFock.RHF.RHFWavefunction,alg::Fermi.CoupledCluster.DPD,int::Fermi.CoupledCluster.Conventional) where T <: AbstractFloat
    error("RCCD with cc_alg $alg and ints_type $int is not implemented. Sorry!")
end


text:
export doubledescription
# Naive implementation of the double description Algorithm
# See JuliaPolyhedra/ConvexHull.jl for an efficient implementation

polytypefor(::Type{T}) where {T <: Integer} = Rational{BigInt}
polytypefor(::Type{Float32}) = Float64
polytypefor(::Type{T}) where {T} = T

polyvectortype(a) = a
# TODO sparse halfspaces does not mean sparse points
polyvectortype(::Type{<:SparseVector{T}}) where T = Vector{T}

dualtype(RepT::Type{<:Representation}) = dualtype(RepT, polyvectortype(vectortype(RepT)))
function dualfullspace(rep::Representation, d::FullDim, ::Type{T}) where T
    dualfullspace(rep, d, T, polyvectortype(similar_type(vectortype(rep), d, T)))
end
function dualfullspace(rep::Representation{N, T}) where {N, T}
    dualfullspace(rep, FullDim{N}(), polytypefor(T))
end

"""
    doubledescription(h::HRepresentation)

Computes the V-representation of the polyhedron represented by `h` using the Double-Description algorithm [1, 2].

    doubledescription(V::VRepresentation)

Computes the H-representation of the polyhedron represented by `v` using the Double-Description algorithm [1, 2].

[1] Motzkin, T. S., Raiffa, H., Thompson, G. L. and Thrall, R. M.
The double description method
*Contribution to the Theory of Games*, *Princeton University Press*, **1953**

[2] Fukuda, K. and Prodon, A.
Double description method revisited
*Combinatorics and computer science*, *Springer*, **1996**, 91-111
"""
function doubledescription end

function doubledescription(h::HRepresentation)
    # The redundancy of V-elements are measured using
    # the number of hyperplane they are in. If there are
    # redundant hyperplanes, it does not matter since no point
    # will be inside them but if there are duplicates it is a problem

    # FIXME Note that removevredundancy, uses `h` which contains all hyperplanes and halfspaces
    # already taken into account but also all the other ones. We should check that this
    # is the right idea.
    h = removeduplicates(h)
    v = dualfullspace(h)
    checkvconsistency(v)
    for hp in hyperplanes(h)
        #v = removeduplicates(v ∩ hel)
        # removeduplicates is cheaper than removevredundancy since the latter
        # needs to go through all the hrep element
        v = removevredundancy(removeduplicates(v ∩ hp), h)
        #v = removevredundancy(v ∩ hel, h)
    end
    for hs in halfspaces(h)
        #v = removeduplicates(v ∩ hel)
        # removeduplicates is cheaper than removevredundancy since the latter
        # needs to go through all the hrep element
        v = removevredundancy(removeduplicates(v ∩ hs), h)
        #v = removevredundancy(v ∩ hel, h)
    end
    v
end

function doubledescription(v::VRepresentation{N, T}) where {N, T}
    checkvconsistency(v)
    lv = LiftedVRepresentation{N, T, Matrix{T}}(v)
    R = -lv.R
    vl = doubledescription(MixedMatHRep{N+1, T}(R, zeros(T, size(R, 1)), lv.linset))
    LiftedHRepresentation{N, T}(vl.R, vl.Rlinset)
end


text:
result_sum = 0
result_prod = 1

for k in 1:10
    global result_sum += 2*k
    global result_prod *= 2*k
end

println("sum=", result_sum)
println("prod=", result_prod)


text:
using DataFrames
using Arrow
using CSV

@testset "dtable" begin
    @testset "constructors - Tables.jl compatibility (NamedTuple)" begin
        size = 1_000
        nt = (a = rand(Int, size), b = rand(Int, size))

        dt = DTable(nt, 100)
        @test fetch(dt) == nt
        @test fetch(dt, NamedTuple) == nt
        @test tabletype(dt) == NamedTuple

        # empty dtable case
        dt = DTable((a = [], b = []), 10)
        @test fetch(dt) == NamedTuple()
    end

    @testset "constructors - Tables.jl compatibility (DataFrames)" begin
        size = 1_000
        nt = (a = rand(Int, size), b = rand(Int, size))
        df = DataFrame(nt)

        dt = DTable(df, 100)
        @test fetch(dt) == df
        @test fetch(dt, DataFrame) == df
        @test tabletype(dt) == DataFrame

        dt = DTable(nt, 100; tabletype=DataFrame)
        @test fetch(dt) == df
        @test fetch(dt, DataFrame) == df
        @test tabletype(dt) == DataFrame

        dt = DTable(df, 100; tabletype=NamedTuple)
        @test fetch(dt) == nt
        @test fetch(dt, DataFrame) == df
        @test tabletype(dt) == NamedTuple
    end

    @testset "constructors - Tables.jl compatibility (Arrow)" begin
        size = 1_000
        nt = (a = rand(Int, size), b = rand(Int, size))
        io = IOBuffer()
        Arrow.write(io, nt)
        arr = Arrow.Table(take!(io))

        dt = DTable(arr, 100)
        @test fetch(dt) == nt
        @test fetch(dt, NamedTuple) == nt
        @test tabletype(dt) == NamedTuple # NamedTuple is the type obtained using the sink function of Arrow.Table
    end

    @testset "constructors - Tables.jl compatibility (CSV)" begin
        size = 1_000
        nt = (a = rand(Int, size), b = rand(Int, size))
        io = IOBuffer()
        CSV.write(io, nt)
        dt1 = CSV.read(take!(io), (csv -> DTable(csv, 100)))
        CSV.write(io, nt)
        # or
        dt2 = DTable(CSV.File(take!(io)), 100)
        @test fetch(dt1) == nt
        @test fetch(dt1, NamedTuple) == nt
        @test tabletype(dt1) == NamedTuple
        @test fetch(dt2) == nt
        @test fetch(dt2, NamedTuple) == nt
        @test tabletype(dt2) == NamedTuple
    end

    @testset "constructors - file input" begin
        n = 20
        size = 1000
        ios = [IOBuffer() for _ in 1:n]
        data = [(a = rand(size), b = rand(size)) for _ in 1:n]
        arr = [Arrow.write(ios[idx], data[idx]) for idx in 1:n]

        da = DTable(x -> Arrow.Table(take!(ios[tryparse(Int64, x)])), [string(i) for i in 1:n])
        db = vcat([DataFrame(d) for d in data]...)
        @test fetch(da, DataFrame) == db
        @test Dagger.resolve_tabletype(da) == NamedTuple
        @test da.tabletype === nothing
        tabletype!(da)
        @test da.tabletype === NamedTuple
    end

    @testset "map" begin
        size = 1_000
        nt = (a = rand(size), b = rand(size))

        df = DataFrame(nt)
        for src in [nt, df]
            dt = DTable(src, 100)

            # Single result
            mt = map(x -> (r = x.a + x.b, ), dt)
            mf = map(x -> x.a + x.b, eachrow(df))
            @test fetch(mt).r == mf

            # Two results
            mt = map(x -> (r1 = x.a + x.b, r2 = x.a - x.b), dt)
            mf = combine(df, [:a, :b] => ((a, b) -> a .+ b) => :r1, [:a, :b] => ((a, b) -> a .- b) => :r2)
            @test fetch(mt).r1 == mf.r1
            @test fetch(mt).r2 == mf.r2
        end

        # Map an empty dtable
        dt = DTable((a = [], b = []), 10)
        m = map(x -> (r = x.a + x.b), dt)
        @test fetch(m) == NamedTuple()
    end

    @testset "filter" begin
        size = 1_000
        nt = (a = rand(size), b = rand(size))

        df = DataFrame(nt)

        for src in [nt, df]
            dt = DTable(src, 100)

            dfr = filter(x -> (x.a .> 0.5) .& (x.b .< 0.5), df)
            dtr = filter(x -> (x.a .> 0.5) .& (x.b .< 0.5), dt)
            @test fetch(dtr, DataFrame) == dfr
        end

        # Filter an empty DTable
        dt = DTable((a = [], b = []), 10)
        f = filter(x -> x.a .> 0.5, dt)
        @test fetch(f) == NamedTuple()
    end

    @testset "reduce" begin
        size = 1_000
        nt = (a=rand(Int, size) .% 100, b=rand(Int, size) .% 100)

        df = DataFrame(nt)
        dtdf = DTable(df, 100)
        dtnt = DTable(nt, 100)

        dtdf1 = reduce(+, dtdf, cols=[:a])
        dtnt1 = reduce(+, dtnt, cols=[:a])
        df1 = reduce((x, y) -> x + y.a, eachrow(df); init=0)

        dtdf2 = reduce(+, map(x -> (r = x.a * x.b,), dtdf), cols=[:r])
        dtnt2 = reduce(+, map(x -> (r = x.a * x.b,), dtnt), cols=[:r])
        df2 = reduce((x, y) -> x + y.a * y.b, eachrow(df);init=0)

        dtdf3 = reduce(+, map(x -> (r = x.a + x.b,), dtdf), cols=[:r])
        dtnt3 = reduce(+, map(x -> (r = x.a + x.b,), dtnt), cols=[:r])
        df3 = reduce((x, y) -> x + y.a + y.b, eachrow(df);init=0)

        dtdf4 = reduce(*, dtdf, cols=[:a])
        dtnt4 = reduce(*, dtnt, cols=[:a])
        df4 = reduce((x, y) -> x * y.a, eachrow(df); init=1)

        @test fetch(dtdf1).a == fetch(dtnt1).a == df1
        @test fetch(dtdf2).r == fetch(dtnt2).r == df2
        @test fetch(dtdf3).r == fetch(dtnt3).r == df3
        @test fetch(dtdf4).a == fetch(dtnt4).a == df4

        all_reduce = reduce(+, dtdf)
        df5 = reduce((x, y) -> x + y.b, eachrow(df); init=0)

        @test fetch(all_reduce).a == df1
        @test fetch(all_reduce).b == df5
    end

    @testset "chaining ops" begin
        nt = (a=1:100, b=(1:100))

        d = DTable(nt, 2)
        f1 = filter(x -> iseven.(x.a), d)
        m1 = map(x -> (r = x.a + x.b,), f1)
        f2 = filter(x -> x.r > 50, m1)
        r = reduce(+, f2)
        @test fetch(r).r == 4788
        r = reduce(*, f2, init=BigInt(1))
        @test fetch(r).r == reduce(*, fetch(f2).r, init=BigInt(1))
    end

    @testset "trim" begin
        nt = (a=1:100, b=(1:100))

        d = DTable(nt, 2)
        fd = filter(x -> x.a > 49, d)

        @test length(d.chunks) == length(fd.chunks)
        @test length(trim(fd).chunks) == 26
        trim!(fd)
        @test length(fd.chunks) == 26

        fd2 = filter(x -> x.a > 100, d)
        @test length(d.chunks) == length(fd2.chunks)
        @test length(trim(fd2).chunks) == 0
        trim!(fd2)
        @test length(fd2.chunks) == 0
    end

    @testset "tabletype" begin
        nt = (a=1:100, b=(1:100))
        d = DTable(nt, 2)
        @test tabletype(d) == NamedTuple
        d.tabletype = nothing
        @test Dagger.resolve_tabletype(d) == NamedTuple
        tabletype!(d)
        @test d.tabletype == NamedTuple

        # empty dtable case
        dt = DTable((a = [], b = []), 10)
        @test tabletype(dt) == NamedTuple # fallback in case it can't be found
    end
end


text:
@testset "save_cells.jl" begin

end


text:
verbose = true
 if verbose && !isdefined(Main,:RvSpectMLBase)   println("# Loading RvSpecMLBase")    end
 using RvSpectMLBase
 if verbose && !isdefined(Main,:EchelleInstruments)   println("# Loading EchelleInstruments")    end
 using EchelleInstruments, EchelleInstruments.EXPRES
 using DataFrames, Query

# USER: You must create a data_paths.jl file in one of the default_paths_to_search listed below. It need only contain one line:
# expres_data_path = "/path/to/EXPRES/data/not/including/target_subdir"
target_subdir = "4628"   # USER: Replace with directory of your choice
 fits_target_str = "4628"
 output_dir = "examples/output/"
 paths_to_search_for_param = [pwd(),"examples",joinpath(pkgdir(RvSpectMLBase),"..","RvSpectML","examples"), "/gpfs/group/ebf11/default/ebf11/expres"]
 # NOTE: make_manifest does not update its paths_to_search when default_paths_to_search is defined here, so if you change the line above, you must also include "paths_to_search=default_paths_to_search" in the make_manifest() function call below
 pipeline_plan = PipelinePlan()
 dont_make_plot!(pipeline_plan, :movie)

reset_all_needs!(pipeline_plan)
if need_to(pipeline_plan,:read_spectra)
   if verbose println("# Finding what data files are avaliable.")  end
   eval(read_data_paths(paths_to_search=paths_to_search_for_param))
   @assert isdefined(Main,:expres_data_path)
   df_files = make_manifest(expres_data_path, target_subdir, EXPRES )

   if verbose println("# Reading in customized parameters from param.jl.")  end
   eval(code_to_include_param_jl(paths_to_search=paths_to_search_for_param))

   if verbose println("# Reading in ", size(df_files_use,1), " FITS files.")  end
   @time all_spectra = map(row->EXPRES.read_data(row,store_min_data=true, store_tellurics=true, normalization=:continuum),eachrow(df_files_use))
   #@time all_spectra = map(row->EXPRES.read_data(row,store_min_data=true, store_tellurics=true, normalization=:raw),eachrow(df_files_use))
   #@time all_spectra = map(row->EXPRES.read_data(row,store_min_data=true, store_tellurics=true, normalization=:blaze),eachrow(df_files_use))
   #RvSpectML.discard_blaze(all_spectra)
   #RvSpectML.discard_continuum(all_spectra)  # Seeing if affects CCFs
   GC.gc()
   dont_need_to!(pipeline_plan,:read_spectra)
   all_spectra
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using Dash
using DashBootstrapComponents
using DataFrames
using Formatting
using HTTP
using HttpCommon
using JSON3
using Logging, LoggingExtras
using Markdown
using Mux
import Mux.WebSockets
using Observables
using PlotlyJS
using PovertyAndInequalityMeasures
using StatsBase
using UUIDs

using ScottishTaxBenefitModel
using .BCCalcs
using .Definitions
using .ExampleHelpers
using .FRSHouseholdGetter
using .GeneralTaxComponents
using .ModelHousehold
using .Monitor
using .Runner
using .RunSettings
using .SimplePovertyCounts: GroupPoverty
using .SingleHouseholdCalculations
using .STBIncomes
using .STBOutput
using .STBParameters
using .Utils

import Base.Threads.@spawn


text:
"""
``MinDistance`` is an abstract type to be used
as the first argument for functions.
"""
abstract MinDistance


text:
import ..UncertainValues: CertainValue

constrain(v::CertainValue) = v
constrain(v::CertainValue, s::SamplingConstraint) = v
constrain(v::CertainValue, s::TruncateLowerQuantile) = v
constrain(v::CertainValue, s::TruncateUpperQuantile) = v
constrain(v::CertainValue, s::TruncateQuantiles) = v
constrain(v::CertainValue, s::TruncateStd) = v

export constrain


text:
module CachedCalls

using FilePathsBase
using FilePathsBase: /
using JLSO

export @cached_call, @hash_call
export cachedcalls_dir

const CACHEDCALLS_PATH = Ref{PosixPath}()

"""
    @cached_call f(args; kwargs)

Caches the result of `f(args; kwargs)` to disk and returns the result. The next time
`f(args; kwargs)` is called with the same values of `args` and `kwargs` the cached result
is returned and `f` is not called again. Splatting is not yet supported.

Restrictions on `f` apply: it must not mutate its arguments or access/mutate globals.
These assumptions will not be checked and if violated could mean that an incorrect
result is returned.

Functions are differentiated by name only, meaning that changing the definition and
rerunning `@cached_call` will return the wrong result.
"""
macro cached_call(ex)
    return quote
        h = @hash_call $(esc(ex))

        fname = $(cachedcalls_dir()) / "$(h).jlso"
        if isfile(fname)
            res = JLSO.load(fname)[:res]
        else
            res = $(esc(ex))
            JLSO.save(fname, :res => res)
        end
        res
    end
end

"""
    @hash_call f(args; kwargs)

Computes the hash of the function call `f(args; kwargs)` by hashing `f`, the values of
`args`, the names of `kwargs`, and the values of `kwargs`.

Different order of kwargs will hash differently. Setting the default values of kwargs
explicitly will hash differently than using the defaults implicitly. Splatting is not yet
supported.
"""
macro hash_call(ex)
    func, args, kw_names, kw_values = _deconstruct(ex)

    return :(hash([
        $(esc(func)), # function name
        $(esc.(args)...), # arg values
        $(kw_names)..., # kwarg names
        $(esc.(kw_values)...) # kwarg values
    ]))
end

"""
    cachedcalls_dir()

Retrieves the path to where the cached files are stored.
"""
function cachedcalls_dir()
    if !isassigned(CACHEDCALLS_PATH)
        CACHEDCALLS_PATH[] = PosixPath(first(Base.DEPOT_PATH)) / ".cachedcalls"
    end

    isdir(CACHEDCALLS_PATH[]) || mkdir(CACHEDCALLS_PATH[])

    return CACHEDCALLS_PATH[]
end

"""
    cachedcalls_dir(p::Union{String, PosixPath})

Sets the path to where the cached files are stored to `p`.
"""
function cachedcalls_dir(p::PosixPath)
    CACHEDCALLS_PATH[] = p
end
function cachedcalls_dir(p::String)
    CACHEDCALLS_PATH[] = PosixPath(p)
end

"""
    _deconstruct(ex)

Deconstruct expression `ex` to a tuple (function name, arguments, kwarg names, kwarg values)
"""
function _deconstruct(ex)
    # escaped expressions need their unescaped subexpression deconstructed
    isesc = Meta.isexpr(ex, :escape)
    nonesc_ex = isesc ? ex.args[1] : ex

    # check assumptions about the call
    Meta.isexpr(nonesc_ex, :call) || error("Only :call expressions are supported, $ex was given.")

    # Extract arguments, kwarg names, and kwarg values.
    func, fargs = Iterators.peel(nonesc_ex.args)

    args = filter(subex -> !Meta.isexpr(subex, [:kw, :parameters]), collect(fargs))
    kwargs = _extract_kwargs(collect(fargs))
    kw_names = first.(kwargs)
    kw_values = last.(kwargs)

    # If the original expression was escaped, we have extracted the non-escaped expression
    # to deconstruct, and must escape the indvidual expression components instead.
    if isesc
        args = esc.(args)
        kw_values = esc.(kw_values)
        func = esc(func)
    end

    return func, args, kw_names, kw_values
end

"""
    _extract_kwargs(callargs::AbstractArray{Any}; keep_args=false)

Extract a tuple of (:kwarg_name, value) from :call expression args.
"""
function _extract_kwargs(x; keep_args=false)
    return keep_args ? [(x, x)] : []
end

function _extract_kwargs(ex::Expr; keep_args=false)
    # kwargs specified without ;
    if Meta.isexpr(ex, :kw)
        return [(ex.args[1], ex.args[2]),]

    # kwargs specified with ;
    elseif Meta.isexpr(ex, :parameters)
        # need to `keep_args` because of `f(;c)` syntactic sugar does not result in :kw
        # expression but in a single arg :c to parameters
        return [(_extract_kwargs.(ex.args; keep_args=true)...)...]

    # container.key, container[index] access, splatting
    # (just pass it on)
    elseif Meta.isexpr(ex, [:., :ref, :(...)])
        return keep_args ? [(ex, ex)] : []

    else
        error("Unexpected input expression to _extract_kwargs: $ex")
    end
end

function _extract_kwargs(a::AbstractArray; keep_args=false)
    return [(_extract_kwargs.(a)...)...]
end

end


text:
include("utils.jl")

@testset "Euclidean" begin
    E = Manifolds.Euclidean(3)
    Ec = Manifolds.Euclidean(3;field=ℂ)
    EM = Manifolds.MetricManifold(E,Manifolds.EuclideanMetric())
    @test is_default_metric(EM) == Val{true}()
    @test is_default_metric(E,Manifolds.EuclideanMetric()) == Val{true}()
    x = zeros(3)
    @test det_local_metric(EM,x) == one(eltype(x))
    @test log_local_metric_density(EM,x) == zero(eltype(x))
    @test project_point!(E,x) == x
    @test manifold_dimension(Ec) == 2*manifold_dimension(E)

    manifolds = [ E, EM, Ec ]
    types = [
        Vector{Float64},
        MVector{3, Float64},
        Vector{Float32},
        Vector{Double64},
    ]

    types_complex = [
        Vector{ComplexF64},
        MVector{3, ComplexF64},
        Vector{ComplexF32},
        Vector{ComplexDF64},
    ]

    for M in manifolds
        basis_types = if M == E
            (ArbitraryOrthonormalBasis(), ProjectedOrthonormalBasis(:svd), DiagonalizingOrthonormalBasis([1.0, 2.0, 3.0]))
        elseif M == Ec
            (ArbitraryOrthonormalBasis(), DiagonalizingOrthonormalBasis([1.0, 2.0, 3.0]))
        else
            ()
        end
        for T in types
            @testset "$M Type $T" begin
                pts = [convert(T, [1.0, 0.0, 0.0]),
                       convert(T, [0.0, 1.0, 0.0]),
                       convert(T, [0.0, 0.0, 1.0])]
                test_manifold(
                    M,
                    pts,
                    test_reverse_diff = isa(T, Vector),
                    test_project_tangent = true,
                    test_musical_isomorphisms = true,
                    test_vector_transport = true,
                    test_mutating_rand = isa(T, Vector),
                    point_distributions = [Manifolds.projected_distribution(M, Distributions.MvNormal(zero(pts[1]), 1.0))],
                    tvector_distributions = [Manifolds.normal_tvector_distribution(M, pts[1], 1.0)],
                    basis_types_vecs = basis_types,
                    basis_types_to_from = basis_types,
                    basis_has_specialized_diagonalizing_get = true
                )
            end
        end
    end
    for T in types_complex
        @testset "Complex Euclidean, type $T" begin
            pts = [convert(T, [1.0im, -1.0im, 1.0]),
                   convert(T, [0.0, 1.0, 1.0im]),
                   convert(T, [0.0, 0.0, 1.0])]
            test_manifold(
                Ec,
                pts,
                test_reverse_diff = isa(T, Vector),
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_vector_transport = true
            )
        end
    end

    @testset "hat/vee" begin
        E = Euclidean(3, 2)
        x = collect(reshape(1.0:6.0, (3, 2)))
        v = collect(reshape(7.0:12.0, (3, 2)))
        @test hat(E, x, vec(v)) ≈ v
        w = similar(v)
        @test hat!(E, w, x, vec(v)) === w
        @test w ≈ v
        @test vee(E, x, v) ≈ vec(v)
        w = similar(vec(v))
        @test vee!(E, w, x, v) === w
        @test w ≈ vec(v)
    end
end


text:
using BenchmarkTools
using DelimitedFiles

q = [0.001,0.002,0.003,0.003,0.004,0.004,0.005,0.007,0.009,0.011]
w = [0.05,0.07,0.08,0.10,0.14,0.20,0.20,0.20,0.10,0.04]
P = 100
S = 25000
r = 0.02

function npv1(q,w,P,S,r)
	inforce = [1.; cumprod(1 .- q .- w)[1:end-1]]
  	ncf = inforce .* P .- inforce .* q .* S
 	d = (1 ./ (1 + r)) .^ (1:length(ncf))
  	return sum(ncf .* d)
end

function npv2(q,w,P,S,r)
	inforce = append!([1.0],@view cumprod(1 .- q .- w)[1:end-1])
  	ncf = inforce .* P .- inforce .* q .* S
 	d = (1 ./ (1 + r)) .^ (1:length(ncf))
  	return sum(ncf .* d)
end

@inline function npv3(q,w,P,S,r,term=nothing)
    term = term === nothing ? length(q) + 1 : term + 1
    inforce = 1.0
    result = 0.0
    v = (1 / ( 1 + r))
    v_t = v

    for (t,(q,w)) in enumerate(zip(q,w))
        deaths = t < term ? inforce * q : 0.0
        lapses = t < term ? inforce * w : 0.0
        premiums = inforce * P
        claims = deaths * S
        ncf = premiums - claims
        result += ncf *  v_t
        v_t *= v
        inforce = inforce - deaths - lapses
    end

    return result
end

@inline function npv4(q,w,P,S,r,term=nothing)
    term = term === nothing ? length(q) : term
    inforce = 1.0
    result = 0.0
    v = (1 / ( 1 + r))
    v_t = v

    for (t,(q,w)) in enumerate(zip(q,w))
        t > term && return result
        result += inforce * (P - S * q) * v_t
        inforce -= inforce * q + inforce * w
        v_t *= v
    end

    return result
end

# suggested by Stefan Karpinski
@inline function npv5(q,w,P,S,r,term=nothing)
    term = term === nothing ? length(q) : term
    inforce = 1.0
    result = 0.0
    v = (1 / ( 1 + r))
    v_t = v
    for (q,w,_) in zip(q,w,1:term)
        result += inforce * (P - S * q) * v_t
        inforce -= inforce * q + inforce * w
        v_t *= v
    end
    return result
end

# suggested by Michael Abbott
@inline function npv6(qs,ws,P,S,r,term=length(qs))
    inforce, result = 1.0, 0.0
    v = 1 / ( 1 + r)
    v_t = v

    for t in 1:min(term, length(qs), length(ws))
        @inbounds q, w = qs[t], ws[t]
        result += inforce * (P - S * q) * v_t
        inforce -= inforce * (q + w)
        v_t *= v
    end

    return result
end

# suggested by Peter J
@inline function npv7(qs,ws,P,S,r,term=length(qs))
    inforce, result = 1.0, 0.0
    v = 1 / ( 1 + r)
    v_t = v

    for t in 1:min(term, length(qs), length(ws))
        @inbounds q, w = qs[t], ws[t]
        result += inforce * (P - S * q) * v_t
        inforce = inforce*(1 - (q + w))
        v_t *= v
    end

    return result
end

# suggested by Peter J
@inline function npv8(qs,ws,P,S,r,term=length(qs))
    inforce, result = 1.0, 0.0
    v = 1 / ( 1 + r)
    v_t = v
    @inbounds @simd for t in 1:min(term, length(qs), length(ws))
        q, w = qs[t], ws[t]
        result += inforce * (P - S * q) * v_t
        inforce = inforce*(1 - (q + w))
        v_t *= v
    end
    return result
end

# suggested by Michael Abbott
@inline function npv9(qs,ws,P,S,r,term=length(qs))
    inforce, result = 1.0, 0.0
    v = 1 / ( 1 + r)
    v_t = v
    @inbounds @fastmath for t in 1:min(term, length(qs), length(ws))
        q, w = qs[t], ws[t]
        result += inforce * (P - S * q) * v_t
        inforce = inforce*(1 - (q + w))
        v_t *= v
    end
    return result
end

results = map([npv1,npv2,npv3,npv4,npv5,npv6,npv7,npv8,npv9]) do f
    b = @benchmark $f($q,$w,$P,$S,$r)
    (f="$f",mean=mean(b),median=median(b))
end

writedlm("julia_results.csv",results, ',')

max_stars_repo_name: JeannotJeannot/Learn

max_issues_repo_name: JeannotJeannot/Learn

max_forks_repo_name: JeannotJeannot/Learn


text:
# tokenizers.jl
#
# author: Wade Shen
# swade@ll.mit.edu
# Copyright © 2014 Massachusetts Institute of Technology, Lincoln Laboratory
# version 0.1
#
# 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.
export english_tokenizer, twenglish_tokenizer, replace_html_entities, twenglish_cleaner

# -------------------------------------------------------------------------------------------------------------------------
# Tokenizers
# -------------------------------------------------------------------------------------------------------------------------
const punct_word      = r"^[\p{P}\p{Po}\p{Sm}]*(.*?)[\p{P}\p{Po}\p{Sm}]*$"
const english_space   = r"[\s\p{Zs}_-]+"
const url_pattern     = r"http://[^\s]*"
const hashtag_pattern = r"^#.*$"
const mention_pattern = r"^@.*$"

function replace_html_entities(s :: AbstractString)
  replace(s, r"&[^;]+?;", s -> s in keys(html_entity_table) ? html_entity_table[s] : s)
end

function pattern_replace(w :: AbstractString)
  if ismatch(r"^[+-]?\p{Sc}\d+([.,]\d+)*$", w) return "--currency--"
  elseif ismatch(r"^[+-]?\d+([.,]\d+)*%$", w) return "--percent--"
  elseif ismatch(r"^[+-]?\d+([.,]\d+)*$", w) return "--number--"
  elseif ismatch(r"^(https?|ftp):.*$"i, w) return "--url--"
  elseif ismatch(r"\d+:\d+(am|pm)$"i, w) return "--time--"
  elseif ismatch(r"\d+(am|pm)$"i, w) return "--time--"
  elseif ismatch(r"\d+[-/]\d+$", w) return "--date--"
  elseif ismatch(r"\d+[-/]\d+[-/]\d+$", w) return "--date--"
  else return w
  end
end

function prereplace(sent :: AbstractString)
  r = replace(sent, r"n't\b", " not")
  r = replace(r, r"'s\b", " s's")
  r = replace(r, r"'d\b", " d'd")
end

function english_tokenizer(s :: AbstractString)
  return [
    begin
      m = match(punct_word, w)
      nw = (m == nothing ? w : m.captures[1])
      pattern_replace(nw)
    end for w in filter(x -> !ismatch(r"^[\p{P}\p{C}\p{S})]+$", x), split(prereplace(strip(normalize_string(s, :NFKC))), english_space))
  ]
end

function twenglish_tokenizer(s :: AbstractString)
  return [
    begin
      m = match(r"^(\p{P}*)(.*?)\p{P}*$", w)
      nw = (m == nothing) ? w : (m.captures[1] == "#" ? (m.captures[1] * m.captures[2]) : m.captures[2])
      patternReplace(nw)
    end for w in filter(x -> !ismatch(r"^(\p{P}|\p{S})+$", x) && !ismatch(r"^@.*$", x), split(strip(s), r"\s+"))
  ]
end

function twenglish_cleaner(tw :: AbstractString; urls = true, hashtags = true, mentions = true)
  ctw = replace(normalize_string(tw, :NFKC), default_space, " ")
  ctw = urls ? replace(ctw, url_pattern, "\u0030\u20E3") : ctw

  while true
    x = replace_html_entities(ctw)
    if x != ctw
      ctw = x
    else
      break
    end
  end

  words = split(ctw)
  words = hashtags ? [ ismatch(hashtag_pattern, w) ? "\u0023\u20E3" : w for w in words ] : words
  words = mentions ? [ ismatch(hashtag_pattern, w) ? "\u0031\u20E3" : w for w in words ] : words

  return join(words, " ")
end


text:
module baffled_piston_examples
using FinEtools
using FinEtoolsAcoustics
using LinearAlgebra
using UnicodePlots

function baffled_piston_H8_transient()
	rho = 1.21*phun("kg/m^3");# mass density
	c  = 343.0*phun("m/s");# sound speed
	bulk =  c^2*rho;
	frequency=2000; # Hz
	omega=2*pi*frequency;
	P_piston = 1.0e-3*phun("MPa"); # amplitude of the piston pressure
	R = 50.0*phun("mm");# radius of the piston
	Ro = 400.0*phun("mm"); # radius of the enclosure
	nref = 2;#number of refinements of the sphere around the piston
	nlayers = 48;                     # number of layers of elements surrounding the piston
	tolerance = R/(2^nref)/10
	dt=1.0/frequency/20;
	tfinal=50*dt;
	nsteps=round(tfinal/dt)+1;

	println("""

		Baffled piston in a half-sphere domain with hard boundary condition.

		Hexahedral mesh. Algorithm version.
		""")

	t0  =  time()

	# Hexahedral mesh
	fens,fes  =  H8sphere(R,nref);
	bfes  =  meshboundary(fes)
	# File  =   "baffledabc_boundary.vtk"
	# vtkexportmesh(File, bfes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
	#  @async run(`"paraview.exe" $File`)

	l = selectelem(fens,bfes,facing = true,direction = [1.0 1.0  1.0])
	ex(xyz, layer) = (R+layer/nlayers*(Ro-R))*xyz/norm(xyz)
	fens1,fes1  =  H8extrudeQ4(fens, subset(bfes,l), nlayers, ex);
	fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
	fes = cat(newfes1,fes2)

	# Piston surface mesh
	bfes  =  meshboundary(fes)
	l1 = selectelem(fens, bfes, facing = true, direction = [-1.0 0.0 0.0])
	l2 = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0], inflate = tolerance)
	piston_fes = subset(bfes,intersect(l1,l2));
	# File  =   "piston_fes.vtk"
	# vtkexportmesh(File, piston_fes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
	#  @async run(`"paraview.exe" $File`)

	# Outer spherical boundary
	louter = selectelem(fens, bfes, facing = true, direction = [1.0 1.0  1.0], dotmin= 0.001)
	outer_fes = subset(bfes,louter);
	# File  =   "outer_fes.vtk"
	# vtkexportmesh(File, outer_fes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
	#  @async run(`"paraview.exe" $File`)

	println("Pre-processing time elapsed  =  ",time() - t0,"s")

	t1  =  time()

	material = MatAcoustFluid(bulk, rho)
	# Region of the fluid
	femm = FEMMAcoust(IntegDomain(fes, GaussRule(3, 2)), material)

	geom = NodalField(fens.xyz);
	P = NodalField(zeros(nnodes(geom),1));
	piston_fenids = connectednodes(piston_fes);
	setebc!(P, piston_fenids, true, 1, 0.0);
	applyebc!(P);
	numberdofs!(P);

	S  =  acousticstiffness(femm, geom, P);
	C  =  acousticmass(femm, geom, P);
	D= (2.0/dt)*S +(dt/2.0)*C;

	P0 = deepcopy(P)
	Pdd0 = deepcopy(P)
	P1 = deepcopy(P)
	Pdd1 = deepcopy(P)
	TMPF = deepcopy(P)
	P = nothing # we don't need this field anymore
	vP0 = zeros(P0.nfreedofs)
	vP1 = deepcopy(vP0)
	vPd0 = deepcopy(vP0)
	vPd1 = deepcopy(vP0)

	t = 0.0
	P0.fixed_values[piston_fenids,1] .= P_piston*sin(omega*t)
	Pdd0.fixed_values[piston_fenids,1] .= P_piston*(-omega^2)*sin(omega*t)
	vP0 = gathersysvec!(P0, vP0)

	nh = selectnode(fens, nearestto = [R+Ro/2, 0.0, 0.0] )
	Pnh = [P1.values[nh, 1][1]]
	ts = [t]
	step =0;
	while t <=tfinal
		step = step  +1;
		t=t+dt;
		P1.fixed_values[piston_fenids,1] .= P_piston*sin(omega*t)
		Pdd1.fixed_values[piston_fenids,1] .= P_piston*(-omega^2)*sin(omega*t)
		TMPF.fixed_values = P0.fixed_values + P1.fixed_values
		F = nzebcloadsacousticmass(femm, geom, TMPF);
		TMPF.fixed_values = Pdd0.fixed_values + Pdd1.fixed_values
		F = F + nzebcloadsacousticstiffness(femm, geom, TMPF);
		# println("$(norm(F))")
		vPd1 = D\((2/dt)*(S*vPd0) - C*(2*vP0+(dt/2)*vPd0) + F);
		vP1 = vP0 + (dt/2)*(vPd0+vPd1);
		scattersysvec!(P1, vP1); # store current pressure
		# Swap variables for the next step
		copyto!(vP0, vP1)
		copyto!(vPd0, vPd1)
		copyto!(P0, P1)
		copyto!(Pdd0, Pdd1)
		# Graphics output
		File  =   "baffled_piston-$(step).vtk"
		vtkexportmesh(File, fes.conn, fens.xyz, FinEtools.MeshExportModule.H8;
			scalars = [(  "P", P1.values+P1.fixed_values)])
		#  @async run(`"paraview.exe" $File`)
		push!(Pnh, P1.values[nh, 1][1])
		push!(ts, t)
		println("step = $( step )")
	end

	# @pgf a = Axis({
	#     xlabel = "Time",
	#     ylabel = "Pnh",
	#     title = "Transient pressure"
	# },
	# Plot(Table([:x => vec(ts), :y => vec(Pnh)])))
	# display(a)
	plt = lineplot(vec(ts), vec(Pnh), canvas = DotCanvas, title = "Transient pressure", name = "P", xlabel = "Time", ylabel = "Pnh")
	display(plt)

	# using Plots
	# plotly()
	# plot(vec(ts), vec(Pnh))
	# gui()
	println("$(Pnh[end])")

	# # Surface for the ABC
	# abc1  =  FDataDict("femm"=>FEMMAcoustSurf(IntegDomain(outer_fes, GaussRule(2, 2)),
	#           material))
	#
	# # Surface of the piston
	# flux1  =  FDataDict("femm"=>FEMMAcoustSurf(IntegDomain(piston_fes, GaussRule(2, 2)),
	#           material),  "normal_flux"=> -rho*a_piston+0.0im);
	#
	# # Make model data
	# modeldata =  FDataDict("fens"=>  fens,
	#                  "omega"=>omega,
	#                  "regions"=>[region1],
	#                  "flux_bcs"=>[flux1], "ABCs"=>[abc1])
	#
	# # Call the solver
	# modeldata = FinEtools.AlgoAcoustModule.steadystate(modeldata)
	#
	# println("Computing time elapsed  =  ",time() - t1,"s")
	# println("Total time elapsed  =  ",time() - t0,"s")
	#
	# geom = modeldata["geom"]
	# P = modeldata["P"]
	#
	# File  =   "baffledabc.vtk"
	# vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H8;
	# scalars = [("absP", abs.(P.values))])
	# @async run(`"paraview.exe" $File`)
	#
	# # using Winston
	# # pl  =  FramedPlot(title = "Matrix",xlabel = "x",ylabel = "Re P, Im P")
	# # setattr(pl.frame, draw_grid = true)
	# # add(pl, Curve([1:length(C[:])],vec(C[:]), color = "blue"))
	#
	# # # pl = plot(geom.values[nLx,1][ix],scalars[nLx][ix])
	# # # xlabel("x")
	# # # ylabel("Pressure")
	# # display(pl)
	#
	# true
end # baffled_piston_H8_transient

function allrun()
    println("#####################################################")
    println("# baffled_piston_H8_transient ")
    baffled_piston_H8_transient()
    return true
end # function allrun

end # module baffled_piston_examples


text:
# ---
# title: 1494. Parallel Courses II
# id: problem1494
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Graph
# link: <https://leetcode.com/problems/parallel-courses-ii/description/>
# hidden: true
# ---
#
# Given the integer `n` representing the number of courses at some university
# labeled from `1` to `n`, and the array `dependencies` where `dependencies[i] =
# [xi, yi]`  represents a prerequisite relationship, that is, the course `xi`
# must be taken before the course `yi`.  Also, you are given the integer `k`.
#
# In one semester you can take **at most** `k` courses as long as you have taken
# all the prerequisites for the courses you are taking.
#
# _Return the minimum number of semesters to take all courses_. It is guaranteed
# that you can take all courses in some way.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/05/22/leetcode_parallel_courses_1.png)**
#
#
#
#     Input: n = 4, dependencies = [[2,1],[3,1],[1,4]], k = 2
#     Output: 3
#     Explanation: The figure above represents the given graph. In this case we can take courses 2 and 3 in the first semester, then take course 1 in the second semester and finally take course 4 in the third semester.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/05/22/leetcode_parallel_courses_2.png)**
#
#
#
#     Input: n = 5, dependencies = [[2,1],[3,1],[4,1],[1,5]], k = 2
#     Output: 4
#     Explanation: The figure above represents the given graph. In this case one optimal way to take all courses is: take courses 2 and 3 in the first semester and take course 4 in the second semester, then take course 1 in the third semester and finally take course 5 in the fourth semester.
#
#
# **Example 3:**
#
#
#
#     Input: n = 11, dependencies = [], k = 2
#     Output: 6
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 15`
#   * `1 <= k <= n`
#   * `0 <= dependencies.length <= n * (n-1) / 2`
#   * `dependencies[i].length == 2`
#   * `1 <= xi, yi <= n`
#   * `xi != yi`
#   * All prerequisite relationships are distinct, that is, `dependencies[i] != dependencies[j]`.
#   * The given graph is a directed acyclic graph.
#
#
## @lc code=start
using LeetCode

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


text:
module SLAMData
using ImageMagick
import CameraModel
import Printf
# import Statistics

export TumData, TumCamROS, read, fastforward!, fastbackward!,
       Trajectory, load_trajectory, load_covs,

       # draw_trajectory!, draw_trajectory,
       getTno2Idx, readid, readk,

       # klg
       KlgReader, readRGBD, readRGBDFrame,

       # traj
       load_redwood_traj, load_redwood_loop, save_redwood_file,
       load_tum_traj, save_tum_traj, load_boundlefusion_traj,
       load_redwood_loop_info, odom2pose,
       # load img to array
       loadpng,

       # sun3d reader
       Sun3DReader

include("types.jl")

include("rgbd_reader.jl")
include("stereo_reader.jl")
include("tum.jl")
include("sun3d.jl")
include("trajectory.jl")
include("klg.jl")

end # module


text:
" Solve a simple Linear least squares problem"
function solve(prob::LinearLeastSquares)
    @unpack K, y = prob
    K\y
end

"Solve a weighted linear least squares problem"
function solve(prob::WeightedLinearLeastSquares)
    @unpack K, y, Sₑ = prob
    Sₑ⁻¹ = inv(Sₑ)
    # Will be re-used:
    KᵀSₑ⁻¹K = K' * Sₑ⁻¹ * K
    # Posterior covariance matrix:
    Ŝ = inv(KᵀSₑ⁻¹K)
    # Averaging kernel matrix:
    A = Ŝ * KᵀSₑ⁻¹K
    # State vector:
    x̂ = Ŝ * K'Sₑ⁻¹ * y
    #sqrtSₑ⁻¹ = sqrt(inv(Sₑ));
    #Ŝ = inv(KᵀSₑ⁻¹K
    #(sqrtSₑ⁻¹ * K) \ (sqrtSₑ⁻¹ * y)
end

max_stars_repo_name: RadiativeTransfer/OptimalEstimation

max_issues_repo_name: RadiativeTransfer/OptimalEstimation

max_forks_repo_name: RadiativeTransfer/OptimalEstimation


text:
using LinearAlgebra, PyPlot, Random, FFTW, Distributions, CircularArrays, Printf, Statistics

close(:all)

ci, co, N, ps, b = 16, 16, 256, 64, 50
plot_I = false

eb = randn(Float32, co*N, ps)

chunk(arr, n) = [arr[i:min(i + n - 1, end)] for i in 1:n:length(arr)]
simil(x, y) = 100*dot(x[:], y[:])/(norm(x)*norm(y))
error(x, y) = norm(x - y)/(norm(x)+norm(y))

function draw_e(ps::Integer, co::Integer, N::Integer)
    n = ps ÷ co
    if n < 8
        n = 8
        inds = chunk(randperm(n*co) .% ps .+ 1, n)
        overlap = 1 #n / (ps ÷ co)
    else
        inds = chunk(randperm(ps), n)
        overlap = 1
    end
    e = zeros(Float32, N*co, ps)
    for i=1:co
        iloc = inds[i]
        ni = length(iloc)
        e[(i-1)*N+1:i*N, iloc] .= Float32(sqrt(ps/(overlap*ni)))*randn(Float32, N, ni)
    end

    return e
end

eortho = draw_e(ps, co, N)

plot_e_scale = ps/50
figure()
subplot(221)
imshow(eb, vmin=-.1, vmax=.1, cmap="seismic", aspect="auto")
title("Z, Random")
subplot(222)
imshow(eortho, vmin=-.1, vmax=.1, cmap="seismic", aspect="auto")
title("Z, Random Orthogonalized")
subplot(223)
imshow(eb*eb', vmin=-plot_e_scale, vmax=plot_e_scale, cmap="seismic", aspect="auto")
title(L"$ZZ^T$ Random")
subplot(224)
imshow(eortho*eortho', vmin=-plot_e_scale, vmax=plot_e_scale, cmap="seismic", aspect="auto")
title(L"$ZZ^T$ Random Orthogonalized")
tight_layout()
savefig("./figures/zortho.png", bbox_inches="tight")

max_stars_repo_name: sheevy/XConv

max_issues_repo_name: sheevy/XConv

max_forks_repo_name: sheevy/XConv


text:
# Test headers written by `unu make`
using NRRD, FileIO, Colors, AxisArrays, ImageAxes, Unitful, SimpleTraits
using Test

headerpath = joinpath(dirname(@__FILE__), "headers")
writepath = mktempdir()
const mm = u"mm"
const μm = u"μm"
const s = u"s"

@traitimpl TimeAxis{Axis{:t}}

@testset "unu headers" begin

for (file, T, axs, perm) in (("test2duchar.nhdr", UInt8, (5,6), ()),
                             ("test2dushort.nhdr", UInt16, (5,6), ()),
                             ("test2dfloat.nhdr", Float32, (5,6), ()),
                             ("test2drgb.nhdr", RGB{N0f8}, (5,6), ()),
                             ("test2drgb_last.nhdr", RGB{N0f8}, (5,6), (3,1,2)),

                             ("test2drgb_spu_spacing.nhdr",
                              RGB{N0f8},
                              (Axis{:space_1}(range(0.0,stop=8.0,length=5)*mm),
                               Axis{:space_2}(range(0.0,stop=17.5,length=6)*μm)),
                              ()),

                             ("test2d_spu_spacing_labels.nhdr",
                              UInt8,
                              (Axis{:first}(range(0.0,stop=8.0,length=5)*mm),
                               Axis{:second}(range(0.0,stop=15.0,length=6)*mm)),
                              ()),

                             ("test2drgb_spu_spacing_labels.nhdr",
                              RGB{N0f8},
                              (Axis{:first}(range(0.0,stop=8.0,length=5)*mm),
                               Axis{:second}(range(0.0,stop=15.0,length=6)*mm)),
                              ()),

                             ("test3d_units_range_axes.nhdr",
                              UInt8,
                              (Axis{:R}(range(10.0,stop=14.0,length=3)*mm),
                               Axis{:A}(range(12.0,stop=16.0,length=3)*mm),
                               Axis{:S}(range(100.0,stop=110.0,length=3)*mm)),
                              ()),

                             ("test3dlist_origin.nhdr",
                              UInt8,
                              (Axis{:R}(1:3),
                               Axis{:A}(2:4),
                               Axis{:S}(3:5),
                               Axis{:dim_4}(0:1)),
                              ()),

                             ("test3d_time.nhdr",
                              UInt16,
                              (Axis{:space_1}(0:7),
                               Axis{:space_2}(0:9),
                               Axis{:space_3}(0:2),
                               Axis{:time}(0:799)),
                              ()),

                             ("test3d_time_labels.nhdr",
                              UInt16,
                              (Axis{:x}(0:7),
                               Axis{:l}(0:9),
                               Axis{:s}(0:2),
                               Axis{:t}(0:799)),
                              ()),

                             ("test3d_time_labels_units.nhdr",
                              UInt16,
                              (Axis{:x}(range(0.0,step=0.25,length=8)*μm),
                               Axis{:l}(range(0.0,step=0.25,length=10)*μm),
                               Axis{:s}(range(0.0,step=2.0,length=3)*μm),
                               Axis{:time}(0:799)),
                              ()),

                             ("test3d_time_labels_units_s.nhdr",
                              UInt16,
                              (Axis{:x}(range(0.0,step=0.25,length=8)*μm),
                               Axis{:l}(range(0.0,step=0.25,length=10)*μm),
                               Axis{:s}(range(0.0,step=2.0,length=3)*μm),
                               Axis{:time}(range(0.0,step=0.8,length=800) * s)),
                              ()))
    try
        # Read tests
        filerd = joinpath(headerpath, file)
        Tr, axsr, permr, need_bswap = NRRD.arraytype(filerd)
        @test Tr == T
        @test axsr === axs
        if !(axsr === axs)
            for (axr, ax) in zip(axsr, axs)
                println("axr == ax: ", axr == ax)
                println(axr)
                println(ax)
            end
        end
        @test permr == perm
        if occursin("time", file)
            @test any(istimeaxis, axsr)
        end
        # Write tests
        filewr = joinpath(writepath, file)
        open(query(filerd)) do iord
            skipmagic(iord)
            version, header, keyvals, comments = NRRD.parse_header(iord)
            open(filewr, "w") do iowr
                print(iowr, "NRRD")  # magic bytes
                NRRD.write_header(iowr, version, header, keyvals, comments)
            end
            Tw, axsw, permw, need_bswapw = NRRD.arraytype(filewr)
            @test Tw == Tr
            @test axsw == axsr
            @test permw == permr
        end
    catch err
        warn("failure on file $file")
        rethrow(err)
    end
end

end # @testset "unu headers"


text:
abstract type AbstractDataset end

# Metadata

struct PLEXOSConfig
    element::String
    value::String
end

PLEXOSConfig(e::Node, ::AbstractDataset) = PLEXOSConfig(
    getchildtext("element", e),
    getchildtext("value", e)
)

struct PLEXOSUnit
    value::String
end

PLEXOSUnit(e::Node, ::AbstractDataset) =
    PLEXOSUnit(getchildtext("value", e))

struct PLEXOSTimeslice
    name::String
end

PLEXOSTimeslice(e::Node, ::AbstractDataset) =
    PLEXOSTimeslice(getchildtext("name", e))

struct PLEXOSModel
    name::String
end

PLEXOSModel(e::Node, ::AbstractDataset) =
    PLEXOSModel(getchildtext("name", e))

struct PLEXOSSample
    name::String
end

PLEXOSSample(e::Node, ::AbstractDataset) =
    PLEXOSSample(getchildtext("sample_name", e))

struct PLEXOSSampleWeight
    sample::PLEXOSSample
    phase::Int # maybe Parametrize on this?
    value::Float64
end

const sample_offset = plexostables_lookup["t_sample"].indexoffset
PLEXOSSampleWeight(e::Node, d::AbstractDataset) =
    PLEXOSSampleWeight(d.samples[getchildint("sample_id", e) + sample_offset],
                       getchildint("phase_id", e),
                       getchildfloat("value", e))

struct PLEXOSClassGroup
    name::String
end

PLEXOSClassGroup(e::Node, ::AbstractDataset) =
    PLEXOSClassGroup(getchildtext("name", e))

struct PLEXOSClass
    name::String
    classgroup::PLEXOSClassGroup
    state::Int
end

PLEXOSClass(e::Node, d::AbstractDataset) =
    PLEXOSClass(getchildtext("name", e),
                d.classgroups[getchildint("class_group_id", e)],
                getchildint("state", e))

struct PLEXOSCategory
    name::String
    rank::Int
    class::PLEXOSClass

    # PLEXOS sometimes reports categories without reporting the classes they
    # refer to: if that happens, just leave the class undefined
    # (category won't have any objects anyways)

    function PLEXOSCategory(e::Node, d::AbstractDataset)

        class_idx = getchildint("class_id", e)

        if isassigned(d.classes, class_idx)
            new(getchildtext("name", e),
                getchildint("rank", e),
                d.classes[class_idx])
        else
            new(getchildtext("name", e),
                getchildint("rank", e))
        end

    end

end

struct PLEXOSAttribute
    name::String
    description::String
    enum::Int
    class::PLEXOSClass

    # PLEXOS sometimes reports attributes without reporting the classes they
    # refer to: if that happens, just leave the classes undefined

    function PLEXOSAttribute(e::Node, d::AbstractDataset)

        name = getchildtext("name", e)
        description = getchildtext("description", e)
        enum = getchildint("enum_id", e)
        class_idx = getchildint("class_id", e)

        if isassigned(d.classes, class_idx)
            new(name, description, enum, d.classes[class_idx])
        else
            new(name, description, enum)
        end

    end

end

struct PLEXOSCollection
    name::String
    complementname::Union{Nothing,String}
    parentclass::PLEXOSClass
    childclass::PLEXOSClass

    # PLEXOS sometimes reports collections without reporting the classes they
    # refer to: if that happens, just leave the classes undefined
    # (collection won't have any members anyways)

    function PLEXOSCollection(e::Node, d::AbstractDataset)

        name = getchildtext("name", e)
        cname = getchildtext("complement_name", e)
        cname == "" && (cname = nothing)

        parent_idx = getchildint("parent_class_id", e)
        child_idx = getchildint("child_class_id", e)

        if isassigned(d.classes, parent_idx) && isassigned(d.classes, child_idx)
            new(name, cname, d.classes[parent_idx], d.classes[child_idx])
        else
            new(name, cname)
        end

    end

end

struct PLEXOSCustomColumn
    name::String
    position::Int
    class::PLEXOSClass

    function PLEXOSCustomColumn(e::Node, d::AbstractDataset)

        name = getchildtext("name", e)
        position = getchildint("position", e)
        class_idx = getchildint("class_id", e)

        if isassigned(d.classes, class_idx)
            new(name, position, d.classes[class_idx])
        else
            new(name, position)
        end

    end

end

struct PLEXOSProperty
    name::String
    summaryname::String
    enum::Int
    unit::PLEXOSUnit
    summaryunit::PLEXOSUnit
    ismultiband::Bool
    isperiod::Bool
    issummary::Bool
    collection::PLEXOSCollection

    # PLEXOS sometimes reports properties without reporting the collections they
    # refer to: if that happens, just leave the collections undefined
    # (property won't have any data anyways)

    function PLEXOSProperty(e::Node, d::AbstractDataset)

        collection_idx = getchildint("collection_id", e)

        name = getchildtext("name", e)
        summary_name = getchildtext("summary_name", e)
        enum = getchildint("enum_id", e)
        unit = d.units[getchildint("unit_id", e) + 1]
        summary_unit = d.units[getchildint("summary_unit_id", e) + 1]
        is_multiband = getchildbool("is_multi_band", e)
        is_period = getchildbool("is_period", e)
        is_summary = getchildbool("is_summary", e)

        if isassigned(d.collections, collection_idx)
            new(name, summary_name, enum,
                unit, summary_unit,
                is_multiband, is_period, is_summary,
                d.collections[collection_idx])
        else
            new(name, summary_name, enum,
                unit, summary_unit,
                is_multiband, is_period, is_summary)
        end

    end

end

# System Data

struct PLEXOSObject
    name::String
    class::PLEXOSClass
    category::PLEXOSCategory
    index::Int
    show::Bool
end

PLEXOSObject(e::Node, d::AbstractDataset) =
    PLEXOSObject(getchildtext("name", e),
                 d.classes[getchildint("class_id", e)],
                 d.categories[getchildint("category_id", e)],
                 getchildint("index", e),
                 getchildbool("show", e))

struct PLEXOSMembership
    parentclass::PLEXOSClass
    childclass::PLEXOSClass
    collection::PLEXOSCollection
    parentobject::PLEXOSObject
    childobject::PLEXOSObject
end

PLEXOSMembership(e::Node, d::AbstractDataset) =
    PLEXOSMembership(d.classes[getchildint("parent_class_id", e)],
                     d.classes[getchildint("child_class_id", e)],
                     d.collections[getchildint("collection_id", e)],
                     d.objects[getchildint("parent_object_id", e)],
                     d.objects[getchildint("child_object_id", e)])

struct PLEXOSAttributeData

    value::Float64
    attribute::PLEXOSAttribute
    object::PLEXOSObject

    function PLEXOSAttributeData(e::Node, d::AbstractDataset)

        value = getchildfloat("value", e)
        attribute = d.attributes[getchildint("attribute_id", e)]
        object_idx = getchildint("object_id", e)

        # PLEXOS sometimes reports attributes without reporting the objects they
        # refer to: if that happens, just leave the objects undefined

        if isassigned(d.objects, object_idx)
            new(value, attribute, d.objects[object_idx])
        else
            new(value, attribute)
        end

    end

end

struct PLEXOSMemoObject
    value::String
    column::PLEXOSCustomColumn
    object::PLEXOSObject

    function PLEXOSMemoObject(e::Node, d::AbstractDataset)

        value = getchildtext("value", e)
        object_idx = getchildint("object_id", e)
        column = d.customcolumns[getchildint("column_id", e)]

        if isassigned(d.objects, object_idx)
            new(value, column, d.objects[object_idx])
        else
            new(value, column)
        end

    end

end

# Results

struct PLEXOSPeriod0
    periodofday::Int
    hour::Int
    day::Int
    week::Int
    month::Int
    quarter::Int
    fiscalyear::Int
    datetime::String
end

PLEXOSPeriod0(e::Node, d::AbstractDataset) =
    PLEXOSPeriod0(getchildint("period_of_day", e),
                  getchildint("hour_id", e),
                  getchildint("day_id", e),
                  getchildint("week_id", e),
                  getchildint("month_id", e),
                  getchildint("quarter_id", e),
                  getchildint("fiscal_year_id", e),
                  getchildtext("datetime", e))

struct PLEXOSPeriod1
    week::Int
    month::Int
    quarter::Int
    fiscalyear::Int
    date::String
end

PLEXOSPeriod1(e::Node, d::AbstractDataset) =
    PLEXOSPeriod1(getchildint("week_id", e),
                  getchildint("month_id", e),
                  getchildint("quarter_id", e),
                  getchildint("fiscal_year_id", e),
                  getchildtext("date", e))

struct PLEXOSPeriod2
    week_ending::String
end

PLEXOSPeriod2(e::Node, d::AbstractDataset) =
    PLEXOSPeriod2(getchildtext("week_ending", e))

struct PLEXOSPeriod3
    month_beginning::String
end

PLEXOSPeriod3(e::Node, d::AbstractDataset) =
    PLEXOSPeriod3(getchildtext("month_beginning", e))

struct PLEXOSPeriod4
    year_ending::String
end

PLEXOSPeriod4(e::Node, d::AbstractDataset) =
    PLEXOSPeriod4(getchildtext("year_ending", e))

struct PLEXOSPeriod6
    day::Int
    datetime::String
end

PLEXOSPeriod6(e::Node, d::AbstractDataset) =
    PLEXOSPeriod6(getchildint("day_id", e),
                  getchildtext("datetime", e))

struct PLEXOSPeriod7
    quarter_beginning::String
end

PLEXOSPeriod7(e::Node, d::AbstractDataset) =
    PLEXOSPeriod7(getchildtext("quarter_beginning", e))

struct PLEXOSPhase1

    period::Int # LT period
    interval::PLEXOSPeriod0

    function PLEXOSPhase1(e::Node, d::AbstractDataset)

        period = getchildint("period_id", e)

        if length(d.intervals) > 0
            new(period, d.intervals[getchildint("interval_id", e)])
        else
            new(period)
        end

    end

end

struct PLEXOSPhase2

    period::Int # PASA period
    interval::PLEXOSPeriod0

    function PLEXOSPhase2(e::Node, d::AbstractDataset)

        period = getchildint("period_id", e)

        if length(d.intervals) > 0
            new(period, d.intervals[getchildint("interval_id", e)])
        else
            new(period)
        end

    end

end

struct PLEXOSPhase3

    period::Int # MT period
    interval::PLEXOSPeriod0

    function PLEXOSPhase3(e::Node, d::AbstractDataset)

        period = getchildint("period_id", e)

        if length(d.intervals) > 0
            new(period, d.intervals[getchildint("interval_id", e)])
        else
            new(period)
        end

    end

end

struct PLEXOSPhase4 # are these values always 1:1 matches?

    period::Int # ST period
    interval::PLEXOSPeriod0

    function PLEXOSPhase4(e::Node, d::AbstractDataset)

        period = getchildint("period_id", e)

        if length(d.intervals) > 0
            new(period, d.intervals[getchildint("interval_id", e)])
        else
            new(period)
        end

    end

end

struct PLEXOSKey
    membership::PLEXOSMembership
    model::PLEXOSModel
    phase::Int
    property::PLEXOSProperty
    periodtype::Int # not accurate, use KeyIndex.periodtype instead
    band::Int
    sample::PLEXOSSample
    timeslice::PLEXOSTimeslice
end

PLEXOSKey(e::Node, d::AbstractDataset) =
    PLEXOSKey(d.memberships[getchildint("membership_id", e)],
              d.models[getchildint("model_id", e)],
              getchildint("phase_id", e),
              d.properties[getchildint("property_id", e)],
              getchildint("period_type_id", e),
              getchildint("band_id", e),
              d.samples[getchildint("sample_id", e) + sample_offset],
              d.timeslices[getchildint("timeslice_id", e) + 1])

struct PLEXOSKeyIndex
    key::PLEXOSKey
    periodtype::Int
    position::Int # bytes from binary file start
    length::Int # 8-byte (Float64) increments
    periodoffset::Int # temporal data offset (if any) in stored times
end

PLEXOSKeyIndex(e::Node, d::AbstractDataset) =
    PLEXOSKeyIndex(d.keys[getchildint("key_id", e)],
                   getchildint("period_type_id", e),
                   getchildint("position", e),
                   getchildint("length", e),
                   getchildint("period_offset", e))


text:
#  Copyright 2019, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See https://github.com/jump-dev/JuMP.jl
#############################################################################

"""
    lp_rhs_perturbation_range(constraint::ConstraintRef;
                              feasibility_tolerance::Float64)
                              ::Tuple{Float64, Float64}

Gives the range by which the rhs coefficient can change and the current LP basis
remains feasible, i.e., where the shadow prices apply.

## Notes
- The rhs coefficient is the value right of the relation, i.e., b for the constraint when of the form a*x □ b, where □ is ≤, =, or ≥.
- The range denotes valid changes, e.g., for a*x <= b + Δ, the LP basis remains feasible for all Δ ∈ [l, u].
- `feasibility_tolerance` is the primal feasibility tolerance, this should
  preferably match the tolerance used by the solver. The default tolerance should
  however apply in most situations (c.f. "Computational Techniques of the
  Simplex Method" by István Maros, section 9.3.4).
"""
function lp_rhs_perturbation_range(
    constraint::ConstraintRef{Model,<:_MOICON};
    feasibility_tolerance::Float64 = 1e-8,
)
    return error(
        "The perturbation range of rhs is not defined or not implemented for this type " *
        "of constraint.",
    )
end

"""
Builds the standard form constraint matrix, variable bounds, and a vector of constraint reference, one for each row.
If Ã is the current constraint matrix, and the corresponding feasible set is of the form
s.t. Row_L <= Ãx <= Row_U,
     Var_L <=  x <= Var_U.
(Equality constraint is interpreted as, e.g., Row_L_i == Row_U_i,
and single sided constraints with an infinite value in the free directions)
Then the function returns a tuple with
    A = [Ã -I]
    var_lower = [Var_L, Row_L]
    var_upper = [Var_U, Row_U]
    affine_constraints::Vector{ConstraintRef} - references corresponding to Row_L <= Ãx <= Row_U
This represents the LP problem feasible set of the form
s.t.            Ax == 0
  var_lower <=   x <= var_upper
"""
function _std_matrix(model::Model)
    con_types = list_of_constraint_types(model)
    con_func_type = first.(con_types)
    if !all(
        broadcast(<:, AffExpr, con_func_type) .|
        broadcast(<:, VariableRef, con_func_type),
    )
        error(
            "The requested operation is not supported because the problem is not a linear optimization problem.",
        )
    end
    con_set_type = last.(con_types)
    if !all(
        broadcast(
            <:,
            con_set_type,
            Union{MOI.LessThan,MOI.GreaterThan,MOI.EqualTo,MOI.Interval},
        ),
    )
        error(
            "The requested operation is not supported because the problem is not a linear optimization problem.",
        )
    end
    vars = all_variables(model)
    var_to_idx = Dict{VariableRef,Int}()
    for (j, var) in enumerate(vars)
        var_to_idx[var] = j
    end
    col_j = Int[]
    row_i = Int[]
    coeffs = Float64[]

    var_lower = fill(-Inf, length(vars))
    var_upper = fill(Inf, length(vars))
    affine_constraints = ConstraintRef[]

    cur_row_idx = 0

    for (F, S) in con_types
        constraints = all_constraints(model, F, S)
        if F <: VariableRef
            for constraint in constraints
                con_obj = constraint_object(constraint)
                j = var_to_idx[con_obj.func]
                range = MOI.Interval(con_obj.set)
                var_lower[j] = max(var_lower[j], range.lower)
                var_upper[j] = min(var_upper[j], range.upper)
            end
            #issue 1892 makes it hard to know how to handle variable constraints, atm. assume unique variable constraints.
        else
            for constraint in constraints
                cur_row_idx += 1
                con_obj = constraint_object(constraint)
                con_terms = con_obj.func.terms
                push!(affine_constraints, constraint)
                for (var, coeff) in con_terms
                    push!(col_j, var_to_idx[var])
                    push!(row_i, cur_row_idx)
                    push!(coeffs, coeff)
                end
                push!(col_j, cur_row_idx + length(vars))
                push!(row_i, cur_row_idx)
                push!(coeffs, -1.0)

                range = MOI.Interval(con_obj.set)
                push!(var_lower, range.lower)
                push!(var_upper, range.upper)
            end
        end
    end

    A = sparse(row_i, col_j, coeffs, cur_row_idx, cur_row_idx + length(vars))
    return A, var_lower, var_upper, affine_constraints
end

"""
Returns the optimal primal values for the problem in standard form, c.f. `_std_matrix(model)`.
"""
function _std_primal_solution(model::Model, constraints::Vector{ConstraintRef})
    vars = all_variables(model)
    return [value.(vars); value.(constraints)]
end

"""
Returns the basis status of all variables of the problem in standard form, c.f. `_std_matrix(model)`.
"""
function _std_basis_status(model::Model)::Vector{MOI.BasisStatusCode}
    con_types = list_of_constraint_types(model)
    vars = all_variables(model)
    var_to_idx = Dict{VariableRef,Int}()
    for (j, var) in enumerate(vars)
        var_to_idx[var] = j
    end
    basis_status = fill(MOI.BASIC, length(vars))

    for (F, S) in con_types
        constraints = all_constraints(model, F, S)
        if F <: VariableRef
            for constraint in constraints
                con_obj = constraint_object(constraint)
                j = var_to_idx[con_obj.func]
                try
                    basis_status[j] =
                        MOI.get(model, MOI.ConstraintBasisStatus(), constraint)
                catch e
                    error(
                        "The perturbation range of rhs is not available because the solver doesn't " *
                        "support ´ConstraintBasisStatus´.",
                    )
                end
                if S <: MOI.LessThan && basis_status[j] == MOI.NONBASIC
                    basis_status[j] = MOI.NONBASIC_AT_UPPER
                elseif S <: MOI.GreaterThan && basis_status[j] == MOI.NONBASIC
                    basis_status[j] = MOI.NONBASIC_AT_LOWER
                end
                # ´_std_matrix´ checks that no invalid sets are used.
            end
            #issue 1892 makes it hard to know how to handle variable constraints, atm. assume unique variable constraints.
        else
            for constraint in constraints
                con_basis_status =
                    MOI.get(model, MOI.ConstraintBasisStatus(), constraint)
                if S <: MOI.LessThan && con_basis_status == MOI.NONBASIC
                    con_basis_status = MOI.NONBASIC_AT_UPPER
                elseif S <: MOI.GreaterThan && con_basis_status == MOI.NONBASIC
                    con_basis_status = MOI.NONBASIC_AT_LOWER
                end
                push!(basis_status, con_basis_status)
            end
        end
    end
    return basis_status
end

function lp_rhs_perturbation_range(
    constraint::ConstraintRef{Model,_MOICON{F,S}};
    feasibility_tolerance::Float64 = 1e-8,
)::Tuple{
    Float64,
    Float64,
} where {
    S<:Union{MOI.LessThan,MOI.GreaterThan,MOI.EqualTo},
    T,
    F<:Union{MOI.ScalarAffineFunction{T},MOI.SingleVariable},
}
    @warn(
        "This function has been deprecated and will be removed in a future " *
        "release. Use `lp_sensitivity_report(model)` instead.",
        maxlog = 1,
    )
    model = owner_model(constraint)
    if termination_status(model) != MOI.OPTIMAL
        error(
            "The perturbation range of rhs is not available because the current solution " *
            "is not optimal.",
        )
    end

    try
        con_status = MOI.get(model, MOI.ConstraintBasisStatus(), constraint)
        # The constraint is inactive.
        if con_status == MOI.BASIC
            if S <: MOI.LessThan
                return value(constraint) -
                       constraint_object(constraint).set.upper,
                Inf
            end
            if S <: MOI.GreaterThan
                return -Inf,
                value(constraint) - constraint_object(constraint).set.lower
            end
            return 0.0, 0.0
        end
    catch e
        error(
            "The perturbation range of rhs is not available because the solver doesn't " *
            "support ´ConstraintBasisStatus´.",
        )
    end
    # The constraint is active.
    # TODO: This could be made more efficient because actually we only need the basic columns.
    A, L, U, constraints = _std_matrix(model)
    var_basis_status = _std_basis_status(model)
    x = _std_primal_solution(model, constraints)
    # Ax = 0 => B*x_B = -N x_N =: b_N
    # Perturb row i by Δ => ̃x_B = B^-1 (b_N + Δ e_i) = x_B + Δ ρ
    # Perturb var bound j by Δ => ̃x_B = B^-1 (b_N - Δ N_j) = x_B + Δ ρ
    # Find min and max values of Δ : L_B <= ̃x_B <= U_B
    # L_B <= x_B + Δ ρ <= U_B  -> L_B - x_B <= Δ ρ <= U_B - x_B
    # Δ <= (U_B - x_B)_+ ./ ρ_+,  (L_B - x_B)_- ./ p_-
    # Δ >= (U_B - x_B)_- ./ ρ_-,  (L_B - x_B)_+ ./ p_+
    basic = var_basis_status .== Ref(MOI.BASIC)
    B = A[:, basic]
    x_B = x[basic]
    local rho::Vector{Float64}
    if F <: MOI.SingleVariable
        var = constraint_object(constraint).func
        j = findfirst(isequal(var), all_variables(model))
        rho = -B \ Vector(A[:, j])
    else
        i = findfirst(isequal(constraint), constraints)
        ei = zeros(size(A, 1))
        ei[i] = 1.0
        rho = B \ ei
    end
    UmX_B = U[basic] - x_B
    LmX_B = L[basic] - x_B
    pos = rho .> 0.0
    neg = rho .< 0.0
    # Find the first basic variable bound that is strictly violated.
    lower_bounds_delta = [-Inf; UmX_B[neg] ./ rho[neg]; LmX_B[pos] ./ rho[pos]]
    lower_bounds_delta_strict =
        lower_bounds_delta + [
            0.0
            feasibility_tolerance ./ rho[neg]
            -feasibility_tolerance ./ rho[pos]
        ]
    upper_bounds_delta = [Inf; UmX_B[pos] ./ rho[pos]; LmX_B[neg] ./ rho[neg]]
    upper_bounds_delta_strict =
        upper_bounds_delta + [
            0.0
            feasibility_tolerance ./ rho[pos]
            -feasibility_tolerance ./ rho[neg]
        ]

    lower_max_idx = last(findmax(lower_bounds_delta_strict))
    upper_min_idx = last(findmin(upper_bounds_delta_strict))
    # Return the unperturbed values
    return lower_bounds_delta[lower_max_idx], upper_bounds_delta[upper_min_idx]
end

"""
Returns the optimal reduced costs for the variables of the problem in standard form, c.f. `_std_matrix(model)`
"""
function _std_reduced_costs(model::Model, constraints::Vector{ConstraintRef})
    con_types = list_of_constraint_types(model)
    vars = all_variables(model)
    var_to_idx = Dict{VariableRef,Int}()
    for (j, var) in enumerate(vars)
        var_to_idx[var] = j
    end
    y = [zeros(length(vars)); dual.(constraints)]

    for (F, S) in con_types
        if F <: VariableRef
            constraints = all_constraints(model, F, S)
            for constraint in constraints
                con_obj = constraint_object(constraint)
                j = var_to_idx[con_obj.func]
                y[j] += dual(constraint)
            end
            #issue 1892 makes is hard to know how to handle variable constraints, atm. assume unique variable constraints
        end
    end

    if objective_sense(model) == MOI.MAX_SENSE
        y .*= -1.0
    end
    return y
end

"""
    lp_objective_perturbation_range(var::VariableRef;
                                    optimality_tolerance::Float64)
                                    ::Tuple{Float64, Float64}

Gives the range by which the cost coefficient can change and the current LP basis
remains optimal, i.e., the reduced costs remain valid.

## Notes
- The range denotes valid changes, Δ ∈ [l, u], for which cost[var] += Δ do not
  violate the current optimality conditions.
- `optimality_tolerance` is the dual feasibility tolerance, this should
  preferably match the tolerance used by the solver. The defualt tolerance should
  however apply in most situations (c.f. "Computational Techniques of the
  Simplex Method" by István Maros, section 9.3.4).
"""
function lp_objective_perturbation_range(
    var::VariableRef;
    optimality_tolerance::Float64 = 1e-8,
)::Tuple{Float64,Float64}
    @warn(
        "This function has been deprecated and will be removed in a future " *
        "release. Use `lp_sensitivity_report(model)` instead.",
        maxlog = 1,
    )
    model = owner_model(var)
    if termination_status(model) != MOI.OPTIMAL
        error(
            "The perturbation range of the objective is not available because the current solution " *
            "is not optimal.",
        )
    end
    if objective_sense(model) == MOI.FEASIBILITY_SENSE
        error(
            "The perturbation range of the objective is not applicable on feasibility problems.",
        )
    end
    if !has_duals(model)
        error(
            "The perturbation range of the objective is not available because no dual result is " *
            "available.",
        )
    end
    # TODO: This could be made more efficient for non-basic variables by checking them
    #       before building calling ´_std_matrix´ and ´_std_reduced_costs´ and only check
    #       the reduced cost of that variable. (issue 1892 makes this somewhat non-trivial)

    # The variable is in the basis, for a minimization problem without variable upper bounds we need:
    # c_N - N^T B^-T (c_B + Δ e_j ) >= 0
    # c_red - Δ N^T (B^-T e_j) = c_red - Δ N^T ρ = c_red - Δ N_red >= 0
    # c_red >= Δ N_red
    # maximum( c_red_- ./ N_red_- ) <= Δ <= minimum( c_red_+ ./ N_red_+ )
    # If upper bounds are present (and active), those inequalities are flipped.
    # If the problem is a maximization problem all inequalities are flipped.
    vars = all_variables(model)
    j = findfirst(isequal(var), vars)

    A, var_lower, var_upper, constraints = _std_matrix(model)
    var_basis_status = _std_basis_status(model)
    c_red = _std_reduced_costs(model, constraints)

    if var_basis_status[j] != MOI.BASIC
        if var_lower[j] == var_upper[j]
            return -Inf, Inf
        end
        if (
            var_basis_status[j] == MOI.NONBASIC_AT_LOWER &&
            objective_sense(model) == MOI.MIN_SENSE
        ) || (
            var_basis_status[j] == MOI.NONBASIC_AT_UPPER &&
            objective_sense(model) == MOI.MAX_SENSE
        )
            return -c_red[j], Inf
        end
        return -Inf, -c_red[j]
    end

    basic = var_basis_status .== Ref(MOI.BASIC)
    upper = var_basis_status[.!basic] .== Ref(MOI.NONBASIC_AT_UPPER)
    ej = zeros(size(A, 1))
    ej[findfirst(isequal(var), vars[basic[1:length(vars)]])] = 1.0
    N_red = A[:, .!basic]' * (A[:, basic]' \ ej)
    c_red = c_red[.!basic]
    pos = N_red .> 0.0
    neg = N_red .< 0.0

    in_lb = (neg .& .!upper) .| (pos .& upper)
    in_ub = (pos .& .!upper) .| (neg .& upper)
    if objective_sense(model) == MOI.MAX_SENSE
        in_lb, in_ub = in_ub, in_lb
    end
    # The reduced cost of variables fixed at equality do not be accounted for (they only change binding direction).
    unfixed_vars = (var_lower.<var_upper)[.!basic]
    in_lb .&= unfixed_vars
    in_ub .&= unfixed_vars
    # Find the first reduced cost that is strictly violated
    lower_bounds_delta = [-Inf; c_red[in_lb] ./ N_red[in_lb]]
    lower_bounds_delta_strict =
        lower_bounds_delta - [0.0; optimality_tolerance ./ abs.(N_red[in_lb])]
    upper_bounds_delta = [Inf; c_red[in_ub] ./ N_red[in_ub]]
    upper_bounds_delta_strict =
        upper_bounds_delta + [0.0; optimality_tolerance ./ abs.(N_red[in_ub])]

    lower_max_idx = last(findmax(lower_bounds_delta_strict))
    upper_min_idx = last(findmin(upper_bounds_delta_strict))
    return lower_bounds_delta[lower_max_idx], upper_bounds_delta[upper_min_idx]
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module LuxCat

using Luxor

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

include("../docs/make.jl")

@testset "JuliaFEM.jl" begin
    @testset "test_dirichlet.jl" begin
        include("test_dirichlet.jl")
    end
    @testset "test_elasticity_1d.jl" begin
        include("test_elasticity_1d.jl")
    end
    @testset "test_elasticity_2d_linear_with_surface_load.jl" begin
        include("test_elasticity_2d_linear_with_surface_load.jl")
    end
    @testset "test_elasticity_2d_nonhomogeneous_boundary_conditions.jl" begin
        include("test_elasticity_2d_nonhomogeneous_boundary_conditions.jl")
    end
    @testset "test_elasticity_2d_nonlinear_with_surface_load.jl" begin
        include("test_elasticity_2d_nonlinear_with_surface_load.jl")
    end
    @testset "test_elasticity_2d_plane_stress_stiffness_matrix.jl" begin
        include("test_elasticity_2d_plane_stress_stiffness_matrix.jl")
    end
    @testset "test_elasticity_2d_residual.jl" begin
        include("test_elasticity_2d_residual.jl")
    end
    @testset "test_elasticity_3d_linear_with_surface_load.jl" begin
        include("test_elasticity_3d_linear_with_surface_load.jl")
    end
    @testset "test_elasticity_3d_nonlinear_with_surface_load.jl" begin
        include("test_elasticity_3d_nonlinear_with_surface_load.jl")
    end
    @testset "test_elasticity_3d_unit_block.jl" begin
        include("test_elasticity_3d_unit_block.jl")
    end
    @testset "test_elasticity_forwarddiff.jl" begin
        include("test_elasticity_forwarddiff.jl")
    end
    @testset "test_elasticity_hollow_sphere_with_surface_pressure.jl" begin
        include("test_elasticity_hollow_sphere_with_surface_pressure.jl")
    end
    @testset "test_elasticity_med_pyr5_point_load.jl" begin
        include("test_elasticity_med_pyr5_point_load.jl")
    end
    @testset "test_elasticity_plane_strain.jl" begin
        include("test_elasticity_plane_strain.jl")
    end
    @testset "test_elasticity_pyr5_point_load.jl" begin
        include("test_elasticity_pyr5_point_load.jl")
    end
    @testset "test_elasticity_tet4_stiffness_matrix.jl" begin
        include("test_elasticity_tet4_stiffness_matrix.jl")
    end
    @testset "test_elasticity_tet10_mass_matrix.jl" begin
        include("test_elasticity_tet10_mass_matrix.jl")
    end
    @testset "test_elasticity_tet10_stiffness_matrix.jl" begin
        include("test_elasticity_tet10_stiffness_matrix.jl")
    end
    @testset "test_elasticity_tetra.jl" begin
        include("test_elasticity_tetra.jl")
    end
    @testset "test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl" begin
        include("test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl")
    end
    @testset "test_elasticplastic_3d_linear_with_surface_load.jl" begin
        include("test_elasticplastic_3d_linear_with_surface_load.jl")
    end
    @testset "test_heat_2d_one_element.jl" begin
        include("test_heat_2d_one_element.jl")
    end
    @testset "test_heat_3d.jl" begin
        include("test_heat_3d.jl")
    end
    @testset "test_heat_tet10_convection.jl" begin
        include("test_heat_tet10_convection.jl")
    end
    @testset "test_heat_3d_2.jl" begin
        include("test_heat_3d_2.jl")
    end
    @testset "test_heat.jl" begin
        include("test_heat.jl")
    end
    @testset "test_heat_2.jl" begin
        include("test_heat_2.jl")
    end
    @testset "test_heat_3.jl" begin
        include("test_heat_3.jl")
    end
    @testset "test_heat_3d_two_rings.jl" begin
        include("test_heat_3d_two_rings.jl")
    end
    @testset "test_heat_4.jl" begin
        include("test_heat_4.jl")
    end
    @testset "test_modal_analysis.jl" begin
        include("test_modal_analysis.jl")
    end
    @testset "test_modal_analysis_elasticity.jl" begin
        include("test_modal_analysis_elasticity.jl")
    end
    @testset "test_modal_analysis_elasticity_2.jl" begin
        include("test_modal_analysis_elasticity_2.jl")
    end
    @testset "test_modal_analysis_zero_eigenmodes.jl" begin
        include("test_modal_analysis_zero_eigenmodes.jl")
    end
    @testset "test_mortar.jl" begin
        include("test_mortar.jl")
    end
    @testset "test_mortar_2d.jl" begin
        include("test_mortar_2d.jl")
    end
    @testset "test_mortar_2d_assembly.jl" begin
        include("test_mortar_2d.jl")
    end
    @testset "test_mortar_2d_contact.jl" begin
        include("test_mortar_2d_contact.jl")
    end
    @testset "test_mortar_2d_mesh_tie.jl" begin
        include("test_mortar_2d_mesh_tie.jl")
    end
    @testset "test_mortar_2d_weighted_gap.jl" begin
        include("test_mortar_2d_weighted_gap.jl")
    end
    @testset "test_mortar_3d_mesh_tie_modal.jl" begin
        include("test_mortar_3d_mesh_tie_modal.jl")
    end
    @testset "test_mortar_3d_mesh_tie_two_rings.jl" begin
        include("test_mortar_3d_mesh_tie_two_rings.jl")
    end
    @testset "test_mortar_3d_polygon_clip.jl" begin
        include("test_mortar_3d_polygon_clip.jl")
    end
    @testset "test_postprocess.jl" begin
        include("test_postprocess.jl")
    end
    @testset "test_potential_energy.jl" begin
        include("test_potential_energy.jl")
    end
    @testset "test_problems_contact_3d.jl" begin
        include("test_problems_contact_3d.jl")
    end
    @testset "test_problems_elasticity.jl" begin
        include("test_problems_elasticity.jl")
    end
    @testset "test_problems_mortar_3d.jl" begin
        include("test_problems_mortar_3d.jl")
    end
    @testset "test_problems_mortar_3d_lowlevel.jl" begin
        include("test_problems_mortar_3d_lowlevel.jl")
    end
    @testset "test_solvers_postprocess.jl" begin
        include("test_solvers_postprocess.jl")
    end
    @testset "test_virtual_work.jl" begin
        include("test_virtual_work.jl")
    end
    @testset "test_von_mises_material.jl" begin
        include("test_von_mises_material.jl")
    end
end

include("../docs/deploy.jl")


text:
using ArgParse

function gt0(x)
    if x> 0
        return true
    else
        return false
    end
end

function leq1gt0(x)
    if x <= 1 && x > 0
        return true
    else
        return false
    end
end

function parse_commandline(a="")
    s = ArgParseSettings()

    @add_arg_table s begin
        ## Parameters
        "--vertices", "-n"
            help = "Total number of vertices"
            arg_type = Int
            default = 10
            range_tester = gt0
        "--edges", "-e"
            help = "Total number of edges"
            arg_type = Int
            default = 9
            range_tester = gt0
        "--iterations", "-i"
            help = "Total number of iterations"
            arg_type = Int
            default = 5000
            range_tester = gt0
        "--edge_rewire", "-r"
            help = "Percent of edges to be rewired"
            arg_type = Float64
            default = 0.75
            range_tester = leq1gt0
        "--temperature", "-t"
            help = "Annealing temperature"
            arg_type = Float64
            default = 1.0
        "--simple_graph"
            help = "use SimpleGraph structure"
            action = :store_true
        "--sparse_rep"
            help = "Use sparse representation"
            action = :store_true

        ## Output config
        "--out_folder"
            help = "name of the folder in which to store output. Will default to res_yyyy_mm_dd_HH_MM:SS"
            arg_type = ASCIIString
            default = ""
        "--init_fname"
            help = "name for the .dot or .png of the original graph to be saved as"
            arg_type = ASCIIString
            default = "initial"
        "--final_fname"
            help = "name for the .dot or .png of the final graph to be saved as"
            arg_type = ASCIIString
            default = "final"
        "--no_figs"
            help = "Should summary figures be produced"
            action = :store_true
        "--plot_graphs"
            help = "Produce GraphViz plots of graphs"
            action = :store_true
        "--save_graphs_plot"
            help = "Produce GraphViz plots of graphs"
            action = :store_true
        "--no_graph_dot"
            help = "Don't write graphs to .dot file"
            action = :store_true
        "--no_console_output"
            help = "not print ot console"
            action = :store_true
    end

    if isequal(a,"")
        return parse_args(s)
    else
        return parse_args(a,s)
    end
end


text:
function _adjust_web_svg_dim(num_str::AbstractString)
    num = parse(Float64, String(num_str[begin:end - 2]))
    return @sprintf("%.6f", num / 10) * "rem"
end

function _adjust_web_svg_display(svg_str::String)
    svgobject = parsexml(svg_str)
    svgroot = root(svgobject)

    height = svgroot["height"]
    svgroot["height"] = _adjust_web_svg_dim(height)
    delete!(svgroot, "width")

    svgroot["class"] = "latexsvg-display"
    svgroot["style"] = "display: block; margin: auto; max-width: 80%"

    return string(svgobject)
end

function _adjust_web_svg_inline(svg_str::String)
    svgobject = parsexml(svg_str)
    svgroot = root(svgobject)

    height = svgroot["height"]
    svgroot["height"] = _adjust_web_svg_dim(height)
    delete!(svgroot, "width")

    svgroot["class"] = "latexsvg-inline"
    svgroot["style"] = "vertical-align: middle"

    return string(svgobject)
end

function Base.show(io::IO, ::MIME"text/html", svg::LaTeXSVG)
    write(io, "
\n")
    write(io, _adjust_web_svg_display(svg.svg))
    write(io, "
")
end

function Base.show(io::IO, ::MIME"juliavscode/html", svg::LaTeXSVG)
    write(io, """

    """)
    write(io, "
\n")
    write(io, _adjust_web_svg_display(svg.svg))
    write(io, "
")
end