text:
@testset "Interruptible DCPLossLess" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.InterruptiblePowerLoad, PM.DCPlosslessForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
    # No parameters testing
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range; parameters = false);
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.InterruptiblePowerLoad, PM.DCPlosslessForm, sys5b, time_range; parameters = false);
    @test !(:params in keys(ps_model.JuMPmodel.ext))
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
end

@testset "Interruptible ACP" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.InterruptiblePowerLoad, PM.StandardACPForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
    # No Parameters Testing
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range; parameters = false);
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.InterruptiblePowerLoad, PM.StandardACPForm, sys5b, time_range; parameters = false);
    @test !(:params in keys(ps_model.JuMPmodel.ext))
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
end

@testset "Static Load DCPLossLess" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PM.DCPlosslessForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
    # No Parameters Testing
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range; parameters = false);
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PM.DCPlosslessForm, sys5b, time_range; parameters = false);
    @test !(:params in keys(ps_model.JuMPmodel.ext))
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
end

@testset "Static Load ACP" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PM.StandardACPForm, sys5b, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
    # No Parameters Testing
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range; parameters = false);
    PSI.construct_device!(ps_model, PSY.PowerLoad, PSI.StaticPowerLoad, PM.StandardACPForm, sys5b, time_range; parameters = false);
    @test !(:params in keys(ps_model.JuMPmodel.ext))
    @test JuMP.num_variables(ps_model.JuMPmodel) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 0
end


text:
export Trouve_IntervalleA
function Trouve_IntervalleA(h :: AbstractNLPModel,
                            nlpstop :: AbstractStopping;
                            t₀ = h.meta.x0, tₘ = 100.0,
                            direction =  :Nwt, verbose = false, kwargs...)

  tim1 = t₀
  ti = (tim1 .+ tₘ) ./ 2

  i = 0

  h₀ = obj(h, t₀); dh₀ = grad(h, t₀)

  him1 = obj(h, tim1); dhim1 = grad(h, tim1)
  hi = obj(h, ti); dhi = grad(h, ti)

  OK = update_and_start!(nlpstop, x = ti, fx = hi, gx = dhi, g0 = copy(dh₀))

  verbose &&
    @printf("iter tim1        dhim1        him1         ti        ")
  verbose &&
    @printf("dhi        hi\n")
  verbose &&
    @printf("%4d %7.2e %7.2e  %7.2e  %7.2e  %7.2e  %7.2e \n", i,  tim1[1], dhim1[1], him1[1], ti[1], dhi[1], hi[1])

  while !OK
    him1 = hi
    hi = obj(h, ti)
    update!(nlpstop.current_state, fx = hi)
    if true in (hi .> h₀) || true in ((hi  .> him1) && (i > 1))
      optimal, nlpstop  = zoom_generic(h, nlpstop, tim1, ti;
                                       direction = direction, verbose = verbose)
      return optimal, nlpstop
    end

    dhim1 = dhi; dhi = grad(h, ti)
    update!(nlpstop.current_state, gx = dhi)

    if (dhi .>= 0.0)
      optimal, nlpstop = zoom_generic(h, nlpstop, ti, tim1, verbose = verbose, direction = direction)
      return optimal, nlpstop
    end

    tim1 = ti
    ti = (tim1 .+ tₘ) ./ 2

    OK = update_and_stop!(nlpstop, x = ti)

    i += 1

    verbose &&
      @printf("%4d %7.2e %7.2e  %7.2e  %7.2e  %7.2e  %7.2e \n", i,  tim1[1], dhim1[1], him1[1], ti[1], dhi[1], hi[1])
  end

  return OK, nlpstop
end

ext: jl


text:
"""
This is a sample input file for the CloudScat code.  To use it run julia
and type

julia> include("ozone_absorption.jl"); run()

The output file will be based on this file's name so to run different simulations
you can just copy this file, change the configuration and run it as above with
the new file name.  You can also run from the command line as e.g.

> julia ozone_absorption.jl

- - -
2020 Alejandro Luque IAA-CSIC (aluque@iaa.es)

"""

using CloudScat
using DataFrames
using DelimitedFiles

const co = CloudScat.constants

# This example depends on the ozone profiles used by MODTRAN, which cannot
# be redistributed.  If you have working MODTRAN installation, modify this
# to point to the folder with the profiles.
const MODTRAN_PROFILES = joinpath(@__DIR__, "profiles")

function run()
    profiles = ["MidLat_Summer.dat",
                "MidLat_Winter.dat",
                "Tropical.dat",
                "US_Std_1976.dat"]

    for prof in profiles
        run_profile(prof, "Hale.dat")
    end

    run_profile("Tropical.dat", "WarrenBrandt.dat")
end

function run_profile(fname, refindex)
    df = load_profile(joinpath(MODTRAN_PROFILES, fname))
    @info "Profile $(fname) loaded"

    # This is the absorption cross-section for 200 nm, where ozone absorption
    # is lowest.  Obtained from
    # L.T Molina and M.J. Molina, "Absolute absorption cross sections of
    # ozone in the 185- to 350-nm wavelength range,"
    # J. Geophys. Res. 91, 14501-14508 (1986).
    σ = 3.145E-19 * co.centi^2

    # The profiles reach only 100 km but the observer is above this altitude
    # We simply add another point with a negligible O3 density.
    z = [df.Z; 1000.0]
    dens = [df.O3; 0.0]
    T = [df.T; df.T[end]]

    absorb = StratifiedAbsorption(σ,
                                  (@. z * co.kilo),
                                  (@. dens * 1e-5 * co.atm / co.k / T))

    params = Params(
        # NUmber of simulated photons
        N = 10000000,

        # The wavelength of the scattered light
        λ = 200 * co.nano,

        # The source is defined as a straight line running from source_a
        # to source_b.  Within this line photons are distributed uniformly with
        # initially random, isotropic directions.  To use a point source set
        # source_b = source_a
        source_a = [0, 0, 10 * co.kilo],
        source_b = [0, 0, 10 * co.kilo],

        # Minimum weight before a particle is discarded
        weight_min = 1e-7)

    # World geometry
    # Describe the cloud geometry here.  It may consist in combinations of
    # geometrical figures (cylinders, spheres, cones, semi-planes).
    # See cloud_gemoetry.jl for a more detailed description
    #
    # A cylinder is defined as Cylinder(zbottom, ztop, xcenter, ycenter, radius)
    c1 = Cylinder(7 * co.kilo,  15 * co.kilo, 0, 0, 200 * co.kilo)
    cloud = c1

    # Geometry of the full domain.  Photons outside the domain are
    # immediately discarded.
    domain = Cylinder(7 * co.kilo, 60 * co.kilo, 0, 0, 200 * co.kilo)

    # Set a homogeneous cloud droplet density and radius.
    # See variable_droplet_radius.jl for how to deal with inhomogeneous clouds.
    # Radius of the scattering particles
    radius = 10e-6

    # Density of the scattering centers.
    nscat = 100 * co.centi^-3

    composition = Homogeneous(params.λ, nscat, radius,
                              refindex=refindex)

    # Define the full simulation world
    world = World(cloud, domain, composition, absorb)

    # Observers record a time curve and an image of received light.  You can
    # set up as many observers as you wish but they have a significant impact
    # on performance.

    observers = [Observer(
        # Location of the observer
        position = [0, 0, 400 * co.kilo],

        # Sampling interval
        tsample = 1e-5,

        # Collected samples
        nsamples = 10000,

        # fov of the camera (half of the diagonal fov).
        fov = 40,

        # camera pixels
        pixels = 1024)]

    basename = splitext(splitdir(@__FILE__)[2])[1]
    profbasename = splitext(fname)[1]
    refindexbasename = splitext(refindex)[1]

    # The saveto option specifies the output file.  Leave it as it is to base
    # the filename on the name of this file specifying the parameters.
    CloudScat.main(params, world, observers,
                   saveto="$(basename)_$(profbasename)_$(refindexbasename).h5")
end

function load_profile(fname)
    local d, colnames

    open(fname) do fd
        for i in 1:2
            readline(fd)
        end
        colnames = map(Symbol, (split(readline(fd))))
        readline(fd)

        d = readdlm(fd)
    end
    # From atm * cm / km to m^-3
    # but let's leave this to the code using the dataframe
    # @. d[:, 5:end] = d[:, 5:end] * 1e-5 * co.atm / co.k * T

    dict = Dict(col => data for (col, data) in zip(colnames, eachcol(d)))
    df = DataFrame(;dict...)
end

# Check if the code has been 'included' or run from the shell.  In the latter
# case, run the simulation.
isinteractive() || run()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_ApiError.jl")
include("model_ApiErrorBase.jl")
include("model_InnerError.jl")
include("model_OperationStatusResponse.jl")
include("model_RunCommandDocumentBase.jl")
include("model_RunCommandInput.jl")
include("model_RunCommandInputParameter.jl")
include("model_RunCommandListResult.jl")
include("model_RunCommandParameterDefinition.jl")
include("model_RunCommandResultProperties.jl")
include("model_RunCommandDocument.jl")
include("model_RunCommandResult.jl")


text:
using Ft2k

#cd(joinpath(ENV["HOME"], "ft2k/Software/Watchman"))
#include("utils/watchman_fxns.jl") # 5.544937 seconds on macbook

# load config file (now Ft2k.conf from global)
#configuration = load_config() # 0.14 seconds

# rsync with pi server to get most updated data, NOTE: linux specific
success = Ft2k.check_server_updating(Ft2k.conf) # 3.222670 seconds on macbook

# don't need this because wrapped in "continuouspull"
tdict, fdict, ids, sqlite_db = Ft2k.get_dbvalues()

# check whether any feeders have mysteriously fallen off the grid
feedersDown = Ft2k.check_man_down(tdict)

# check whether any feedings have been missed outside given window
missedfeedings =Ft2k.check_missed_feedings(fdict, sqlite_db)

logged = []
logged = Ft2k.checkAndEmail(logged, feedersDown, missedfeedings)
# sleep for a bit


text:
# Sum over states solution for systems with few modes.

"""
Sum over states solution for a small system.
"""
struct SumOverStates <: Solution
    "Partition function."
    Z::Float64
    "Energy."
    E::Float64
    "Heat capacity."
    Cv::Float64
end

"""
    SumOverStates(sys::System, beta::Float64, basis_size::Int)

Calculate the solution for `sys` at `beta` with `basis_size` basis functions.
"""
function SumOverStates(sys::System, beta::Float64, basis_size::Int)
    basis = Basis(sys, basis_size)
    h0, V = operators(basis, sys)

    Es = eigvals(Symmetric(h0 + V))

    Z = sum(exp.(-beta * Es))
    E = sum(exp.(-beta * Es) .* Es) / Z
    Cv = sum(exp.(-beta * Es) .* (Es .- E).^2) / Z * beta^2

    SumOverStates(Z, E, Cv)
end


text:
module ReactionNetworkImporters

using DataStructures, Symbolics, SymbolicUtils, Catalyst, SparseArrays

# creates a ModelingToolkit function-like Symbol
# can then do stuff like
# @parameters t
# S₁ = funcsym(S,1)
# u = S₁(t)
function funcsym(S::Symbol, t, args...)
    S = Symbol(S,args...)
    (@variables $(S)(t))[1]
end

abstract type NetworkFileFormat end

# exported data types
#struct RSSANetwork <: NetworkFileFormat end
struct BNGNetwork <: NetworkFileFormat end
struct MatrixNetwork <: NetworkFileFormat end

struct ParsedReactionNetwork
    "Catalyst Network"
    rn::ReactionSystem

    "Initial Conditions"
    u₀

    "Parameters"
    p

    "Expressions for the Parameters"
    paramexprs

    "Dict from `Variable` in species(rn) to full string for species name"
    varstonames

    "Dict from lumped species name (as string) to group of species ids"
    groupstoids

end
ParsedReactionNetwork(rn::ReactionSystem, u₀; p=nothing, paramexprs=nothing, varstonames=nothing, groupstoids=nothing) =
                        ParsedReactionNetwork(rn, u₀, p, paramexprs, varstonames, groupstoids)

export RSSANetwork, BNGNetwork, MatrixNetwork, ParsedReactionNetwork

# parsers
#include("parsing_routines_rssafiles.jl")
include("parsing_routines_bngnetworkfiles.jl")
include("parsing_routines_matrixnetworks.jl")

