text:
##### Beginning of file

import Dates
import HTTP
import TimeZones

function default_repo_description(
        ;
        env::AbstractDict = ENV,
        from::Any = "",
        when::Any = Dates.now(TimeZones.localzone(),),
        time_zone::Dates.TimeZone = TimeZones.TimeZone("America/New_York"),
        by::Any = "",
        )::String

    _from::String = strip(string(from))
    from_string::String = ""
    if length(_from) == 0
        from_string = ""
    else
        from_string = string(
            " ",
            strip(string("from $(_from)",)),
            )
    end

    _when::String = ""
    date_time_string::String = ""
    if isa(when, TimeZones.ZonedDateTime)
        _when = strip(
            string(TimeZones.astimezone(when,time_zone,))
            )
    else
        _when = strip(
            string(when)
            )
    end
    if length(_when) == 0
        date_time_string = ""
    else
        date_time_string = string(
            " ",
            strip(string("on ",_when,)),
            )
    end

    by_string::String = ""
    _by::String = strip(
        string(by)
        )
    if length(_by) == 0
        by_string = ""
    else
        by_string = string(
            " ",
            strip(string("by ",_by,)),
            )
    end

    travis_string::String = ""
    if _is_travis_ci(env)
        TRAVIS_BUILD_NUMBER::String = strip(
            get(env, "TRAVIS_BUILD_NUMBER", "")
            )
        TRAVIS_JOB_NUMBER::String = strip(
            get(env, "TRAVIS_JOB_NUMBER", "")
            )
        TRAVIS_EVENT_TYPE::String = strip(
            get(env, "TRAVIS_EVENT_TYPE", "unknown-travis-event")
            )
        TRAVIS_BRANCH = strip(
            get(env, "TRAVIS_BRANCH", "unknown-branch")
            )
        TRAVIS_COMMIT = strip(
            get(env, "TRAVIS_COMMIT", "unknown-commit")
            )
        TRAVIS_PULL_REQUEST = strip(
            get(env,"TRAVIS_PULL_REQUEST","unknown-pull-request-number")
            )

        job_or_build_string::String = ""
        if length(TRAVIS_JOB_NUMBER) > 0
            job_or_build_string = string(
                " ",
                strip(string("job $(TRAVIS_JOB_NUMBER)")),
                )
        elseif length(TRAVIS_BUILD_NUMBER) > 0
            job_or_build_string = string(
                " ",
                strip(string("build $(TRAVIS_BUILD_NUMBER)")),
                )
        else
            job_or_build_string = ""
        end

        triggered_by_string::String = ""
        if lowercase(TRAVIS_EVENT_TYPE) == "push"
            triggered_by_string = string(
                " ",
                strip(
                    string(
                        ", triggered by the push of",
                        " commit \"$(TRAVIS_COMMIT)\"",
                        " to branch \"$(TRAVIS_BRANCH)\"",
                        )
                    ),
                )
        elseif lowercase(TRAVIS_EVENT_TYPE) == "pull_request"
            triggered_by_string = string(
                " ",
                strip(
                    string(
                        ", triggered by",
                        " pull request #$(TRAVIS_PULL_REQUEST)",
                        )
                    ),
                )
        elseif lowercase(TRAVIS_EVENT_TYPE) == "cron"
            triggered_by_string = string(
                " ",
                strip(
                    string(
                        ", triggered by Travis",
                        " cron job on",
                        " branch \"$(TRAVIS_BRANCH)\"",
                        )
                    ),
                )
        else
            triggered_by_string = string(
                " ",
                strip(
                    string(
                        ", triggered by Travis \"",
                        strip(TRAVIS_EVENT_TYPE),
                        "\" event on",
                        " branch \"$(TRAVIS_BRANCH)\"",
                        )
                    ),
                )
        end
        travis_string = string(
            " ",
            strip(
                string(
                    "via Travis",
                    job_or_build_string,
                    triggered_by_string,
                    )
                ),
            )
    else
        travis_string = ""
    end

    new_description::String = strip(
        string(
            "Last mirrored",
            from_string,
            date_time_string,
            by_string,
            travis_string,
            )
        )

    return new_description
end

function _url_exists(url::AbstractString)::Bool
    _url::String = strip(convert(String, url))
    result::Bool = try
        r = HTTP.request("GET", _url)
        @debug("HTTP GET request result: ", _url, r.status,)
        r.status == 200
    catch exception
        @debug(string("Ignoring exception"), exception,)
        false
    end
    if result
    else
        @debug(string("URL does not exist"), _url,)
    end
    return result
end

mutable struct _DummyOutputWrapperStruct{I, F, S, O}
    previous_time_seconds::I
    f::F
    interval_seconds::I
    dummy_output::S
    io::O
end

function _DummyOutputWrapperStruct(
        ;
        interval_seconds::I = 60,
        initial_offset_seconds::I = interval_seconds,
        f::F,
        dummy_output::S = "This is a dummy line of output",
        io::O = Base.stdout,
        )::_DummyOutputWrapperStruct{I, F, S, O} where
            I <: Integer where
            F <: Function where
            S <: AbstractString where
            O <: IO
    current_time_seconds::I = floor(I, time())
    initial_time_seconds::I = current_time_seconds + initial_offset_seconds
    wrapper_struct::_DummyOutputWrapperStruct{I, F, S} =
        _DummyOutputWrapperStruct(
            initial_time_seconds,
            f,
            interval_seconds,
            dummy_output,
            io,
            )
    return wrapper_struct
end

function (x::_DummyOutputWrapperStruct{I, F, S, O})() where
        I <: Integer where
        F <: Function where
        S <: AbstractString where
        O <: IO
    current_time_seconds::I = floor(I, time())
    previous_time_seconds::I = x.previous_time_seconds
    f::F = x.f
    interval_seconds::I = x.interval_seconds
    dummy_output::S = x.dummy_output
    io::O = x.io
    elapsed_seconds::Int = current_time_seconds - previous_time_seconds
    print_dummy_output::Bool = elapsed_seconds > interval_seconds
    if print_dummy_output
        println(io, dummy_output)
        x.previous_time_seconds = current_time_seconds
    end
    f_result = f()
    return f_result
end

function dummy_output_wrapper(
        ;
        f::F,
        interval_seconds::I = 60,
        initial_offset_seconds::I = interval_seconds,
        dummy_output::S = "This is a dummy line of output",
        io::O = Base.stdout,
        ) where
            I <: Integer where
            F <: Function where
            S <: AbstractString where
            O <: IO
    wrapper_struct::_DummyOutputWrapperStruct{I, F, S, O} =
        _DummyOutputWrapperStruct(
            ;
            f = f,
            interval_seconds = interval_seconds,
            initial_offset_seconds = initial_offset_seconds,
            dummy_output = dummy_output,
            )
    function my_wrapper_function()
        result = wrapper_struct()
        return result
    end
    return my_wrapper_function
end

function command_ran_successfully!!(
        cmd::Base.AbstractCmd;
        max_attempts::Integer = 10,
        max_seconds_per_attempt::Real = 540,
        seconds_to_wait_between_attempts::Real = 30,
        error_on_failure::Bool = true,
        last_resort_run::Bool = true,
        )::Bool
    success_bool::Bool = false

    my_false = dummy_output_wrapper(
        ;
        f = () -> false,
        interval_seconds = 60,
        initial_offset_seconds = 60,
        dummy_output = "Still waiting between attempts...",
        io = Base.stdout,
        )

    for attempt = 1:max_attempts
        if success_bool
        else
            @debug(string("Attempt $(attempt)"))
            if attempt > 1
                timedwait(
                    () -> my_false(),
                    float(seconds_to_wait_between_attempts);
                    pollint = float(1.0),
                    )
            end
            p = run(cmd; wait = false,)
            my_process_exited = dummy_output_wrapper(
                ;
                f = () -> process_exited(p),
                interval_seconds = 60,
                initial_offset_seconds = 60,
                dummy_output = "The process is still running...",
                io = Base.stdout,
                )
            timedwait(
                () -> my_process_exited(),
                float(max_seconds_per_attempt),
                pollint = float(1.0),
                )
            if process_running(p)
                success_bool = false
                try
                    kill(p, Base.SIGTERM)
                catch exception
                    @warn("Ignoring exception: ", exception)
                end
                try
                    kill(p, Base.SIGKILL)
                catch exception
                    @warn("Ignoring exception: ", exception)
                end
            else
                success_bool = try
                    success(p)
                catch exception
                    @warn("Ignoring exception: ", exception)
                    false
                end
            end
        end
    end
    if !success_bool && last_resort_run
        @debug(string("Attempting the last resort run..."))
        run(cmd)
        success_bool = true
    end
    if success_bool
        @debug(string("Command ran successfully."),)
    else
        if error_on_failure
            error(string("Command did not run successfully."),)
        else
            @warn(string("Command did not run successfully."),)
        end
    end
    return success_bool
end

function retry_function_until_success(
        f::Function;
        max_attempts::Integer = 10,
        seconds_to_wait_between_attempts::Real = 30,
        )
    success_bool::Bool = false
    f_result = nothing

    my_false = dummy_output_wrapper(
        ;
        f = () -> false,
        interval_seconds = 60,
        initial_offset_seconds = 60,
        dummy_output = "Still waiting between attempts...",
        io = Base.stdout,
        )

    for attempt = 1:max_attempts
        if success_bool
        else
            @debug(string("Attempt $(attempt)"))
            if attempt > 1
                timedwait(
                    () -> my_false(),
                    float(seconds_to_wait_between_attempts);
                    pollint = float(1.0),
                    )
            end
            @debug(string("Running the provided function..."))
            success_bool = true
            f_result = try
                f()
            catch exception
                success_bool = false
                @warn("Ignoring exception: ", exception)
                nothing
            end
        end
    end

    if success_bool
        @debug(string("Function ran successfully."),)
        return f_result
    else
        error(string("Function did not run successfully."),)
    end
end

function _is_travis_ci(
        a::AbstractDict = ENV,
        )::Bool
    ci::String = lowercase(
        strip(get(a, "CI", "false"))
        )
    travis::String = lowercase(
        strip(get(a, "TRAVIS", "false"))
        )
    continuous_integration::String = lowercase(
        strip(get(a, "CONTINUOUS_INTEGRATION", "false"))
        )
    ci_is_true::Bool = ci == "true"
    travis_is_true::Bool = travis == "true"
    continuous_integration_is_true::Bool = continuous_integration == "true"
    answer::Bool = ci_is_true &&
        travis_is_true &&
        continuous_integration_is_true
    return answer
end

##### End of file


text:
module ParetoRecipes

using RecipesBase

"""
Compute the Pareto front of a set of points.

Uses an inefficient brute-force exact method.

The Pareto front is also known as the Pareto frontier or Pareto set.
"""
function pareto_front(points, minimize = true)
    cmp_strict = minimize ? (<) : (>);
    cmp = minimize ? (<=) : (>=);
    dims = size(points,2)
    strictly_dominates(u, v) = all(cmp(u[i], v[i]) for i = 1:dims+1) && any(cmp_strict(u[i], v[i]) for i = 1:dims+1)
    undominated(p) = !any(strictly_dominates(p2, p) for p2 ∈ points if p2 != p)
    front = filter(undominated, points)
    sort(front)
end

@userplot Pareto

@recipe function f(h::Pareto)
    if length(h.args) != 1 # TODO: add a typecheck if feasible
        error("Pareto plots require a single array of vectors.  Got: $(typeof(h.args))")
    end
    points = h.args[1]

    legend := false

    # all points without annotated front
    @series begin
        seriestype := :scatter
        [p[1] for p in points], [p[2] for p in points]
    end

    # add front
    @series begin
        seriestype := :path
        front = pareto_front(points)
        ([p[1] for p in front], [p[2] for p in front])
    end

    ()
end

export pareto_front
export pareto

end

ext: jl


text:
baremodule BenchPerf

function run end
# function save end
# function load end

function perfpath end
function loadperf end
function wrap end

module Internal

import JSON
import PrettyTables
import Tables
using BenchmarkTools

if !isdefined(@__MODULE__, Symbol("@something"))
    using Compat: @something
    const USE_COMPAT = true
else
    const USE_COMPAT = false
end

using ..BenchPerf: BenchPerf

include("perf.jl")
include("execution.jl")
include("wrap.jl")
include("accessors.jl")
include("show.jl")
include("tables.jl")

end  # module Internal

end  # baremodule BenchPerf


text:
"""
    GreensMeasurement(mc::DQMC, model)

Measures the equal time Greens function of the given DQMC simulation and model.
"""
struct GreensMeasurement{OT <: AbstractObservable} <: AbstractMeasurement
    obs::OT
end
function GreensMeasurement(mc::DQMC, model)
    o = LightObservable(
        LogBinner(zeros(eltype(mc.s.greens), size(mc.s.greens))),
        "Equal-times Green's function",
        "Observables.jld",
        "Equal-times Green's function"
    )
    GreensMeasurement{typeof(o)}(o)
end
function measure!(m::GreensMeasurement, mc::DQMC, model, i::Int64)
    push!(m.obs, greens(mc))
end

"""
    BosonEnergyMeasurement(mc::DQMC, model)

Measures the bosnic energy of the given DQMC simulation and model.

Note that this measurement requires `energy_boson(mc, model, conf)` to be
implemented for the specific `model`.
"""
struct BosonEnergyMeasurement{OT <: AbstractObservable} <: AbstractMeasurement
    obs::OT
