text:
module AnalyticalWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractAnalyticalPropagator = @jimport org.orekit.propagation.analytical.AbstractAnalyticalPropagator
const AbstractPropagator = @jimport org.orekit.propagation.AbstractPropagator
const AdapterPropagator = @jimport org.orekit.propagation.analytical.AdapterPropagator
const AdapterPropagator_DifferentialEffect = @jimport org.orekit.propagation.analytical.AdapterPropagator$DifferentialEffect
const AdditionalStateProvider = @jimport org.orekit.propagation.AdditionalStateProvider
const AggregateBoundedPropagator = @jimport org.orekit.propagation.analytical.AggregateBoundedPropagator
const AttitudeProvider = @jimport org.orekit.attitudes.AttitudeProvider
const BoundedPropagator = @jimport org.orekit.propagation.BoundedPropagator
const CartesianOrbit = @jimport org.orekit.orbits.CartesianOrbit
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const EcksteinHechlerPropagator = @jimport org.orekit.propagation.analytical.EcksteinHechlerPropagator
const Ephemeris = @jimport org.orekit.propagation.analytical.Ephemeris
const EventDetector = @jimport org.orekit.propagation.events.EventDetector
const Field = @jimport org.hipparchus.Field
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldAbstractAnalyticalPropagator = @jimport org.orekit.propagation.analytical.FieldAbstractAnalyticalPropagator
const FieldAbstractPropagator = @jimport org.orekit.propagation.FieldAbstractPropagator
const FieldAdditionalStateProvider = @jimport org.orekit.propagation.FieldAdditionalStateProvider
const FieldBoundedPropagator = @jimport org.orekit.propagation.FieldBoundedPropagator
const FieldCartesianOrbit = @jimport org.orekit.orbits.FieldCartesianOrbit
const FieldEcksteinHechlerPropagator = @jimport org.orekit.propagation.analytical.FieldEcksteinHechlerPropagator
const FieldEventDetector = @jimport org.orekit.propagation.events.FieldEventDetector
const FieldKeplerianPropagator = @jimport org.orekit.propagation.analytical.FieldKeplerianPropagator
const FieldOrbit = @jimport org.orekit.orbits.FieldOrbit
const FieldOrekitFixedStepHandler = @jimport org.orekit.propagation.sampling.FieldOrekitFixedStepHandler
const FieldOrekitStepHandler = @jimport org.orekit.propagation.sampling.FieldOrekitStepHandler
const FieldPVCoordinatesProvider = @jimport org.orekit.utils.FieldPVCoordinatesProvider
const FieldSpacecraftState = @jimport org.orekit.propagation.FieldSpacecraftState
const Frame = @jimport org.orekit.frames.Frame
const J2DifferentialEffect = @jimport org.orekit.propagation.analytical.J2DifferentialEffect
const JString = @jimport java.lang.String
const KeplerianPropagator = @jimport org.orekit.propagation.analytical.KeplerianPropagator
const List = @jimport java.util.List
const Object = @jimport java.lang.Object
const Orbit = @jimport org.orekit.orbits.Orbit
const OrekitFixedStepHandler = @jimport org.orekit.propagation.sampling.OrekitFixedStepHandler
const OrekitStepHandler = @jimport org.orekit.propagation.sampling.OrekitStepHandler
const PVCoordinatesProvider = @jimport org.orekit.utils.PVCoordinatesProvider
const PropagationType = @jimport org.orekit.propagation.PropagationType
const Propagator = @jimport org.orekit.propagation.Propagator
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const TimeStampedFieldPVCoordinates = @jimport org.orekit.utils.TimeStampedFieldPVCoordinates
const TimeStampedPVCoordinates = @jimport org.orekit.utils.TimeStampedPVCoordinates
const UnnormalizedSphericalHarmonicsProvider = @jimport org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider
const UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics = @jimport org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider$UnnormalizedSphericalHarmonics

export AbstractAnalyticalPropagator
export AdapterPropagator
export AggregateBoundedPropagator
export EcksteinHechlerPropagator
export Ephemeris
export FieldAbstractAnalyticalPropagator
export FieldEcksteinHechlerPropagator
export FieldKeplerianPropagator
export J2DifferentialEffect
export KeplerianPropagator
export add_additional_state_provider
export add_effect
export add_event_detector
export apply
export basic_propagate
export clear_events_detectors
export equals
export get_additional_state_providers
export get_attitude_provider
export get_class
export get_effects
export get_events_detectors
export get_extrapolation_threshold
export get_field
export get_frame
export get_generated_ephemeris
export get_initial_state
export get_managed_additional_states
export get_max_date
export get_min_date
export get_mode
export get_propagator
export get_pv_coordinates
export get_pv_provider
export hash_code
export is_additional_state_managed
export notify
export notify_all
export propagate
export propagate_orbit
export reset_initial_state
export set_attitude_provider
export set_ephemeris_mode
export set_master_mode
export set_slave_mode
export to_string
export wait

include("GnssWrapper/GnssWrapper.jl")
include("TleWrapper/TleWrapper.jl")
include("abstract_analytical_propagator.jl")
include("adapter_propagator.jl")
include("aggregate_bounded_propagator.jl")
include("eckstein_hechler_propagator.jl")
include("ephemeris.jl")
include("field_abstract_analytical_propagator.jl")
include("field_eckstein_hechler_propagator.jl")
include("field_keplerian_propagator.jl")
include("j2_differential_effect.jl")
include("keplerian_propagator.jl")

end


text:
function integer_regex(fm::CSVFormat)
    if fm.thousands_separator == nothing
        return r"^\s*-?\d+\s*$"
    else
        if fm.thousands_separator == '.'
            return r"^\s*-?((\d{1,3}\.)?(\d{3}\.)*\d{3}|\d{1,3})\s*$"
        elseif fm.thousands_separator == ','
            return r"^\s*-?((\d{1,3},)?(\d{3},)*\d{3}|\d{1,3})\s*$"
        else
            error("thousands_separator ($(fm.thousands_separator)) not supported.")
        end
    end
end

function float_regex(fm::CSVFormat)
    if fm.thousands_separator == nothing
        if fm.decimal_separator == '.'
            return r"^\s*-?(\d+(\.\d*)?|\d*\.\d+)\s*$"
        elseif fm.decimal_separator == ','
            return r"^\s*-?(\d+(,\d*)?|\d*,\d+)\s*$"
        else
            error("Decimal separator not supported: $(fm.decimal_separator)")
        end
    else
        if fm.thousands_separator == '.' && fm.decimal_separator == ','
            return r"^\s*-?((\d{1,3}\.)?(\d{3}\.)*\d{3}|\d{1,3})(,\d*)?\s*$"
        elseif fm.thousands_separator == ',' && fm.decimal_separator == '.'
            return r"^\s*-?((\d{1,3},)?(\d{3},)*\d{3}|\d{1,3})(\.\d*)?\s*$"
        else
            error("Combination of thousands_separator '$(fm.thousands_separator)' and decimal_separator '$(fm.decimal_separator)' not supported.")
        end
    end
end

mutable struct InferenceState
    datatype::DataType
    allow_missing::Bool
end

InferenceState() = InferenceState(Any, false)

function datatype(s::InferenceState)
    if s.allow_missing
        return Union{Missing, s.datatype}
    else
        return s.datatype
    end
end

Base.:(==)(s1::InferenceState, s2::InferenceState) = s1.datatype == s2.datatype && s1.allow_missing == s2.allow_missing

function infer_type(value::String, format::CSVFormat, state::InferenceState=InferenceState(), fr::Regex=float_regex(format), ir::Regex=integer_regex(format)) :: InferenceState
    inferred_type = Any
    inferred_missing = false

    if ismissingstr(value, format)
        inferred_missing = true
    else
        # If it's already a String, there's nothing to improve with inference
        if state.datatype == String
            inferred_type = String
        else

            if occursin(ir, value)
                inferred_type = Int
            elseif occursin(fr, value)
                inferred_type = Float64
            else
                try
                    Date(value, format.date_format)
                    inferred_type = Date
                catch e
                    inferred_type = String
                end
            end
        end
    end

    # Promote type
    if type_order(inferred_type) > type_order(state)
        state.datatype = inferred_type
    end

    # Promote Missing
    state.allow_missing = state.allow_missing || inferred_missing
    return state
end

type_order(state::InferenceState) = type_order(state.datatype)
type_order(::Type{Any}) = 0
type_order(::Type{Int}) = 1
type_order(::Type{Float64}) = 2
type_order(::Type{Date}) = 3
type_order(::Type{String}) = 4

function infer_schema(raw::Array{String, 2}, format::CSVFormat=CSVFormat(), header::Bool=true) :: Schema
    FST_DATAROW = header ? 2 : 1
    fr = float_regex(format)
    ir = integer_regex(format)
    local s::InferenceState
    (rows, cols) = size(raw)
    @assert rows >= FST_DATAROW "Empty data"

    schema_header = Vector{Symbol}(undef, cols)
    if header
        @inbounds for c in 1:cols
            schema_header[c] = Symbol(raw[1,c])
        end
    else
        @inbounds for c in 1:cols
            schema_header[c] = Symbol(c)
        end
    end

    schema_datatypes = Vector{Type}(undef, cols)

    @inbounds for c in 1:cols
        s = InferenceState()

        for r in FST_DATAROW:rows
            infer_type(raw[r,c], format, s, fr, ir)
        end

        schema_datatypes[c] = datatype(s)
    end

    return Schema(schema_header, schema_datatypes)
end

function infer_schema(input, format::CSVFormat=CSVFormat(); header::Bool=true, use_mmap::Bool=false)
    raw = readdlm(input, format.dlm, String; use_mmap=use_mmap)
    return infer_schema(raw, format, header)
end

ext: jl


text:
using SpectralKit
using Test
import ForwardDiff

####
#### utility functions for tests
####

function test_roots(family, N; atol = 1e-13)
    r = roots(family, N)
    @test r isa Vector{Float64}
    @test length(r) == N
    @test all(abs.(evaluate.(family, N, r, Val(0))) .≤ atol)
end

function test_augmented_extrema(family, N; atol = 1e-13)
    a = augmented_extrema(family, N)
    @test a isa Vector{Float64}
    @test length(a) == N
    @test all(abs.(evaluate.(family, N - 1, a[2:(end-1)], Val(1))) .≤ atol)
    mi, ma = domain_extrema(family)
    â = sort(a)                 # non-increasing transformations switch signs
    @test mi ≈ â[1] atol = 1e-16
    @test ma ≈ â[end] atol = 1e-16
end

function test_endpoint_continuity(family, expected_extrema, ns;
                                  eps_step = eps(), inf_proxy = 1e6, atol = eps()^0.25)
    mi, ma = domain_extrema(family)
    @test mi == expected_extrema[1]
    @test ma == expected_extrema[2]
    for n in ns
        ev(x, order) = evaluate(family, n, x, order)
        # at domain minimum
        f_mi, f′_mi = @inferred ev(mi, Val(0:1))
        mi_approx = isfinite(mi) ? mi + eps() : -inf_proxy
        @test f_mi isa Float64
        @test f′_mi isa Float64
        @test f_mi ≈ ev(mi, Val(0))
        @test f_mi ≈ ev(mi_approx, Val(0)) atol = atol
        @test f′_mi ≈ ev(mi, Val(1))
        @test f′_mi ≈ ev(mi_approx, Val(1)) atol = atol

        # at domain maximum
        f_ma, f′_ma = @inferred ev(ma, Val(0:1))
        ma_approx = isfinite(ma) ? ma - eps() : inf_proxy
        @test f_ma isa Float64
        @test f′_ma isa Float64
        @test f_ma ≈ ev(ma, Val(0))
        @test f_ma ≈ ev(ma_approx, Val(0)) atol = atol
        @test f′_ma ≈ ev(ma, Val(1))
        @test f′_ma ≈ ev(ma_approx, Val(1)) atol = atol
    end