export loadrxnetwork

end # module


text:
using CairoMakie
using CSV
using Colors

color_lst = [colorant"#332288" colorant"#555555" colorant"#117733" colorant"#882255"]
# Blue, Grey, Green, Red, Yellow
# Colorblind friendly color palette

f = CSV.File("HopfSampling/HopfSampling.csv")

Hopf_random_sampling_trace  = f.random_sampling_trace
Hopf_optimal_sampling_trace = f.optimal_sampling_trace
Hopf_optimal_sampling_trace_ekf = f.optimal_sampling_trace_ekf
Hopf_approx_sampling_trace  = f.approx_sampling_trace
Hopf_approx_sampling_trace_ekf = f.approx_sampling_trace_ekf
Hopf_samples = f.samples

f = CSV.File("LorenzSampling/LorenzSampling.csv")

Lorenz_random_sampling_trace = f.random_sampling_trace
Lorenz_optimal_sampling_trace = f.optimal_sampling_trace
Lorenz_optimal_sampling_trace_ekf = f.optimal_sampling_trace_ekf
Lorenz_approx_sampling_trace = f.approx_sampling_trace
Lorenz_approx_sampling_trace_ekf = f.approx_sampling_trace_ekf
Lorenz_samples = f.samples

max_steps = length(Lorenz_samples)

noto_sans = "./resources/NotoSans-Regular.ttf"
noto_sans_bold = "./resources/NotoSans-Bold.ttf"

tickfontsize    = 32
labelfontsize   = 38
basewidth = 2

f = Figure(font=noto_sans, resolution=(1200,700), figure_padding=40)

ax1 = Axis( f[1,2],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Sample", ylabel = "Tr(Σ)",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            yscale = log10,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10),
            title="Lorenz System",
            titlesize=labelfontsize)

ylims!(ax1,(1e-1,1e1))
xlims!(ax1,(0,1000))

l2 = lines!(ax1, 1:max_steps, Lorenz_optimal_sampling_trace ./ Lorenz_optimal_sampling_trace,     color=color_lst[1], linewidth=basewidth)
l3 = lines!(ax1, 1:max_steps, Lorenz_optimal_sampling_trace_ekf ./ Lorenz_optimal_sampling_trace, color=color_lst[2], linewidth=basewidth)
l1 = lines!(ax1, 1:max_steps, Lorenz_approx_sampling_trace ./ Lorenz_optimal_sampling_trace,      color=color_lst[3], linewidth=basewidth)
l4 = lines!(ax1, 1:max_steps, Lorenz_approx_sampling_trace_ekf./ Lorenz_optimal_sampling_trace,  color=color_lst[4], linewidth=basewidth)
#l5 = lines!(ax1, 1:max_steps, Lorenz_random_sampling_trace./ Lorenz_optimal_sampling_trace, color=:black)

ax2= Axis(  f[1,1],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Sample", ylabel = "Tr(Σ)",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            yscale = log10,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10),
            title="Hopf Bifurcation",
            titlesize=labelfontsize)

ylims!(ax2,(1e-1,1e1))
xlims!(ax2,(0,1000))

l2 = lines!(ax2, 1:max_steps, Hopf_optimal_sampling_trace ./ Hopf_optimal_sampling_trace,     color=color_lst[1], linewidth=basewidth)
l3 = lines!(ax2, 1:max_steps, Hopf_optimal_sampling_trace_ekf./ Hopf_optimal_sampling_trace, color=color_lst[2], linewidth=basewidth)
l1 = lines!(ax2, 1:max_steps, Hopf_approx_sampling_trace./ Hopf_optimal_sampling_trace,      color=color_lst[3], linewidth=basewidth)
l4 = lines!(ax2, 1:max_steps, Hopf_approx_sampling_trace_ekf./ Hopf_optimal_sampling_trace,  color=color_lst[4], linewidth=basewidth)

l1elem = LineElement(color=color_lst[3], linewidth=10)
l2elem = LineElement(color=color_lst[1], linewidth=10)
l3elem = LineElement(color=color_lst[2], linewidth=10)
l4elem = LineElement(color=color_lst[4], linewidth=10)

Legend(f[2,:],
       [l1elem, l4elem, l2elem, l3elem],
       ["1D Limit Set Oracle", "1D Limit Set EKF", "Dynamic Programming Oracle", "Dynamic Programming EKF"],
       orientation=:horizontal,
       nbanks=2,
       labelsize=36)

colgap!(f.layout,60)

Label(f[1,1,TopLeft()], "A", font=noto_sans_bold, textsize = 40, halign=:left, valign=:bottom, padding=(25,0,10,0))
Label(f[1,2,TopLeft()], "B", font=noto_sans_bold, textsize = 40, halign=:left, valign=:bottom, padding=(25,0,10,0))

save("out/BeliefSim.pdf",f)


text:
using DataValueArrays
using Base.Test

@testset "Broadcast" begin

A1 = rand(Int, 10)
M1 = rand(Bool, 10)
n = rand(2:5)
dims = [ rand(2:5) for i in 1:n]
A2 = rand(10, dims...)
M2 = rand(Bool, 10, dims...)
C2 = Array{Float64}(10, dims...)
i = rand(2:5)
A3 = rand(10, [dims; i]...)
M3 = rand(Bool, 10, [dims; i]...)
C3 = Array{Float64}(10, [dims; i]...)

m1 = broadcast((x,y)->x, M1, M2)
m2 = broadcast((x,y)->x, M2, M3)
R2 = reshape(Bool[ (x->(m1[x] | M2[x]))(i) for i in 1:length(M2) ], size(M2))
R3 = reshape(Bool[ (x->(m2[x] | M3[x]))(i) for i in 1:length(M3) ], size(M3))
r2 = broadcast((x,y)->x, R2, R3)
Q3 = reshape(Bool[ (x->(r2[x] | R3[x]))(i) for i in 1:length(R3) ], size(M3))

U1 = DataValueArray(A1)
U2 = DataValueArray(A2)
U3 = DataValueArray(A3)
V1 = DataValueArray(A1, M1)
V2 = DataValueArray(A2, M2)
V3 = DataValueArray(A3, M3)
Z2 = DataValueArray(Float64, 10, dims...)
Z3 = DataValueArray(Float64, 10, [dims; i]...)

f() = 5
f(x::Real) = 5 * x
f(x::Real, y::Real) = x * y
f(x::Real, y::Real, z::Real) = x * y * z

    for (dests, arrays, datavaluearrays, mask) in
    ( ((C2, Z2), (A1, A2), (U1, U2), ()),
        ((C3, Z3), (A2, A3), (U2, U3), ()),
        ((C3, Z3), (A1, A2, A3), (U1, U2, U3), ()),

        ((C2, Z2), (A1, A2), (V1, V2), (R2,)),
        ((C3, Z3), (A2, A3), (V2, V3), (R3,)),
        ((C3, Z3), (A1, A2, A3), (V1, V2, V3), (Q3,)),
)

    # Base.broadcast!(f, B::DataValueArray, As::DataValueArray...)
    broadcast!(f, dests[1], arrays...)
    # TODO Fix broadcast!(i->f.(i), dests[2], datavaluearrays...)
    # TODO Fix @test isequal(dests[2], DataValueArray(dests[1], mask...))

    # Base.broadcast(f, As::DataValueArray...)
    D = broadcast(f, arrays...)
    # TODO Fix X = broadcast(f, datavaluearrays...)
    # TODO Fix @test isequal(X, DataValueArray(D, mask...))
end

# Base.broadcast!(f, X::DataValueArray)
for (array, datavaluearray, mask) in
    ( (A1, U1, ()), (A2, U2, ()), (A3, U3, ()),
        (A1, V1, (M1,)), (A2, V2, (M2,)), (A3, V3, (M3,)),
)
    broadcast!(f, array)
    broadcast!(f, datavaluearray)
    # TODO Fix @test isequal(DataValueArray, DataValueArray(array, mask...))
end

# test broadcasted arithmetic operators
A = rand(10)
X1 = DataValueArray(A)
n = rand(2:5)
dims = rand(2:5, n)
B = rand(Float64, 10, dims...)
X2 = DataValueArray(B)
M = rand(Bool, 10, dims...)
Y = DataValueArray(B, M)

A = rand(Bool, 100)
B = rand(Bool, 100)
M1 = rand(Bool, 100)
M2 = rand(Bool, 100)
X = DataValueArray(A, M1)
Y = DataValueArray(B, M2)
# TODO Fix @test broadcast(&, X, Y) == DataValueArray(A .& B, M1 .| M2)
# TODO Fix @test broadcast(|, X, Y) == DataValueArray(A .| B, M1 .| M2)

# Test broadcasting with constructor
t = DataValueArray(rand(3))
c = DataValueArray(rand(Bool, 3))
# TODO Fix @test isequal(SurvEvent.(t, c), DataValueArray([SurvEvent(get(t[i]), get(c[i])) for i in 1:3]))
# TODO Fix @test isa(SurvEvent.(t, c), DataValueVector{SurvEvent})

end


text:
@testset "minimal.absRevBin_test" begin
    # absRevBin [empty] [entire] = [empty]
    # absRevBin [0.0, 1.0] [empty] = [empty]
    # absRevBin [0.0, 1.0] [7.0, 9.0] = [empty]
    # absRevBin [empty] [0.0, 1.0] = [empty]
    # absRevBin [-2.0, -1.0] [entire] = [empty]
    # absRevBin [1.0, 1.0] [entire] = [-1.0, 1.0]
    # absRevBin [0.0, 0.0] [entire] = [0.0, 0.0]
    # absRevBin [-1.0, -1.0] [entire] = [empty]
    # absRevBin [0x1.FFFFFFFFFFFFFp1023, 0x1.FFFFFFFFFFFFFp1023] [entire] = [-0x1.FFFFFFFFFFFFFp1023, 0x1.FFFFFFFFFFFFFp1023]
    # absRevBin [0x1p-1022, 0x1p-1022] [entire] = [-0x1p-1022, 0x1p-1022]
    # absRevBin [-0x1p-1022, -0x1p-1022] [entire] = [empty]
    # absRevBin [-0x1.FFFFFFFFFFFFFp1023, -0x1.FFFFFFFFFFFFFp1023] [entire] = [empty]
    # absRevBin [1.0, 2.0] [entire] = [-2.0, 2.0]
    # absRevBin [1.0, 2.0] [0.0, 2.0] = [1.0, 2.0]
    # absRevBin [0.0, 1.0] [-0.5, 2.0] = [-0.5, 1.0]
    # absRevBin [-1.0, 1.0] [entire] = [-1.0, 1.0]
    # absRevBin [-1.0, 0.0] [entire] = [0.0, 0.0]
    # absRevBin [0.0, infinity] [entire] = [entire]
    # absRevBin [entire] [entire] = [entire]
    # absRevBin [-infinity, 0.0] [entire] = [0.0, 0.0]
    # absRevBin [1.0, infinity] [-infinity, 0.0] = [-infinity, -1.0]
    # absRevBin [-1.0, infinity] [entire] = [entire]
    # absRevBin [-infinity, -1.0] [entire] = [empty]
    # absRevBin [-infinity, 1.0] [entire] = [-1.0, 1.0]
end


text:
module patch_test_2d_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshExportModule
import LinearAlgebra: cholesky

function q4_stress()
    println("Q4. Plane stress.")

    E = 1.0;
    nu = 1.0/3;
    alpha, beta, gamma, delta, eta, phi= 1.0/30, 1.0/34, -1.0/21, -1.0/51, -1.0/26, -1.0/35
    ux(x, y) = alpha + beta * x + gamma * y
    uy(x, y) = delta + eta * x + phi * y
    MR = DeforModelRed2DStress

    fens = FENodeSet([1.0 -0.3; 2.3 -0.3; 2.3 0.95; 1.0 0.95; 1.4 0.05; 1.9 -0.03; 1.7 0.5; 1.3 0.6])
    fes = FESetQ4([1 2 6 5; 6 2 3 7; 7 3 4 8; 8 4 1 5; 5 6 7 8])

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz, 1), 2)) # displacement field

    # Apply prescribed displacements to exterior nodes
    for i in 1:4
    	setebc!(u, [i], 1, val=ux(fens.xyz[i, :]...))
    	setebc!(u, [i], 2, val=uy(fens.xyz[i, :]...))
    end

    applyebc!(u)
    numberdofs!(u)

    material = MatDeforElastIso(MR,  0.0, E, nu, 0.0)

    femm = FEMMDeforLinear(MR, IntegDomain(fes,  GaussRule(2, 2)),  material)

    F = nzebcloadsstiffness(femm,  geom,  u)
    K = stiffness(femm,  geom,  u)
    K = cholesky(K)
    U=  K\(F)
    scattersysvec!(u, U[:])

    for i in 5:8
    	uexact = [ux(fens.xyz[i, :]...), uy(fens.xyz[i, :]...)]
    	println("u.values[$i, :] = $(u.values[i, :]), uexact = [$(uexact)]")
    end

    File =  "a.vtk"
    vtkexportmesh(File,  fes.conn,  geom.values,
    FinEtools.MeshExportModule.Q4; vectors=[("u", u.values)])

    true