end
function BosonEnergyMeasurement(mc::DQMC, model)
    o = LightObservable(Float64, name="Bosonic Energy", alloc=1_000_000)
    BosonEnergyMeasurement{typeof(o)}(o)
end
function measure!(m::BosonEnergyMeasurement, mc::DQMC, model, i::Int64)
    push!(m.obs, energy_boson(mc, model, conf(mc)))
end

function default_measurements(mc::DQMC, model)
    Dict(
        :conf => ConfigurationMeasurement(mc, model),
        :Greens => GreensMeasurement(mc, model),
        :BosonEnergy => BosonEnergyMeasurement(mc, model)
    )
end


text:
# p26.m - eigenvalues of 2nd-order Chebyshev diff. matrix

N = 60; (D,x) = cheb(N); D2 = D^2; D2 = D2[2:N,2:N];
(lam,V) = eig(D2);
ii = sortperm(-lam); e = lam[ii]; V = V[:,ii];

# Plot eigenvalues:
clf(); axes([.1, .62, .8, .3]);
loglog(-e,".",markersize=4); ylabel("eigenvalue");
title("N = $N       max |λ| = $(signif(maximum(-e)/N^4,5)) \$N^4\$");
semilogy(2*N/pi*[1,1],[1,1e6],"--r");
text(2.1*N/pi,24,"2π / N",fontsize=12);

# Plot eigenmodes N/4 (physical) and N (nonphysical):
vN4 = [0; V[:,Int(N/4-1)]; 0];
xx = -1:.01:1; vv = polyval(polyfit(x,vN4),xx);
axes([.1, .36, .8, .15]); plot(xx,vv);
plot(x,vN4,"k.",markersize=4); title("eigenmode N/4");
vN = V[:,N-1];
axes([.1, .1, .8, .15]);
semilogy(x[2:N],abs.(vN)); axis([-1,1,5e-6,1]);
plot(x[2:N],abs.(vN),"k.",markersize=4);
title("absolute value of eigenmode N    (log scale)");


text:
for F in [:expect, :fidelity, :operator_fidelity]
    @eval Base.adjoint(::typeof($F)) = Adjoint($F)
    @eval Base.show(io::IO, ::Adjoint{Any,typeof($F)}) = print(io, "$($F)'")
    @eval Base.show(io::IO, ::MIME"text/plain", ::Adjoint{Any,typeof($F)}) = print(io, "$($F)'")
end

function (::Adjoint{Any,typeof(expect)})(op::AbstractBlock, reg_or_circuit)
    expect_g(op, reg_or_circuit)
end

function expect_g(op::AbstractBlock, circuit::Pair{<:ArrayReg,<:AbstractBlock})
    reg, c = circuit
    out = copy(reg) |> c
    outδ = copy(out) |> op
    (in, inδ), paramsδ = apply_back((out, outδ), c)
    return inδ => paramsδ .* 2
end

function expect_g(op::AbstractBlock, reg::ArrayReg)
    copy(reg) |> op
end

_eval(p::Pair{<:AbstractRegister,<:AbstractBlock}) = copy(p.first) |> p.second
_eval(reg::AbstractRegister) = reg
YaoBase.fidelity(p1, p2) = fidelity(_eval(p1), _eval(p2))

function (::Adjoint{Any,typeof(fidelity)})(
    reg1::Union{ArrayReg,Pair{<:ArrayReg,<:AbstractBlock}},
    reg2::Union{ArrayReg,Pair{<:ArrayReg,<:AbstractBlock}},
)
    fidelity_g(reg1, reg2)
end