end

function test_derivatives(family, generator, ns; M = 100)
    for n in ns
        for _ in 1:M
            x = generator()
            f, f′ = @inferred evaluate(family, n, x, Val(0:1))
            @test f ≈ @inferred evaluate(family, n, x, Val(0))
            @test f′ ≈ @inferred evaluate(family, n, x, Val(1))
            @test f′ ≈ ForwardDiff.derivative(x -> evaluate(family, n, x, Val(0)), x)
        end
    end
end

####
#### tests
####

@testset "Chebyshev" begin
    F = Chebyshev()

    test_roots(F, 11)
    @test roots(F, 11)[6] == 0  # precise 0

    test_augmented_extrema(F, 11)
    @test augmented_extrema(F, 11)[6] == 0 # precise 0

    test_endpoint_continuity(F, (-1, 1), 0:10)

    test_derivatives(F, () -> rand() * 2 - 1, 0:10)
end

@testset "ChebyshevSemiInf" begin
    F = ChebyshevSemiInf(2.0, 4.7)

    test_roots(F, 9)

    test_augmented_extrema(F, 10)

    test_endpoint_continuity(F, (2.0, Inf), 0:10; atol = 1e-3)

    test_derivatives(F, () -> 2.0 + abs2(randn()), 0:10)

    F = ChebyshevSemiInf(3.0, -1.9)

    test_roots(F, 11; atol = 1e-13)

    test_augmented_extrema(F, 7; atol = 1e-10)

    test_endpoint_continuity(F, (-Inf, 3.0), 0:10; atol = 1e-3)

    test_derivatives(F, () -> 3.0 - abs2(randn()), 0:10)

    @test_throws ArgumentError ChebyshevSemiInf(0.0, 0.0)
end

@testset "ChebyshevInf" begin
    F = ChebyshevInf(0.0, 1.0)

    test_roots(F, 11)
    @test roots(F, 11)[6] == 0  # precise 0

    test_augmented_extrema(F, 11)
    @test augmented_extrema(F, 11)[6] == 0 # precise 0

    test_endpoint_continuity(F, (-Inf, Inf), 0:10)

    test_derivatives(F, randn, 0:10)

    @test_throws ArgumentError ChebyshevInf(0.0, -3.0)
    @test_throws ArgumentError ChebyshevInf(0.0, 0.0)
end

@testset "ChebyshevInterval" begin
    F = ChebyshevInterval(2.0, 5)

    test_roots(F, 11)

    test_augmented_extrema(F, 11)

    test_endpoint_continuity(F, (2.0, 5.0), 0:10)

    test_derivatives(F, () -> rand() * 3 + 2, 0:10)

    @test_throws ArgumentError ChebyshevInterval(2.0, 1.0)
    @test_throws ArgumentError ChebyshevInterval(2.0, 2)
    @test_throws ArgumentError ChebyshevInterval(-Inf, Inf)
end


text:
using Knet,JLD,ArgParse
using Knet: readvocab, lookup

function main(args=ARGS)
    s = ArgParseSettings()
    @add_arg_table s begin
        ("--modelfile"; required=true)
        ("--dictfiles"; nargs=2; required=true)
        ("--beamsize"; arg_type=Int; default=16)
        ("--nbest"; arg_type=Int; default=1)
    end
    isa(args, AbstractString) && (args=split(args))
    opts = parse_args(args, s)
    model = gpucopy(load(opts["modelfile"],"model"))
    dict1 = readvocab(opts["dictfiles"][1])
    dict2 = readvocab(opts["dictfiles"][2])
    for l in eachline(STDIN)
        word = map(c->dict1[c], split(l))
        pred = predict(model, Any[word]; beamsize=opts["beamsize"], nbest=opts["nbest"], dense=true)
        pstr = lookup(pred[1], dict2)
        println(join(map(x->join(x," "), pstr), "\t"))
    end
end

!isinteractive() && !isdefined(:load_only) && main(ARGS)


text:
module simple_debugger
using MLStyle
using Crayons
using PrettyPrinting
using ReplMaker
using REPL

export @debug!, @nodebug!, @watch!, @report, @watch_locals!, @inject

DebugScope = Dict{String, Tuple{LineNumberNode, Any}}

mutable struct MetaData
    current_line :: LineNumberNode
    to_watch     :: DebugScope
    committed    :: Vector{DebugScope}
    debug        :: Bool
end

const line = LineNumberNode(32, :fake_line)
const empty_dbg_scope = DebugScope()
const meta = MetaData(line, empty_dbg_scope, [], false)

Base.show(io::IO, ::MetaData) = print(io, "meta")

macro debug!()
    meta.debug = true
end

macro nodebug!()
    meta.debug = false
end

function commit!(meta)
    push!(meta.committed, meta.to_watch)
end

function rollback!(meta)
    meta.to_watch = pop!(meta.committed)
end

function watch!(meta, name, var, line)
    meta.to_watch[name] = (line, var)
    var
end

const fake_exps    = Expr[]
const dbg_scoperef = Ref{Vector{Expr}}(fake_exps)

function parse_expr(s)
    ex = Meta.parse(s)
    let_bindings = dbg_scoperef.x
    :(let $(let_bindings...); $ex end)
end

const dbg_repl = Base.active_repl

dbg_mode = initrepl(
    parse_expr,
    prompt_text="dbg> ",
    prompt_color = :green,
    start_key="|",
    repl = dbg_repl,
    mode_name="debug_mode"
)

function report(meta, msg)
    bold = Crayon(bold=true)
    style_title = Crayon(foreground=:light_magenta)
    style_line = Crayon(foreground=:green, italics=true, underline=true)
    style_var = Crayon(foreground=:light_blue, bold=true)
    println(args...) = begin
        Base.print(args...)
        Base.print(Crayon(reset=true), "\n")
    end
    println(style_title, "Debugging ", bold, msg)
    println(style_line, "[where] ", meta.current_line)

    for (k, (line, v)) in meta.to_watch
        print(style_var, "  ", k, " ", style_line, bold, "[", line, "]:", Crayon(reset=true))
        pprint(v)
        println()
    end
    print(style_title, "Open a REPL here and check manually?[Y/N]")
    option = stdin |> readline |> strip |> lowercase |> Symbol
    @match option begin
        :n => ()
        :y =>
            begin
                dbg_scoperef.x = map(meta.to_watch |> collect) do (k, (line, v))
                    sym = Symbol(k)
                    :($sym = $v)
                end
                println(
                    style_var,
                    "Press | and switch to dbg-mode, enter Ctrl+D to leave sub-REPL."
                )
                REPL.run_repl(dbg_repl)
            end
        _ => throw(error("invalid option $option"))
    end
    nothing
end

function inject_func(msg, expr, meta)
    quote
        $commit!($meta)
        $meta.to_watch = $DebugScope()
        try
            $(inject(expr, meta))
        finally
            if $meta.debug && $(!isempty)($meta.to_watch)
                $report($meta, $msg)
            end
            $rollback!($meta)
        end
    end
end

is_call(ex) =
    @match ex begin
        Expr(:call, _...) => true
        :($a :: $_)       => is_call(a)
        :($a where $_)    => is_call(a)
        _                 => false
    end

function inject(expr, meta)
    rec(ex) = inject(ex, meta)

    @match expr begin
        Expr(:(=), func_head, body) && if is_call(func_head) end ||
        Expr(:function || :(->), func_head, body) =>
            let body = inject_func(string(func_head), body, meta)
                Expr(:function, func_head, body)
            end

        line::LineNumberNode =>
            Expr(:block, line, :($meta.current_line = $(QuoteNode(line))))

        :(@report $line) =>
            quote
                if $meta.debug && $(!isempty)($meta.to_watch)
                    $report($meta, $msg)
                end
            end

        :(@watch_locals! $(_...)) =>
            quote
                for (k, v) in $Base.@locals
                    $watch!($meta, $string(k), v, $(QuoteNode(line)))
                end
            end

        :(@watch! $line $var) || :(@watch $var) =>
            let name = string(var)
                :($watch!($meta, $name, $var, $(QuoteNode(line))))
            end

        Expr(hd, tl...) =>
            Expr(hd, map(rec, tl)...)

        a => a
    end
end

macro watch!()
    throw(error("only used as tokens of @inject"))
end

macro report!()
    throw(error("only used as tokens of @inject"))
end

macro watch_locals!()
    throw(error("only used as tokens of @inject"))
end

macro inject(ex)
    inject(ex, meta) |> esc
end

end

max_stars_repo_name: thautwarm/MLStyle-Playground

max_issues_repo_name: thautwarm/MLStyle-Playground

max_forks_repo_name: thautwarm/MLStyle-Playground


text:
using Pathfinder
using Test

@testset "Pathfinder.jl" begin
    include("woodbury.jl")
    include("maximize.jl")
    include("inverse_hessian.jl")
    include("mvnormal.jl")
    include("elbo.jl")
    include("resample.jl")
    include("singlepath.jl")
    include("multipath.jl")
end


text:
struct PoincareInvariant1stCanonical{ET,DT,TT}
    equ::ET
    Δt::TT
    nloop::Int
    ntime::Int
    nsave::Int
    nt::Int
    I::OffsetArray{DT,1,Vector{DT}}
end

function PoincareInvariant1stCanonical(f_equ::Function, f_loop_q::Function, f_loop_p::Function, Δt::TT, d::Int, nloop::Int, ntime::Int, nsave::Int, DT=Float64) where {TT}

    if get_config(:verbosity) > 1
        println()
        println("First Euler-Poincaré Integral Invariant")
        println("=======================================")
        println()
        println(" nloop = ", nloop)
        println(" ntime = ", ntime)
        println(" nsave = ", nsave)
        println(" Δt    = ", Δt)
        println()
    end

    # compute initial conditions
    q₀ = [f_loop_q(i/nloop) for i in 1:nloop]
    p₀ = [f_loop_p(i/nloop) for i in 1:nloop]

    # initialise euation
    equ = f_equ(q₀, p₀)

    # create arrays for results
    nt = div(ntime, nsave)

    I = OffsetArray(zeros(DT, nt+1), 0:nt)

    PoincareInvariant1stCanonical{typeof(equ),DT,TT}(equ, Δt, nloop, ntime, nsave, nt, I)
end

function evaluate_poincare_invariant(pinv::PoincareInvariant1stCanonical, sol::SolutionPODE)
    v = zeros(size(sol.q[begin],1), size(sol.q,2))

    for i in axes(sol.q,1)
        compute_velocity(hcat(sol.q[i,:]...), v)
        pinv.I[i] = compute_loop_integral(hcat(sol.p[i,:]...), v)
    end

    pinv.I
end

function write_to_hdf5(pinv::PoincareInvariant1stCanonical, sol::Solution, output_file::String)
    # h5open(output_file, isfile(output_file) ? "r+" : "w") do h5
    h5open(output_file, "w") do h5
        write(h5, "t", sol.t)
        write(h5, "I", pinv.I)
    end