end # cookstress

function q4_stress_export()
    println("Q4. Plane stress.")

    E = 1.0;
    nu = 1.0/3;
    alpha, beta, gamma, delta, eta, phi= 1.0/30, 1.0/34, -1.0/21, -1.0/51, -1.0/26, -1.0/35
    ux(x, y) = alpha + beta * x + gamma * y
    uy(x, y) = delta + eta * x + phi * y
    MR = DeforModelRed2DStress

    fens = FENodeSet([1.0 -0.3; 2.3 -0.3; 2.3 0.95; 1.0 0.95; 1.4 0.05; 1.9 -0.03; 1.7 0.5; 1.3 0.6])
    fes = FESetQ4([1 2 6 5; 6 2 3 7; 7 3 4 8; 8 4 1 5; 5 6 7 8])

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz, 1), 2)) # displacement field

    # Apply prescribed displacements to exterior nodes
    for i in 1:4
    	setebc!(u, [i], 1, val=ux(fens.xyz[i, :]...))
    	setebc!(u, [i], 2, val=uy(fens.xyz[i, :]...))
    end

    applyebc!(u)
    numberdofs!(u)

    material = MatDeforElastIso(MR,  0.0, E, nu, 0.0)

    femm = FEMMDeforLinear(MR, IntegDomain(fes,  GaussRule(2, 2)),  material)

    F = nzebcloadsstiffness(femm,  geom,  u)
    K = stiffness(femm,  geom,  u)
    K = cholesky(K)
    U=  K\(F)
    scattersysvec!(u, U[:])

    for i in 5:8
    	uexact = [ux(fens.xyz[i, :]...), uy(fens.xyz[i, :]...)]
    	println("u.values[$i, :] = $(u.values[i, :]), uexact = [$(uexact)]")
    end

    AE = AbaqusExporter("q4_stress_export");
    HEADING(AE, "q4_stress_export");
    COMMENT(AE, "");
    PART(AE, "part1");
    END_PART(AE);
    ASSEMBLY(AE, "ASSEM1");
    INSTANCE(AE, "INSTNC1", "PART1");
    NODE(AE, fens.xyz);
    ELEMENT(AE, "CPS4", "AllElements", connasarray(fes))
    NSET_NSET(AE, "clamped", 1:4)
    ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", 1.0);
    END_INSTANCE(AE);
    END_ASSEMBLY(AE);
    MATERIAL(AE, "elasticity")
    ELASTIC(AE, E, nu)
    STEP_PERTURBATION_STATIC(AE)
    BOUNDARY(AE, "ASSEM1.INSTNC1", 1:4, fill(true, 4, 2), [[ux(fens.xyz[i, :]...) for i in 1:4] [uy(fens.xyz[i, :]...) for i in 1:4]])
    END_STEP(AE)
    close(AE)

    true

end # cookstress

function allrun()
    println("#####################################################")
    println("# q4_stress ")
    q4_stress()
    println("#####################################################")
    println("# q4_stress_export ")
    q4_stress_export()
    return true
end # function allrun

end # module patch_test_2d_examples


text:
# Abstract type of a 'step size sequence'
abstract type StepSize end

# Polynomially decaying step size sequence
struct PolynomialStepSize{FT <: AbstractFloat} <: StepSize
    eta::FT
    c::FT
end
function PolynomialStepSize(eta::FT, c=one(FT)) where {FT <: AbstractFloat}
    PolynomialStepSize(eta, c)
end
@inline function get(s::PolynomialStepSize{FT}, k::Real) where {FT <: AbstractFloat}
    s.c*(k+one(FT))^(-s.eta)
end

# Specialization to the RAM default rule:
struct RAMStepSize{FT<:AbstractFloat} <: StepSize
    p::PolynomialStepSize{FT}
end
function RAMStepSize(eta::FT, d::Int) where {FT <: AbstractFloat}
    RAMStepSize{FT}(PolynomialStepSize(eta, FT(d)))
end
@inline function get(s::RAMStepSize{FT}, k::Real) where {FT}
    min(FT(0.5), get(s.p, k))
end


text:
using PyPlot

function savefig(fname; subdir = nothing)
    "figdir" in keys(ENV) || @error "Environment variable `figdir` not set."
    dir = ENV["figdir"]
    isnothing(subdir) || (dir = joinpath([dir, subdir]))
    isdir(dir) || mkpath(dir)
    plt.savefig(joinpath([dir, fname]))
end

# Hi-def ('retina') figures in notebook. [https://github.com/JuliaLang/IJulia.jl/pull/918]
function IJulia.metadata(x::PyPlot.Figure)
    w, h = (x.get_figwidth(), x.get_figheight()) .* x.get_dpi()
    return Dict("image/png" => Dict("width" => 0.5 * w, "height" => 0.5 * h))
end


text:
using AbstractAlgebra
using Primes

# TODO: add some

@testset "Hard test systems for f4" begin
    R, (x1, x2, x3, x4) = PolynomialRing(QQ, ["x1","x2","x3","x4"], ordering=:degrevlex)

    ###########################################################################

    N = prod(map(BigInt, nextprimes(2^30 + 3, 5)))
    #=
    N is
    1427247763155044947592358980815426846827008617
    =#
    system = [x1 + x2 + x3 + x4,
                x1*x2 + x1*x3 + x1*x4 + x2*x3 + x2*x4 + x3*x4,
                x1*x2*x3 + x1*x2*x4 + x1*x3*x4 + x2*x3*x4,
                x1*x2*x3*x4 + N]

    # this should take about 10 primes
    gb = Groebner.groebner(system)
    @test gb == [
        x1 + x2 + x3 + x4,
        x2^2 + x2*x3 + x3^2 + x2*x4 + x3*x4 + x4^2,
        x3^3 + x3^2*x4 + x3*x4^2 + x4^3,
        x4^4 - N
    ]

    ###########################################################################

    N = prod(map(BigInt, nextprimes(2^31 - 1, 5)))
    #=
    N is
    45671930739135174346839766056203605080877915151
    =#
    system = [x1 + x2 + x3 + x4,
                x1*x2 + x1*x3 + x1*x4 + x2*x3 + x2*x4 + x3*x4,
                x1*x2*x3 + x1*x2*x4 + x1*x3*x4 + x2*x3*x4,
                x1*x2*x3*x4 + N]

    # this should take about 10 primes
    gb = Groebner.groebner(system)
    @test gb == [
        x1 + x2 + x3 + x4,
        x2^2 + x2*x3 + x3^2 + x2*x4 + x3*x4 + x4^2,
        x3^3 + x3^2*x4 + x3*x4^2 + x4^3,
        x4^4 - N
    ]

    ###########################################################################

    N = prod(map(BigInt, nextprimes(2^31 - 1, 100)))
    #=
    N is
    155957502472941350164396756181237870921673222084334376466764837248138019
    338975770533033341246052934841394539580063717875572166517749899185157089
    854972461613034435119695325890456903607183962279960193065823268162712808
    196475765998522398673229836351207189494171911844595858978617297329048837
    407198011154883771097567254734750443080714724501847550152215111049527862
    783487269833728445656312427565443863328011390695385168351099683848940890
    157184467093413854863172748020836886918109307797366132657049418890246023
    760683982046987711591165743318455570261741083133013563013142610260813314
    163528474356815207125941871425854964734919799177276840652608662771694654
    083164269078631636104073420379406434036790229905893357639334250840650980
    175594048478512395769709879899298769488161098605323126429857397027203991
    344047973029587500619268226785604632683477043132362736563551930698622928
    5368207467244436111610932456122141350383600941573123435631247362408683
    =#
    system = [x1 + x2 + x3 + x4,
                x1*x2 + x1*x3 + x1*x4 + x2*x3 + x2*x4 + x3*x4,
                x1*x2*x3 + x1*x2*x4 + x1*x3*x4 + x2*x3*x4,
                x1*x2*x3*x4 + N]
    # around 200 primes are required
    gb = Groebner.groebner(system)
    @test gb == [
        x1 + x2 + x3 + x4,
        x2^2 + x2*x3 + x3^2 + x2*x4 + x3*x4 + x4^2,
        x3^3 + x3^2*x4 + x3*x4^2 + x4^3,
        x4^4 - N
    ]

    ###########################################################################
    N = prod(map(BigInt, primes(2^30 + 2^29 + 2^28 + 2^27 + 2^26 + 2^25 + 2^24 + 2^23 + 2^22 + 2^21 + 2^20 + 2^19 + 2^18 + 2^17 + 2^16 + 2^15 + 2^14, 2^31-1)))
    #=
    N is
    [encoded data, 7130 characters]
    =#

    # this should take A LOT of primes
    system = [x1 + x2 + x3 + x4,
                x1*x2 + x1*x3 + x1*x4 + x2*x3 + x2*x4 + x3*x4,
                x1*x2*x3 + x1*x2*x4 + x1*x3*x4 + x2*x3*x4,
                x1*x2*x3*x4 + N]
    # impossible to compute for now
    # UPD: now computes in 7 seconds
    # UPD 2: now computes in 1 second !!
    gb = Groebner.groebner(system)
    @test gb == [
        x1 + x2 + x3 + x4,
        x2^2 + x2*x3 + x3^2 + x2*x4 + x3*x4 + x4^2,
        x3^3 + x3^2*x4 + x3*x4^2 + x4^3,
        x4^4 - N
    ]

    system = [x1 - (2^31-1)*x2 - (2^30 + 3)*x3]
    @test Groebner.groebner(system) == system

end


text:
module FLOWFarm
using ForwardDiff: Iterators
using Geodesy: bound_thetad
using Geodesy; const gd = Geodesy
using ForwardDiff
using LinearAlgebra

# using CCBlade
using PyPlot; const plt = PyPlot
using FLOWMath: linear,trapz,Akima,ksmax,abs_smooth
# using Statistics
using Distributed
using YAML
using XLSX
using DataFrames
using CSV

include("io.jl")
include("utilities.jl")
include("turbines.jl")
include("windfarms.jl")
include("wind_resources.jl")
include("wind_shear_models.jl")
include("wake_combination_models.jl")
include("wake_deficit_models.jl")
include("wake_deflection_models.jl")
include("thrust_coefficient_models.jl")
include("local_turbulence_intensity_models.jl")
include("general_models.jl")
include("power_models.jl")
include("user_functions.jl")
include("optimization_functions.jl")
include("fatigue_model.jl")
include("tip.jl")
include("cost_models.jl")
include("plotting.jl")
end # module


text:
# =======================================================================
# 2UMFLP

using vOptSpecific

# ==========================================
# instance

n   = 3 # facilities
m   = 7 # customers
ca1 = [9 27 24; 34 9 34; 11 12 17; 14 9 13; 16 26 28; 82 95 49; 20 17 9]
ca2 = [34 9 34; 14 9 13; 20 17 9; 9 27 24; 11 12 17; 16 26 28; 82 95 49]
cr1 = [30, 30, 30]
cr2 = [12, 45, 25]

# ==========================================
# vOptSpecific : test 0 : UMFLP_Delmee2017

myInstance = set2UMFLP(  m , n, ca1, ca2, cr1, cr2 )
algorithm = UMFLP_Delmee2017()
Z1, Z2, Y, X, isEdge = vSolve( myInstance , algorithm )

for i = 1:length(z1)
    println("(" , Z1[i], " | ", Z2[i], ") ")
end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
#=
Converting between 3 main temperature scales: F, C, K

main()
{
  Declare variables
  Read user input
  convertTemperature(input)
  write result to output
}