function fidelity_g(
    reg1::Union{ArrayReg,Pair{<:ArrayReg,<:AbstractBlock}},
    reg2::Union{ArrayReg,Pair{<:ArrayReg,<:AbstractBlock}},
)
    if reg1 isa Pair
        in1, c1 = reg1
        out1 = copy(in1) |> c1
    else
        out1 = reg1
    end

    if reg2 isa Pair
        in2, c2 = reg2
        out2 = copy(in2) |> c2
    else
        out2 = reg2
    end
    if nremain(out1) != 0
        throw(ArgumentError("The gradient of registers with environment is not implemented yet.
        However, back propagating over a focused register is possible,
        please file an issue if you really need this feature."))
    end
    overlap = out1' * out2

    out1δ = copy(out2)
    regscale!.(out1δ, conj.(overlap) ./ 2 ./ abs.(overlap))
    out2δ = copy(out1)
    regscale!.(out2δ, overlap ./ 2 ./ abs.(overlap))

    if reg1 isa Pair
        (_, in1δ), params1δ = apply_back((out1, out1δ), c1)
        res1 = in1δ => params1δ .* 2
    else
        res1 = out1δ
    end

    if reg2 isa Pair
        (_, in2δ), params2δ = apply_back((out2, out2δ), c2)
        res2 = in2δ => params2δ .* 2
    else
        res2 = out2δ
    end

    return res1, res2
end

function (::Adjoint{Any,typeof(operator_fidelity)})(b1::AbstractBlock, b2::AbstractBlock)
    operator_fidelity_g(b1, b2)
end

function operator_fidelity_g(b1::AbstractBlock, b2::AbstractBlock)
    U1 = mat(b1)
    U2 = mat(b2)
    s = sum(conj(U1) .* U2)
    adjs = conj(s) / abs(s) / size(U1, 1)
    adjm1 = U2 * adjs
    adjm2 = U1 * conj(adjs)
    mat_back(b1, adjm1), mat_back(b2, adjm2)
end


text:
module cemgym

# this is a port of the example found at
#	http://rl-gym-doc.s3-website-us-west-2.amazonaws.com/mlss/lab1.html#starter

using Reinforce
# using OpenAIGym
# using Distributions
using Transformations
using StochasticOptimization
using LearningStrategies

using Plots; gr(size=(500,200))

# ----------------------------------------------------------------

"""
Wraps a LearnBase.Transformation which converts an input vector to action values.
For discrete actions, it chooses the action which produces the highest value.
For continuous (interval) actions, it squashes actions to [0,1].
"""
struct TransformPolicy{T} <: AbstractPolicy
	trans::T
end
Transformations.params(tp::TransformPolicy) = params(tp.trans)
Reinforce.action(π::TransformPolicy, r, s, A::DiscreteSet) = A[indmax(transform!(π.trans, s))]

# # continuous: return the transform value, squashed to [0,1]
# # TODO: remove the squashing when a "Affine + Sigmoid" transformation is available
# function Reinforce.action(π::TransformPolicy, r, s, A::IntervalSet)
# 	Transformations.sigmoid(transform(π.trans, s)[1]) * (A.amax-A.amin) + A.amin
# end

# ----------------------------------------------------------------

# initialize a policy, do the learning, then return the policy
function do_cem_test(sublearners...; env = GymEnv("CartPole-v0"),
                       maxsteps = 200, # max number of steps in one episode
                       maxiter = 1000, # max learning iterations
					   # noise_max = 1.0,
					   # noise_steps = 20,
					   # noise_func = t -> max(noise_max - t/noise_steps, 0.0),
					   kw...)
    # generic query of state and action size
    s = state(env)
    A = actions(env,s)
    @assert isa(A, DiscreteSet) # this is the only kind that will work right now
	nS, nA = map(length, (s, A))

    # create a simple policy: action = argmax(wx+b)
    policy = TransformPolicy(Affine(nS, nA))

    # initialize the CEM, which will learn Θ = {w,b}
    strat = CrossEntropyMethod(;maxsteps=maxsteps, kw...)

    # keep a vector of test episode returns.  after each iteration, run (and plot)
    # an episode using the current CEM μ
    Rs = zeros(0)
    theme(:dark)
    function iterfunc(model, i)
        copy!(params(policy), strat.μ)
        R = run_episode(env, policy, maxsteps=maxsteps) do
            plot(plot(env), plot(Rs, label="Test Reward")) |> display
        end
        push!(Rs, R)
    end

    # create a MetaLearner driven by the CEM strategy
    learner = strategy(strat, sublearners...;
        maxiter=maxiter,
        oniter=iterfunc,
        kw...)

    # do the learning.  our iterator just repeatedly gives us the environment
    learn!(policy, learner, repeated(env))

    @show policy strat
    policy, strat
end

end #module


text:
module Test

include("util_macros.jl")
include("runtime.jl")
include("agent.jl")
include("signal.jl")
include("queue.jl")
include("aql.jl")
include("memory.jl")
include("ext_finalization.jl")
include("code.jl")
include("isa.jl")
include("intrinsics.jl")
include("emulation.jl")
if HSA.has_hsa_codegen()
include("execution.jl")
else
warn("codegen is unavailable: skipping 'Execution'")
end

end


text:
function LofOffsetPointing(arg0::Frame, arg1::BodyShape, arg2::AttitudeProvider, arg3::Vector3D)
    return LofOffsetPointing((Frame, BodyShape, AttitudeProvider, Vector3D), arg0, arg1, arg2, arg3)
end

function get_attitude(obj::LofOffsetPointing, arg0::FieldPVCoordinatesProvider, arg1::FieldAbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", FieldAttitude, (FieldPVCoordinatesProvider, FieldAbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_attitude(obj::LofOffsetPointing, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", Attitude, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_target_pv(obj::LofOffsetPointing, arg0::FieldPVCoordinatesProvider, arg1::FieldAbsoluteDate, arg2::Frame)
    return jcall(obj, "getTargetPV", TimeStampedFieldPVCoordinates, (FieldPVCoordinatesProvider, FieldAbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_target_pv(obj::LofOffsetPointing, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getTargetPV", TimeStampedPVCoordinates, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end


text:
import Base: +, -, *, /, \

#--------------------------------------------------
# Vector space algebra

# Unary ops
-(a::StaticArrayLite) = map(-, a)

# Binary ops
# Between arrays
+(a::StaticArrayLite, b::StaticArrayLite) = map(+, a, b)
+(a::AbstractArray, b::StaticArrayLite) = map(+, a, b)
+(a::StaticArrayLite, b::AbstractArray) = map(+, a, b)

-(a::StaticArrayLite, b::StaticArrayLite) = map(-, a, b)
-(a::AbstractArray, b::StaticArrayLite) = map(-, a, b)
-(a::StaticArrayLite, b::AbstractArray) = map(-, a, b)

# Scalar-array
*(a::Number, b::StaticArrayLite) = map(x -> a*x, b)
*(a::StaticArrayLite, b::Number) = map(x -> x*b, a)

/(a::StaticArrayLite, b::Number) = map(x -> x/b, a)
\(a::Number, b::StaticArrayLite) = map(x -> a\x, b)

#--------------------------------------------------
# Matrix algebra

# Transpose, conjugate, etc
conj(a::StaticArrayLite) = map(conj, a)

dot(a::StaticArrayLite, b::LinearAlgebra.Adjoint{<:StaticArrayLite}) = a*adjoint(b)
bilinear_vecdot(a::StaticArrayLite, b::StaticArrayLite) = a*transpose(b)

*(a::LinearAlgebra.Transpose{<:StaticArrayLite{<:Number, 1}}, b::StaticArrayLite{<:Number,1}) = mapreduce(*, +, a, b)
*(a::LinearAlgebra.Adjoint{<:StaticArrayLite{<:Real, 1}}, b::StaticArrayLite{<:Number,1}) = mapreduce(*, +, a, b)
*(a::LinearAlgebra.Adjoint{<:StaticArrayLite{<:Complex, 1}}, b::StaticArrayLite{<:Number,1}) = mapreduce((x, y) -> conj(x)*y, +, a, b)

#--------------------------------------------------
# Norms
LinearAlgebra.norm_sqr(a::StaticArrayLite) = mapreduce(abs2, +, a; init=zero(real(eltype(a))))
LinearAlgebra.norm(a::StaticArrayLite) = sqrt(LinearAlgebra.norm_sqr(a))

LinearAlgebra.normalize(a::StaticArrayLite) = inv(norm(a))*a

function LinearAlgebra.normalize!(a::StaticArrayLite)
    c = inv(norm(a))
    map!(x -> c*x, a, a)
    return a
end


text:
using TransformVariables
using NamedTupleTools

@inline function asflat(d::Dists.RealInterval)
    lb = isfinite(d.lb) ? d.lb : -∞
    ub = isfinite(d.lb) ? d.ub : ∞
    return as(Real, lb, ub)
end

@inline asflat(m::Dists.ContinuousUnivariateDistribution) = asflat(Dists.support(m))

@inline asflat(d::Dists.Dirichlet) = TransformVariables.UnitSimplex(length(d.alpha))

@inline asflat(d::Dists.MvNormal) = as(Vector, length(d.μ))

@inline function asflat(d::Dists.Product{T}) where {T<:Dists.ContinuousUnivariateDistribution}
    @assert !Base.isabstracttype(T) "$d is abstract type this isn't a homogenous product dist which isn't currently supported"
    as(Vector, asflat(d[1]), length(d.pars))
end

@inline asflat(d::MT.ProductMeasure) = as(Vector, asflat(first(MT.marginals(d))), length(d.pars))

@inline asflat(d::NamedTuple) = as(prototype(d)(asflat.(fieldvalues(d))))
@inline asflat(d::Tuple) = as(asflat.(d))


text:
# Conversions between Julia and Python types for the PyCall module.

#########################################################################
# Conversions of simple types (numbers and nothing)

# conversions from Julia types to PyObject:

@static if pyversion < v"3"
    PyObject(i::Unsigned) = PyObject(@pycheckn ccall(@pysym(:PyInt_FromSize_t),
                                                    PyPtr, (UInt,), i))
    PyObject(i::Integer) = PyObject(@pycheckn ccall(@pysym(:PyInt_FromSsize_t),
                                                    PyPtr, (Int,), i))
else
    PyObject(i::Unsigned) = PyObject(@pycheckn ccall(@pysym(:PyLong_FromUnsignedLongLong),
                                                    PyPtr, (Culonglong,), i))
    PyObject(i::Integer) = PyObject(@pycheckn ccall(@pysym(:PyLong_FromLongLong),
                                                    PyPtr, (Clonglong,), i))
end

PyObject(b::Bool) = PyObject(@pycheckn ccall((@pysym :PyBool_FromLong),
                                             PyPtr, (Clong,), b))

PyObject(r::Real) = PyObject(@pycheckn ccall((@pysym :PyFloat_FromDouble),
                                             PyPtr, (Cdouble,), r))

PyObject(c::Complex) = PyObject(@pycheckn ccall((@pysym :PyComplex_FromDoubles),
                                                PyPtr, (Cdouble,Cdouble),
                                                real(c), imag(c)))

PyObject(n::Nothing) = pyerr_check("PyObject(nothing)", pyincref(pynothing[]))

# conversions to Julia types from PyObject

@static if pyversion < v"3"
    convert(::Type{T}, po::PyObject) where {T<:Integer} =
        T(@pycheck ccall(@pysym(:PyInt_AsSsize_t), Int, (PyPtr,), po))
elseif pyversion < v"3.2"
    convert(::Type{T}, po::PyObject) where {T<:Integer} =
        T(@pycheck ccall(@pysym(:PyLong_AsLongLong), Clonglong, (PyPtr,), po))
else
    function convert(::Type{T}, po::PyObject) where {T<:Integer}
        overflow = Ref{Cint}()
        val = T(@pycheck ccall(@pysym(:PyLong_AsLongLongAndOverflow), Clonglong, (PyPtr, Ref{Cint}), po, overflow))
        iszero(overflow[]) || throw(InexactError(:convert, T, po))
        return val
    end
    function convert(::Type{Integer}, po::PyObject)
        overflow = Ref{Cint}()
        val = @pycheck ccall(@pysym(:PyLong_AsLongLongAndOverflow), Clonglong, (PyPtr, Ref{Cint}), po, overflow)
        iszero(overflow[]) || return convert(BigInt, po)
        return val
    end
end

convert(::Type{Bool}, po::PyObject) =
    0 != @pycheck ccall(@pysym(:PyObject_IsTrue), Cint, (PyPtr,), po)

convert(::Type{T}, po::PyObject) where {T<:Real} =
    T(@pycheck ccall(@pysym(:PyFloat_AsDouble), Cdouble, (PyPtr,), po))

convert(::Type{T}, po::PyObject) where T<:Complex =
    T(@pycheck ccall(@pysym(:PyComplex_AsCComplex), Complex{Cdouble}, (PyPtr,), po))

convert(::Type{Nothing}, po::PyObject) = nothing

function Base.float(o::PyObject)
    a = PyAny(o)
    if a isa PyObject
        hasproperty(o, :__float__) && return o.__float__()
        throw(ArgumentError("don't know how convert $o to a Julia floating-point value"))
    end
    return float(a)
end

#########################################################################
# String conversions (both bytes arrays and unicode strings)

function PyObject(s::AbstractString)
    sb = String(s)
    if pyunicode_literals || !isascii(sb)
        PyObject(@pycheckn ccall(@pysym(PyUnicode_DecodeUTF8),
                                 PyPtr, (Ptr{UInt8}, Int, Ptr{UInt8}),
                                 sb, sizeof(sb), C_NULL))
    else
        pybytes(sb)
    end
end

const _ps_ptr= Ptr{UInt8}[C_NULL]
const _ps_len = Int[0]
function convert(::Type{T}, po::PyObject) where T<:AbstractString
    if pyisinstance(po, @pyglobalobj :PyUnicode_Type)
        convert(T, PyObject(@pycheckn ccall(@pysym(PyUnicode_AsUTF8String),
                                             PyPtr, (PyPtr,), po)))
    else
        @pycheckz ccall(@pysym(PyString_AsStringAndSize),
                        Cint, (PyPtr, Ptr{Ptr{UInt8}}, Ptr{Int}),
                        po, _ps_ptr, _ps_len)
        convert(T, unsafe_string(_ps_ptr[1], _ps_len[1]))
    end
end

# TODO: should symbols be converted to a subclass of Python strings/bytes,
#       so that PyAny conversion can convert it back to a Julia symbol?
PyObject(s::Symbol) = PyObject(string(s))
convert(::Type{Symbol}, po::PyObject) = Symbol(convert(AbstractString, po))

#########################################################################
# ByteArray conversions

function PyObject(a::DenseVector{UInt8})
  if stride(a,1) != 1
    try
        return NpyArray(a, true)
    catch
        return array2py(a) # fallback to non-NumPy version
    end
  end
  PyObject(@pycheckn ccall((@pysym :PyByteArray_FromStringAndSize),
                           PyPtr, (Ptr{UInt8}, Int), a, length(a)))
end

ispybytearray(po::PyObject) =
  pyisinstance(po, @pyglobalobj :PyByteArray_Type)

function convert(::Type{Vector{UInt8}}, po::PyObject)
    b = PyBuffer(po)
    iscontiguous(b) || error("a contiguous buffer is required")
    return copy(unsafe_wrap(Array, Ptr{UInt8}(pointer(b)), sizeof(b)))
end

# TODO: support zero-copy PyByteArray <: AbstractVector{UInt8} object

#########################################################################
# Pointer conversions, using ctypes or PyCapsule

PyObject(p::Ptr) = pycall(c_void_p_Type, PyObject, UInt(p))

function convert(::Type{Ptr{Cvoid}}, po::PyObject)
    if pyisinstance(po, c_void_p_Type)
        v = po."value"
        # ctypes stores the NULL pointer specially, grrr
        pynothing_query(v) == Nothing ? C_NULL :
          convert(Ptr{Cvoid}, convert(UInt, po."value"))
    elseif pyisinstance(po, @pyglobalobj(:PyCapsule_Type))
        @pycheck ccall((@pysym :PyCapsule_GetPointer),
                       Ptr{Cvoid}, (PyPtr,Ptr{UInt8}),
                       po, ccall((@pysym :PyCapsule_GetName),
                                 Ptr{UInt8}, (PyPtr,), po))
    else
        convert(Ptr{Cvoid}, convert(UInt, po))
    end
end

pyptr_query(po::PyObject) = pyisinstance(po, c_void_p_Type) || pyisinstance(po, @pyglobalobj(:PyCapsule_Type)) ? Ptr{Cvoid} : Union{}

#########################################################################
# for automatic conversions, I pass Vector{PyAny}, NTuple{N, PyAny}, etc.,
# but since PyAny is an abstract type I need to convert this to Any
# before actually creating the Julia object

# I want to use a union, but this seems to confuse Julia's method
# dispatch for the convert function in some circumstances
# const PyAny = Union{PyObject, Int, Bool, Float64, ComplexF64, AbstractString, Function, Dict, Tuple, Array}
abstract type PyAny end

function pyany_toany(T::Type)
    T === Vararg{PyAny} ? Vararg{Any} : T
end
pyany_toany(::Type{PyAny}) = Any
pyany_toany(t::Type{T}) where {T<:Tuple} = Tuple{map(pyany_toany, t.types)...}

# PyAny acts like Any for conversions, except for converting PyObject (below)
convert(::Type{PyAny}, x) = x

#########################################################################
# Function conversion (see callback.jl for conversion the other way)
# (rarely needed given call overloading in Julia 0.4)

convert(::Type{Function}, po::PyObject) =
    function fn(args...; kwargs...)
        pycall(po, PyAny, args...; kwargs...)
    end

#########################################################################
# Tuple conversion.  Julia Pairs are treated as Python tuples.

function PyObject(t::Union{Tuple,Pair})
    len = lastindex(t) # lastindex, not length, because of julia#14924
    o = PyObject(@pycheckn ccall((@pysym :PyTuple_New), PyPtr, (Int,), len))
    for i = 1:len
        oi = PyObject(t[i])
        @pycheckz ccall((@pysym :PyTuple_SetItem), Cint, (PyPtr,Int,PyPtr),
                         o, i-1, oi)
        pyincref(oi) # PyTuple_SetItem steals the reference
    end
    return o
end

# somewhat annoying to get the length and types in a tuple type
# ... would be better not to have to use undocumented internals!
istuplen(T,isva,n) = isva ? n ≥ length(T.parameters)-1 : n == length(T.parameters)
function tuptype(T::DataType,isva,i)
    if isva && i ≥ length(T.parameters)
        return Base.unwrapva(T.parameters[end])
    else
        return T.parameters[i]
    end
end
tuptype(T::UnionAll,isva,i) = tuptype(T.body,isva,i)
isvatuple(T::UnionAll) = isvatuple(T.body)
isvatuple(T::DataType) = !isempty(T.parameters) && Base.isvarargtype(T.parameters[end])

function convert(tt::Type{T}, o::PyObject) where T<:Tuple
    isva = isvatuple(T)
    len = @pycheckz ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
    if !istuplen(tt, isva, len)
        throw(BoundsError())
    end
    ntuple((i ->
            convert(tuptype(T, isva, i),
                    PyObject(ccall((@pysym :PySequence_GetItem), PyPtr,
                                   (PyPtr, Int), o, i-1)))),
           len)
end

function convert(::Type{Pair{K,V}}, o::PyObject) where {K,V}
    k, v = convert(Tuple{K,V}, o)
    return Pair(k, v)
end

#########################################################################
# PyVector: no-copy wrapping of a Julia object around a Python sequence

"""
    PyVector(o::PyObject)

This returns a PyVector object, which is a wrapper around an arbitrary Python list or sequence object.

Alternatively, `PyVector` can be used as the return type for a `pycall` that returns a sequence object (including tuples).
"""
mutable struct PyVector{T} <: AbstractVector{T}
    o::PyObject
    function PyVector{T}(o::PyObject) where T
        if ispynull(o)
            throw(ArgumentError("cannot make PyVector from NULL PyObject"))
        end
        new{T}(o)
    end
end

PyVector(o::PyObject) = PyVector{PyAny}(o)
PyObject(a::PyVector) = a.o
convert(::Type{PyVector}, o::PyObject) = PyVector(o)
convert(::Type{PyVector{T}}, o::PyObject) where {T} = PyVector{T}(o)
unsafe_convert(::Type{PyPtr}, a::PyVector) = PyPtr(a.o)
PyVector(a::PyVector) = a
PyVector(a::AbstractVector{T}) where {T} = PyVector{T}(array2py(a))

# when a PyVector is copied it is converted into an ordinary Julia Vector
similar(a::PyVector, T, dims::Dims) = Array{T}(dims)
similar(a::PyVector{T}) where {T} = similar(a, pyany_toany(T), size(a))
similar(a::PyVector{T}, dims::Dims) where {T} = similar(a, pyany_toany(T), dims)
similar(a::PyVector{T}, dims::Int...) where {T} = similar(a, pyany_toany(T), dims)
eltype(::PyVector{T}) where {T} = pyany_toany(T)
eltype(::Type{PyVector{T}}) where {T} = pyany_toany(T)

size(a::PyVector) = (length(a.o),)

getindex(a::PyVector) = getindex(a, 1)
getindex(a::PyVector{T}, i::Integer) where {T} = convert(T, PyObject(@pycheckn ccall((@pysym :PySequence_GetItem), PyPtr, (PyPtr, Int), a, i-1)))

setindex!(a::PyVector, v) = setindex!(a, v, 1)
function setindex!(a::PyVector, v, i::Integer)
    @pycheckz ccall((@pysym :PySequence_SetItem), Cint, (PyPtr, Int, PyPtr), a, i-1, PyObject(v))
    v
end

summary(a::PyVector{T}) where {T} = string(Base.dims2string(size(a)), " ",
                                          string(pyany_toany(T)), " PyVector")

splice!(a::PyVector, i::Integer) = splice!(a.o, i)
function splice!(a::PyVector{T}, indices::AbstractVector{I}) where {T,I<:Integer}
    v = pyany_toany(T)[a[i] for i in indices]
    for i in sort(indices, rev=true)
        @pycheckz ccall((@pysym :PySequence_DelItem), Cint, (PyPtr, Int), a, i-1)
    end
    v
end
pop!(a::PyVector) = pop!(a.o)
popfirst!(a::PyVector) = popfirst!(a.o)
empty!(a::PyVector) = empty!(a.o)

# only works for List subtypes:
push!(a::PyVector, item) = push!(a.o, item)
insert!(a::PyVector, i::Integer, item) = insert!(a.o, i, item)
pushfirst!(a::PyVector, item) = pushfirst!(a.o, item)
prepend!(a::PyVector, items) = prepend!(a.o, items)
append!(a::PyVector{T}, items) where {T} = PyVector{T}(append!(a.o, items))

#########################################################################
# Lists and 1d arrays.

if VERSION < v"1.1.0-DEV.392" # #29440
    cirange(I,J) = CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))
else
    cirange(I,J) = I:J
end

# recursive conversion of A to a list of list of lists... starting
# with dimension dim and Cartesian index i in A.
function array2py(A::AbstractArray{<:Any, N}, dim::Integer, i::CartesianIndex{N}) where {N}
    if dim > N # base case
        return PyObject(A[i])
    else # recursively store multidimensional array as list of lists
        ilast = CartesianIndex(ntuple(j -> j == dim ? lastindex(A, dim) : i[j], Val{N}()))
        o = PyObject(@pycheckn ccall((@pysym :PyList_New), PyPtr, (Int,), size(A, dim)))
        for icur in cirange(i,ilast)
            oi = array2py(A, dim+1, icur)
            @pycheckz ccall((@pysym :PyList_SetItem), Cint, (PyPtr,Int,PyPtr),
                             o, icur[dim]-i[dim], oi)
            pyincref(oi) # PyList_SetItem steals the reference
        end
        return o
    end
end

array2py(A::AbstractArray) = array2py(A, 1, first(CartesianIndices(A)))

PyObject(A::AbstractArray) =
   ndims(A) <= 1 || hasmethod(stride, Tuple{typeof(A),Int}) ? array2py(A) :
   pyjlwrap_new(A)

function py2array(T, A::Array{TA,N}, o::PyObject,
                  dim::Integer, i::Integer) where {TA,N}
    if dim > N
        A[i] = convert(T, o)
        return A
    elseif dim == N
        len = @pycheckz ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
        if len != size(A, dim)
            error("dimension mismatch in py2array")
        end
        s = stride(A, dim)
        for j = 0:len-1
            A[i+j*s] = convert(T, PyObject(ccall((@pysym :PySequence_GetItem),
                                                 PyPtr, (PyPtr, Int), o, j)))
        end
        return A
    else # dim < N: recursively extract list of lists into A
        len = @pycheckz ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
        if len != size(A, dim)
            error("dimension mismatch in py2array")
        end
        s = stride(A, dim)
        for j = 0:len-1
            py2array(T, A, PyObject(ccall((@pysym :PySequence_GetItem),
                                       PyPtr, (PyPtr, Int), o, j)),
                     dim+1, i+j*s)
        end
        return A
    end
end

# figure out if we can treat o as a multidimensional array, and return
# the dimensions
function pyarray_dims(o::PyObject, forcelist=true)
    if !(forcelist || pyisinstance(o, @pyglobalobj :PyList_Type))
        return () # too many non-List types can pretend to be sequences
    end
    len = ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
    len < 0 && error("not a PySequence object")
    if len == 0
        return (0,)
    end
    dims0 = pyarray_dims(PyObject(ccall((@pysym :PySequence_GetItem),
                                        PyPtr, (PyPtr, Int), o, 0)),
                         false)
    if isempty(dims0) # not a nested sequence
        return (len,)
    end
    for j = 1:len-1
        dims = pyarray_dims(PyObject(ccall((@pysym :PySequence_GetItem),
                                           PyPtr, (PyPtr, Int), o, j)),
                            false)
        if dims != dims0
            # elements don't have equal lengths, cannot
            # treat as multidimensional array
            return (len,)
        end
    end
    return tuple(len, dims0...)
end

function py2array(T, o::PyObject)
    b = PyBuffer()
    if isbuftype!(o, b)
        dims = size(b)
    else
        dims = pyarray_dims(o)
    end
    pydecref(b) # safe for immediate release
    A = Array{pyany_toany(T)}(undef, dims)
    py2array(T, A, o, 1, 1) # fixme: faster conversion for supported buffer types?
end

function py2vector(T, o::PyObject)
    len = ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
    if len < 0 || # not a sequence
       len+1 < 0  # object pretending to be a sequence of infinite length
        pyerr_clear()
        throw(ArgumentError("expected Python sequence"))
    end
    py2array(T, Array{pyany_toany(T)}(undef, len), o, 1, 1)
end
convert(::Type{Vector{T}}, o::PyObject) where T = py2vector(T, o)

convert(::Type{Array}, o::PyObject) = map(identity, py2array(PyAny, o))
convert(::Type{Array{T}}, o::PyObject) where {T} = py2array(T, o)

PyObject(a::BitArray) = PyObject(Array(a))

# NumPy conversions (multidimensional arrays)
include("numpy.jl")

#########################################################################
# PyDict: no-copy wrapping of a Julia object around a Python dictionary

# we check for "items" attr since PyMapping_Check doesn't do this (it only
# checks for __getitem__) and PyMapping_Check returns true for some
# scipy scalar array members, grrr.
function is_mapping_object(o::PyObject)
    pyisinstance(o, @pyglobalobj :PyDict_Type) ||
    (pyquery((@pyglobal :PyMapping_Check), o) &&
      ccall((@pysym :PyObject_HasAttrString), Cint, (PyPtr,Ptr{UInt8}), o, "items") == 1)
end

"""
    PyDict(o::PyObject)
    PyDict(d::Dict{K,V})

This returns a PyDict, which is a no-copy wrapper around a Python dictionary.

Alternatively, you can specify the return type of a `pycall` as PyDict.
"""
mutable struct PyDict{K,V,isdict} <: AbstractDict{K,V}
    o::PyObject
    # isdict = true for python dict, otherwise is a generic Mapping object

    function PyDict{K,V,isdict}(o::PyObject) where {K,V,isdict}
        if !isdict && !ispynull(o) && !is_mapping_object(o)
            throw(ArgumentError("only Dict and Mapping objects can be converted to PyDict"))
        end
        return new{K,V,isdict}(o)
    end
end

PyDict{K,V}(o::PyObject) where {K,V} = PyDict{K,V,pyisinstance(o, @pyglobalobj :PyDict_Type)}(o)
PyDict{K,V}() where {K,V} = PyDict{K,V,true}(PyObject(@pycheckn ccall((@pysym :PyDict_New), PyPtr, ())))

PyDict(o::PyObject) = PyDict{PyAny,PyAny}(o)
PyObject(d::PyDict) = d.o
PyDict() = PyDict{PyAny,PyAny}()
PyDict(d::AbstractDict{K,V}) where {K,V} = PyDict{K,V}(PyObject(d))
PyDict(d::AbstractDict{Any,Any}) = PyDict{PyAny,PyAny}(PyObject(d))
PyDict(d::AbstractDict{Any,V}) where {V} = PyDict{PyAny,V}(PyObject(d))
PyDict(d::AbstractDict{K,Any}) where {K} = PyDict{K,PyAny}(PyObject(d))
convert(::Type{PyDict}, o::PyObject) = PyDict(o)
convert(::Type{PyDict{K,V}}, o::PyObject) where {K,V} = PyDict{K,V}(o)
unsafe_convert(::Type{PyPtr}, d::PyDict) = PyPtr(d.o)

haskey(d::PyDict{K,V,true}, key) where {K,V} = 1 == ccall(@pysym(:PyDict_Contains), Cint, (PyPtr, PyPtr), d, PyObject(key))
keys(::Type{T}, d::PyDict{K,V,true}) where {T,K,V} = convert(Vector{T}, PyObject(@pycheckn ccall((@pysym :PyDict_Keys), PyPtr, (PyPtr,), d)))
values(::Type{T}, d::PyDict{K,V,true}) where {T,K,V} = convert(Vector{T}, PyObject(@pycheckn ccall((@pysym :PyDict_Values), PyPtr, (PyPtr,), d)))

keys(::Type{T}, d::PyDict{K,V,false}) where {T,K,V} = convert(Vector{T}, pycall(d.o["keys"], PyObject))
values(::Type{T}, d::PyDict{K,V,false}) where {T,K,V} = convert(Vector{T}, pycall(d.o["values"], PyObject))
haskey(d::PyDict{K,V,false}, key) where {K,V} = 1 == ccall(@pysym(:PyMapping_HasKey), Cint, (PyPtr, PyPtr), d, PyObject(key))

similar(d::PyDict{K,V}) where {K,V} = Dict{pyany_toany(K),pyany_toany(V)}()
eltype(::Type{PyDict{K,V}}) where {K,V} = Pair{pyany_toany(K),pyany_toany(V)}
Base.keytype(::PyDict{K,V}) where {K,V} = pyany_toany(K)
Base.valtype(::PyDict{K,V}) where {K,V} = pyany_toany(V)
Base.keytype(::Type{PyDict{K,V}}) where {K,V} = pyany_toany(K)
Base.valtype(::Type{PyDict{K,V}}) where {K,V} = pyany_toany(V)

function setindex!(d::PyDict, v, k)
    @pycheckz ccall((@pysym :PyObject_SetItem), Cint, (PyPtr, PyPtr, PyPtr),
                     d, PyObject(k), PyObject(v))
    v
end

get(d::PyDict{K,V}, k, default) where {K,V} = get(d.o, V, k, default)

function pop!(d::PyDict{K,V,true}, k) where {K,V}
    v = d[k]
    @pycheckz ccall(@pysym(:PyDict_DelItem), Cint, (PyPtr, PyPtr), d, PyObject(k))
    return v
end
function pop!(d::PyDict{K,V,false}, k) where {K,V}
    v = d[k]
    @pycheckz ccall(@pysym(:PyObject_DelItem), Cint, (PyPtr, PyPtr), d, PyObject(k))
    return v
end

function pop!(d::PyDict, k, default)
    try
        return pop!(d, k)
    catch
        return default
    end
end

function delete!(d::PyDict{K,V,true}, k) where {K,V}
    e = ccall(@pysym(:PyDict_DelItem), Cint, (PyPtr, PyPtr), d, PyObject(k))
    e == -1 && pyerr_clear() # delete! ignores errors in Julia
    return d
end
function delete!(d::PyDict{K,V,false}, k) where {K,V}
    e = ccall(@pysym(:PyObject_DelItem), Cint, (PyPtr, PyPtr), d, PyObject(k))
    e == -1 && pyerr_clear() # delete! ignores errors in Julia
    return d
end

function empty!(d::PyDict{K,V,true}) where {K,V}
    @pycheck ccall((@pysym :PyDict_Clear), Cvoid, (PyPtr,), d)
    return d
end
function empty!(d::PyDict{K,V,false}) where {K,V}
    # for generic Mapping items we must delete keys one by one
    for k in keys(d)
        delete!(d, k)
    end
    return d
end

length(d::PyDict{K,V,true}) where {K,V} = @pycheckz ccall(@pysym(:PyDict_Size), Int, (PyPtr,), d)
length(d::PyDict{K,V,false}) where {K,V} = @pycheckz ccall(@pysym(:PyObject_Size), Int, (PyPtr,), d)
isempty(d::PyDict) = length(d) == 0

struct PyDict_Iterator
    # arrays to pass key, value, and pos pointers to PyDict_Next
    ka::Ref{PyPtr}
    va::Ref{PyPtr}
    pa::Ref{Int}
    i::Int # current position in items list (0-based)
    len::Int # length of items list
end

function Base.iterate(d::PyDict{K,V,true}, itr=PyDict_Iterator(Ref{PyPtr}(), Ref{PyPtr}(), Ref(0), 0, length(d))) where {K,V}
    itr.i >= itr.len && return nothing
    if 0 == ccall((@pysym :PyDict_Next), Cint,
                    (PyPtr, Ref{Int}, Ref{PyPtr}, Ref{PyPtr}),
                    d, itr.pa, itr.ka, itr.va)
        error("unexpected end of PyDict_Next")
    end
    ko = pyincref(itr.ka[]) # PyDict_Next returns
    vo = pyincref(itr.va[]) #   borrowed ref, so incref
    (Pair(convert(K,ko), convert(V,vo)),
        PyDict_Iterator(itr.ka, itr.va, itr.pa, itr.i+1, itr.len))
end

# Iterator for generic mapping, using Python items iterator.
# Our approach is to wrap an iterator over d.o["items"]
# which necessitates including d.o["items"] in the state.
function _start(d::PyDict{K,V,false}) where {K,V}
    d_items = pycall(d.o."items", PyObject)
    (d_items, iterate(d_items))
end
function Base.iterate(d::PyDict{K,V,false}, itr=_start(d)) where {K,V}
    d_items, iter_result = itr
    iter_result === nothing && return nothing
    item, state = iter_result
    iter_result = iterate(d_items, state)
    (item[1] => item[2], (d_items, iter_result))
end

#########################################################################
# Dictionary conversions (copies)

function PyObject(d::AbstractDict)
    o = PyObject(@pycheckn ccall((@pysym :PyDict_New), PyPtr, ()))
    for k in keys(d)
        @pycheckz ccall((@pysym :PyDict_SetItem), Cint, (PyPtr,PyPtr,PyPtr),
                         o, PyObject(k), PyObject(d[k]))
    end
    return o
end

function convert(::Type{Dict{K,V}}, o::PyObject) where {K,V}
    copy(PyDict{K,V}(o))
end

#########################################################################
# AbstractRange: integer ranges are converted to xrange,
#                while other ranges (<: AbstractVector) are converted to lists

xrange(start, stop, step) = pycall(pyxrange[], PyObject,
                                   start, stop, step)

function PyObject(r::AbstractRange{T}) where T<:Integer
    s = step(r)
    f = first(r)
    l = last(r) + s
    if max(f,l) > typemax(Clong) || min(f,l) < typemin(Clong)
        # in Python 2.x, xrange is limited to Clong
        PyObject(T[r...])
    else
        xrange(f, l, s)
    end
end

function convert(::Type{T}, o::PyObject) where T<:AbstractRange
    v = PyVector(o)
    len = length(v)
    if len == 0
        return 1:0 # no way to get more info from an xrange
    elseif len == 1
        start = v[1]
        return start:start
    else
        start = v[1]
        stop = v[len]
        step = v[2] - start
        return step == 1 ? (start:stop) : (start:step:stop)
    end
end

#########################################################################
# BigFloat and Complex{BigFloat}: convert to/from Python mpmath types

# load mpmath module & initialize.  Currently, this is done
# the first time a BigFloat is converted to Python.  Alternatively,
# we could do it when PyCall is initialized (if mpmath is available),
# at the cost of slowing down initialization in the common case where
# BigFloat conversion is not needed.
const mpprec = [0]
const mpmath = PyNULL()
const mpf = PyNULL()
const mpc = PyNULL()
function mpmath_init()
    if ispynull(mpmath)
        copy!(mpmath, pyimport("mpmath"))
        copy!(mpf, mpmath."mpf")
        copy!(mpc, mpmath."mpc")
    end
    curprec = precision(BigFloat)
    if mpprec[1] != curprec
        mpprec[1] = curprec
        mpmath."mp"."prec" = mpprec[1]
    end
end

# TODO: When mpmath uses MPFR internally, can we avoid the string conversions?
# Using strings will work regardless of the mpmath backend, but is annoying
# both from a performance perspective and because it is a lossy conversion
# (since strings use a decimal representation, while MPFR is binary).

function PyObject(x::BigFloat)
    mpmath_init()
    pycall(mpf, PyObject, string(x))
end

function PyObject(x::Complex{BigFloat})
    mpmath_init()
    pycall(mpc, PyObject, string(real(x)), string(imag(x)))
end

convert(::Type{BigFloat}, o::PyObject) = parse(BigFloat, pystr(o))

function convert(::Type{Complex{BigFloat}}, o::PyObject)
    try
        Complex{BigFloat}(convert(BigFloat, o."real"),
                          convert(BigFloat, o."imag"))
    catch
        convert(Complex{BigFloat}, convert(Complex{Float64}, o))
    end
end

pymp_query(o::PyObject) = pyisinstance(o, mpf) ? BigFloat : pyisinstance(o, mpc) ? Complex{BigFloat} : Union{}

#########################################################################
# (Int64), Int128 and BigInt conversion to Python "long" integers

const LongInt = @static (Sys.WORD_SIZE==32) ? Union{Int64,UInt64,Int128,UInt128,BigInt} : Union{Int128,UInt128,BigInt}

function PyObject(i::LongInt)
    PyObject(@pycheckn ccall((@pysym :PyLong_FromString), PyPtr,
                             (Ptr{UInt8}, Ptr{Cvoid}, Cint),
                             String(string(i)), C_NULL, 10))
end

convert(::Type{BigInt}, o::PyObject) = parse(BigInt, pystr(o))

#########################################################################
# Dates (Calendar time)

include("pydates.jl")
#init_datetime() = nothing
#pydate_query(o) = Union{}

#########################################################################
# Inferring Julia types at runtime from Python objects:
#
# [Note that we sometimes use the PyFoo_Check API and sometimes we use
#  PyObject_IsInstance(o, PyFoo_Type), since sometimes the former API
#  is a macro (hence inaccessible in Julia).]

# A type-query function f(o::PyObject) returns the Julia type
# for use with the convert function, or Union{} if there isn't one.

@static if pyversion < v"3"
    pyint_query(o::PyObject) = pyisinstance(o, @pyglobalobj :PyInt_Type) ?
        (pyisinstance(o, @pyglobalobj :PyBool_Type) ? Bool : Int) :
        pyisinstance(o, @pyglobalobj :PyLong_Type) ? BigInt :
        pyisinstance(o, npy_integer) ? Int : Union{}
else
    pyint_query(o::PyObject) = pyisinstance(o, @pyglobalobj :PyLong_Type) ?
        (pyisinstance(o, @pyglobalobj :PyBool_Type) ? Bool : Integer) :
        pyisinstance(o, npy_integer) ? Integer : Union{}
end

pyfloat_query(o::PyObject) = pyisinstance(o, @pyglobalobj :PyFloat_Type) ||  pyisinstance(o, npy_floating) ? Float64 : Union{}

pycomplex_query(o::PyObject) =
    pyisinstance(o, @pyglobalobj :PyComplex_Type) ||  pyisinstance(o, npy_complexfloating) ? ComplexF64 : Union{}

pystring_query(o::PyObject) = pyisinstance(o, @pyglobalobj PyString_Type) ? AbstractString : pyisinstance(o, @pyglobalobj :PyUnicode_Type) ? String : Union{}

# Given call overloading, all PyObjects are callable already, so
# we never automatically convert to Function.
pyfunction_query(o::PyObject) = Union{}

pynothing_query(o::PyObject) = o ≛ pynothing[] ? Nothing : Union{}

# We refrain from converting all objects that support the mapping protocol (PyMapping_Check)
# to avoid converting types like Pandas `DataFrame` that are only lossily
# representable as a Julia dictionary (issue #376).
pydict_query(o::PyObject) = pyisinstance(o, @pyglobalobj :PyDict_Type) ? Dict{PyAny,PyAny} : Union{}

typetuple(Ts) = Tuple{Ts...}

function pysequence_query(o::PyObject)
    # pyquery(:PySequence_Check, o) always succeeds according to the docs,
    # but it seems we need to be careful; I've noticed that things like
    # scipy define "fake" sequence types with intmax lengths and other
    # problems
    if pyisinstance(o, @pyglobalobj :PyTuple_Type)
        len = @pycheckz ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)
        return typetuple(pytype_query(PyObject(ccall((@pysym :PySequence_GetItem), PyPtr, (PyPtr,Int), o,i-1)), PyAny) for i = 1:len)
    elseif pyisinstance(o, pyxrange[])
        return AbstractRange
    elseif ispybytearray(o)
        return Vector{UInt8}
    elseif !isbuftype(o)
        # only handle PyList for now
        return pyisinstance(o, @pyglobalobj :PyList_Type) ? Array : Union{}
    else
        T, native_byteorder = array_format(o)
        if T == PyPtr
            T = PyObject
        end
        return Array{T}
    end
end

macro return_not_None(ex)
    quote
        T = $(esc(ex))
        if T != Union{}
            return T
        end
    end
end

const pytype_queries = Tuple{PyObject,Type}[]
"""
    pytype_mapping(pytype, jltype)

Given a Python type object `pytype`, tell PyCall to convert it to
`jltype` in `PyAny(object)` conversions.
"""
function pytype_mapping(py::PyObject, jl::Type)
    for (i,(p,j)) in enumerate(pytype_queries)
        if p == py
            pytype_queries[i] = (py,jl)
            return pytype_queries
        end
    end
    push!(pytype_queries, (py,jl))
end
"""
    pytype_query(o::PyObject, default=PyObject)

Given a Python object `o`, return the corresponding
native Julia type (defaulting to `default`) that we convert
`o` to in `PyAny(o)` conversions.
"""
function pytype_query(o::PyObject, default::TypeTuple=PyObject)
    # TODO: Use some kind of hashtable (e.g. based on PyObject_Type(o)).
    #       (A bit tricky to correctly handle Tuple and other containers.)
    @return_not_None pyint_query(o)
    pyisinstance(o, npy_bool) && return Bool
    @return_not_None pyfloat_query(o)
    @return_not_None pycomplex_query(o)
    @return_not_None pystring_query(o)
    @return_not_None pyfunction_query(o)
    @return_not_None pydate_query(o)
    @return_not_None pydict_query(o)
    @return_not_None pyptr_query(o)
    @return_not_None pysequence_query(o)
    @return_not_None pynothing_query(o)
    @return_not_None pymp_query(o)
    for (py,jl) in pytype_queries
        if pyisinstance(o, py)
            return jl
        end
    end
    return default
end

function convert(::Type{PyAny}, o::PyObject)
    if ispynull(o)
        return o
    end
    try
        T = pytype_query(o)
        if T == PyObject && is_pyjlwrap(o)
            return unsafe_pyjlwrap_to_objref(o)
        end
        convert(T, o)
    catch
        pyerr_clear() # just in case
        o
    end
end


text:
module Models

using MLJModelInterface
import MLJBase # needed for UnivariateFinite in ConstantClassifier

const MMI = MLJModelInterface

include("Constant.jl")
include("DecisionTree.jl")
include("NearestNeighbors.jl")
include("MultivariateStats.jl")
include("Transformers.jl")
include("foobarmodel.jl")
include("simple_composite_model.jl")

end


text:
macro S_str(terms)
    _lookup(terms)
end


text:
function segway_rk4(xt::AbstractVector, ut::AbstractVector, dt::Number, seg::Segway = Segway_unmutable())
    xt1 = copy(xt)
    segway_rk4!(xt1, ut, dt, seg)
    return xt1
end

function segway_rk4!(xt::AbstractVector, ut::AbstractVector, dt::Number, seg::Segway = Segway_unmutable())
    temp = similar(xt)

    #get partial results dy1-dy4
    dy1 = dxdt_segway(xt, ut, seg);
    temp .= dy1.*(dt/2.0) .+ xt;

    dy2 = dxdt_segway(temp, ut, seg);
    temp .= dy2.*(dt/2.0) .+ xt;

    dy3 = dxdt_segway(temp, ut, seg);
    temp .= dy3.*dt .+ xt;

    dy4 = dxdt_segway(temp, ut, seg);

    # update xt
    @. temp = dt*(dy1+2.0*(dy2+dy3)+dy4)/6.0

    #Check if segway has fallen #TODO limits anpassen
    #if -pi/2>=xt[4]+temp[4]
    #    xt[4] = -pi/2
	#	xt[5:7] = 0.0
    #elseif pi/2<=xt[4]+temp[4]
    #    xt[4] = pi/2
	#	xt[5:7] = 0.0
    #else
        xt .+= temp
    #end
    return nothing
end

"""
Segway dynamics, taken from Pathak2005.
State vector is [x0, y0, phi, alpha, dalpha, v, dphi]
"""
@inline function dxdt_segway(x::AbstractVector, u::AbstractVector, seg::Segway = Segway_unmutable())
    g = 9.81 #N/kg

    # Rename variables
    phi = x[3]
    alpha = x[4]
    dalpha = x[5]
    v = x[6]
    dphi = x[7]

    # Rigid Body Dynamics
    Dalpha = seg.Mb^2*cos(alpha)^2*seg.cz^2*seg.R^2+((-seg.Mb^2-2*seg.Mw*seg.Mb)*seg.cz^2-2*seg.Iyy*seg.Mw-seg.Iyy*seg.Mb)*seg.R^2-2*seg.Mb*seg.cz^2*seg.Iwa-2*seg.Iyy*seg.Iwa
    Galpha = (-seg.Mb*seg.cz^2+seg.Izz-seg.Ixx)*seg.R^2*cos(alpha)^2 + (seg.Mb*seg.cz^2+seg.Ixx+2*seg.Iwd+2*seg.b^2*seg.Mw)*seg.R^2+2*seg.b^2*seg.Iwa
    H = 1/2*seg.Mb*seg.R^2*seg.Izz+seg.Iwa*seg.Izz-seg.Mw*seg.R^2*seg.Ixx-seg.Iwa*seg.Ixx-seg.Mb*seg.cz^2*seg.Mw*seg.R^2 - seg.Mb*seg.cz^2*seg.Iwa -1/2*seg.Mb*seg.R^2*seg.Ixx+seg.Mw*seg.R^2*seg.Izz
    Kalpha = (-4*seg.Iyy*seg.Mb*seg.R^2*seg.cz-3*seg.R^2*seg.Mb^2*seg.cz^3+seg.Mb*seg.R^2*seg.cz*(seg.Ixx-seg.Izz))*sin(alpha) + (seg.Mb*seg.R^2*seg.cz*(seg.Ixx-seg.Izz)+seg.R^2*seg.Mb^2*seg.cz^3)*sin(3*alpha)
    f21 = sin(2*alpha)*dphi^2*H/Dalpha + seg.Mb^2*seg.cz^2*seg.R^2*sin(2*alpha)*dalpha^2/(2*Dalpha) + (-2*seg.Mb^2*seg.R^2*seg.cz-4*seg.Iwa*seg.Mb*seg.cz-4*seg.Mw*seg.R^2*seg.Mb*seg.cz)*g*sin(alpha)/(2*Dalpha)
    f22 = Kalpha*dphi^2+(seg.Mb^2*seg.cz^2*seg.R^2*g*sin(2*alpha))/(2*Dalpha)+(-4*seg.Iyy*seg.Mb*seg.R^2*seg.cz-4*seg.R^2*seg.Mb^2*seg.cz^3)*sin(alpha)*dalpha^2/(4*Dalpha)
    f23 = (-(seg.Ixx-seg.Izz)*seg.R^2-seg.Mb*seg.cz^2*seg.R^2)*sin(2*alpha)*dalpha*dphi/Galpha-sin(alpha)*seg.R^2*seg.Mb*seg.cz*v*dphi/Galpha
    g21 = (u[1] + u[2])*(seg.Mb*seg.R^2+2*seg.Mw*seg.R^2+2*seg.Iwa+seg.Mb*cos(alpha)*seg.cz*seg.R)/Dalpha
    g22 = -(u[1] + u[2])*seg.R*(seg.Mb*cos(alpha)*seg.cz*seg.R+seg.Iyy+seg.Mb*seg.cz^2)/Dalpha
    g23 = (u[1] - u[2])*seg.R*seg.b/Galpha

    xdot = similar(x)
    xdot[1] = cos(phi)*v
    xdot[2] = sin(phi)*v
    xdot[3] = x[7]
    xdot[4] = x[5]
    xdot[5] = f21 + g21
    xdot[6] = f22 + g22
    xdot[7] = f23 + g23

    return xdot
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
module AbsTypes

using DataFrames

export fit!, transform!, fit_transform!
export Machine, Learner, Transformer, Workflow, Computer

abstract type Machine end
abstract type Computer <: Machine end # does computation: learner and transformer
abstract type Workflow <: Machine end # types: Linear vs Combine
abstract type Learner <: Computer end
abstract type Transformer <: Computer end

# multiple dispatch for fit
function fit!(mc::Machine, input::DataFrame, output::Vector)
	error(typeof(mc),"not implemented")
end

function transform!(mc::Machine, input::DataFrame)
	error(typeof(mc),"not implemented")
end

# dynamic dispatch based Machine subtypes
function fit_transform!(mc::Machine, input::DataFrame, output::Vector=Vector())
	fit!(mc,input,output)
	transform!(mc,input)
end

end


text:
module Jitterplot
using CairoMakie
# A jitter function to randomly assign the x-axis positions for each x-parameter
function jitter(a, width_jit)
    a = Base.convert.(Float64, a)
    a .+= Base.rand(size(a)...) .* width_jit .- width_jit .* 0.5
    return a
end
# jitterplot! functions; kwargs are the same as the arguments in the scatter fuction of Makie.jl; dodge_group and dodge_color are the attributes for multiple dodged jitterplots; when dodge_color is assigned as one of the kwargs, the color argument in the kwargs does not apply.
function jitterplot!(x, y, width_jit :: Float64 = 0.2; kwargs...)
    if haskey(kwargs, :dodge_gap)
        dodge_gap = kwargs[:dodge_gap]
    else
        dodge_gap = 1.0
    end
    if (haskey(kwargs, :dodge_group) == true) & (haskey(kwargs, :dodge_color) == true)
        dodge_group = kwargs[:dodge_group]
        dodge_color = kwargs[:dodge_color]
        array0 = Array{Union{Float64, Int64}}(undef, 0 , 3)
        for (x, y, z) in zip(x, y, dodge_group)
            array1 = [x y z]
            array0 = vcat(array0, array1)
        end
        for (grp, grpl, dcolor) in zip(unique(dodge_group), 1:length(unique(dodge_group)), dodge_color)
            row_ind = findall(q -> (q .== grp), array0[:, 3])
            x = jitter((array0[row_ind, 1] .- 0.25 * dodge_gap) .+ 0.5 * dodge_gap * (grpl - 1)/(length(unique(dodge_group)) - 1), width_jit)
            y = array0[row_ind, 2]
            if haskey(kwargs, :direction)
                if kwargs[:direction] == :horizontal
                    Makie.scatter!(y, x, color = dcolor; kwargs...)
                else
                    Makie.scatter!(x, y, color = dcolor; kwargs...)
                end
            else
                Makie.scatter!(x, y, color = dcolor; kwargs...)
            end
        end
    elseif (haskey(kwargs, :dodge_group) == true) & (haskey(kwargs, :dodge_color) == false)
        dodge_group = kwargs[:dodge_group]
        array0 = Array{Union{Float64, Int64}}(undef, 0 , 3)
        for (x, y, z) in zip(x, y, dodge_group)
            array1 = [x y z]
            array0 = vcat(array0, array1)
        end
        for (grp, grpl) in zip(unique(dodge_group), 1:length(unique(dodge_group)))
            row_ind = findall(q -> (q .== grp), array0[:, 3])
            x = jitter((array0[row_ind, 1] .- 0.25 * dodge_gap) .+ 0.5 * dodge_gap * (grpl - 1)/(length(unique(dodge_group)) - 1), width_jit)
            y = array0[row_ind, 2]
            if haskey(kwargs, :direction)
                if kwargs[:direction] == :horizontal
                    Makie.scatter!(y, x; kwargs...)
                else
                    Makie.scatter!(x, y; kwargs...)
                end
            else
                Makie.scatter!(x, y; kwargs...)
            end
        end
    else
        x = jitter(x, width_jit)
        if haskey(kwargs, :direction)
            if kwargs[:direction] == :horizontal
                Makie.scatter!(y, x; kwargs...)
            else
                Makie.scatter!(x, y; kwargs...)
            end
        else
            Makie.scatter!(x, y; kwargs...)
        end
    end
end
# jitterplot function; not applicable for multiple dodged jitterplots.
function jitterplot(x, y, width_jit :: Float64; kwargs...)
    x = jitter(x, width_jit)
    if haskey(kwargs, :direction)
        if kwargs[:direction] == :horizontal
            Makie.scatter(y, x; kwargs...)
        else
            Makie.scatter(x, y; kwargs...)
        end
    else
        Makie.scatter(x, y; kwargs...)
    end
end
end


text:
@doc raw"""
    map2salm(map, spin, ℓmax)
    map2salm(map, plan)

Transform `map` values sampled on the sphere to ``{}_sa_{\ell, m}`` modes.

The `map` array should have size Nφ along its first dimension and Nϑ along its
second; any number of dimensions may follow.  The `spin` must be entered
explicitly, and `ℓmax` is the highest ℓ value you want in the output.

For repeated applications of this function with different values of `map`, it
is more efficient to pre-compute `plan` using [`plan_map2salm`](@ref).  These
functions will create a new `salm` array on each call.  To operate in place on
a pre-allocated `salm` array, use [`map2salm!`](@ref).

The core of this function follows the method described by [Reinecke and
Seljebotn](https://dx.doi.org/10.1051/0004-6361/201321494).

"""
function map2salm(map::AbstractArray{Complex{T}}, spin::Int, ℓmax::Int, show_progress=false) where {T<:Real}
    Nφ, Nϑ, Nextra... = size(map)
    salm = zeros(complex(T), (Ysize(ℓmax), Nextra...))
    map2salm!(salm, map, spin, ℓmax, show_progress)
    return salm
end

@doc raw"""
    map2salm!(salm, map, spin, ℓmax)
    map2salm!(salm, map, plan)

Transform `map` values sampled on the sphere to ``{}_sa_{\ell, m}`` modes in
place.

For details, see [`map2salm`](@ref).

"""
function map2salm!(
    salm::AbstractArray{Complex{T}},
    map::AbstractArray{Complex{T}},
    spin::Int, ℓmax::Int, show_progress=false
) where {T<:Real}
    plan = plan_map2salm(map, spin, ℓmax)
    map2salm!(salm, map, plan, show_progress)
end

"""
    plan_map2salm(map, spin, ℓmax)

Precompute values to use in executing [`map2salm`](@ref) or
[`map2salm!`](@ref).

The arguments to this function exactly mirror those of the first form of
[`map2salm`](@ref), and all but the first argument in the first form of
[`map2salm!`](@ref).  The `plan` returned by this function can be passed to the
second forms of those functions to avoid some computation and allocation costs.

Note that the `plan` object is not thread safe; a separate `plan` should be
created for each thread that will use one, or locks should be used to ensure
that a single `plan` is not used at the same time on different threads.

"""
function plan_map2salm(map_data::AbstractArray{Complex{T}}, spin::Int, ℓmax::Int) where {T<:Real}
    Nφ, Nϑ, Nextra... = size(map_data)
    Gs = [Array{complex(T)}(undef, (Nφ,)) for i = 1:nthreads()]
    m′max = abs(spin)
    Hwedge = Array{T}(undef, WignerHsize(ℓmax, m′max))
    H_rec_coeffs = H_recursion_coefficients(ℓmax, T)
    weight = clenshaw_curtis(Nϑ, T)
    expiθ = complex_powers(cis(π / T(Nϑ-1)), Nϑ-1)
    ϵs = SphericalFunctions.ϵ(-spin)
    extra_dims = Base.Iterators.product((1:e for e in Nextra)...)
    fftplan = T<:MachineFloat ? plan_fft(map_data[:, 1, first(extra_dims)...]) : nothing

    return (spin, ℓmax, Nφ, Nϑ, Nextra, Gs, m′max, Hwedge, H_rec_coeffs, weight, expiθ, ϵs, extra_dims, fftplan)
end

function computeG!(
    G::AbstractArray{Complex{T}},
    map::AbstractArray{Complex{T}},
    weight::T, fftplan
) where {T<:Real}
    G[:] = weight * fft(map)
end

function computeG!(
    G::AbstractArray{Complex{T}},
    map::AbstractArray{Complex{T}},
    weight::T, fftplan
) where {T<:MachineFloat}
    @views mul!(G[:], fftplan, map)
    @views G[:] *= weight
end

function map2salm!(
    salm::AbstractArray{Complex{T}}, map::AbstractArray{Complex{T}},
    (spin, ℓmax, Nφ, Nϑ, Nextra, Gs, m′max, Hwedge, H_rec_coeffs, weight, expiθ, ϵs, extra_dims, fftplan),
    show_progress=false
) where {T<:Real}
    s1 = size(salm)
    s2 = (Ysize(ℓmax), Nextra...)
    @assert s1==s2 "size(salm)=$s1  !=  (Ysize(ℓmax), Nextra...)=$s2"

    absspin = abs(spin)
    progress = Progress(Nϑ * prod(Nextra); showspeed=true, enabled=show_progress)

    @inbounds for ϑ ∈ 1:Nϑ
        H!(Hwedge, expiθ[ϑ], ℓmax, m′max, H_rec_coeffs, WignerHindex)
        @threads for extra ∈ collect(extra_dims)
            # NOTE: We can't thread at a higher level because each thread could access the
            # same element of `salm` simultaneously below; by threading at this level, we
            # are assured that the `extra...` index used below is unique to each thread.
            iₜ = Threads.threadid()
            G = Gs[iₜ]
            computeG!(G, map[:, ϑ, extra...], weight[ϑ], fftplan)
            for ℓ ∈ absspin:ℓmax
                λ_factor = ϵs * √((2ℓ+1)*T(π)) / Nφ

                i₀ = WignerHindex(ℓ, spin, 0, m′max)

                let m=0
                    salm[Yindex(ℓ, m), extra...] +=
                        G[m+1] * λ_factor * Hwedge[i₀]
                end

                i₊ = i₀
                i₋ = i₀
                if !signbit(spin)
                    for m ∈ 1:min(ℓ, absspin)
                        i₊ -= ℓ-m+2
                        i₋ += ℓ-m+1
                        salm[Yindex(ℓ, m), extra...] +=
                            G[m+1] * ϵ(m) * λ_factor * Hwedge[i₊]
                        salm[Yindex(ℓ, -m), extra...] +=
                            G[Nφ-m+1] * λ_factor * Hwedge[i₋]
                    end
                else
                    for m ∈ 1:min(ℓ, absspin)
                        i₊ += ℓ-m+1
                        i₋ -= ℓ-m+2
                        salm[Yindex(ℓ, m), extra...] +=
                            G[m+1] * ϵ(m) * λ_factor * Hwedge[i₊]
                        salm[Yindex(ℓ, -m), extra...] +=
                            G[Nφ-m+1] * λ_factor * Hwedge[i₋]
                    end
                end
                for m ∈ absspin+1:ℓ
                    i₊ += 1
                    i₋ += 1
                    salm[Yindex(ℓ, m), extra...] +=
                        G[m+1] * ϵ(m) * λ_factor * Hwedge[i₊]
                    salm[Yindex(ℓ, -m), extra...] +=
                        G[Nφ-m+1] * λ_factor * Hwedge[i₋]
                end
            end
            next!(progress)
        end
    end
end

function map2salm(map::AbstractArray{Complex{T}}, plan, show_progress=false) where {T<:Real}
    salm = zeros(complex(T), (Ysize(plan[2]), plan[5]...))
    map2salm!(salm, map, plan, show_progress)
    return salm
end


text:
#
# cobyla.jl --
#
# Julia interface to Mike Powell's COBYLA method.
#
# -----------------------------------------------------------------------------
#
# This file is part of OptimPackNextGen.jl which is licensed under the MIT
# "Expat" License:
#
# Copyright (C) 2015-2019, Éric Thiébaut
# <https://github.com/emmt/OptimPackNextGen.jl>.
#

module Cobyla

export
    cobyla,
    cobyla!

using Printf

import
    ..AbstractContext,
    ..AbstractStatus,
    ..getncalls,
    ..getradius,
    ..getreason,
    ..getstatus,
    ..grow!,
    ..iterate,
    ..restart

# The dynamic library implementing the method.
import ..libcobyla

# Status returned by most functions of the library.
struct Status <: AbstractStatus
    _code::Cint
end

# Possible status values returned by COBYLA.
const INITIAL_ITERATE      = Status( 2)
const ITERATE              = Status( 1)
const SUCCESS              = Status( 0)
const BAD_NVARS            = Status(-1)
const BAD_NCONS            = Status(-2)
const BAD_RHO_RANGE        = Status(-3)
const BAD_SCALING          = Status(-4)
const ROUNDING_ERRORS      = Status(-5)
const TOO_MANY_EVALUATIONS = Status(-6)
const BAD_ADDRESS          = Status(-7)
const CORRUPTED            = Status(-8)

# Get a textual explanation of the status returned by COBYLA.
function getreason(status::Status)
    ptr = ccall((:cobyla_reason, libcobyla), Ptr{UInt8}, (Cint,), status._code)
    if ptr == C_NULL
        error("unknown COBYLA status: ", status._code)
    end
    unsafe_string(ptr)
end

"""
# Minimizing a function of many variables subject to inequality constraints

Mike Powell's **COBYLA** algorithm attempts to find the variables `x` which
solve the problem:

    min f(x)    s.t.   c(x) <= 0

where `x` is a vector of variables that has `n ≥ 1` components, `f(x)` is an
objective function and `c(x)` implement `m` inequality constraints.  The
algorithm employs linear approximations to the objective and constraint
functions, the approximations being formed by linear interpolation at `n+1`
points in the space of the variables.  We regard these interpolation points as
vertices of a simplex.  The parameter `rho` controls the size of the simplex
and it is reduced automatically from `rhobeg` to `rhoend`.  For each `rho`,
COBYLA tries to achieve a good vector of variables for the current size, and
then `rho` is reduced until the value `rhoend` is reached.  Therefore `rhobeg`
and `rhoend` should be set to reasonable initial changes to and the required
accuracy in the variables respectively, but this accuracy should be viewed as a
subject for experimentation because it is not guaranteed.  The subroutine has
an advantage over many of its competitors, however, which is that it treats
each constraint individually when calculating a change to the variables,
instead of lumping the constraints together into a single penalty function.
The name of the subroutine is derived from the phrase "Constrained Optimization
BY Linear Approximations".

The in-place version of the algorithm is called as:

    Cobyla.minimize!(fc, x, m, rhobeg, rhoend) -> (status, x, fx)

where `x` is a vector with the initial and final variables, `m`, `rhobeg`, and
`rhoend` have been defined already, while `fc` is a Julia function which is
called as:

    fc(x, cx) -> fx

to store in `cx` the values of the constraints at `x` and to return `fx` the
value of the objective function at `x`.  If there are no constraints
(i.e. `m=0`), then `fc` is called without the `cx` argument as:

    fc(x) -> fx

The method returns a tuple with `status` the termination condition (should be
`Cobyla.SUCCESS` unless keyword `check` is set `false`, see below), `x` the
solution found by the algorithm and `fx` the corresponding function value.

The method:

    Cobyla.minimize(fc, x0, m, rhobeg, rhoend) -> (status, x, fx)

is identical but to `Cobyla.minimize!` does not modify the vector `x0` of
initial variables.

## Scaling of variables

The proper scaling of the variables is important for the success of the
algorithm and the optional `scale` keyword should be specified if the typical
precision is not the same for all variables.  If specified, `scale` is an array
of strictly nonnegative values and of same size as the variables `x`, such that
`scale[i]*rho` (with `rho` the trust region radius) is the size of the trust
region for the `i`-th variable.  If keyword `scale` is not specified, a unit
scaling for all the variables is assumed.

## Keywords

The following keywords are available:

* `scale` specifies the typical magnitudes of the variables.  If specified, it
  must have as many elements as `x`, all strictly positive.

* `check` (`true` by default) specifies whether to throw an exception if the
  algorithm is not fully successful.

* `verbose` (`0` by default) set the amount of printing.

* `maxeval` set the maximum number of calls to the objective function.  The
  default setting is `maxeval = 30n` with `n = length(x)` the number of
  variables.

## References

The algorithm is described in:

* M.J.D. Powell, "A direct search optimization method that models the objective
  and constraint functions by linear interpolation," in Advances in
  Optimization and Numerical Analysis Mathematics and Its Applications,
  vol. 275 (eds. Susana Gomez and Jean-Pierre Hennart), Kluwer Academic
  Publishers, pp. 51-67 (1994).

"""
minimize(args...; kwds...) = optimize(args...; maximize=false, kwds...)
minimize!(args...; kwds...) = optimize!(args...; maximize=false, kwds...)
@doc @doc(minimize) minimize!

"""

    Cobyla.maximize(fc, x0, m, rhobeg, rhoend) -> (status, x, fx)
    Cobyla.maximize!(fc, x, m, rhobeg, rhoend) -> (status, x, fx)

are similar to `Cobyla.minimize` and `Cobyla.minimize!` respectively but
solve the contrained maximization problem:

    max f(x)    s.t.   c(x) <= 0

"""
maximize(args...; kwds...) = optimize(args...; maximize=true, kwds...)
maximize!(args...; kwds...) = optimize!(args...; maximize=true, kwds...)
@doc @doc(maximize) maximize!

# `_wrklen(...)` yields the number of elements in COBYLA workspace.
_wrklen(n::Integer, m::Integer) = _wrklen(Int(n), Int(m))
_wrklen(n::Int, m::Int) = n*(3*n + 2*m + 11) + 4*m + 6

# `_work(...)` yields a large enough workspace for NEWUOA.
_work(::Type{T}, len::Integer) where {T} = Vector{T}(undef, len)

# Wrapper for the objective function in COBYLA, the actual objective
# function is provided by the client data as a `jl_value_t*` pointer.
function _objfun(n::Cptrdiff_t, m::Cptrdiff_t, xptr::Ptr{Cdouble},
                 _c::Ptr{Cdouble}, fptr::Ptr{Cvoid})::Cdouble
    x = unsafe_wrap(Array, xptr, n)
    f = unsafe_pointer_to_objref(fptr)
    return (m > 0 ? Cdouble(f(x, unsafe_wrap(Array, _c, m))) : Cdouble(f(x)))
end

# With precompilation, `__init__()` carries on initializations that must occur
# at runtime like `@cfunction` which returns a raw pointer.
const _objfun_c = Ref{Ptr{Cvoid}}()
function __init__()
    _objfun_c[] = @cfunction(_objfun, Cdouble,
                             (Cptrdiff_t, Cptrdiff_t, Ptr{Cdouble},
                              Ptr{Cdouble}, Ptr{Cvoid}))
end

"""
The methods:

    Cobyla.optimize(fc, x0, m, rhobeg, rhoend) -> (status, x, fx)
    Cobyla.optimize!(fc, x, m, rhobeg, rhoend) -> (status, x, fx)

are identical to `Cobyla.minimize` and `Cobyla.minimize!` respectively but have
an additional `maximize` keyword which is `false` by default and which
specifies whether to maximize the objective function; otherwise, the method
attempts to minimize the objective function.

"""
optimize(fc::Function, x0::AbstractVector{<:Real}, args...; kwds...) =
    optimize!(fc, copyto!(Array{Cdouble}(undef, length(x0)), x0),
              args...; kwds...)

function optimize!(fc::Function, x::DenseVector{Cdouble},
                   m::Integer, rhobeg::Real, rhoend::Real;
                   scale::DenseVector{Cdouble} = Array{Cdouble}(undef, 0),
                   maximize::Bool = false,
                   check::Bool = false,
                   verbose::Integer = 0,
                   maxeval::Integer = 30*length(x),
                   work::Vector{Cdouble} = _work(Cdouble, _wrklen(length(x), m)),
                   iact::Vector{Cptrdiff_t} = _work(Cptrdiff_t, m + 1))
    n = length(x)
    nscl = length(scale)
    if nscl == 0
        sclptr = Ptr{Cdouble}(0)
    elseif nscl == n
        sclptr = pointer(scale)
    else
        error("bad number of scaling factors")
    end
    grow!(work, _wrklen(n, m))
    grow!(iact, m + 1)
    status = Status(ccall((:cobyla_optimize, libcobyla), Cint,
                          (Cptrdiff_t, Cptrdiff_t, Cint, Ptr{Cvoid}, Any,
                           Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble,
                           Cptrdiff_t, Cptrdiff_t, Ptr{Cdouble},
                           Ptr{Cptrdiff_t}),
                          n, m, maximize, _objfun_c[], fc,
                          x, sclptr, rhobeg, rhoend, verbose, maxeval,
                          work, iact))
    if check && status != SUCCESS
        error(getreason(status))
    end
    return (status, x, work[1])
end

@doc @doc(optimize) optimize!

# Simpler version, mostly for testing.

function cobyla!(f::Function, x::DenseVector{Cdouble},
                 m::Integer, rhobeg::Real, rhoend::Real;
                 check::Bool = true,
                 verbose::Integer = 0,
                 maxeval::Integer = 30*length(x),
                 work::Vector{Cdouble} = _work(Cdouble, _wrklen(length(x), m)),
                 iact::Vector{Cptrdiff_t} = _work(Cptrdiff_t, m + 1))
    n = length(x)
    grow!(work, _wrklen(n, m))
    grow!(iact, m + 1)
    status = Status(ccall((:cobyla, libcobyla), Cint,
                          (Cptrdiff_t, Cptrdiff_t, Ptr{Cvoid}, Any,
                           Ptr{Cdouble}, Cdouble, Cdouble, Cptrdiff_t,
                           Cptrdiff_t, Ptr{Cdouble}, Ptr{Cptrdiff_t}),
                          n, m, _objfun_c[], f, x, rhobeg, rhoend,
                          verbose, maxeval, work, iact))
    if check && status != SUCCESS
        error(getreason(status))
    end
    return (status, x, work[1])
end

cobyla(f::Function, x0::DenseVector{Cdouble}, args...; kwds...) =
    cobyla!(f, copy(x0), args...; kwds...)

"""

```julia
using OptimPackNextGen.Powell
ctx = Cobyla.Context(n, m, rhobeg, rhoend; verbose=0, maxeval=500)
```

creates a new reverse communication workspace for COBYLA algorithm.  A typical
usage is:

```julia
x = Array{Cdouble}(undef, n)
c = Array{Cdouble}(undef, m)
x[...] = ... # initial solution
ctx = Cobyla.Context(n, m, rhobeg, rhoend, verbose=1, maxeval=500)
status = getstatus(ctx)
while status == Cobyla.ITERATE
    fx = ...       # compute function value at X
    c[...] = ...   # compute constraints at X
    status = iterate(ctx, fx, x, c)
end
if status != Cobyla.SUCCESS
    println("Something wrong occured in COBYLA: ", getreason(status))
end
```

""" Context

# Context for reverse communication variant of COBYLA.
# Must be mutable to be finalized.
mutable struct Context <: AbstractContext
    ptr::Ptr{Cvoid}
    n::Int
    m::Int
    rhobeg::Cdouble
    rhoend::Cdouble
    verbose::Int
    maxeval::Int
    function Context(n::Integer, m::Integer,
                     rhobeg::Real, rhoend::Real;
                     verbose::Integer=0, maxeval::Integer=500)
        n ≥ 2 || throw(ArgumentError("bad number of variables"))
        m ≥ 0 || throw(ArgumentError("bad number of constraints"))
        0 ≤ rhoend ≤ rhobeg ||
            throw(ArgumentError("bad trust region radius parameters"))
        ptr = ccall((:cobyla_create, libcobyla), Ptr{Cvoid},
                    (Cptrdiff_t, Cptrdiff_t, Cdouble, Cdouble,
                     Cptrdiff_t, Cptrdiff_t),
                    n, m, rhobeg, rhoend, verbose, maxeval)
        ptr != C_NULL || error(errno() == Base.Errno.ENOMEM
                               ? "insufficient memory"
                               : "unexpected error")
        return finalizer(ctx -> ccall((:cobyla_delete, libcobyla), Cvoid,
                                      (Ptr{Cvoid},), ctx.ptr),
                         new(ptr, n, m, rhobeg, rhoend, verbose, maxeval))
    end
end

@deprecate create(args...; kwds...) Context(args...; kwds...)

function iterate(ctx::Context, f::Real, x::DenseVector{Cdouble},
                 c::DenseVector{Cdouble})
    length(x) == ctx.n || error("bad number of variables")
    length(c) == ctx.m || error("bad number of constraints")
    Status(ccall((:cobyla_iterate, libcobyla), Cint,
                 (Ptr{Cvoid}, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
                 ctx.ptr, f, x, c))
end

function iterate(ctx::Context, f::Real, x::DenseVector{Cdouble})
    length(x) == ctx.n || error("bad number of variables")
    ctx.m == 0 || error("bad number of constraints")
    Status(ccall((:cobyla_iterate, libcobyla), Cint,
                 (Ptr{Cvoid}, Cdouble, Ptr{Cdouble}, Ptr{Cvoid}),
                 ctx.ptr, f, x, C_NULL))
end

restart(ctx::Context) =
    Status(ccall((:cobyla_restart, libcobyla), Cint, (Ptr{Cvoid},), ctx.ptr))

getstatus(ctx::Context) =
    Status(ccall((:cobyla_get_status, libcobyla), Cint, (Ptr{Cvoid},), ctx.ptr))

# Get the current number of function evaluations.  Result is -1 if
# something is wrong (e.g. CTX is NULL), nonnegative otherwise.
getncalls(ctx::Context) =
    Int(ccall((:cobyla_get_nevals, libcobyla), Cptrdiff_t, (Ptr{Cvoid},), ctx.ptr))

getradius(ctx::Context) =
    ccall((:cobyla_get_rho, libcobyla), Cdouble, (Ptr{Cvoid},), ctx.ptr)

getlastf(ctx::Context) =
    ccall((:cobyla_get_last_f, libcobyla), Cdouble, (Ptr{Cvoid},), ctx.ptr)

end # module Cobyla


text:
#
# Test the specific vector trnasport along implementations that do iterative transport,
# also the Schild and pole special cases
#
using ManifoldsBase, Test
import ManifoldsBase: parallel_transport_to!, parallel_transport_along!

struct NonDefaultEuclidean <: AbstractManifold{ManifoldsBase.ℝ} end
ManifoldsBase.log!(::NonDefaultEuclidean, v, x, y) = (v .= y .- x)
ManifoldsBase.exp!(::NonDefaultEuclidean, y, x, v) = (y .= x .+ v)
function ManifoldsBase.parallel_transport_to!(::NonDefaultEuclidean, Y, p, X, q)
    return copyto!(Y, X)
end
function ManifoldsBase.parallel_transport_along!(
    ::NonDefaultEuclidean,
    Y,
    p,
    X,
    c::AbstractVector,
)
    return copyto!(Y, X)
end

@testset "vector_transport_along" begin
    M = NonDefaultEuclidean()
    types = [Vector{Float64}, Vector{Float32}]
    for T in types
        @testset "Type $T" begin
            pts = convert.(Ref(T), [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])
            X2 = log(M, pts[1], pts[3])
            c = [0.5 * (pts[1] + pts[2]), pts[2], 0.5 * (pts[2] + pts[3]), pts[3]]
            @test vector_transport_along(M, pts[1], X2, c, SchildsLadderTransport()) == X2
            @test vector_transport_along(M, pts[1], X2, c, PoleLadderTransport()) == X2
            @test vector_transport_along(M, pts[1], X2, c, ParallelTransport()) == X2
            @test vector_transport_along(M, pts[1], X2, [], SchildsLadderTransport()) == X2
            @test vector_transport_along(M, pts[1], X2, [], PoleLadderTransport()) == X2
            @test vector_transport_along(M, pts[1], X2, [], ParallelTransport()) == X2
            # check mutating ones with defaults
            p = allocate(pts[1])
            ManifoldsBase.pole_ladder!(M, p, pts[1], pts[2], pts[3])
            # -log_p3 p == log_p1 p2
            @test isapprox(M, -log(M, pts[3], p), log(M, pts[1], pts[2]))
            ManifoldsBase.schilds_ladder!(M, p, pts[1], pts[2], pts[3])
            @test isapprox(M, log(M, pts[3], p), log(M, pts[1], pts[2]))
        end
    end
end

@testset "vector-transport fallback types" begin
    VT = VectorTransportDirection()
    M = NonDefaultEuclidean()
    p = [1.0, 0.0, 0.0]
    q = [0.0, 1.0, 0.0]
    X = [0.0, 0.0, 1.0]
    Y = similar(X)
    @test vector_transport_direction(M, p, X, p - q, VT) == X
    @test vector_transport_direction!(M, Y, p, X, p - q, VT) == X
    VT2 = VectorTransportTo()
    @test vector_transport_to(M, p, X, q, VT2) == X
    @test vector_transport_to!(M, Y, p, X, q, VT2) == X
end


text:
using MaxLFSR
using Test
using Random

# Set the seed for consistent trials
Random.seed!(123)

function _fill!(A, itr)
    A .= false
    count = 0
    for i in itr
        A[i] = true
        count += 1
    end
    return count
end

@testset "MaxLFSR.jl" begin
    # Test `nbits`
    @test MaxLFSR.nbits(1) == 4
    @test MaxLFSR.nbits(2) == 4
    @test MaxLFSR.nbits(15) == 4
    @test MaxLFSR.nbits(16) == 5
    @test MaxLFSR.nbits(32) == 6

    @test_throws ArgumentError MaxLFSR.nbits(0)
    @test_throws ArgumentError MaxLFSR.nbits(typemax(UInt64))

    # Try messing with the `seed` argument
    @test_throws ArgumentError LFSR(10; seed = 0)
    @test_throws ArgumentError LFSR(10; seed = 20)

    # Test the 1 length of 1 works
    itr = LFSR(1)
    @test collect(itr) == [1]

    # Some more argument handling erros
    @test_throws ArgumentError LFSR(0)
    @test_throws ArgumentError LFSR(typemax(UInt))

    # Iterate through powers of 2 from 2 to 30 with random seeds.
    # Then generate some random lengths in the same rangs.
    nloops = 3
    for pow in 2:30
        println("Testing Length: 2 ^ $pow")
        len = 2^pow
        A = falses(len)

        for i in 1:nloops
            if i == 1
                seed = 1
            else
                seed = rand(1:len)
            end

            # Clear the tracking array
            itr = LFSR(len; seed = seed)
            runtime = @elapsed count = _fill!(A, itr)
            if i == nloops
                @show runtime
            end

            @test count == length(itr) == len
            @test all(isequal(true), A)
        end
    end

    # Test 10000 different lengths between 2 and 2^20
    A = falses(2^20)
    total = 10000
    for i in 1:total
        if iszero(mod(i, div(total, 10)))
            println("Performing Fuzz Test $i of $total")
        end

        len = rand(1:2^20)
        seed = rand(1:len)
        resize!(A, len)

        # Normal LFSR
        itr = LFSR(len; seed = seed)
        count = _fill!(A, itr)

        @test count == length(itr) == len
        @test all(isequal(true), A)
    end

    #####
    ##### FastLFSR
    #####

    @test_throws ArgumentError MaxLFSR.FastLFSR(32)
    @test_throws ArgumentError MaxLFSR.FastLFSR(7)

    # Iterate through powers of 2 from 2 to 30 with random seeds.
    # Then generate some random lengths in the same rangs.
    nloops = 3
    for pow in 4:30
        println("Testing Length: 2 ^ $pow")
        len = 2^pow - 1
        A = falses(len)

        for i in 1:nloops
            if i == 1
                seed = 1
            else
                seed = rand(1:len)
            end

            # Clear the tracking array
            itr = FastLFSR(len; seed = seed)
            runtime = @elapsed count = _fill!(A, itr)

            # Report the time for the last iteration.
            if i == nloops
                @show runtime
            end

            @test count == length(itr) == len
            @test all(isequal(true), A)
        end
    end

end


text:
function Prew(protocol::Float64,poke::Int64)
    protocol*(ℯ^(-(poke-1)/5))
end

function Prew(protocol,poke::AbstractVector{Int64})
    Prew.(protocol,poke)
end

function Prew(protocol,poke::UnitRange{Int64})
    Prew.(protocol,collect(poke))
end

function Prew(poke::UnitRange{Int64})
    df = DataFrame(Poke = Int64[], Prew = Float64[], Protocol = String[])
    for protocol in [1.0,0.75,0.5]
        dd = DataFrame(Poke = collect(poke), Prew = Prew.(protocol,collect(poke)))
        dd[!,:Protocol] .= string(protocol)
        append!(df,dd)
    end
    return df
end
"""
    Poutcome(outcome, poke, protocol)

Compute the probability of a positive or negative outcome at the nth poke
for a given protocol
"""
function Poutcome(outcome::Bool,poke::Int64,protocol::Float64)
    if outcome
        return Prew(protocol,poke)
    elseif !outcome
        return 1 - Prew(protocol,poke)
    end
end

"""
    Pobservation(observations,protocol)

Compute the probability of a specific sequence of observation for a given
protocol
"""
function Pobservations(observations::AbstractVector{Bool},protocol::Float64)
    cumprod([Poutcome(o,p,protocol) for (p,o) in enumerate(observations)])
end

function Pobservations_lastval(observations::AbstractVector{Bool},protocol::Float64)
    cumprod([Poutcome(o,p,protocol) for (p,o) in enumerate(observations)])[end]
end

"""
    Pprotocol(protocol,observations)

Compute the probability of the current trial being equal to protocol
given a series of (observations). Compute the probability of a series of
observation to be genrated by the specified scaling factor (protocol)
given the 3 possible protocol type (0.5,0.75,1.0)

"""
function Pprotocol(protocol::Float64,observations::AbstractVector{Bool})
    possibleprotocols = Pobservations(observations,0.5) +
        Pobservations(observations,0.75) +
        Pobservations(observations,1.0)
    Pobservations(observations,protocol)./possibleprotocols
end

function Pprotocol_lastval(protocol::Float64,observations::AbstractVector{Bool})
    possibleprotocols = Pobservations(observations,0.5) +
        Pobservations(observations,0.75) +
        Pobservations(observations,1.0)
    Pobservations(observations,protocol)/possibleprotocols
end

"""
    Pnext(observations)

Compute the probability of the next poke to be rewarded after a series of
observations and given the 3 possible protocol type (0.5,0.75,1.0)

"""
function Pnext(observations::AbstractVector{Bool})
    res = zeros(length(observations))
    for prot in [0.5,0.75,1.00]
        res += Prew(prot,2:length(observations)+1) .* Pprotocol(prot,observations)
    end
    return res
end

function Pnext_lastval(observations::AbstractVector{Bool})
    res = 0
    for prot in [0.5,0.75,1.00]
        res += Prew(prot,length(observations)+1) * Pprotocol(prot,observations)
    end
    return res
end


text:
# calculate residuals for straight line fits (of R2 vs B0²)
line(B0², β, c) = @. β*B0² + c
function resid(B0², R2, β, c)
    y = Measurements.value.(R2)
    σ = Measurements.uncertainty.(R2)
    yhat = line(B0², β, c)
    return vec((y - yhat) ./ σ)
end

# define shortcuts for passing parameter vectors
line(x, p) = line(x, p[1], p[2])
resid(B0², R, p) = resid(B0², R, p[1], p[2])

# calculate χ² score using all ZQ, DQ, DQ' and QQ measurements
# performs least squares fitting to determine intercepts for each line
function calcχ2(ξH, ξC, resi)
    βZQ = Measurements.value((ξCcsa[resi]-ξHcsa[resi])^2 + (ξC-ξH)^2)
    βDQ = Measurements.value((ξCcsa[resi]+ξHcsa[resi])^2 + (ξC+ξH)^2)
    βDQQ = Measurements.value((ξCcsa[resi]-3ξHcsa[resi])^2 + (ξC-3ξH)^2)
    βQQ = Measurements.value((ξCcsa[resi]+3ξHcsa[resi])^2 + (ξC+3ξH)^2)

    c0 = [10.0]
    χ2 = 0.0
    fit = LsqFit.lmfit(c -> resid(B0², zq[resi,:], βZQ, c), c0, Float64[])
    χ2 += sum(fit.resid.^2)
    fit = LsqFit.lmfit(c -> resid(B0², dq[resi,:], βDQ, c), c0, Float64[])
    χ2 += sum(fit.resid.^2)
    fit = LsqFit.lmfit(c -> resid(B0², dqq[resi,:], βDQQ, c), c0, Float64[])
    χ2 += sum(fit.resid.^2)
    fit = LsqFit.lmfit(c -> resid(B0², qq[resi,:], βQQ, c), c0, Float64[])
    χ2 += sum(fit.resid.^2)
end

function fitβ(B0², R2)
    p0 = [.1, 10] # slope, intercept
    fit = LsqFit.lmfit(p -> resid(B0², R2, p), p0, Float64[])
    β = fit.param[1] ± stderror(fit)[1]
end

max_stars_repo_name: chriswaudby/qq-hahnecho

max_issues_repo_name: chriswaudby/qq-hahnecho

max_forks_repo_name: chriswaudby/qq-hahnecho


text:
version = v"11.0.1"

include("../common.jl")

# Build the tarballs
build_tarballs(ARGS, configure(version; experimental=true)...;
               preferred_gcc_version=v"6.1.0", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Copyright(c) Rafael Dias
# MIT Licensed
include("analise-estruturas.jl")

# Ordem deve ser horária
# coordenadas = [
#     PropriedadesSecoes.Coordenada(53f0, 0f0)
#     PropriedadesSecoes.Coordenada(53f0, 43f0)
#     PropriedadesSecoes.Coordenada(0f0, 43f0)
#     PropriedadesSecoes.Coordenada(0f0, 60f0)
#     PropriedadesSecoes.Coordenada(70f0, 60f0)
#     PropriedadesSecoes.Coordenada(70f0, 0f0)
# ]

coordenadas = [
    PropriedadesSecoes.Coordenada(0f0, 0f0)
    PropriedadesSecoes.Coordenada(0f0, 2f0)
    PropriedadesSecoes.Coordenada(5f0, 2f0)
    PropriedadesSecoes.Coordenada(5f0, 12f0)
    PropriedadesSecoes.Coordenada(0f0, 12f0)
    PropriedadesSecoes.Coordenada(0f0, 14f0)
    PropriedadesSecoes.Coordenada(12f0, 14f0)
    PropriedadesSecoes.Coordenada(12f0, 12f0)
    PropriedadesSecoes.Coordenada(7f0, 12f0)
    PropriedadesSecoes.Coordenada(7f0, 2f0)
    PropriedadesSecoes.Coordenada(12f0, 2f0)
    PropriedadesSecoes.Coordenada(12f0, 0f0)
]

@show PropriedadesSecoes.areaPoligono(coordenadas)
@show PropriedadesSecoes.centroGravidadePoligonoX(coordenadas)
@show PropriedadesSecoes.centroGravidadePoligonoY(coordenadas)
@show PropriedadesSecoes.perimetroPoligono(coordenadas)
@show PropriedadesSecoes.comprimentoEntrePontos(coordenadas)

elementos = [AnaliseEstruturas.buildElemento(0f0, 0f0, 1f0, 1f0) AnaliseEstruturas.buildElemento(1f0, 1f0, 1f0, 3f0)]
estrutura = AnaliseEstruturas.Estrutura(elementos)

print("\n\nMétodo dos deslocamentos")
for el::AnaliseEstruturas.Elemento in elementos
    AnaliseEstruturas.metodoDeslocamentos(el)
end

print("\nMétodo das forças")
for el::AnaliseEstruturas.Elemento in elementos
    AnaliseEstruturas.metodoForcas(el)
end

max_stars_repo_name: rsd-consultorias/scripts-engenharia

max_issues_repo_name: rsd-consultorias/scripts-engenharia

max_forks_repo_name: rsd-consultorias/scripts-engenharia