end


text:
#
# Subject
#

a = Subject(10)
b = Subject(30)
@test get(a) == 10
@test _map(a, x -> x + 2) == Subject(12)
@test _filter(a, x -> x > 2) == a
@test _filter(a, x -> x > 20) == Subject(nothing)
@test _reduce(a, b, (x, y) -> (x + y)/2) == Subject(20.0)

#
# init_task
#

t = init_task(() -> 5)
@test consume(t) == 5

t = init_task(x -> x + 5, 1)
@test consume(t) == 6

t = init_task((x, y) -> x*y + 5, 2, 3)
@test consume(t) == 11


text:
using POMDPModels
using POMDPs

let
    problem = BabyPOMDP(-5, -10, 0.1, 0.8, 0.1, 0.9)

    n = 100000
    results = Vector{Float64}(n)

    i = 1
    rng = MersenneTwister(i)
    init_state = false
    od = observation(problem, init_state, false, init_state)
    obs = rand(rng, od)
    sim = RolloutSimulator(rng=rng, initialstate=init_state, eps=0.0001)
    policy = FeedWhenCrying()
    results[i] = simulate(sim, problem, policy, updater(policy), PreviousObservation(obs))

    rngs = [MersenneTwister(i) for i in 1:n]

    @time for i in 1:n
        rng = rngs[i]
        init_state = false
        od = observation(problem, init_state, false, init_state)
        obs = rand(rng, od)
        sim = RolloutSimulator(rng=rng, initialstate=init_state, eps=0.0001)
        policy = FeedWhenCrying()
        results[i] = simulate(sim, problem, policy, updater(policy), PreviousObservation(obs))
    end
end


text:
# Pkg.clone("git://github.com/lejon/TSne.jl.git")
# Pkg.add("Gadfly")


text:
include("css.jl")
include("formats.jl")
using Lathe.stats: mean
using Dates
import Base: getindex, show, Meta
#=============
OddFrame Type
=============#
mutable struct OddFrame <: AbstractOddFrame
        labels::Array{Symbol}
        columns::Array{Any}
        coldata::Array{Pair}
        head::Function
        drop::Function
        dropna::Function
        #==
        Constructors
        ==#
        function OddFrame(p::Pair ...)
                # Labels/Columns
                labels = [pair[1] for pair in p]
                columns = [pair[2] for pair in p]
                length_check(columns)
                name_check(labels)
                types = [typeof(x[1]) for x in columns]
                # coldata
                coldata = generate_coldata(columns, types)
                # Head
                head(x::Int64) = _head(labels, columns, coldata, x)
                head() = _head(labels, columns, coldata, 5)
                # Drop
                drop(x) = _drop(x, columns)
                drop(x::Symbol) = _drop(x, labels, columns, coldata)
                drop(x::String) = _drop(Symbol(x), labels, columns, coldata)
                dropna() = _dropna(columns)
                # type
                new(labels, columns, coldata, head, drop, dropna);
        end
        function OddFrame(file_path::String)
                # Labels/Columns
                extensions = Dict("csv" => read_csv)
                extension = split(file_path, '.')[2]
                labels, columns = extensions[extension](file_path)
                length_check(columns)
                name_check(labels)
                types, columns = read_types(columns)
                # Coldata
                coldata = generate_coldata(columns, types)
                # Head
                head(x::Int64) = _head(labels, columns, coldata, x)
                head() = _head(labels, columns, coldata, 5)
                # Drop
                drop(x) = _drop(x, columns)
                drop(x::Symbol) = _drop(x, labels, columns, coldata)
                drop(x::String) = _drop(Symbol(x), labels, columns, coldata)
                dropna() = _dropna(columns)
                # type
                new(labels, columns, coldata, head, drop, dropna);
        end
        #==
        Supporting
                Functions
        ( Support constructors )
        ==#
        function generate_coldata(columns::Array, types::Array)
                pairs = []
                for (i, T) in enumerate(types)
                        if T == String
                                push!(pairs, T => string("Data-type: ",
                                T, "\nFeature Type: Categorical\n",
                                "Categories: ", length(Set(columns[i]))))
                        elseif T == Bool
                                push!(pairs, T => string("Data-type: ",
                                T, "\nFeature Type: Categorical\n",
                                "Categories: ", length(Set(columns[i]))))
                        elseif length(columns[i]) / length(Set(columns[i])) <= 1.8
                                push!(pairs, T => string("Data-type: ",
                                T, "\nFeature Type: Continuous\n", "Mean: ",
                                mean(columns[i])))
                        else
                                push!(pairs, T => string("Data-type: ",
                                T, "\nFeature Type: Categorical\n",
                                "Categories: ", length(Set(columns[i]))))
                        end
                end
                pairs
        end

        #==
        THROWS
        ==#
        function length_check(ps)
                ourlen = length(ps[1])
[if length(x) != ourlen throw(DimensionMismatch("Columns must be the same size")) end for x in ps]
        end

        function name_check(labels)
                if length(Set(labels)) != length(labels)
                throw(ErrorException("Column names may not be duplicated!"))
                end
        end
        #==
        Child
            methods
            ==#
        function _head(labels::Array{Symbol},
                columns::Array, coldata::Array, count::Int64)
                # Create t-header and t-body tags
                thead = "

"
                tbody = "