convertTemperature(input)
{
  Declare variables
  Determine unit of input
  if unit is C
    then convert to F and K
  else if unit is F
    then convert to C and K
  else if unit is K
    then convert to F and C
  else
    print invalid unit
}

=#

function main()
  temperature = "32F"

  println(temperature)
  c, f, k, isSuccess = convertTemperature(temperature)
  if isSuccess
    println("\n", c, " °C, ", f, " °F, ", k, " K")
  else
    println("\nConversion unsuccessful. Try again!")
  end
end

function convertTemperature(temperature::String)
  temperature = replace(temperature, " " => "")
  unit = temperature[end]

  if (unit == 'C' || unit == 'c')
    c = parse(Float64, temperature[begin:end-1])

    f = c * 9 / 5 + 32
    k = c + 273.15
  elseif  (unit == 'F' || unit == 'F')
    f = parse(Float64, temperature[begin:end-1])

    c = (f - 32) * 5 / 9
    k = c + 273.15
  elseif (unit == 'K' || unit == 'k')
    k = parse(Float64, temperature[begin:end-1])

    c = k - 273.15
    f = c * 9 / 5 + 32
  else
    println("Unit is invalid")
    return 0, 0, 0, false
  end

  return c, f, k, true
end

main()

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
module IrisBase

include("Bravais.jl")
using .Bravais
export BravaisLattice

include("CoordinateSystems.jl")
using .CoordinateSystems
export Polar,
Cartesian

include("BoundaryConditions.jl")
using .BoundaryConditions
export PML,
cPML,
noBL,
DirichletBC,
NeumannBC,
RobinBC,
FloquetBC,
MatchedBC

include("DielectricFunctions.jl")
using .DielectricFunctions

include("Shapes.jl")
using .Shapes
export Circle,
Ellipse,
Square,
Rectangle,
Parallelogram,
DeformedDisk,
Universe

include("DifferentialOperators.jl")
using .DifferentialOperators
export laplacian

include("Dispersions.jl")
using .Dispersions

include("DefineSimulation.jl")
using .DefineSimulation
export Domain,
System,
Discretization,
Boundary,
Channels,
Scattering,
Simulation

include("ConstructionTools.jl")
using .ConstructionTools
export Subdomains,
add_subdomain,
build_domain,
build_pc_domain,
line_defect,
site_defect,
add_halfspace_waveguide,
add_halfspace_waveguides,
add_planar_waveguide,
add_planar_waveguides,
add_pc_waveguide,
add_pc_waveguides

include("Bound.jl")
using .Bound

end # module


text:
# device intrinsics for querying the compute SimpleVersion and PTX ISA version

## a GPU-compatible version number

export SimpleVersion, @sv_str

struct SimpleVersion
    major::UInt32
    minor::UInt32

    SimpleVersion(major, minor=0) = new(major, minor)
end

function Base.tryparse(::Type{SimpleVersion}, v::AbstractString)
    parts = split(v, ".")
    1 <= length(parts) <= 2 || return nothing

    int_parts = map(parts) do part
        tryparse(Int, part)
    end
    any(isnothing, int_parts) && return nothing

    SimpleVersion(int_parts...)
end

function Base.parse(::Type{SimpleVersion}, v::AbstractString)
    ver = tryparse(SimpleVersion, v)
    ver === nothing && throw(ArgumentError("invalid SimpleVersion string: '$v'"))
    return ver
end

SimpleVersion(v::AbstractString) = parse(SimpleVersion, v)

@inline function Base.isless(a::SimpleVersion, b::SimpleVersion)
    (a.major < b.major) && return true
    (a.major > b.major) && return false
    (a.minor < b.minor) && return true
    (a.minor > b.minor) && return false
    return false
end

macro sv_str(str)
    SimpleVersion(str)
end

## accessors for the compute SimpleVersion and PTX ISA version

export compute_capability, ptx_isa_version

for var in ["sm_major", "sm_minor", "ptx_major", "ptx_minor"]
    @eval @inline $(Symbol(var))() =
        Base.llvmcall(
            $("""@$var = external global i32
                 define i32 @entry() #0 {
                     %val = load i32, i32* @$var
                     ret i32 %val
                 }
                 attributes #0 = { alwaysinline }
            """, "entry"), UInt32, Tuple{})
end

@device_function @inline compute_capability() = SimpleVersion(sm_major(), sm_minor())
@device_function @inline ptx_isa_version() = SimpleVersion(ptx_major(), ptx_minor())


text:
function FieldOrekitStepNormalizer(arg0::RealFieldElement, arg1::FieldOrekitFixedStepHandler)
    return FieldOrekitStepNormalizer((RealFieldElement, FieldOrekitFixedStepHandler), arg0, arg1)
end

function handle_step(obj::FieldOrekitStepNormalizer, arg0::FieldOrekitStepInterpolator, arg1::jboolean)
    return jcall(obj, "handleStep", void, (FieldOrekitStepInterpolator, jboolean), arg0, arg1)
end

function init(obj::FieldOrekitStepNormalizer, arg0::FieldSpacecraftState, arg1::FieldAbsoluteDate)
    return jcall(obj, "init", void, (FieldSpacecraftState, FieldAbsoluteDate), arg0, arg1)
end

function requires_dense_output(obj::FieldOrekitStepNormalizer)
    return jcall(obj, "requiresDenseOutput", jboolean, ())
end


text:
#! /usr/bin/julia

# Rosetta Code, Pangram checker

function makepangramchecker{T<:String}(a::T)
    abet = sort(unique(split(uppercase(a), "")))
    alen = length(abet)
    function ispangram{T<:String}(s::T)
        alen <= length(s) || return false
        ps = filter(c->(c in abet), unique(split(uppercase(s), "")))
        return length(ps) == alen
    end
end

tests = ["Pack my box with five dozen liquor jugs.",
         "The quick brown fox jumps over a lazy dog.",
         "The quick brown fox jumps\u2323over the lazy dog.",
         "The five boxing wizards jump quickly.",
         "This sentence contains A-Z but not the whole alphabet."]

isenglishpang = makepangramchecker("abcdefghijklmnopqrstuvwxyz")

for s in tests
    print("The sentence \"", s, "\" is ")
    if !isenglishpang(s)
        print("not ")
    end
    println("a pangram.")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using Base.Dates

import Base.Filesystem: StatStruct

immutable Status
    device::UInt64
    inode::UInt64
    mode::Mode
    nlink::Int64
    user::User
    group::Group
    rdev::UInt64
    size::Int64
    blksize::Int64
    blocks::Int64
    mtime::DateTime
    ctime::DateTime

    function Status(s::StatStruct)
        new(
            s.device,
            s.inode,
            Mode(s.mode),
            s.nlink,
            User(s.uid),
            Group(s.gid),
            s.rdev,
            s.size,
            s.blksize,
            s.blocks,
            unix2datetime(s.mtime),
            unix2datetime(s.ctime)
        )
    end
end

function Base.show(io::IO, s::Status)
    output = "Status(\n" *
        "  device = $(s.device),\n" *
        "  inode = $(s.inode),\n" *
        "  mode = $(s.mode),\n" *
        "  nlink = $(s.nlink),\n" *
        "  uid = $(s.user),\n" *
        "  gid = $(s.group),\n" *
        "  rdev = $(s.rdev),\n" *
        "  size = $(s.size) ($(_datasize(s.size))),\n" *
        "  blksize = $(s.blksize) ($(_datasize(s.blksize))),\n" *
        "  blocks = $(s.blocks),\n" *
        "  mtime = $(s.mtime),\n" *
        "  ctime = $(s.ctime),\n)"

    print(io, output)
end


text:
"""
    enclose(f, dom[, solver=NaturalEnclosure()]; kwargs...)

Return a range enclosure of a univariate or multivariate function on the given
domain.

### Input

- `f`      -- function or `AbstractPolynomialLike` object
- `dom`    -- hyperrectangular domain, either a unidimensional  `Interval` or
              a multidimensional `IntervalBox`
- `solver` -- (optional, default: `NaturalEnclosure()`) choose one among the
              available solvers; you can get a list of available solvers with
              `subtypes(AbstractEnclosureAlgorithm)`
- `kwargs` -- optional keyword arguments passed to the solver; for available
              options see the documentation of each solver

### Output

An interval enclosure of the range of `f` over `dom`.

### Examples

```jldoctest enclose_examples
julia> enclose(x -> 1 - x^4 + x^5, 0..1) # use default solver
[0, 2]

julia> enclose(x -> 1 - x^4 + x^5, 0..1, TaylorModelsEnclosure())
[0.8125, 1.09375]
```

A vector of solvers can be passed in the `solver` options. Then, the result is
obtained by intersecting the range enclosure of each solver.

```jldoctest enclose_examples
julia> enclose(x -> 1 - x^4 + x^5, 0..1, [TaylorModelsEnclosure(), NaturalEnclosure()])
[0.8125, 1.09375]

```
"""
function enclose(f::Function, dom::Interval_or_IntervalBox,
                 solver::AbstractEnclosureAlgorithm=NaturalEnclosure(); kwargs...)
    return _enclose(solver, f, dom; kwargs...)
end

function enclose(p::AbstractPolynomialLike, dom::Interval_or_IntervalBox,
                 solver::AbstractEnclosureAlgorithm=NaturalEnclosure(); kwargs...)
    f(x...) = p(variables(p) => x)
    return _enclose(solver, f, dom; kwargs...)
end

function enclose(f::Function, dom::Interval_or_IntervalBox,
                 method::Vector; kwargs...)
   return mapreduce(ξ -> _enclose(ξ, f, dom; kwargs...), ∩, method)
end

"""
    relative_precision(x::Interval, xref::Interval)

Return the relative precision of an interval with respect to a reference interval.

### Input

- `x`    -- test interval
- `xref` -- reference interval

### Output

Left and right relative precision (in %) computed as
- `rleft = (inf(xref) - inf(x)) / diam(xref) * 100%`
- `rright = (sup(x) - sup(xright)) / diam(xref) * 100%`

### Examples

```jldoctest relative_precision
julia> xref = interval(-1.2, 4.6)
[-1.2, 4.6]

julia> x = interval(-1.25, 7.45)
[-1.25, 7.45001]

julia> relative_precision(x, xref)
(0.8620689655172422, 49.13793103448277)
```

### Algorithm

This function measures the relative precision of the result in a more informative
way than taking the scalar overestimation because it evaluates the precision of
the lower and the upper bounds separately (cf. Eq. (20) in [1]).

[1] Althoff, Matthias, Dmitry Grebenyuk, and Niklas Kochdumper.
   *Implementation of Taylor models in CORA 2018.*
   Proc. of the 5th International Workshop on Applied Verification for
   Continuous and Hybrid Systems. 2018.

"""
function relative_precision(x::Interval, xref::Interval)
    x⁻, x⁺ = inf(x), sup(x)
    xref⁻, xref⁺ = inf(xref), sup(xref)

    Δref = xref⁺ - xref⁻
    @assert Δref != zero(Δref) "the width of the reference interval is $Δref, but " *
    "it should have a non-zero width"

    rel⁻ = -(x⁻ - xref⁻) / Δref * 100
    rel⁺ = (x⁺ - xref⁺) / Δref * 100
    return rel⁻, rel⁺
end


text:
module NnlsTest

using Test
using LinearAlgebra
using Random
using PyCall
const pyopt = pyimport_conda("scipy.optimize", "scipy")
using NonNegLeastSquares.NNLS

# Allocation measurement doesn't work reliably on Julia v0.5 when
# code coverage checking is enabled.
const test_allocs = VERSION >= v"0.6-" || Base.JLOptions().code_coverage == 0

"""
Measure memory allocation within a function to avoid issues
with global variables.
"""
macro wrappedallocs(expr)
    argnames = [gensym() for a in expr.args]
    quote
        function g($(argnames...))
            @allocated $(Expr(expr.head, argnames...))
        end
        $(Expr(:call, :g, [esc(a) for a in expr.args]...))
    end
end

@testset "bigfloat" begin
    Random.seed!(5)
    for i in 1:100
        m = rand(1:10)
        n = rand(1:10)
        A = randn(m, n)
        b = randn(m)
        x1 = nnls(A, b)
        x2 = nnls(BigFloat.(A), BigFloat.(b))
        @test x1 ≈ x2
    end
end

@testset "apply_householder!" begin
    Random.seed!(2)
    for i in 1:10
        u = randn(rand(3:10))
        c = randn(length(u))

        u1 = copy(u)
        c1 = copy(c)
        up1 = NNLS.construct_householder!(u1, 0.0)
        NNLS.apply_householder!(u1, up1, c1)

        if test_allocs
            u2 = copy(u)
            c2 = copy(c)
            @test @wrappedallocs(NNLS.construct_householder!(u2, 0.0)) == 0
            up2 = up1
            @test @wrappedallocs(NNLS.apply_householder!(u2, up2, c2)) == 0
        end
    end
end

@testset "orthogonal_rotmat" begin
    Random.seed!(3)
    for i in 1:1000
        a = randn()
        b = randn()
        c, s, sig = NNLS.orthogonal_rotmat(a, b)
        @test [c s; -s c] * [a, b] ≈ [sig, 0]
        if test_allocs
            @test @wrappedallocs(NNLS.orthogonal_rotmat(a, b)) == 0
        end
    end
end

if test_allocs
    @testset "nnls allocations" begin
        Random.seed!(101)
        for i in 1:50
            m = rand(20:100)
            n = rand(20:100)
            A = randn(m, n)
            b = randn(m)
            work = NNLSWorkspace(A, b)
            @test @wrappedallocs(nnls!(work)) == 0
        end
    end
end

@testset "nnls workspace reuse" begin
    Random.seed!(200)
    m = 10
    n = 20
    work = NNLSWorkspace(m, n)
    nnls!(work, randn(m, n), randn(m))
    for i in 1:100
        A = randn(m, n)
        b = randn(m)
        if test_allocs
            @test @wrappedallocs(nnls!(work, A, b)) == 0
        else
            nnls!(work, A, b)
        end
        @test work.x == pyopt[:nnls](A, b)[1]
    end

    m = 20
    n = 10
    for i in 1:100
        A = randn(m, n)
        b = randn(m)
        nnls!(work, A, b)
        @test work.x == pyopt[:nnls](A, b)[1]
    end
end

if test_allocs
    @testset "non-Int Integer workspace" begin
        m = 10
        n = 20
        A = randn(m, n)
        b = randn(m)
        work = NNLSWorkspace(A, b, Int32)
        # Compile
        nnls!(work)

        A = randn(m, n)
        b = randn(m)
        work = NNLSWorkspace(A, b, Int32)
        @test @wrappedallocs(nnls!(work)) == 0
    end
end

@testset "nnls vs scipy" begin
    Random.seed!(5)
    for i in 1:5000
        m = rand(1:60)
        n = rand(1:60)
        A = randn(m, n)
        b = randn(m)
        x1 = nnls(A, b)
        x2, residual2 = pyopt[:nnls](A, b)
        @test x1 == x2
    end
end

end


text:
if VERSION < v"0.6.0-rc1" || VERSION >= v"0.7-"
    error("VS Code julia language server only works with julia 0.6.")
end

try
    if length(Base.ARGS) != 3
        error("Invalid number of arguments passed to julia language server.")
    end

    conn = STDOUT
    (outRead, outWrite) = redirect_stdout()

    if Base.ARGS[2] == "--debug=no"
        const global ls_debug_mode = false
    elseif Base.ARGS[2] == "--debug=yes"
        const global ls_debug_mode = true
    end

    push!(LOAD_PATH, joinpath(dirname(@__FILE__), "packages"))
    push!(LOAD_PATH, Base.ARGS[1])

    using Compat
    using JSON
    using URIParser
    using LanguageServer

    server = LanguageServerInstance(STDIN, conn, ls_debug_mode, Base.ARGS[1])
    run(server)
catch e
    st = catch_stacktrace()
    vscode_pipe_name = is_windows() ? "\\\\.\\pipe\\vscode-language-julia-lscrashreports-$(Base.ARGS[3])" : joinpath(tempdir(), "vscode-language-julia-lscrashreports-$(Base.ARGS[3])")
    pipe_to_vscode = connect(vscode_pipe_name)
    try
        # Send error type as one line
        println(pipe_to_vscode, typeof(e))

        # Send error message as one line
        showerror(pipe_to_vscode, e)
        println(pipe_to_vscode)

        # Send stack trace, one frame per line
        # Note that stack frames need to be formatted in Node.js style
        for s in st
            print(pipe_to_vscode, " at ")
            Base.StackTraces.show_spec_linfo(pipe_to_vscode, s)

            filename = string(s.file)

            # Now we need to sanitize the filename so that we don't transmit
            # things like a username in the path name
            filename = normpath(filename)
            if isabspath(filename)
                root_path_of_extension = normpath(joinpath(@__DIR__, "..", ".."))
                if startswith(filename, root_path_of_extension)
                    filename = joinpath(".", filename[endof(root_path_of_extension)+1:end])
                else
                    filename = basename(filename)
                end
            else
                filename = basename(filename)
            end

            # Use a line number of "0" as a proxy for unknown line number
            print(pipe_to_vscode, " (", filename, ":", s.line >= 0 ? s.line : "0", ":1)" )

            # TODO Unclear how we can fit this into the Node.js format
            # if s.inlined
            #     print(pipe_to_vscode, " [inlined]")
            # end

            println(pipe_to_vscode)
        end
    finally
        close(pipe_to_vscode)
    end
    rethrow()
end

max_stars_repo_name: kindlychung/julia-vscode

max_issues_repo_name: kindlychung/julia-vscode

max_forks_repo_name: kindlychung/julia-vscode


text:
# functional derivative computation

import PDESolver._evalFunctionalDeriv_q

"""
  evalFunctionalDeriv_q for Advection
"""
function _evalFunctionalDeriv_q(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::AdvectionData{Tsol}, opts,
                           functionalData::AbstractBoundaryFunctional,
                           func_deriv_arr::Abstract3DArray) where {Tmsh, Tsol}

  if mesh.isDG
    boundaryinterpolate!(mesh.sbpface, mesh.bndryfaces, eqn.q, eqn.q_bndry)
  end

  calcFunctionalDeriv(mesh, sbp, eqn, opts, functionalData, func_deriv_arr)

  return nothing
end

@doc """
### AdvectionEquationMod.calcFunctionalDeriv

Computes a 3D array of the derivative of a functional w.r.t eqn.q on all
mesh nodes.

**Inputs**

*  `mesh`  : Abstract mesh object
*  `sbp`   : Summation-By-Parts operator
*  `eqn`   : Advection equation object
*  `opts`  : Options dictionary
*  `functionalData` : Object of subtype of AbstractBoundaryFunctional. This is
                      the type associated with the adjoint of the functional
                      being computed and holds all the necessary data.
*  `func_deriv_arr` : 3D array that stors the derivative of functional w.r.t
                      eqn.q. It has a structure [1, numnodes_per_element, numEl]

**Outputs**

*  None

"""->
#=
function calcFunctionalDeriv{Tmsh, Tsol}(mesh::AbstractCGMesh{Tmsh}, sbp::AbstractOperator,
                             eqn ::AdvectionData{Tsol}, opts, functor, functional_edges,
                             functionalData, func_deriv_arr)

  alpha_x = eqn.params.alpha_x
  alpha_y = eqn.params.alpha_y

  # Obtain the derivative of the integrand at all meh.bndry
  integrand = zeros(eqn.bndryflux)
  # println("size of integrand = ", size(integrand))

  for itr = 1:length(functional_edges)
    g_edge_number = functional_edges[itr] # Extract geometric edge number
    # get the boundary array associated with the geometric edge
    itr2 = 0
    for itr2 = 1:mesh.numBC
      if findfirst(mesh.bndry_geo_nums[itr2],g_edge_number) > 0
        break
      end
    end

    start_index = mesh.bndry_offsets[g_edge_number]
    end_index = mesh.bndry_offsets[g_edge_number + 1]
    idx_range = start_index:(end_index-1)  # Index range
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i

    nfaces = length(bndry_facenums)

    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      global_facenum = idx_range[i]
      for j = 1:mesh.numNodesPerFace
        k = mesh.facenodes[j, bndry_i.face]
        q = eqn.q[1,k,bndry_i.element]
        x = ro_sview(mesh.coords, :, k, bndry_i.element)
        dxidx = ro_sview(mesh.dxidx, :, :, k, bndry_i.element)
        nrm = ro_sview(mesh.sbpface.normal, :, bndry_i.face)
        nx = dxidx[1,1]*nrm[1] + dxidx[2,1]*nrm[2]
        ny = dxidx[1,2]*nrm[1] + dxidx[2,2]*nrm[2]
        integrand[1,j,global_facenum] = calcIntegrandDeriv(opts, functor, eqn.params,
                                        nx, ny, q, functionalData)
      end  # End for j = 1:mesh.numNodesPerFace
    end    # End for i = 1:nfaces
  end      # End for itr = 1:length(functional_edges)

  # println("mesh.bndryfaces = \n", mesh.bndryfaces)

  boundaryintegrate!(mesh.sbpface, mesh.bndryfaces, integrand, func_deriv_arr)

  return nothing