"
                # populate row headers
                [thead = string(thead, "",  string(name),
                 "	") for (n, name) in enumerate(labels)]
                 # finish t-head
                 thead = string(thead, "

")
                 # populate each row iteratively.
                 for i in 1:count
                         obs = [row[i] for row in columns]
                         tbody = string(tbody, "
")
                         [tbody = string(tbody, ""
                         , observ,
        "	") for (count, observ) in enumerate(obs)]
                        tbody = string(tbody, "
")
                 end
                 # Finish tags:
                 tbody = string(tbody, "
")
                 final = string("
", thead, tbody,
                  "
", _css)
                 # Display
                 # TODO: Figure out how to determine whether one is in
                 #    the REPL, prefereably before any of this function is
                 # called.
                 display("text/html", final)
        end

        function _drop(column::Symbol, labels::Array{Symbol}, columns::Array,
                coldata::Array)
                pos = findall(x->x==column, labels)[1]
                deleteat!(labels, pos)
                deleteat!(coldata, pos)
                deleteat!(columns, pos)
        end

        function _drop(row::Int64, columns::Array)
                [deleteat!(col, row) for col in columns]
        end

        function _drop(row::Array, columns::Array)
                [deleteat!(col, row) for col in columns]
        end

        function _dropna(columns::Array)
                for col in columns
                        mask = [ismissing(x) for x in col]
                        pos = findall(x->x==0, mask)
                        _drop(pos, columns)
                end
        end

end

#===
Indexing
===#
function getindex(od::AbstractOddFrame, col::Symbol)
        pos = findall(x->x==col, od.labels)[1]
        return(od.columns[pos])
end
getindex(od::AbstractOddFrame, col::String) = od[Symbol(col)]
getindex(od::AbstractOddFrame, axis::Int64) = od.columns[axis]
function getindex(od::OddFrame, mask::BitArray)
        pos = findall(x->x==0, mask)
        od.drop(pos)
end

#===
Iterators
===#
columns(od::OddFrame) = od.columns
labels(od::OddFrame) = od.labels
names(od::OddFrame) = od.labels
pairs(od::OddFrame) = [od.labels[i] => od.columns[i] for i in 1:length(od.labels)]


text:
"""
Data structure for a mesh material that is used to desribe the light components of a [Phong reflection](https://en.wikipedia.org/wiki/Phong_reflection_model)
type model. All data is stored as RGBα for Red, Green, Blue and transparency.

# Exemples
> a= Component_light(0.5,0.2,0.1,0.0)
> mat= material(a,a,a,a,0.5)
> mat.shininess
0.5
"""
struct material
    emission::Array{AbstractFloat,1}
    ambiant::Array{AbstractFloat,1}
    diffuse::Array{AbstractFloat,1}
    specular::Array{AbstractFloat,1}
    shininess::AbstractFloat
end


text:
module DifferentialEvolution
using ProgressMeter

# Write your package code here.
export Evolver, Options, Controller, DEEvolver, DEOptions

export pop, fitness, senses, bounds, options, nindividuals, nobjectives, best_fitness, best_individual,
    StdDE, StdOptions, mutate_rand_1, mutate_rand_1!, crossover_binomial!, StdSoController, evolve!,
    nondominated_sort, get_nondomination_rank, crowded_compare, dominate, assign_crowding_distance,
    bounce_back!, demo!, StdSoDE, StdMoDE

include("common.jl")
include("utils.jl")
include("std_de.jl")
include("mutation.jl")
include("crossover.jl")
include("constraint.jl")
include("mo/nsga_ii.jl")
include("mo/mo.jl")
include("mo/demo.jl")
include("mo/debi.jl")
end


text:
module Optimization

using ..Diagram
using ..Fortune
using ..Intersect

pointsₖ = nothing # temporary points array for backtracking

function norm2(v)
   sum = 0

   for vᵢ in v
      sum += vᵢ[1]^2 + vᵢ[2]^2
   end

   return sum
end

function normInf(v)
   max = 0

   for vᵢ in v
      if (vᵢ[1]^2) > max
         max = vᵢ[1]^2
      end
      if (vᵢ[2]^2) > max
         max = vᵢ[2]^2
      end
   end

   return max
end

function errorInTriangle(z, a, b)
   p = [a[1] - z[1], a[2] - z[2]]
   q = [b[1] - z[1], b[2] - z[2]]

   jacobian = abs(p[1]q[2] - p[2]q[1])

   # error = c∫∫x² + y² dudv, see the monograph appendix for a derivation
   # u and v are variables such that Triangle(x,y) -> Triangle(u,v) right triangle with catheti of size 1
   # c is just so that I can take the common factor 1/12 out of both integrals and multiply once

   c∫x²dudv = p[1]^2 + q[1]^2 + p[1]*q[1]
   c∫y²dudv = p[2]^2 + q[2]^2 + p[2]*q[2]

   return jacobian * ((c∫x²dudv + c∫y²dudv) / 12) # testing with area
end

function errorInRegion(region::Diagram.Region)
   z = region.generator

   edge = region.borderHead

   sum = errorInTriangle(z, edge.origin, edge.next.origin)
   edge = edge.next
   while edge != region.borderHead
      sum += errorInTriangle(z, edge.origin, edge.next.origin)
      edge = edge.next
   end

   return sum
end

function defaultf(V::Diagram.DCEL)
   sum = 0
   for region in V.regions
      sum += errorInRegion(region)
   end

   return sum
end

function default∇f(V::Diagram.DCEL, points::Array{Tuple{Real, Real}, 1})
   function weightedDifference(a::Array{Tuple{Real, Real}, 1}, b::Array{Tuple{Real, Real}, 1}, c::Array{Real, 1})
      r = Array{Tuple{Real, Real}}([])

      for i in 1:size(a)[1]
         push!(r, (c[i]*(a[i][1] - b[i][1]), c[i]*(a[i][2] - b[i][2])))
      end

      return r
   end

   centroids, areas = Diagram.centroidsAndAreas(V)

   return weightedDifference(points, centroids, areas)
end

function modf(V::Diagram.DCEL)
   function g(V)
   end

   return defaultf(V) + g(V)
end

function mod∇f(V::Diagram.DCEL, points::Array{Tuple{Real, Real}, 1})
   function ∇g(V, points)
   end

   return default∇f(V, points) + ∇g(V, points)
end

function f(V::Diagram.DCEL)
   return defaultf(V)
end

function ∇f(V::Diagram.DCEL, points::Array{Tuple{Real, Real}, 1})
   return default∇f(V, points)
end

function init(points::Array{Tuple{Real, Real}, 1}, w::Number, h::Number)
   global WIDTH, HEIGHT, pointsₖ

   WIDTH = w
   HEIGHT = h
   pointsₖ = deepcopy(points)

   V = Fortune.compute(points, w, h)
   Intersect.intersect(V, Intersect.Rectangle(WIDTH, HEIGHT))

   f = Optimization.f(V)
   ∇f = Optimization.∇f(V, points)
   ξ = Optimization.normInf(∇f)

   return V, f, ∇f, ξ
end

function lloydIteration(points::Array{Tuple{Real, Real}, 1}, V::Diagram.DCEL)
   points, areas = Diagram.centroidsAndAreas(V)

   V = Fortune.compute(points)
   Intersect.intersect(V, Intersect.Rectangle(WIDTH, HEIGHT))

   f = Optimization.f(V)
   ∇f = Optimization.∇f(V, points)
   ξ = Optimization.normInf(∇f)

   return V, f, ξ
end

function lineSearch(points::Array{Tuple{Real, Real}, 1}, d::Array{Tuple{Real, Real}, 1}, α::Number, λ₀::Number, μ::Number, ϵ::Number, V::Diagram.DCEL, f::Number, ∇f::Array{Tuple{Real, Real}, 1}, ξ)
   global WIDTH, HEIGHT, pointsₖ

   evaluations = 0

   n = size(points)[1]

   # backtracking
   λ = λ₀
   while true
      outOfBounds = false

      for i in 1:n
         x = points[i][1] - λ*d[i][1]
         y = points[i][2] - λ*d[i][2]

         if !(0 <= x <= WIDTH && 0 <= y <= HEIGHT)
            outOfBounds = true

            break
         end

         pointsₖ[i] = (x, y)
      end

      dᵗ∇f = 0
      if outOfBounds
         fₖ = Inf
         ∇fₖ = Inf
         ξₖ = Inf
         V = nothing
      else
         Vₖ = Fortune.compute(pointsₖ)
         Intersect.intersect(Vₖ, Intersect.Rectangle(WIDTH, HEIGHT))
         evaluations += 1

         fₖ = Optimization.f(Vₖ)
         ∇fₖ = Optimization.∇f(Vₖ, pointsₖ)
         ξₖ = Optimization.normInf(∇fₖ)

         for i in 1:n
            dᵗ∇f += d[i][1]*∇f[i][1]
            dᵗ∇f += d[i][2]*∇f[i][2]
         end
      end

      #println("fₖ = ", fₖ, "\t\tf = ", f, "\t\tλ = ", λ, "\t\tξₖ = ", ξₖ)
      # armijo's condition || failsafes
      if fₖ < f - α*λ*dᵗ∇f || ξₖ <= ϵ || λ == 0
         for i in 1:n
            points[i] = pointsₖ[i]
         end

         V = Vₖ
         f = fₖ
         ∇f = ∇fₖ
         ξ = ξₖ
         break
      else
         old = λ
         λ *= μ
         if old == λ
            ξ = ϵ
            break
         end
      end
   end

   return V, f, ∇f, ξ, evaluations
end

end # module

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
"""
SPQR_RANK_FORM_BASIS forms the basis for the null space of a matrix.

 Called from spqr_basic and spqr_cod after these routines call spqr and
 spqr_rank_ssi.  The input parameters are as used in spqr_basic and spqr_cod.
 The parameter call_from indicates the type of call:
      call_from = 1 -- a call from spqr_basic, form null space basis for A'
      call_from = 2 -- a call from spqr_cod, form null space basis for A
      call_from = 3 -- a call from spqr_cod, form null space basis for A'
 Output:
   N -- basis for the null space in implicit form
   stats -- an update of the stats vector in spqr_basic or spqr_cod
   stats_ssp -- information about the call to spqr_ssp
   est_sval_upper_bounds -- an update to the estimated singular value upper
      bounds
"""
function spqr_rank_form_basis(call_from, A::AbstractMatrix{T},
                              U::AbstractMatrix{T}, V::AbstractMatrix{T},
                              Q1::AbstractMatrix{T}, prow::AbstractVector{Int},
                              rank_spqr::Int, numerical_rank::Int,
                              stats::Statistics, opts,
                              est_sval_upper_bounds::AbstractVector{S},
                              nsvals_small::Int, nsvals_large::Int,
                              p1::AbstractVector{Int} = 1:0,
                              Q2::AbstractMatrix{T} = zeros(0,0),
                              p2::AbstractVector{Int} = 1:0
                             ) where {T<:Number, S<:Real}

# Copyright 2012, Leslie Foster and Timothy A Davis.

    opts, get_details, start_with_A_transpose, k =
    get_opts(opts, :get_details, :start_with_A_transpose, :k)

    if get_details == 1
        t = time_ns()
    end
    m , n = size(A)
    @assert rank_spqr == size(U, 1)
    @assert rank_spqr == size(V, 1)
    nullity_R11 = rank_spqr - numerical_rank

    #-------------------------------------------------------------------------------
    # form null space basis for A or A' using X and Q from QR factorization
    #-------------------------------------------------------------------------------

    N = if call_from == 1
        ImplicitProd(invperm(prow), Q1, U)
    elseif start_with_A_transpose == (call_from == 2)
        ImplicitProd(Q1, invperm(p2), V)
    else
        ImplicitProd(invperm(p1), Q2, U)
    end

    if get_details == 1
        stats.time_basis = time_ns() - t
    end

    #-------------------------------------------------------------------------------
    # call spqr_ssp to enhance, potentially, est_sval_upper_bounds, the estimated
    #    upper bounds on the singular values
    # and / or
    #    to estimate ||A * N || or || A' * N ||
    #-------------------------------------------------------------------------------

    # Note: nullity = m - numerical_rank;   # the number of columns in X and N

    kk = max(nsvals_small, opts.data.k)

    if call_from == 1    # call from spqr_basic (for null space basis for A')

        # Note: opts.k is not the same as k in the algorithm description above.

        # Note that, by the interleave theorem for singular values, for
        # i=1:nullity, singular value i of A'*N will be an upper bound on singular
        # value numerical_rank + i of A.  S[i] is an estimate of singular value i
        # of A'*N with an estimated accuracy of stats_ssp.est_error_bounds[i].
        # Therefore let
        s_ssp, stats_ssp = spqr_ssp(A', N; nargout=2, opts..., k = kk)

    elseif call_from == 2    # call from spqr_cod (for null space basis of A)

        # Note that, by the interleave theorem for singular
        # values, for i = 1, ..., nullity, singular value i of A*N will be
        # an upper bound on singular value numerical_rank + i of A.
        # S[i] is an estimate of singular value i of A*N with an estimated
        # accuracy of stats_ssp.est_error_bounds[i].  Therefore let
        s_ssp, stats_ssp = spqr_ssp(A, N; nargout=2, opts..., k = kk)
    end

    if call_from == 1 || call_from == 2

        # By the comments prior to the call to spqr_ssp we have
        # s_ssp + stats_ssp.est_error_bounds are estimated upper bounds
        # for singular values (numerical_rank+1):(numerical_rank+nsvals_small)
        # of A.  We have two estimates for upper bounds on these singular
        # values of A.  Choose the smaller of the two:
        if nsvals_small > 0
            est_sval_upper_bounds[nsvals_large+1:end] .=
                min.(est_sval_upper_bounds[nsvals_large+1:end],
                vec(s_ssp[1:nsvals_small]) + stats_ssp.est_error_bounds[1:nsvals_small])
        end

    elseif call_from == 3  # call from spqr_cod (for null space basis for A')

        # call ssp to estimate nsvals_small sing. values of A' * N
        #    (useful to estimate the norm of A' * N)
        _, stats_ssp = spqr_ssp(A', N; nargout=2, opts..., nsvals_small = nsvals_small)   ##ok

    end

    if (get_details == 1)
        if call_from == 1 || call_from == 3
            stats.stats_ssp_NT = stats_ssp
        elseif call_from == 2
            stats.stats_ssp_N = stats_ssp
        end
        stats.time_svd = stats.time_svd + stats_ssp.time_svd
    end

    N , stats, stats_ssp, est_sval_upper_bounds
end


text:
# function cutrod_naive(p::Array{Int, 1}, n::Int)::Int
# 	if n == 0
# 		return 0
# 	end

# 	q = typemin(Int)    # q = -∞
# 	for i ∈ 1:n
# 		q = max(q, p[i] + cutrod(p, n - i))
# 	end

# 	return q
# end

# function cutrod_memoized(p::Array{Int, 1}, n::Int)::Int
# 	r = fill(typemin(Int), length(p))

# 	function cutrod_memoized_aux(p::Array{Int, 1}, n::Int, r::Array{Int, 1})::Int
# 		if n == 0
# 			return 0
# 		end

# 		if r[n] ≥ 0
# 			return r[n]
# 		end

# 		q = typemin(Int)    # q = -∞
# 		for i ∈ 1:n
# 			q = max(q, p[i] + cutrod_memoized_aux(p, n - i, r))
# 		end

# 		return q
# 	end

# 	return cutrod_memoized_aux(p, n, r)
# end

# prices = [1, 5, 8, 9, 10, 17, 17, 20, 24, 30]
# for j ∈ 1:10
# 	println("Best reveneu for length $j is $(cutrod_memoized(prices, j)) \$.")
# end

max_stars_repo_name: SudipSinha/edu

max_issues_repo_name: SudipSinha/edu

max_forks_repo_name: SudipSinha/edu


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using  Base.MultiplicativeInverses: SignedMultiplicativeInverse

struct ReshapedArray{T,N,P<:AbstractArray,MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int}}}} <: AbstractArray{T,N}
    parent::P
    dims::NTuple{N,Int}
    mi::MI
end
ReshapedArray(parent::AbstractArray{T}, dims::NTuple{N,Int}, mi) where {T,N} = ReshapedArray{T,N,typeof(parent),typeof(mi)}(parent, dims, mi)

# IndexLinear ReshapedArray
const ReshapedArrayLF{T,N,P<:AbstractArray} = ReshapedArray{T,N,P,Tuple{}}

# Fast iteration on ReshapedArrays: use the parent iterator
struct ReshapedArrayIterator{I,M}
    iter::I
    mi::NTuple{M,SignedMultiplicativeInverse{Int}}
end
ReshapedArrayIterator(A::ReshapedArray) = _rs_iterator(parent(A), A.mi)
function _rs_iterator(P, mi::NTuple{M}) where M
    iter = eachindex(P)
    ReshapedArrayIterator{typeof(iter),M}(iter, mi)
end

struct ReshapedIndex{T}
    parentindex::T
end

# eachindex(A::ReshapedArray) = ReshapedArrayIterator(A)  # TODO: uncomment this line
start(R::ReshapedArrayIterator) = start(R.iter)
@inline done(R::ReshapedArrayIterator, i) = done(R.iter, i)
@inline function next(R::ReshapedArrayIterator, i)
    item, inext = next(R.iter, i)
    ReshapedIndex(item), inext
end
length(R::ReshapedArrayIterator) = length(R.iter)

"""
    reshape(A, dims...) -> R
    reshape(A, dims) -> R

Return an array `R` with the same data as `A`, but with different
dimension sizes or number of dimensions. The two arrays share the same
underlying data, so that setting elements of `R` alters the values of
`A` and vice versa.

The new dimensions may be specified either as a list of arguments or
as a shape tuple. At most one dimension may be specified with a `:`,
in which case its length is computed such that its product with all
the specified dimensions is equal to the length of the original array
`A`. The total number of elements must not change.

```jldoctest
julia> A = collect(1:16)
16-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16

julia> reshape(A, (4, 4))
4×4 Array{Int64,2}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> reshape(A, 2, :)
2×8 Array{Int64,2}:
 1  3  5  7   9  11  13  15
 2  4  6  8  10  12  14  16
```

"""
reshape

reshape(parent::AbstractArray, dims::IntOrInd...) = reshape(parent, dims)
reshape(parent::AbstractArray, shp::NeedsShaping) = reshape(parent, to_shape(shp))
reshape(parent::AbstractArray, dims::Dims)        = _reshape(parent, dims)

# Allow missing dimensions with Colon():
reshape(parent::AbstractArray, dims::Int...) = reshape(parent, dims)
reshape(parent::AbstractArray, dims::Union{Int,Colon}...) = reshape(parent, dims)
reshape(parent::AbstractArray, dims::Tuple{Vararg{Union{Int,Colon}}}) = _reshape(parent, _reshape_uncolon(parent, dims))
@inline function _reshape_uncolon(A, dims)
    @noinline throw1(dims) = throw(DimensionMismatch(string("new dimensions $(dims) ",
        "may have at most one omitted dimension specified by `Colon()`")))
    @noinline throw2(A, dims) = throw(DimensionMismatch(string("array size $(length(A)) ",
        "must be divisible by the product of the new dimensions $dims")))
    pre = _before_colon(dims...)
    post = _after_colon(dims...)
    any(d -> d isa Colon, post) && throw1(dims)
    sz, remainder = divrem(length(A), prod(pre)*prod(post))
    remainder == 0 || throw2(A, dims)
    (pre..., sz, post...)
end
@inline _before_colon(dim::Any, tail...) =  (dim, _before_colon(tail...)...)
@inline _before_colon(dim::Colon, tail...) = ()
@inline _after_colon(dim::Any, tail...) =  _after_colon(tail...)
@inline _after_colon(dim::Colon, tail...) = tail

reshape(parent::AbstractArray{T,N}, ndims::Val{N}) where {T,N} = parent
function reshape(parent::AbstractArray, ndims::Val{N}) where N
    reshape(parent, rdims(Val(N), indices(parent)))
end

# Move elements from inds to out until out reaches the desired
# dimensionality N, either filling with OneTo(1) or collapsing the
# product of trailing dims into the last element
rdims_trailing(l, inds...) = length(l) * rdims_trailing(inds...)
rdims_trailing(l) = length(l)
rdims(out::Val{N}, inds::Tuple) where {N} = rdims(ntuple(i -> OneTo(1), Val(N)), inds)
rdims(out::Tuple{}, inds::Tuple{}) = () # N == 0, M == 0
@noinline rdims(out::Tuple{}, inds::Tuple{Any}) = throw(ArgumentError("new dimensions cannot be empty")) # N == 0
@noinline rdims(out::Tuple{}, inds::NTuple{M,Any}) where {M} = throw(ArgumentError("new dimensions cannot be empty")) # N == 0
rdims(out::Tuple{Any}, inds::Tuple{}) = out # N == 1, M == 0
rdims(out::NTuple{N,Any}, inds::Tuple{}) where {N} = out # N > 1, M == 0
rdims(out::Tuple{Any}, inds::Tuple{Any}) = inds # N == 1, M == 1
rdims(out::Tuple{Any}, inds::NTuple{M,Any}) where {M} = (OneTo(rdims_trailing(inds...)),) # N == 1, M > 1
rdims(out::NTuple{N,Any}, inds::NTuple{N,Any}) where {N} = inds # N > 1, M == N
rdims(out::NTuple{N,Any}, inds::NTuple{M,Any}) where {N,M} = (first(inds), rdims(tail(out), tail(inds))...) # N > 1, M > 1, M != N

# _reshape on Array returns an Array
_reshape(parent::Vector, dims::Dims{1}) = parent
_reshape(parent::Array, dims::Dims{1}) = reshape(parent, dims)
_reshape(parent::Array, dims::Dims) = reshape(parent, dims)

# When reshaping Vector->Vector, don't wrap with a ReshapedArray
function _reshape(v::AbstractVector, dims::Dims{1})
    len = dims[1]
    len == length(v) || _throw_dmrs(n, "length", len)
    v
end
# General reshape
function _reshape(parent::AbstractArray, dims::Dims)
    n = _length(parent)
    prod(dims) == n || _throw_dmrs(n, "size", dims)
    __reshape((parent, IndexStyle(parent)), dims)
end

@noinline function _throw_dmrs(n, str, dims)
    throw(DimensionMismatch("parent has $n elements, which is incompatible with $str $dims"))
end

# Reshaping a ReshapedArray
_reshape(v::ReshapedArray{<:Any,1}, dims::Dims{1}) = _reshape(v.parent, dims)
_reshape(R::ReshapedArray, dims::Dims) = _reshape(R.parent, dims)

function __reshape(p::Tuple{AbstractArray,IndexCartesian}, dims::Dims)
    parent = p[1]
    strds = front(size_to_strides(size(parent)..., 1))
    strds1 = map(s->max(1,s), strds)  # for resizing empty arrays
    mi = map(SignedMultiplicativeInverse, strds1)
    ReshapedArray(parent, dims, reverse(mi))
end

function __reshape(p::Tuple{AbstractArray,IndexLinear}, dims::Dims)
    parent = p[1]
    ReshapedArray(parent, dims, ())
end

size(A::ReshapedArray) = A.dims
similar(A::ReshapedArray, eltype::Type, dims::Dims) = similar(parent(A), eltype, dims)
IndexStyle(::Type{<:ReshapedArrayLF}) = IndexLinear()
parent(A::ReshapedArray) = A.parent
parentindexes(A::ReshapedArray) = map(s->1:s, size(parent(A)))
reinterpret(::Type{T}, A::ReshapedArray, dims::Dims) where {T} = reinterpret(T, parent(A), dims)

@inline ind2sub_rs(::Tuple{}, i::Int) = i
@inline ind2sub_rs(strds, i) = _ind2sub_rs(strds, i - 1)
@inline _ind2sub_rs(::Tuple{}, ind) = (ind + 1,)
@inline function _ind2sub_rs(strds, ind)
    d, r = divrem(ind, strds[1])
    (_ind2sub_rs(tail(strds), r)..., d + 1)
end

@inline function getindex(A::ReshapedArrayLF, index::Int)
    @boundscheck checkbounds(A, index)
    @inbounds ret = parent(A)[index]
    ret
end
@inline function getindex(A::ReshapedArray{T,N}, indexes::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, indexes...)
    _unsafe_getindex(A, indexes...)
end
@inline function getindex(A::ReshapedArray, index::ReshapedIndex)
    @boundscheck checkbounds(parent(A), index.parentindex)
    @inbounds ret = parent(A)[index.parentindex]
    ret
end

@inline function _unsafe_getindex(A::ReshapedArray{T,N}, indexes::Vararg{Int,N}) where {T,N}
    i = sub2ind(size(A), indexes...)
    I = ind2sub_rs(A.mi, i)
    _unsafe_getindex_rs(parent(A), I)
end
_unsafe_getindex_rs(A, i::Integer) = (@inbounds ret = A[i]; ret)
@inline _unsafe_getindex_rs(A, I) = (@inbounds ret = A[I...]; ret)

@inline function setindex!(A::ReshapedArrayLF, val, index::Int)
    @boundscheck checkbounds(A, index)
    @inbounds parent(A)[index] = val
    val
end
@inline function setindex!(A::ReshapedArray{T,N}, val, indexes::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, indexes...)
    _unsafe_setindex!(A, val, indexes...)
end
@inline function setindex!(A::ReshapedArray, val, index::ReshapedIndex)
    @boundscheck checkbounds(parent(A), index.parentindex)
    @inbounds parent(A)[index.parentindex] = val
    val
end

@inline function _unsafe_setindex!(A::ReshapedArray{T,N}, val, indexes::Vararg{Int,N}) where {T,N}
    @inbounds parent(A)[ind2sub_rs(A.mi, sub2ind(size(A), indexes...))...] = val
    val
end

# helpful error message for a common failure case
const ReshapedRange{T,N,A<:Range} = ReshapedArray{T,N,A,Tuple{}}
setindex!(A::ReshapedRange, val, index::Int) = _rs_setindex!_err()
setindex!(A::ReshapedRange{T,N}, val, indexes::Vararg{Int,N}) where {T,N} = _rs_setindex!_err()
setindex!(A::ReshapedRange, val, index::ReshapedIndex) = _rs_setindex!_err()

@noinline _rs_setindex!_err() = error("indexed assignment fails for a reshaped range; consider calling collect")

unsafe_convert(::Type{Ptr{T}}, a::ReshapedArray{T}) where {T} = unsafe_convert(Ptr{T}, parent(a))

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
@testset "JuliaDemoCard" begin
    # default behavior
    simplest = democard("simplest.jl")
    @test simplest.cover === nothing
    @test simplest.id == "Simplest"
    @test simplest.path == "simplest.jl"
    @test simplest.title == "Simplest"
    @test simplest.description == "This is the content"

    @testset "parse" begin
        @testset "title, id and description" begin
            # JuliaDemoCard doesn't parse title from the markdown contents
            title_1 = JuliaDemoCard("title_1.jl")
            @test title_1.id == "Custom-Title"
            @test title_1.title == "Custom Title"
            @test title_1.description == "This is the content"
            @test title_1.hidden == false
            @test title_1.author == ""
            @test title_1.date == DateTime(0)

            title_2 = JuliaDemoCard("title_2.jl")
            @test title_2.id == "custom_id_2"
            @test title_2.title == "Custom Title 2"
            @test title_2.description == "This is the content"

            title_3 = JuliaDemoCard("title_3.jl")
            @test title_3.id == "Custom-Title-3-1"
            @test title_3.title == "Custom Title 3-1"
            @test title_3.description == "This is the content"

            title_4 = JuliaDemoCard("title_4.jl")
            @test title_4.id == "Custom-Title"
            @test title_4.title == "Custom Title"
            @test title_4.description == "This is the content"

            title_5 = JuliaDemoCard("title_5.jl")
            @test title_5.id == "Custom-Title"
            @test title_5.title == "Custom Title"
            @test title_5.description == "Custom Description"

            title_6 = JuliaDemoCard("title_6.jl")
            @test title_6.id == "custom_id"
            @test title_6.title == "Custom Title"
            @test title_6.description == "This is the content"

            title_7 = JuliaDemoCard("title_7.jl")
            @test title_7.id == "custom_id"
            @test title_7.title == "Custom Title"
            @test title_7.description == "This is the content"
            @test title_7.hidden == false
            @test title_7.author == "Jane Doe; John Roe"
            @test title_7.date == DateTime("2020-01-31")
            @test title_7.cover == "https://juliaimages.org/latest/assets/logo.png"

            title_8 = JuliaDemoCard("title_8.jl")
            @test title_8.id == "custom_id"
            @test title_8.title == "Custom Title"
            @test title_8.description == "This is the content"

            description_1 = JuliaDemoCard("description_1.jl")
            @test description_1.id == "Custom-Title"
            @test description_1.title == "Custom Title"
            @test description_1.description == "this is a single line destiption that spans over multiple lines\n"

            description_2 = JuliaDemoCard("description_2.jl")
            @test description_2.id == "Custom-Title"
            @test description_2.title == "Custom Title"
            @test description_2.description == "this is a multi line\ndestiption that spans\nover multiple lines\n"
        end

        @testset "cover" begin
            cover_1 = JuliaDemoCard("cover_1.jl")
            @test cover_1.cover == nothing

            cover_2 = JuliaDemoCard("cover_2.jl")
            @test cover_2.cover ==  joinpath("..", "logo.png")

            cover_3 = JuliaDemoCard("cover_3.jl")
            @test cover_3.cover == "nonexistence.jpg"

            cover_4 = JuliaDemoCard("cover_4.jl")
            @test cover_4.cover == joinpath("..", "logo.png")
        end

        @testset "version" begin
            card = JuliaDemoCard("version_1.jl")
            card.julia == v"1.2.3"

            card = JuliaDemoCard("version_2.jl")
            card.julia == v"999.0.0"
            warn_msg = @capture_err preview_demos("version_2.jl"; require_html=false)
            @test occursin("The running Julia version `$(VERSION)` is older than the declared compatible version `$(card.julia)`.", warn_msg)
        end

        @testset "notebook" begin
            card = JuliaDemoCard("notebook_1.jl")
            @test card.notebook == true

            card = JuliaDemoCard("notebook_2.jl")
            @test card.notebook == false

            card = JuliaDemoCard("notebook_3.jl")
            @test card.notebook == false

            card = @suppress_err JuliaDemoCard("notebook_4.jl")
            @test card.notebook == nothing
            warn_msg = @capture_err JuliaDemoCard("notebook_4.jl")
            @test occursin("`notebook` option should be either `\"true\"` or `\"false\"`, instead it is: nothing. Fallback to unconfigured.", warn_msg)
        end

        @testset "generate" begin
            page_dir = @suppress_err preview_demos("title_7.jl", theme="grid", require_html=false)
            card_dir = joinpath(page_dir, "julia")

            function strip_notebook(contents)
                # the notebook contents are different in CI and non-CI cases
                # [![notebook](badge_url)](notebook_url)
                replace(contents, r"\[!\[notebook\]\(.*\)\]\(.*\.ipynb\) " => "")
            end
            isempty(DemoCards.get_nbviewer_root_url("master")) || println(read(joinpath(card_dir, "title_7.md"), String)) # for debug usage
            @test_reference joinpath(test_root, "references", "cards", "julia_md.md") strip_notebook(read(joinpath(card_dir, "title_7.md"), String)) by=ignore_CR
            @test_reference joinpath(test_root, "references", "cards", "julia_src.jl") read(joinpath(card_dir, "title_7.jl"), String) by=ignore_CR
            @test isfile(joinpath(card_dir, "title_7.ipynb"))

            # check if notebook keyword actually affects the generation
            page_dir = @suppress_err preview_demos("notebook_2.jl", theme="grid", require_html=false)
            card_dir = joinpath(page_dir, "julia")
            contents = read(joinpath(card_dir, "notebook_2.md"), String)
            @test !occursin("[notebook]", contents)
        end
    end
end


text:
using PyCall
using Test

function testSerial(port)
	serial = pyimport("serial")

	try
		uc = serial.Serial(port, 115200, timeout=0.2)
		return true
	catch e
		return false
	end
end

function testLEDs(port, leds)
	serial = pyimport("serial")
	uc::PyObject = undef

	try
		uc = serial.Serial(port, 115200, timeout=0.2)
	catch e
		return false
	end

	for i in 1:(leds + 1)
		uc.write([UInt8(i - 1), 0x00, 0x00, 0x00])
	end

	uc.write([0x00, 0xFF, 0x00, 0x00, 0xFF])

	for i in 1:(leds + 1)
		if (i != 1)
			uc.write([UInt8(i - 2), 0x00, 0x00, 0x00])
		end

		uc.write([UInt8(i - 1), 0xFF, 0x00, 0x00, 0xFF])
		sleep(0.1)
	end

	return true
end

@test testSerial("/dev/ttyACM0")
@test testLEDs("/dev/ttyACM0", 68)


text:
function dnsol(n, L::TronSparseMatrixCSC, r)
    # Solve L*x = r and store the result in r.

    @inbounds for j=1:n
        temp = r[j]/L.diag_vals[j]
        for k=L.colptr[j]:L.colptr[j+1]-1
            r[L.rowval[k]] = r[L.rowval[k]] - L.tril_vals[k]*temp
        end
        r[j] = temp
    end

    return
end

function dnsol(n, L::TronDenseMatrix, r)
    # Solve L*x = r and store the result in r.

    @inbounds for j=1:n
        temp = r[j]/L.vals[j,j]
        for k=j+1:n
            r[k] = r[k] - L.vals[k,j]*temp
        end
        r[j] = temp
    end
end


text:
using BenchmarkTools
using Gadfly
using StatsBase
BLAS.set_num_threads(1)

using MLKernels

isequal(GaussianKernel(3.0), GaussianKernel(3.0))
GaussianKernel(3.0) == GaussianKernel(3.0)

l = GaussianKernel(3.0)
k = GaussianKernel(3.0)

l == k                         # false
l.alpha == k.alpha             # false
l.alpha.value == k.alpha.value # true

isimmutable(GaussianKernel(3.0)) # true
isimmutable(l.alpha)             # false
isimmutable(l.alpha.bounds)      # true
isimmutable(l.alpha.value)       # true

N = 5000
x = rand(1, N) * 4π - 2π
yy = sinc.(x) # vec(sinc.(4 .* x) .+ 0.2 .* sin.(30 .* x))
y = squeeze(yy + 0.1randn(1, N), 1)
xnew = collect(-2.5π:0.01:2.5π)'

reload("KernelRidgeRegression")

using KernelRidgeRegression
using StatsBase
using MLKernels

@time mykrr = fit(KRR, x, y, 1e-3/5000, GaussianKernel(1.0))
@time ynew = predict(mykrr, xnew);

show(mykrr)

@time mynystkrr = fit(KernelRidgeRegression.NystromKRR, x, y, 1e-3/5000, 280, MLKernels.GaussianKernel(100.0))
@time ynystnew = predict(mynystkrr, xnew)

@time myfastkrr = fit(KernelRidgeRegression.FastKRR, x, y, 4/5000, 11, MLKernels.GaussianKernel(100.0))
@time yfastnew = predict(myfastkrr, xnew)

@time mytnkrr = fit(KernelRidgeRegression.TruncatedNewtonKRR, x, y, 4/5000, MLKernels.GaussianKernel(100.0), 0.5, 200)
@time ytnnew = predict(mytnkrr, xnew)

@time myrandkrr = fit(KernelRidgeRegression.RandomFourierFeatures, x, y, 1e-3/5000, 200 , 1.0)
@time yrandnew = predict(myrandkrr, xnew)

# tanh makes the whole thing rotation-symetric and go through the origin
@time myrandkrr2 = fit(KernelRidgeRegression.RandomKRR, x, y,;
                                             4/5000, 2000, 1.0, (X, W) -> tanh(X' * W));
@time yrandnew2 = predict(myrandkrr2, xnew);

@time myrandkrr3 = fit(KernelRidgeRegression.RandomFourierFeatures, x, y,
                                             4/5000, 2000, 1.0, (X, W) -> 1 ./ ((X' * W) .^ 2 + 1))
@time yrandnew3 = predict(myrandkrr3 , xnew);
@time myrandkrr4 = fit(KernelRidgeRegression.RandomFourierFeatures, x, y,
                                             4/5000, 2000, 1.0, (X, W) -> ((X' * W) ./ maximum(X'*W)))
@time yrandnew4 = predict(myrandkrr4 , xnew);
@time myrandkrr5 = fit(KernelRidgeRegression.RandomFourierFeatures, x, y,
                                             4/5000, 2000, 1.0, (X, W) -> begin
                                             z = X' * W
                                             f(x) = x > 0 ? x / (x + 1) : x / (x - 1)
                                             for i in eachindex(z)
                                                 z[i] = f(z[i])
                                             end
                                             z
                                             end)
@time yrandnew5 = predict(myrandkrr5 , xnew);

extrema(ynew - yfastnew)
extrema(ynew - yrandnew)
extrema(ynew - ytnnew)
extrema(ynew - ynystnew)

sqrt(mean((ynew - ynystnew) .^ 2))
sqrt(mean((ynew - yfastnew) .^ 2))
sqrt(mean((ynew - yrandnew) .^ 2))

plot(
    # layer(x = xnew, y = yfastnew,  Geom.line, Theme(default_color = colorant"yellow")),
    # layer(x = xnew, y = ynystnew,  Geom.line, Theme(default_color = colorant"red")),
    # layer(x = xnew, y = ytnnew,    Geom.line, Theme(default_color = colorant"purple")),
    layer(x = xnew, y = yrandnew,  Geom.line, Theme(default_color = colorant"purple")),
    # layer(x = xnew, y = yrandnew2,  Geom.line, Theme(default_color = colorant"yellow")),
    # layer(x = xnew, y = ynew,      Geom.line, Theme(default_color = colorant"green")),
    layer(x = x,    y = y,        Geom.line, Theme(default_color = colorant"blue")),
    Coord.cartesian(ymin = -1.5, ymax = 1.5),
    Guide.manual_color_key(
        "", ["Fast", "Nystrom", "RFF", "Data"],
        [colorant"red", colorant"purple", colorant"green", colorant"blue"]
    )
)

using MLDataUtils
using LossFunctions
using Combinatorics
using DataStructures
reload("KernelRidgeRegression")

N = 500
x = rand(1, N) * 4π - 2π
yy = sinc.(x) # vec(sinc.(4 .* x) .+ 0.2 .* sin.(30 .* x))
y = squeeze(yy + 0.1randn(1, N), 1)
xnew = collect(-2.5π:0.01:2.5π)'

reload("KernelRidgeRegression")
KernelRidgeRegression.crossvalidate_parameters(
    KernelRidgeRegression.KRR,
    x, y, 6,
    [10.0^x for x in -5:3],
    [MLKernels.GaussianKernel(10.0^x) for x in -6:3]
)

tmp = Iterators.product([10.0^x for x in -3:3], [10.0^x for x in -3:3])
collect(tmp)
for (k, v) in tmp
    println(k)
    println(v)
end

kk = collect(keys(tmp))
for i in Iterators.product(values(tmp)...)
    @show OrderedDict(kk[j] => i[j] for j = 1:length(i))
end

lossₘᵢₙ
λₘᵢₙ
σₘᵢₙ
ynew = predict(mₘᵢₙ, xnew)

plot(
    # layer(x = xnew, y = ytnnew,    Geom.line, Theme(default_color = colorant"purple")),
    # layer(x = xnew, y = yrandnew,  Geom.line, Theme(default_color = colorant"red")),
    # layer(x = xnew, y = yrandnew2,  Geom.line, Theme(default_color = colorant"yellow")),
    # layer(x = xnew, y = yfastnew,  Geom.line, Theme(default_color = colorant"yellow")),
    layer(x = xnew, y = ynew,      Geom.line, Theme(default_color = colorant"green")),
    # layer(x = x,    y = y,        Geom.line, Theme(default_color = colorant"blue")),
    Coord.cartesian(ymin = -1.5, ymax = 1.5),
    Guide.manual_color_key("", ["RFF", "KRR", "Data"], [colorant"red", colorant"green", colorant"blue"])
)


text:
"""
    AbstractNearestNeighbors

Abstract type for nearest neighbor plug-in implementations.
"""
abstract type AbstractNearestNeighbors end

"""
    knn(NN::AbstractNearestNeighbors, X::AbstractVecOrMat{T}; kwargs...) -> (D,E)

Perform construction of the distance matrix `D` and neares neighbor weighted graph `E` from the `NN` object
"""
function knn(NN::AbstractNearestNeighbors, X::AbstractVecOrMat{T}; kwargs...) where T<:Real end

"""
    BruteForce

Calculate NN using pairwise distance matrix.
"""
struct BruteForce{T<:Real} <: AbstractNearestNeighbors
    k::Integer
    fitted::AbstractMatrix{T}
end
fit(::Type{BruteForce}, X::AbstractMatrix{T}, k::Integer) where {T<:Real} = BruteForce(k, X)
show(io::IO, NN::BruteForce) = print(io, "BruteForce(k=$(NN.k))")

function knn(NN::BruteForce{T}, X::AbstractVecOrMat{T}; self=false) where T<:Real
    m, n = size(X)
    k = NN.k
    @assert n > k "Number of observations must be more then $(k)"

    D = pairwise((x,y)->norm(x-y), eachcol(NN.fitted), eachcol(X))

    d = Array{T}(undef, k, n)
    e = Array{Int}(undef, k, n)
    idxs = self ? (1:k) : (2:k+1)

    @inbounds for j = 1 : n
        e[:, j] = sortperm(D[:,j])[idxs]
        d[:, j] = D[e[:, j],j]
    end

    return (d, e)
end


text:
# conversion methods for DataArrays and DataFrames

function rcopy{T,S<:VectorSxp}(::Type{DataArray{T}}, s::Ptr{S})
    DataArray(rcopy(Array{T},s), isna(s))
end
function rcopy{S<:VectorSxp}(::Type{DataArray}, s::Ptr{S})
    DataArray(rcopy(Array,s), isna(s))
end

function rcopy(::Type{DataArray}, s::Ptr{IntSxp})
    isFactor(s) && error("$s is a R factor")
    DataArray(rcopy(Array,s), isna(s))
end
function rcopy(::Type{PooledDataArray}, s::Ptr{IntSxp})
    isFactor(s) || error("$s is not a R factor")
    refs = DataArrays.RefArray([isna(x) ? zero(Int32) : x for x in s])
    compact(PooledDataArray(refs,rcopy(getattrib(s,Const.LevelsSymbol))))
end

function rcopy(::Type{DataFrame}, s::Ptr{VecSxp})
    isFrame(s) || error("s is not a R data frame")
    DataFrame(Any[isFactor(c)? rcopy(PooledDataArray, c) : rcopy(DataArray, c) for c in s],
              rcopy(Array{Symbol},getnames(s)))
end

## DataArray to sexp conversion.
function sexp(v::DataArray)
    rv = protect(sexp(v.data))
    try
        for (i,isna) = enumerate(v.na)
            if isna
                rv[i] = naeltype(eltype(rv))
            end
        end
    finally
        unprotect(1)
    end
    rv
end

## PooledDataArray to sexp conversion.
function sexp{T<:Compat.String,R<:Integer}(v::PooledDataArray{T,R})
    rv = sexp(v.refs)
    setattrib!(rv, Const.LevelsSymbol, sexp(v.pool))
    setattrib!(rv, Const.ClassSymbol, sexp("factor"))
    rv
end

## DataFrame to sexp conversion.
function sexp(d::DataFrame)
    nr,nc = size(d)
    rd = protect(allocArray(VecSxp, nc))
    try
        for i in 1:nc
            rd[i] = sexp(d[d.colindex.names[i]])
        end
        setattrib!(rd,Const.NamesSymbol, sexp([string(n) for n in d.colindex.names]))
        setattrib!(rd,Const.ClassSymbol, sexp("data.frame"))
        setattrib!(rd,Const.RowNamesSymbol, sexp(1:nr))
    finally
        unprotect(1)
    end
    rd
end

# R formula objects
function sexp(f::Formula)
    s = protect(rlang_p(:~,rlang_formula(f.lhs),rlang_formula(f.rhs)))
    try
        setattrib!(s,Const.ClassSymbol,sexp("formula"))
        setattrib!(s,".Environment",Const.GlobalEnv)
    finally
        unprotect(1)
    end
    s
end

function rlang_formula(e::Expr)
    e.head == :call || error("invalid formula object")
    op = e.args[1]
    if op == :&
        op = :(:)
    end
    if length(e.args) > 3 && op in (:+,:*,:(:))
        rlang_p(op,
                rlang_formula(Expr(e.head,e.args[1:end-1]...)),
                rlang_formula(e.args[end]))
    else
        rlang_p(op,map(rlang_formula,e.args[2:end])...)
    end
end
rlang_formula(e::Symbol) = e
rlang_formula(n::Number) = n


text:
using ReachabilityAnalysis
#=low = [0.2; -0.1; zeros(n1)]
high = [0.3; 0.1; zeros(n1)]
if one_loop_iteration
    low[end] = 1.0
    high[end] = 1.0
end
X0 = Hyperrectangle(low=low, high=high)=#
Dict([
    "opts" => Dict(:T=>99.0, :alg=>BOX(δ=0.01)),
    #"X0"   => [(3, X0)],
    "info" => Dict(
        "name" => "filtered_oscillator",
        "dim" => 3,
        "linear" => true,
        "hybrid" => true
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using FactCheck
using MarketData

include("financialtimeseries.jl")

#include("instruments.jl")

exitstatus()


text:
export audible
import Base: eltype, isinf

abstract type AbstractStream{R,T} end
const Audible{R,T} = Union{Sound{R,T},AbstractStream{R,T}}
stream(x::AbstractStream) = x

samplerate{R}(s::AbstractStream{R}) = R*Hz
eltype{R,T}(x::AbstractStream{R,T}) = T
isinf(x::AbstractStream) = isinf(duration(x))
done(x::AbstractStream) = duration(x) <= 0s

struct Stream{R,T} <: AbstractStream{R,T}
  fn::Function
  index::Array{Int,0}
end
stream(x::Sound) = Stream(x)
Stream{T}(R::Int,::Type{T},fn::Function) = Stream{R,T}(fn,Array{Int}() .= 1)

duration(x::Stream) = Inf*s
nsamples(x::Stream) = typemax(Int)
index(x::Stream) = x.index[]
index_ids(x::Stream) = [object_id(x.index)]

Stream{R,T}(x::Sound{R,T,1}) = limit(Stream(R,T,i -> x[i]),nsamples(x))
function Stream{R,T}(x::Sound{R,T,2})
  if nchannels(x) == 2
    leftright(Stream(x[:,left]),Stream(x[:,right]))
  else
    Stream(x[:])
  end
end

Stream(x::Array) = Stream(sound(x))

"""
    sound(stream,[len])

Consume some amount of the stream, converting it to a finite `sound`.

If left unspecified, the entire stream is consumed.  Infinite streams throw an
error.
"""
sound{R}(fs::AbstractStream{R},len::Quantity) = sound(fs,insamples(len,R*Hz))
function sound(fs::AbstractStream)
  if isinf(fs)
    error("Cannot turn an infinite stream into a sound. Specify a length as a ",
          "second argument.")
  else
    sound(fs,nsamples(fs))
  end
end

function sound{R,T}(fs::Stream{R,T},len::Int)
  i = fs.index[]
  fs.index[] .+= len
  Sound{R,T,1}(fs.fn(i:(i+len-1)))
end

"""
    audible(fn,len=Inf,asseconds=true;[sample_rate=samplerate(),eltype=Float64])

Creates monaural sound where `fn(t)` returns the amplitudes for a given `Range`
of time points, value ranging between -1 and 1.

If `asseconds` is false, `audible` creates a monaural sound where `fn(i)`
returns the amplitudes for a given `Range` of sample indices.

The function `fn` should always return elements of type `eltype`.

If an infinite length is specified, a stream is created rather than a sound.

The function `fn` need not be pure and it can be safely assumed that `fn` will
only be called for a given range of indices once. While indices and times passed
to `fn` normally begin from 1 and 0, respectively, this is not always the case.
"""
function audible(fn::Function,len=Inf,asseconds=true;
                 sample_rate=samplerate(),eltype=Float64)
  sample_rate_Hz = inHz(sample_rate)
  if ustrip(len) < Inf
    n = ustrip(insamples(len,sample_rate_Hz))
    R = ustrip(sample_rate_Hz)
    Sound{R,eltype,1}(!asseconds ? fn(1:n) : fn(((1:n)-1)/R))
  else
    R = ustrip(sample_rate_Hz)
    Stream(R,eltype,!asseconds ? fn : i -> fn((i-1)/R))
  end
end

mutable struct LimitStream{R,T} <: AbstractStream{R,T}
  data::AbstractStream{R,T}
  limit::Int
  limit_start::Int
end
duration{R}(x::LimitStream{R}) = inseconds(nsamples(x)*samples,R)
nsamples(x::LimitStream) =
  clamp(min(x.limit,(x.limit_start + x.limit) - index(x.data)),
        0,nsamples(x.data))
index(x::LimitStream) = index(x.data)
index_ids(x::LimitStream) = index_ids(x.data)
function sound(ls::LimitStream,len::Int = typemax(Int))
  ls.limit_start = index(ls)
  sound(ls.data,min(len,nsamples(ls)))
end
right(x::LimitStream) = limit(right(x.data),n)
left(x::LimitStream) = limit(left(x.data),n)

limit{R}(stream::AbstractStream{R},len::Time) = limit(stream,insamples(len,R*Hz))
limit{R,T}(s::AbstractStream{R,T},n::Int) = LimitStream{R,T}(s,n,1)

mutable struct CatStream{R,T} <: AbstractStream{R,T}
  a::AbstractStream{R,T}
  b::AbstractStream{R,T}
end

vcat(xs::Audible...) = error("Sample rates differ, fix with `resample`")
vcat{R}(xs::Audible{R}...) = reduce(vcat_helper,map(stream,xs))
function vcat_helper{R,T}(a::AbstractStream{R,T},b::AbstractStream{R,T})
  if isinf(a)
    error("Connot concatenate to the end of an infinite stream")
  else
    CatStream{R,T}(a,b)
  end
end
function vcat_helper{R,T,S}(a::AbstractStream{R,T},b::AbstractStream{R,S})
  error("Cannot concatenate different element types, convert to common ",
        "element type.")
  # doesn't quite work...
  # P = promote_type(T,S)
  # helper(s) = map(e -> convert(P,e),s)
  # vcat(helper(a),helper(b))
end

duration(x::CatStream) = duration(x.a) + duration(x.b)
nsamples(x::CatStream) =
  isinf(x.a) || isinf(x.b) ? typemax(Int) : nsamples(x.a) + nsamples(x.b)
index(x::CatStream) = index(x.a)
index_ids(x::CatStream) = [index_ids(x.a); index_ids(x.b)]
left(x::CatStream) = CatStream(left(x.a),left(x.b))
right(x::CatStream) = CatStream(right(x.a),right(x.b))

struct EmptyStream{R,T} <: AbstractStream{R,T}
end

duration(x::EmptyStream) = 0s
nsamples(x::EmptyStream) = 0
sound{R,T}(x::EmptyStream{R,T},len::Int) = Sound{R,T,1}(T[])
index(x::EmptyStream) = 0
index_ids(x::EmptyStream) = UInt64[]
left(x::EmptyStream) = x
right(x::EmptyStream) = y

function step_next{R,T}(x::CatStream{R,T},b::CatStream{R,T})
  x.a = b
  x.b = b.b
end

function step_next{R,T}(x::CatStream{R,T},b::AbstractStream{R,T})
  x.a = b
  x.b = EmptyStream{R,T}()
end

function sound(cs::CatStream,len::Int)
  first_part = sound(cs.a,len)
  if nsamples(first_part) < len
    step_next(cs,cs.b)
    result = [first_part; sound(cs,len - nsamples(first_part))]

    result
  else
    first_part
  end
end

function playable{R}(s::AbstractStream{R},cache=true,sample_rate=samplerate())
  if ustrip(sample_rate) != R
    error("Cannot convert a stream at sample rate $(R*Hz) to $sample_rate. ",
          "Create a stream at $sample_rate instead.")
  end
  s
end

"""
For a stream, `fn` will be applied to each unit of sound as it is requested
from the stream.
"""
@inline function audiofn{R,T}(fn::Function,stream::Stream{R,T})
  Stream{R,T}(i -> fn(stream.fn(i)),stream.index)
end

# TODO: we need to lookup the guidelines for implementing broadcast
# before implementing this
# # TODO: this is mostly useful Julia 0.6 where ./ etc... invokes broadcast
# function broadcast{R,T}(f::Function,xs::Stream{R,T}...)
#   let offsets = Int[]
#     Stream{R,T}(xs[1].index) do i
#       if empty(offsets)
#         indices = map(x -> x.index,xs)
#         append!(offsets,[0; indices[1] - indices[2:end]])

#         # TODO: find max, as per soundop
#       end
#       f.(map((x,offset) -> x.fn(i-offset),xs,offsets)...)
#     end
#   end
# end

const emptyfn = () -> nothing
mutable struct OpStream{R,T} <: AbstractStream{R,T}
  op::Function
  a::AbstractStream{R,T}
  b::AbstractStream{R,T}
end

asstream(x::AbstractArray) = Stream(x)
asstream(x::AbstractStream) = x

function *(x::Number,stream::Stream)
  audiofn(stream) do stream
    x * stream
  end
end

function *(stream::Stream,x::Number)
  audiofn(stream) do stream
    x * stream
  end
end

soundop(op,xs::AbstractStream...) =
  error("Incompatible sample rates. Use `resample`.")
soundop{R}(op,xs::Union{Audible{R},Array}...) = soundop(op,map(asstream,xs)...)
soundop{R}(op,xs::AbstractStream{R}...) = reduce((x,y) -> soundop(op,x,y),xs)
soundop(op,x::AbstractStream) = x
function soundop{R,T,S}(op,x::AbstractStream{R,T},y::AbstractStream{R,S})
  error("The streams have different element types, convert to a common, ",
        "element type first.")
  #P = promote_type(T,S)
  #helper(x,_) = map(e -> convert(P,e),x)
  #soundop(op,...
end
function soundop{R,T}(op,x::AbstractStream{R,T},y::AbstractStream{R,T})
  if x !== y && any(i ∈ index_ids(x) for i ∈ index_ids(y))
    error("Cannot combine two streams derrived from the same underlying ",
          " stream. Please use `deepcopy` on one of the streams you are trying",
          " to combine.")
  end
  OpStream{R,T}(op,x,y)
end

duration(x::OpStream) = max(duration(x.a),duration(x.b))
nsamples(x::OpStream) = max(nsamples(x.a),nsamples(x.b))
index(x::OpStream) = index(x.a)
index_ids(x::OpStream) = [index_ids(x.a); index_ids(x.b)]

function sound{R,T}(os::OpStream{R,T},len::Int)
  if os.op === emptyfn
    sound(os.a,len)
  elseif os.a === os.b
    a = sound(os.a,len)
    os.op(a,a)
  else
    na = nsamples(os.a)
    nb = nsamples(os.b)
    N_all = min(na,nb,len)
    x = sound(os.op(sound(os.a,N_all),sound(os.b,N_all)))

    if nb < len && na > nb
      x = [x; sound(os.a,min(len - nb,na - nb))]
      os.op = emptyfn
      os.b = os.a
    elseif na < len && nb > na
      x = [x; sound(os.b,min(len - na,nb - na))]
      os.op = emptyfn
      os.a = os.b
    end

    x
  end
end

function audiofn(fn::Function,stream::AbstractStream)
  soundop(stream,stream) do x,_
    fn(x)
  end
end

leftright{R,T}(x::AbstractStream{R,T},y::AbstractStream{R,T}) =
  soundop(leftright,x,y)

left(x::AbstractStream) = audiofn(left,x)
right(x::AbstractStream) = audiofn(right,x)


text:
import PeriodicTable

# Alias to avoid similarity of elements and Element in DFTK module namespace
periodic_table = PeriodicTable.elements

# Data structure for chemical element and the potential model via which
# they interact with electrons. A compensating charge background is
# always assumed. It is assumed that each implementing struct
# defines at least the functions `local_potential_fourier` and `local_potential_real`.
# Very likely `charge_nuclear` and `charge_ionic` need to be defined as well.
abstract type Element end

"""Return the total nuclear charge of an atom type"""
charge_nuclear(::Element) = 0
# This is a fallback implementation that should be altered as needed.

"""Return the total ionic charge of an atom type (nuclear charge - core electrons)"""
charge_ionic(el::Element) = charge_nuclear(el)

"""Return the number of valence electrons"""
n_elec_valence(el::Element) = charge_ionic(el)

"""Return the number of core electrons"""
n_elec_core(el::Element) = charge_nuclear(el) - charge_ionic(el)

"""Radial local potential, in Fourier space: V(q) = int_{R^3} V(x) e^{-iqx} dx."""
function local_potential_fourier(el::Element, q::AbstractVector)
    local_potential_fourier(el, norm(q))
end

"""Radial local potential, in real space."""
function local_potential_real(el::Element, q::AbstractVector)
    local_potential_real(el, norm(q))
end

struct ElementCoulomb <: Element
    Z::Int  # Nuclear charge
    symbol  # Element symbol
end
charge_ionic(el::ElementCoulomb)   = el.Z
charge_nuclear(el::ElementCoulomb) = el.Z

"""
Element interacting with electrons via a bare Coulomb potential
(for all-electron calculations)
`key` may be an element symbol (like `:Si`), an atomic number (e.g. `14`)
or an element name (e.g. `"silicon"`)
"""
ElementCoulomb(key) = ElementCoulomb(periodic_table[key].number, Symbol(periodic_table[key].symbol))

function local_potential_fourier(el::ElementCoulomb, q::T) where {T <: Real}
    q == 0 && return zero(T)  # Compensating charge background
    # General atom => Use default Coulomb potential
    # We use int_{R^3} -Z/r e^{-i q⋅x} = 4π / |q|^2
    return -4T(π) * el.Z / q^2
end

local_potential_real(el::ElementCoulomb, r::Real) = -el.Z / r

struct ElementPsp <: Element
    Z::Int  # Nuclear charge
    symbol  # Element symbol
    psp     # Pseudopotential data structure
end

"""
Element interacting with electrons via a pseudopotential model.
`key` may be an element symbol (like `:Si`), an atomic number (e.g. `14`)
or an element name (e.g. `"silicon"`)
"""
function ElementPsp(key; psp)
    ElementPsp(periodic_table[key].number, Symbol(periodic_table[key].symbol), psp)
end
charge_ionic(el::ElementPsp)   = el.psp.Zion
charge_nuclear(el::ElementPsp) = el.Z

function local_potential_fourier(el::ElementPsp, q::T) where {T <: Real}
    q == 0 && return zero(T)  # Compensating charge background
    eval_psp_local_fourier(el.psp, q)
end

function local_potential_real(el::ElementPsp, r::Real)
    # Use local part of pseudopotential defined in Element object
    return eval_psp_local_real(el.psp, r)
end

struct ElementCohenBergstresser <: Element
    Z::Int  # Nuclear charge
    symbol  # Element symbol
    V_sym   # Map |G|^2 (in units of (2π / lattice_constant)^2) to form factors
    lattice_constant  # Lattice constant (in Bohr) which is assumed
end
charge_ionic(el::ElementCohenBergstresser)   = 2
charge_nuclear(el::ElementCohenBergstresser) = el.Z

"""
Element where the interaction with electrons is modelled
as in [CohenBergstresser1966](https://doi.org/10.1103/PhysRev.141.789).
Only the homonuclear lattices of the diamond structure
are implemented (i.e. Si, Ge, Sn).

`key` may be an element symbol (like `:Si`), an atomic number (e.g. `14`)
or an element name (e.g. `"silicon"`)
"""
function ElementCohenBergstresser(key; lattice_constant=nothing)
    # Form factors from Cohen-Bergstresser paper Table 2, converted to Bohr
    # Lattice constants from Table 1, converted to Bohr
    data = Dict(:Si => (form_factors=Dict( 3 => -0.21u"Ry",
                                           8 =>  0.04u"Ry",
                                          11 =>  0.08u"Ry"),
                        lattice_constant=5.43u"Å"),
                :Ge => (form_factors=Dict( 3 => -0.23u"Ry",
                                           8 =>  0.01u"Ry",
                                          11 =>  0.06u"Ry"),
                        lattice_constant=5.66u"Å"),
                :Sn => (form_factors=Dict( 3 => -0.20u"Ry",
                                           8 =>  0.00u"Ry",
                                          11 =>  0.04u"Ry"),
                        lattice_constant=6.49u"Å"),
            )

    symbol = Symbol(periodic_table[key].symbol)
    if !(symbol in keys(data))
        error("Cohen-Bergstresser potential not implemented for element $symbol.")
    end
    isnothing(lattice_constant) && (lattice_constant = data[symbol].lattice_constant)
    lattice_constant = austrip(lattice_constant)

    # Unit-cell volume of the primitive lattice (used in DFTK):
    unit_cell_volume = det(lattice_constant / 2 .* [[0 1 1]; [1 0 1]; [1 1 0]])

    # The form factors in the Cohen-Bergstresser paper Table 2 are
    # with respect to normalized planewaves (i.e. not plain Fourier coefficients)
    # and are already symmetrized into a sin-cos basis (see derivation p. 141)
    # => Scale by Ω / 2 to get them into the DFTK convention
    V_sym = Dict(key => austrip(value) * unit_cell_volume / 2
                 for (key, value) in pairs(data[symbol].form_factors))

    ElementCohenBergstresser(periodic_table[key].number, symbol, V_sym, lattice_constant)
end

function local_potential_fourier(el::ElementCohenBergstresser, q::T) where {T <: Real}
    q == 0 && return zero(T)  # Compensating charge background

    # Get |q|^2 in units of (2π / lattice_constant)^2
    qsq_pi = Int(round(q^2 / (2π / el.lattice_constant)^2, digits=2))
    T(get(el.V_sym, qsq_pi, 0.0))
end