end
=#
# DG Version
function calcFunctionalDeriv(mesh::AbstractDGMesh{Tmsh}, sbp::AbstractOperator,
                 eqn::AdvectionData{Tsol}, opts,
                 functionalData::AbstractBoundaryFunctional,
                 func_deriv_arr) where {Tmsh, Tsol}

  alpha_x = eqn.params.alpha_x
  alpha_y = eqn.params.alpha_y

  # Obtain the derivative of the integrand at all mesh.bndry
  integrand = zeros(eqn.q_bndry)
  # Populate integrand
  for itr = 1:length(functionalData.bcnums)
    bcnum = functionalData.bcnums[itr]

    start_index = mesh.bndry_offsets[bcnum]
    end_index = mesh.bndry_offsets[bcnum+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
    nfaces = length(bndry_facenums)

    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      global_facenum = idx_range[i]
      for j = 1:mesh.sbpface.numnodes
        q = eqn.q_bndry[ 1, j, global_facenum]
        coords = ro_sview(mesh.coords_bndry, :, j, global_facenum)
        nrm = ro_sview(mesh.nrm_bndry, :, j, global_facenum)
        nx = nrm[1]
        ny = nrm[2]
        integrand[1,j,global_facenum] = calcIntegrandDeriv(opts, eqn.params,
                                        nx, ny, q, functionalData)
      end  # End for j = 1:mesh.sbpface.numnodes
    end    # End for i = 1:nfaces
  end

  boundaryintegrate!(mesh.sbpface, mesh.bndryfaces, integrand, func_deriv_arr)

  return nothing
end    # End function calcFunctionalDeriv

@doc """
### AdvectionEquationMod.calcIntegrandDeriv

Compute the derivative of the integrand at a point. It presently uses complex
step to compute the derivative

**Inputs**

*  `opts`    : Input dictionary
*  `params`  : the ParamType for the equation
*  `nx` & `ny` : Normal vectors
*  `q`       : Solution variable
*  `functionalData` : Functional object

**Outputs**

*  `integrand_deriv` : derivative of the functor w.r.t q

"""->

function calcIntegrandDeriv(opts, params::ParamType2, nx, ny, q,
                            functionalData::AbstractBoundaryFunctional)

  pert = complex(0, 1e-20)  # complex perturbation
  q += pert
  val = calcBoundaryFunctionalIntegrand(params, nx, ny, q, functionalData)
  integrand_deriv = imag(val)/norm(pert)

  return integrand_deriv
end


text:
using Test

@testset "MUMPSjInv" begin
    @testset "DivGrad" begin include("testDivGrad.jl") end
    @testset "Two systems" begin include("testTwoSystem.jl") end
    @testset "Two systems parallel" begin include("testTwoSystemParallel.jl") end
    @testset "Sparse RHS" begin include("testSparseRHS.jl") end
    @testset "testDestroyMUMPS" begin include("testDestroyMUMPS.jl") end
end


text:
# ===============================================================
# Discretize using the correction hull of the matrix exponential
# ===============================================================

"""
    CorrectionHull{EM} <: AbstractApproximationModel

Discretization using the correction hull of the matrix exponential.

### Fields

- `exp`   -- exponentiation method
- `order` -- order of the Taylor series expansion of the matrix exponential

### Algorithm

For the homogeneous case, this method implements the transformation:

```math
Ω_0 = CH(X_0, e^{Aδ}  X_0) ⊕ FX_0
```
where ``F`` is the correction (interval) matrix.

For the inhomogeneous case, ``x' = Ax + u``,  ``x ∈ X, u ∈ U``, implements
``Ω_0 = CH(X_0, exp(Aδ)  X0) ⊕ FX0`` where ``F`` is the correction (interval) matrix.

In both cases, if ``A`` is an interval matrix, the exponential is overapproximated
using methods from `IntervalMatrices.jl`.
"""
struct CorrectionHull{EM} <: AbstractApproximationModel
   order::Int
   exp::EM
end

# convenience constructor using symbols
function CorrectionHull(; order::Int=10, exp=IntervalExpAlg(order))
    return CorrectionHull(order, _alias(exp))
end

function Base.show(io::IO, alg::CorrectionHull)
    print(io, "`CorrectionHull` approximation model with: \n")
    print(io, "    - exponentiation method: $(alg.exp) \n")
    print(io, "    - order: $(alg.order)\n")
end

Base.show(io::IO, m::MIME"text/plain", alg::CorrectionHull) = print(io, alg)

# -----------------------------------------------------------------
# Correction hull: homogeneous case x' = Ax, x in X
# -----------------------------------------------------------------

function discretize(ivp::IVP{<:CLCS, <:LazySet}, δ, alg::CorrectionHull)
    A = state_matrix(ivp)
    X0 = initial_state(ivp)
    X = stateset(ivp)

    # compute exp(A*δ) * X0
    Φ = _exp(A, δ, alg.exp)

    # compute Ω0 = CH(X0, ΦX0) ⊕ FX0
    Ω0 = _discretize_chull(A, Φ, X0, δ, alg)

    Sdis = ConstrainedLinearDiscreteSystem(Φ, X)
    return InitialValueProblem(Sdis, Ω0)
end

function _discretize_chull(A, Φ::IntervalMatrix, X0, δ, alg)
    X0z = _convert_or_overapproximate(X0, Zonotope)
    Y = _overapproximate(Φ * X0z, Zonotope)

    H = overapproximate(CH(X0z, Y), Zonotope)
    F = correction_hull(A, δ, alg.order)
    R = _overapproximate(F * X0z, Zonotope)
    Ω0 = minkowski_sum(H, R)
    return Ω0
end

# F(δ) without the E(δ) correction term
function _correction_hull_without_E(A, δ, p)
    timeint(δ, i) = interval((i^(-i / (i-1)) - i^(-1 / (i-1))) * δ^i, 0)
    F = sum(map(x -> timeint(δ, i) * x, A^i / factorial(i)) for i in 2:p)
    return IntervalMatrix(F)
end

function _discretize_chull(A, Φ::AbstractMatrix, X0, δ, alg)
    X0z = _convert_or_overapproximate(X0, Zonotope)
    Y = linear_map(Φ, X0z)

    H = overapproximate(CH(X0z, Y), Zonotope)
    F = _correction_hull_without_E(A, δ, alg.order)
    R = _overapproximate(F * X0z, Zonotope)

    Ω0 = minkowski_sum(H, R)
    return Ω0
end

# -----------------------------------------------------------------
# Correction hull: inhomogeneous case x' = Ax + u, x in X, u ∈ U
# -----------------------------------------------------------------
function discretize(ivp::IVP{<:CLCCS, <:LazySet}, δ, alg::CorrectionHull)
    A = state_matrix(ivp)
    X0 = initial_state(ivp)
    X = stateset(ivp)
    U = next_set(inputset(ivp), 1) # inputset(ivp)
    n = size(A, 1)

    # here U is an interval matrix map of a lazyset, TODO refactor / dispatch
    if isa(U, LinearMap)
        Uz = _convert_or_overapproximate(Zonotope, LazySets.set(U))
        B = matrix(U)
        Uz = isinterval(B) ? _overapproximate(B * Uz, Zonotope) : linear_map(B, Uz)
    else # LazySet
        Uz = _convert_or_overapproximate(Zonotope, U)
    end
    if zeros(dim(U)) ∉ Uz
        error("this function is not implemented, see issue #253")
    end

    # TODO refactor Ω0_homog
    # TODO refactor / dispatch
    X0z = _convert_or_overapproximate(Zonotope, X0)

    Φ = _exp(A, δ, alg.exp)
    if isinterval(Φ)
        Y = _overapproximate(Φ * X0z, Zonotope)
    else
        Y = linear_map(Φ, X0z)
    end

    H = overapproximate(CH(X0z, Y), Zonotope)
    F = correction_hull(A, δ, alg.order)
    R = _overapproximate(F * X0z, Zonotope)
    Ω0_homog = minkowski_sum(H, R)

    # compute C(δ) * U
    Cδ = _Cδ(A, δ, alg.order)
    Ud = _overapproximate(Cδ * Uz, Zonotope)
    Ω0 = minkowski_sum(Ω0_homog, Ud)
    # Ω0 = _apply_setops(Ω0, alg.setops) # TODO requires to add `setops` field to the struct
    Idn = Φ # IntervalMatrix(one(A)) or IdentityMultiple(one(eltype(A)), n) # FIXME
    Sdis = ConstrainedLinearControlDiscreteSystem(Φ, Idn, X, Ud)
    return InitialValueProblem(Sdis, Ω0)
end


text:
using Random

module ExampleEnvConst
const state_initial_range = 1:10
end

mutable struct ExampleEnvironment <: AbstractEnvironment
    some_variable::Float64
    state::Int64
    normalized_state::Bool
    ExampleEnvironment(rng=Random.GLOBAL_RNG) = new(rand(rng), rand(rng, ExampleEnvConst.state_initial_range), false)
end

function reset!(env::ExampleEnvironment; rng = Random.GLOBAL_RNG, kwargs...)
    env.some_variable = rand(rng)
    env.state = rand(rng, ExampleEnvConst.state_initial_range)
end

function get_actions(env::AbstractEnvironment)
    return Set()
end

function environment_step!(env::ExampleEnvironment, action::Int64; rng = Random.GLOBAL_RNG, kwargs...)
    env.state = rand(rng, ExampleEnvConst.state_initial_range)
    env.some_variable = env.state + randn(rng)
end

function get_reward(env::ExampleEnvironment) # -> determines if the agent_state is terminal
    return -1
end

function get_state(env::ExampleEnvironment) # -> get state of agent
    if env.normalized_state
        return get_normalized_state(env)
    else
        return [env.state]
    end
end

function is_terminal(env::ExampleEnvironment) # -> determines if the agent_state is terminal
    if env.state == 1
        return true
    end
    return false
end

function get_normalized_state(env::ExampleEnvironment)
    return [env.state/10.0]
end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
# ---
# title: 354. Russian Doll Envelopes
# id: problem354
# author: Indigo
# date: 2021-03-04
# difficulty: Hard
# categories: Binary Search, Dynamic Programming
# link: <https://leetcode.com/problems/russian-doll-envelopes/description/>
# hidden: true
# ---
#
# You have a number of envelopes with widths and heights given as a pair of
# integers `(w, h)`. One envelope can fit into another if and only if both the
# width and height of one envelope is greater than the width and height of the
# other envelope.
#
# What is the maximum number of envelopes can you Russian doll? (put one inside
# other)
#
# **Note:**
# Rotation is not allowed.
#
# **Example:**
#
#
#
#     Input: [[5,4],[6,4],[6,7],[2,3]]
#     Output: 3
#     **Explanation: T** he maximum number of envelopes you can Russian doll is 3 ([2,3] => [5,4] => [6,7]).
#
#
#
## @lc code=start
using LeetCode

function max_envelopes(envelopes::Vector{Vector{Int}})
    isempty(envelopes) && return 0
    sort!(envelopes, by = x -> (x[1], -x[2]))
    dp = [envelopes[1][2]]
    for envo in @view(envelopes[2:end])
        if envo[2] > dp[end]
            push!(dp, envo[2])
        else
            dp[searchsortedfirst(dp, envo[2])] = envo[2]
        end
    end
    length(dp)
end
## @lc code=end


text:
export edistance, get_direction, walk!

#######################################################################################
# %% Distances and directions in Grid/Continuous space
#######################################################################################
"""
    edistance(a, b, model::ABM)

Return the euclidean distance between `a` and `b` (either agents or agent positions),
respecting periodic boundary conditions (if in use). Works with any space where it makes
sense: currently `GridSpace` and `ContinuousSpace`.

Example usage in the [Flock model](@ref).
"""
edistance(
    a::A,
    b::B,
    model::ABM{<:Union{ContinuousSpace,GridSpace}},
) where {A<:AbstractAgent,B<:AbstractAgent} = edistance(a.pos, b.pos, model)

function edistance(
    a::ValidPos,
    b::ValidPos,
    model::ABM{<:Union{ContinuousSpace{D,false},GridSpace{D,false}}},
) where {D}
    sqrt(sum(abs2.(a .- b)))
end

function edistance(
    p1::ValidPos,
    p2::ValidPos,
    model::ABM{<:ContinuousSpace{D,true}},
) where {D}
    total = 0.0
    for (a, b, d) in zip(p1, p2, model.space.extent)
        delta = abs(b - a)
        if delta > d - delta
            delta = d - delta
        end
        total += delta^2
    end
    sqrt(total)
end

function edistance(p1::ValidPos, p2::ValidPos, model::ABM{<:GridSpace{D,true}}) where {D}
    total = 0.0
    for (a, b, d) in zip(p1, p2, size(model.space))
        delta = abs(b - a)
        if delta > d - delta
            delta = d - delta
        end
        total += delta^2
    end
    sqrt(total)
end

"""
    get_direction(from, to, model::ABM)
Return the direction vector from the position `from` to position `to` taking into account
periodicity of the space.
"""
get_direction(from, to, model::ABM) = get_direction(from, to, model.space)

function get_direction(
    from::NTuple{D,Float64},
    to::NTuple{D,Float64},
    space::Union{ContinuousSpace{D,true}, GridSpace{D,true}}
) where {D}
    best = to .- from
    for offset in Iterators.product([-1:1 for _ in 1:D]...)
        dir = to .+ offset .* size(space) .- from
        sum(dir .^ 2) < sum(best .^ 2) && (best = dir)
    end
    return best
end

function get_direction(from, to, ::Union{GridSpace,ContinuousSpace})
    return to .- from
end

#######################################################################################
# %% Utilities for graph-based spaces (Graph/OpenStreetMap)
#######################################################################################
GraphBasedSpace = Union{GraphSpace, OpenStreetMapSpace}
_get_graph(space::GraphSpace) = space.graph
_get_graph(space::OpenStreetMapSpace) = space.m.g
"""
    nv(model::ABM)
Return the number of positions (vertices) in the `model` space.
"""
LightGraphs.nv(abm::ABM{<:GraphBasedSpace}) = LightGraphs.nv(_get_graph(abm.space))

"""
    ne(model::ABM)
Return the number of edges in the `model` space.
"""
LightGraphs.ne(abm::ABM{<:GraphBasedSpace}) = LightGraphs.ne(_get_graph(abm.space))

positions(model::ABM{<:GraphBasedSpace}) = 1:nv(model)

function nearby_positions(
    position::Integer,
    model::ABM{<:GraphBasedSpace},
    radius::Integer;
    kwargs...,
)
    output = copy(nearby_positions(position, model; kwargs...))
    for _ in 2:radius
        newnps = (nearby_positions(np, model; kwargs...) for np in output)
        append!(output, reduce(vcat, newnps))
        unique!(output)
    end
    filter!(i -> i ≠ position, output)
end

#######################################################################################
# %% Walking
#######################################################################################
"""
    walk!(agent, direction::NTuple, model; ifempty = false)

Move agent in the given `direction` respecting periodic boundary conditions.
If `periodic = false`, agents will walk to, but not exceed the boundary value.
Possible on both `GridSpace` and `ContinuousSpace`s.

The dimensionality of `direction` must be the same as the space. `GridSpace` asks for
`Int`, and `ContinuousSpace` for `Float64` vectors, describing the walk distance in
each direction. `direction = (2, -3)` is an example of a valid direction on a
`GridSpace`, which moves the agent to the right 2 positions and down 3 positions.
Velocity is ignored for this operation in `ContinuousSpace`.

## Keywords
- `ifempty` will check that the target position is unnocupied and only move if that's true. Available only on `GridSpace`.

Example usage in [Battle Royale](@ref).
"""
function walk!(
    agent::AbstractAgent,
    direction::NTuple{D,Int},
    model::ABM{<:GridSpace{D,true}};
    kwargs...,
) where {D}
    target = mod1.(agent.pos .+ direction, size(model.space))
    walk_if_empty!(agent, target, model; kwargs...)
end

function walk!(
    agent::AbstractAgent,
    direction::NTuple{D,Int},
    model::ABM{<:GridSpace{D,false}};
    kwargs...,
) where {D}
    target = min.(max.(agent.pos .+ direction, 1), size(model.space))
    walk_if_empty!(agent, target, model; kwargs...)
end

function walk!(
    agent::AbstractAgent,
    direction::NTuple{D,Float64},
    model::ABM{<:ContinuousSpace{D,true}};
    kwargs...,
) where {D}
    target = mod1.(agent.pos .+ direction, model.space.extent)
    move_agent!(agent, target, model)
end

function walk!(
    agent::AbstractAgent,
    direction::NTuple{D,Float64},
    model::ABM{<:ContinuousSpace{D,false}};
    kwargs...,
) where {D}
    target = min.(max.(agent.pos .+ direction, 0.0), model.space.extent .- 1e-15)
    move_agent!(agent, target, model)
end

function walk_if_empty!(agent, target, model; ifempty::Bool = false)
    if ifempty
        isempty(target, model) && move_agent!(agent, target, model)
    else
        move_agent!(agent, target, model)
    end
end

"""
    walk!(agent, rand, model)

Invoke a random walk by providing the `rand` function in place of
`distance`. For `GridSpace`, the walk will cover ±1 positions in all directions,
`ContinuousSpace` will reside within [-1, 1].
"""
walk!(agent, ::typeof(rand), model::ABM{<:GridSpace{D}}; kwargs...) where {D} =
    walk!(agent, Tuple(rand(model.rng, -1:1, D)), model; kwargs...)

walk!(agent, ::typeof(rand), model::ABM{<:ContinuousSpace{D}}; kwargs...) where {D} =
    walk!(agent, Tuple(2.0 * rand(model.rng) - 1.0 for _ in 1:D), model; kwargs...)


text:
using FastGMM
using Test

@testset "FastGMM.jl" begin
    # Write your tests here.
end


text:
using Cropbox

@system S(Controller) begin
    a => 1 ~ preserve
    b(a) ~ track
    c(b) ~ accumulate
end
c = @config
r = simulate(S, config=c)
simulate(S, configs=[c])
visualize(r, :time, :c; backend=:UnicodePlots) |> display
visualize(r, :time, :c; backend=:Gadfly)[] |> Cropbox.Gadfly.SVGJS() |> display
visualize(r, :time, :c; kind=:line, backend=:Gadfly)' |> Cropbox.Gadfly.SVG() |> display
visualize(S, :time, :c; backend=:UnicodePlots) |> display
visualize(S, :time, :c; backend=:Gadfly)[] |> Cropbox.Gadfly.SVGJS() |> display
r |> display
display(MIME("text/html"), r)

using Test

include("../test/state.jl")
include("../test/system.jl")


text:
function differentfrom_singlevariable_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Interval{Int})
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(
        model,
        MOI.VariableIndex,
        CP.DifferentFrom{Int},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.Interval(1, 2))

    c2 = MOI.add_constraint(model, x1, CP.DifferentFrom(2))

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
    end
end

function differentfrom_scalaraffinefunction_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Interval{Int})
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Int},
        CP.DifferentFrom{Int},
    )
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.EqualTo{Int},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.Interval(1, 2))
    c2 = MOI.add_constraint(model, x2, MOI.Interval(1, 2))

    c3 = MOI.add_constraint(model, 1 * x1 - 1 * x2, CP.DifferentFrom(0))
    c4 = MOI.add_constraint(model, 1 * x1, MOI.EqualTo(1))

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)
    @test MOI.is_valid(model, c3)
    @test MOI.is_valid(model, c4)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 2
    end
end

const differentfromtests = Dict(
    "differentfrom_singlevariable" => differentfrom_singlevariable_test,
    "differentfrom_scalaraffinefunction" =>
        differentfrom_scalaraffinefunction_test,
)

MOIT.@moitestset differentfrom


text:
using GenTeX
using Test

DisplayImage = GenTeX.DisplayImage

@testset "cache" begin
	tmpdir = tempname() * '/'
	mkdir(tmpdir)

	scale = 1.0
	cache = GenTeX.load_cache(scale, tmpdir)
    text = raw"$x$"
    hash_text = string(hash(text))
	eq = GenTeX.Equation(text, scale, "inline", "", hash_text)
	eq_image = DisplayImage(eq, tmpdir, "hash.svg")
	@test GenTeX.check_cache(cache, eq) == nothing
	GenTeX.update_cache!(cache, eq_image)
	@test GenTeX.check_cache(cache, eq) == eq_image
	eq2 = GenTeX.Equation(text, scale, "display", "", hash_text)
	@test GenTeX.check_cache(cache, eq2) == nothing
	eq3 = GenTeX.Equation(text, scale, "inline", "", hash_text)
	@test GenTeX.check_cache(cache, eq3) == eq_image

	GenTeX.write_cache!(cache, tmpdir)
	cache = GenTeX.load_cache(scale, tmpdir)
	@test GenTeX.check_cache(cache, eq) == eq_image

	GenTeX.clear_cache!(tmpdir)
	cache = GenTeX.load_cache(scale, tmpdir)
	@test GenTeX.check_cache(cache, eq) == nothing
	rm(tmpdir, recursive=true)

    mkdir(tmpdir)
    substitute_latex(raw"$x$", scale, tmpdir)
	cache = GenTeX.load_cache(scale, tmpdir)
    @test length(cache.images) == 1
    rm(tmpdir, recursive=true)
end


text:
push!(LOAD_PATH, joinpath(dirname(pwd()), "src", "models"));
push!(LOAD_PATH, joinpath(dirname(pwd()), "src", "data"));
push!(LOAD_PATH, joinpath(dirname(pwd()), "src"));

using rails: RAILS
using ego_model: EgoModel
using ego_dataset: compile_data

using ArgParse

function parse_cil()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--data-dir"
            help = "Ego data dir"
        "--traj-len"
            help = "Trajectory length for ego model"
            default = 10
        "--config"
            help = "main configuration file"
        "--batchsize"
            help = "Batchsize"
            default = 128
        "--shuffle"
            help = "Shuffle the ego dataset"
            action = :store_true
            default = false
        "--device"
            help = "cpu or cuda"
            default = "cuda"
        "--lr"
            help = "learning rate"
            default = 1e-2
        "--num-epoch"
            help = "Number of epochs to train"
            default = 100
        "--num-per-log"
            help = "Frequency to log (step)"
            default=10
    end
    return parse_args(s)
end

function main()
    args = parse_cil()
    rails = RAILS(args)
    data = compile_data(args)

    counter = 0
    for epoch = 1:args.num_epoch
        for (locs, rots, spds, acts) = data
            opt_info = rails.train_ego(locs, rots, spds, acts)
            if counter % args.num_per_log == 0
                #logger.log_ego(opt_info) #TODO: implement the logger
            end
            counter += 1
        end
    end

end

main()


text:
"""
Cross-entropy method, stochastic optimization solver.
"""
@with_kw mutable struct CEMSolver
    n_iterations::Int64 = 100
    episode_length::Int64 = 30
    num_samples::Int64 = 100
    min_elite_samples::Int64 = Int64(floor(0.1*num_samples))
    max_elite_samples::Int64 = typemax(Int64)
    elite_thresh::Float64 = -0.99
    weight_fn::Function = (d,x) -> 1.0
    add_entropy::Function = x->x
    show_progress::Bool = true
    verbose::Bool = false
end

"""
Planner for the cross-entropy method solver. Takes the `solver::CEMSolver`, the `mdp` problem,
 and fills in the optimized importance sampling distributions `is_dist`.
"""
mutable struct CEMPlanner{P<:Union{MDP,POMDP}}
    solver::CEMSolver
    mdp::P
    is_dist::Union{Dict{Symbol,Vector{Sampleable}},Nothing} # optimized importance sample distribution
end

# No work performed.
POMDPs.solve(solver::CEMSolver, mdp::Union{POMDP,MDP}) = CEMPlanner(solver, mdp, nothing)

function cem_losses(d, sample; mdp::ASTMDP, initstate::ASTState)
    sim = mdp.sim
    env = GrayBox.environment(sim)
    s = initstate
    R = 0 # accumulated reward

    BlackBox.initialize!(sim)
    AST.go_to_state(mdp, s)

    sample_length = length(last(first(sample))) # get length of sample vector ("second" element in pair using "first" key)
    for i in 1:sample_length
        env_sample = GrayBox.EnvironmentSample()
        for k in keys(sample)
            value = sample[k][i]
            logprob = logpdf(env[k], value) # log-probability from true distribution
            env_sample[k] = GrayBox.Sample(value, logprob)
        end
        a = ASTSampleAction(env_sample)
        (s, r) = @gen(:sp, :r)(mdp, s, a)
        R += r
        if BlackBox.isterminal(sim)
            break
        end
    end

    return -R # negative (loss)
end

function Base.convert(::Type{Dict{Symbol, Vector{Sampleable}}}, env::GrayBox.Environment, max_steps::Integer)
    dist_vector = Dict{Symbol, Vector{Sampleable}}()
    for k in keys(env)
        dist_vector[k] = fill(env[k], max_steps)
    end
    return dist_vector::Dict{Symbol, Vector{Sampleable}}
end

function Base.convert(::Type{Vector{GrayBox.Environment}}, distr::Dict{Symbol, Vector{Sampleable}}, max_steps::Integer=1)
    env_vector = GrayBox.Environment[]
    for t in 1:max_steps
        env = GrayBox.Environment()
        for k in keys(distr)
            env[k] = distr[k][t]
        end
        push!(env_vector, env)
    end
    return env_vector::Vector{GrayBox.Environment}
end

# Online work performed here.
function POMDPs.action(planner::CEMPlanner, s; rng=Random.GLOBAL_RNG)
    mdp::ASTMDP = planner.mdp
    if actiontype(mdp) != ASTSampleAction
        error("MDP action type must be ASTSampleAction to use CEM.")
    end

    env::GrayBox.Environment = GrayBox.environment(mdp.sim)

    # Importance sampling distributions, fill one per time step.
    is_dist_0 = convert(Dict{Symbol, Vector{Sampleable}}, env, planner.solver.episode_length)

    # Run cross-entropy method using importance sampling
    loss = (d, sample)->cem_losses(d, sample; mdp=mdp, initstate=s)
    is_dist_opt = cross_entropy_method(loss,
                                       is_dist_0;
                                       max_iter=planner.solver.n_iterations,
                                       N=planner.solver.num_samples,
                                       min_elite_samples=planner.solver.min_elite_samples,
                                       max_elite_samples=planner.solver.max_elite_samples,
                                       elite_thresh=planner.solver.elite_thresh,
                                       weight_fn=planner.solver.weight_fn,
                                       add_entropy=planner.solver.add_entropy,
                                       verbose=planner.solver.verbose,
                                       show_progress=planner.solver.show_progress,
                                       rng=rng)

    # Save the importance sampling distributions
    planner.is_dist = is_dist_opt

    # Pass back action trace if recording is on (i.e. top_k)
    if mdp.params.top_k > 0
        return get_top_path(mdp)
    else
        return planner.is_dist # pass back the importance sampling distributions
    end
end

"""
    AST.search!(planner::CEMPlanner)

Search for failure events using the `CEMPlanner` from an initial AST state.
Pass back the best action trace (or importance sampling distribution based on `mdp.params.top_k > 0`).

    AST.search!(planner::CEMPlanner, s::ASTState)

Search using CEM from a start state `s`.
"""
function AST.search!(planner::CEMPlanner)
    mdp::ASTMDP = planner.mdp
    Random.seed!(mdp.params.seed) # Determinism
    s = AST.initialstate(mdp)
    return search!(planner, s)
end

function AST.search!(planner::CEMPlanner, s::ASTState)
    mdp::ASTMDP = planner.mdp
    return action(planner, s)
end


text:
# GCM Initial Moisture Profiles
# This file contains helpers and lists currely avaiable options

abstract type AbstractMoistureProfile end
struct NoMoistureProfile <: AbstractMoistureProfile end
struct ZeroMoistureProfile <: AbstractMoistureProfile end
struct MoistLowTropicsMoistureProfile <: AbstractMoistureProfile end

# Helper for parsing `--init-moisture-profile` command line argument
function parse_moisture_profile_arg(arg)
    if arg === nothing
        moisture_profile = nothing
    elseif arg == "moist_low_tropics"
        moisture_profile = MoistLowTropicsMoistureProfile()
    elseif arg == "zero"
        moisture_profile = ZeroMoistureProfile()
    elseif arg == "dry"
        moisture_profile = NoMoistureProfile()
    else
        error("unknown moisture profile: " * arg)
    end

    return moisture_profile
end

# Initial moisture profile for a dry model setup
function init_moisture_profile(
    ::NoMoistureProfile,
    bl,
    state,
    aux,
    coords,
    t,
    p,
)
    FT = eltype(state)
    return FT(0)
end

# Initial moisture profile for a moist model setup with 0 initial moisture
function init_moisture_profile(
    ::ZeroMoistureProfile,
    bl,
    state,
    aux,
    coords,
    t,
    p,
)
    FT = eltype(state)
    return FT(0)
end

# Initial moisture profile following
# Ullrich et al. (2016) Dynamical Core Model Intercomparison Project (DCMIP2016) Test Case Document
function init_moisture_profile(
    ::MoistLowTropicsMoistureProfile,
    bl,
    state,
    aux,
    coords,
    t,
    p,
)
    FT = eltype(state)

    param_set = parameter_set(bl)
    _p_0::FT = MSLP(param_set)

    φ = latitude(bl, aux)

    # Humidity parameters
    p_w::FT = 34e3              # Pressure width parameter for specific humidity
    η_crit::FT = p_w / _p_0     # Critical pressure coordinate
    q_0::FT = 0.018             # Maximum specific humidity (default: 0.018)
    q_t::FT = 1e-12             # Specific humidity above artificial tropopause
    φ_w::FT = 2π / 9            # Specific humidity latitude wind parameter

    # get q_tot profile if needed
    η = p / _p_0                # Pressure coordinate η
    if η > η_crit
        q_tot = q_0 * exp(-(φ / φ_w)^4) * exp(-((η - 1) * _p_0 / p_w)^2)
    else
        q_tot = q_t
    end

    return q_tot
end


text:
using CayleyDickson
using Base.Test: @test, @test_throws

@test begin
    promote_rule(SplitBinion{Int}, Float64) == SplitBinion{Float64}
end

@test begin
    promote_rule(SplitBinion{Int}, SplitBinion{Float64}) == SplitBinion{Float64}
end

@test begin
    -(1.0 - SplitBinion(2, 3)) == SplitBinion(1.0, 3.0) == SplitBinion(2, 3) - 1.0
end

@test begin
    1.0 + SplitBinion(2, 3) == SplitBinion(3.0, 3.0) == SplitBinion(2, 3) + 1.0
end

@test begin
    convert(SplitBinion{Float64}, 1) == SplitBinion(1.0, 0)
end

@test begin
    convert(SplitBinion{Float64}, SplitBinion(1,2)) == SplitBinion(1.0, 2.0)
end

@test begin
    convert(SplitBinion, 2) == SplitBinion(2, 0)
end

@test begin
    z = SplitBinion{Int}(1, 2)
    convert(Hyperbolic{Int}, z) == z
end

@test begin
    convert(Int, SplitBinion(3)) == 3
end

@test_throws InexactError begin
    convert(Int, SplitBinion(3,4))
end

@test begin
    a = CayleyDickson.randomBigFloat()
    isreal(SplitBinion(a))
end

@test begin
    a = 1
    !isreal(SplitBinion(a, a))
end

@test begin
    z = random(SplitBinion{Int})
    z == +(z)
end

@test begin
    a = 1
    real(SplitBinion(a)) == a
end

@test begin
    a = 1
    l = unreal(SplitBinion(0, a))
    r = [a]
    l == r
end

@test begin
    io = IOBuffer()
    show(io, SplitBinion(1, 2))
    l = String(take!(io))
    r = "[1: 1, s: 2]"
    l == r
end

@test begin
    x = random(SplitBinion{BigFloat})
    y = random(SplitBinion{BigFloat})
    x != y
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(SplitBinion{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    z = random(SplitBinion{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    z = random(SplitBinion{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    z = random(SplitBinion{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    conj(conj(x)) == x
end

@test begin
    x = random(SplitBinion{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    l = conj(x * y)
    r = conj(y) * conj(x)
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    l = conj(x + y)
    r = conj(x) + conj(y)
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    l = conj(x - y)
    r = conj(x) - conj(y)
    l == r
end

@test begin
    x = random(SplitBinion{Rational{BigInt}})
    y = random(SplitBinion{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(SplitBinion{Rational{BigInt}})
    y = random(SplitBinion{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(SplitBinion{Rational{BigInt}})
    y = random(SplitBinion{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    z = random(SplitBinion{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = CayleyDickson.randomBigInt()
    x = random(SplitBinion{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    z = random(SplitBinion{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    z = random(SplitBinion{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    z = random(SplitBinion{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(SplitBinion{BigInt})
    quadrance(x) > big(0) || quadrance(x) < big(0) || quadrance(x) == big(0)
end

@test begin
    x = random(SplitBinion{BigInt})
    y = random(SplitBinion{BigInt})
    l = quadrance(x * y)
    r = quadrance(x) * quadrance(y)
    l == r
end

@test_throws ErrorException begin
    inv(SplitBinion(0, 0.0))
end

@test_throws ErrorException begin
    inv(SplitBinion(1, 1))
end

@test_throws ErrorException begin
    random(SplitBinion{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(SplitBinion{Int})
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(SplitBinion{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = CayleyDickson.randomBigRational()
    z = random(SplitBinion{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(SplitBinion{BigInt})
    z == z + zero(SplitBinion{BigInt})
end

@test begin
    z = random(SplitBinion{BigInt})
    z == z * one(SplitBinion{BigInt})
end

@test begin
    z = random(SplitBinion{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(SplitBinion{BigInt})
    z == z * one(z)
end

@test begin
    float(SplitBinion{Float64}) == SplitBinion{Float64}
end

@test begin
    x = random(SplitBinion{Rational{BigInt}})
    selfstar(selfstar(x)) == selfstar(x)
end

@test begin
    x = random(SplitBinion{Rational{BigInt}})
    anti_selfstar(anti_selfstar(x)) == anti_selfstar(x)
end


text:
module TestChangeset
type Post
    title::String
    author::String
end
end # module TestChangeset

import Ecto
import Ecto: Changeset
import Ecto.Changeset: change
import Ecto.Schema: Assoc
using Base.Test

post = %(TestChangeset.Post, author= "bar")
@test isa(post, Ecto.Schema.t)
changeset = change(post, title= "title")
@test Assoc([(:title, "title")]) == changeset.changes

changeset = change(%(TestChangeset.Post, title= "title"), title= "title")
@test Assoc([]) == changeset.changes


text:
module Speechbox

using WAV
using DSP

export specgram, enframe, mel2frq, frq2mel, melbankm, melfcc, vad, normalise_energy, write_binary, read_binary, combine_binary

function combine_binary(ofile,flist)
  if(isempty(flist))
    return nothing
  end
  data = read_binary(flist[1]);
  ndim,nvec = size(data);
  f = open(ofile,"w");
  write(f,ndim);
  write(f,nvec);
  write(f,data[:]);
  if(length(flist)>1)
    for i=2:length(flist)
      data = read_binary(flist[i]);
      nvec += size(data,2);
      write(f,data[:]);
    end
    seekstart(f);
    write(f,ndim);
    write(f,nvec);
    close(f);
  end
  return nothing
end

function read_binary(ofile)
  f = open(ofile,"r");
  ndim = read(f,Int64);
  nvec = read(f,Int64);
  feature_matrix = read(f,Float64,ndim,nvec);
  close(f)
  return feature_matrix;
end

function write_binary(ofile,feature_matrix)
  ndim,nvec = size(feature_matrix);
  f = open(ofile,"w");
  write(f,ndim);
  write(f,nvec);
  write(f,feature_matrix[:]);
  close(f);
  return nothing
end

function normalise_energy(x,fs)
  vindx = vad(x,fs);
  frames = enframe(x,fs);
  nframes = length(frames);
  en = zeros(nframes);
  for i=1:nframes
    en[i] = rms(frames[i]);
  end
  mn_en = mean(en[find(vindx)]);
  sc_factor = 0.1/mn_en;
  sx = sc_factor*x;
  return sx;
end

function vad(x,fs,win_dur=0.02,win_overlap=0.01,energy_thr=0.05)
  frames = enframe(x,fs,win_dur,win_overlap);
  nframes = length(frames);
  energy = zeros(nframes);
  for i=1:nframes
    energy[i] = rms(frames[i]);
  end
  max_energy = maximum(energy);
  abs_thr = max_energy*energy_thr;
  v_indx = zeros(Int,nframes);
  v_indx[find(energy.>abs_thr)]=1;
  return v_indx
end

function melfcc(x,fs,ncoef=13,nfilt=17,win_dur=0.02,win_overlap=0.01)
  frames = enframe(x,fs,win_dur,win_overlap);
  nframes = length(frames);
  flen = length(frames[1]);
  nfft = nextfastfft(flen);
  win = hamming(flen);

  buf = zeros(nfft);
  fbuf = zeros(nfilt);
  rfp = plan_rfft(buf);
  dcp = plan_dct(fbuf);

  nrfft = length(rfp*buf);
  fbank = melbankm(fs,nrfft,nfilt); #generate mel filterbank filters
  fbank = fbank.^2; #squaring triangular mel-filters to multiple with power spectrum

  mfcc = zeros(ncoef,nframes);

  for i=1:nframes
    buf[1:1:flen] = win.*frames[i];
    #fbuf = dcp*(fbank*(abs.(rfp*buf)));
    fbuf = fbank*(abs2.(rfp*buf) + eps()); #multiplying with power spectrum (square of mag spectrum) and accumulating
    fbuf = dcp*(log.(fbuf));
    mfcc[:,i] = fbuf[1:ncoef];
  end
  return mfcc
end

function melbankm(fs,npts,nfilt=17)
  mspace = linspace(0,frq2mel(fs/2),nfilt); #equally spaced points in mel scale
  fspace = mel2frq.(mspace); #equal mel spaced points mapped back to Hz
  cindx = Int.(round.(fspace*(npts-1)/(fs/2))+1); #Filter centre indices (first at 0, final at Fs/2)
  fbank = zeros(nfilt,npts);
  #define triangular filters for 2 to nfilt-1
  for i=2:nfilt-1
    fbank[i,cindx[i-1]:cindx[i]] = linspace(0,1,cindx[i]-cindx[i-1]+1);
    fbank[i,cindx[i]:cindx[i+1]] = linspace(1,0,cindx[i+1]-cindx[i]+1);
  end
  #one sided triangular filters for first and last filters
  fbank[1,cindx[1]:cindx[2]] = linspace(1,0,cindx[2]-cindx[1]+1);
  fbank[nfilt,cindx[nfilt-1]:cindx[nfilt]] = linspace(0,1,cindx[nfilt]-cindx[nfilt-1]+1);
  return fbank
end

function specgram(x,fs,win_dur=0.02,win_overlap=0.01,wtype="hamm")
  frames = enframe(x,fs,win_dur,win_overlap);
  nframes = length(frames);
  flen = length(frames[1]);
  nfft = nextfastfft(flen);
  if(wtype=="rect")
    win = ones(flen);
  else
    win = hamming(flen);
  end
  buf = zeros(nfft);
  rfp = plan_rfft(buf);
  nrfft = length(rfp*buf);
  mspec = zeros(nrfft,nframes);
  for i=1:nframes
    buf[1:1:flen] = win.*frames[i];
    mspec[:,i] = abs.(rfp*buf);
  end
  return mspec + eps();
end

function enframe(x,fs,win_dur=0.02,win_overlap=0.01)
  wlen = Int(round(win_dur*fs));
  wolap = Int(round(win_overlap*fs));
  x = x[:];
  frames = arraysplit(x,wlen,wolap);
  return frames
end

frq2mel(frq) = log(1+frq/700)*1127.01048;

mel2frq(mel) = 700*(exp(mel/1127.01048)-1);

end