text:
## container of points into vectors n vectors of length N
## N points, each of size n

## Lesson learned -- this is a very bad idea!
## better to handle the T a different way
evec(T,n) =  Tuple(T[] for _ in 1:n)
evec(T, N, n) =  Tuple(Vector{T}(undef, N) for _ in 1:n)

## julia> @btime xs_ys1(vs) setup=(vs=[randn(1000) for i in 1:3]);
##   83.308 μs (1013 allocations: 172.67 KiB)

## julia> @btime xs_ys2(vs) setup=(vs=[randn(1000) for i in 1:3]);
##   222.371 μs (2016 allocations: 180.72 KiB)

## julia> @btime xs_ys3(vs) setup=(vs=[randn(1000) for i in 1:3]);
##   1.003 ms (1019 allocations: 165.20 KiB)

## julia> @btime xs_ys4(vs) setup=(vs=[randn(1000) for i in 1:3]);
##   1.115 ms (5474 allocations: 210.95 KiB)

## julia> @btime xs_ys5(vs) setup=(vs=[randn(1000) for i in 1:3]);
##   1.120 ms (5474 allocations: 210.95 KiB)

## julia> @btime xs_ys6(vs) setup=(vs=[randn(1000) for i in 1:3]);
##   76.604 μs (1008 allocations: 164.63 KiB)

## julia> @btime xs_ys7(vs) setup=(vs=[randn(1000) for i in 1:3]);
##   74.306 μs (1008 allocations: 164.63 KiB)

## julia> @btime xs_ys8(vs) setup=(vs=[randn(1000) for i in 1:3]);
##   36.098 μs (2006 allocations: 94.25 KiB)

## julia> @btime xs_ys9(vs) setup=(vs=[randn(1000) for i in 1:3]);
##   85.732 μs (3006 allocations: 203.63 KiB)

## ....

    ## THE WINNER, but we would use one with keywords
## julia> @btime xs_ys13a(vs) setup=(vs=[randn(1000) for i in 1:3]);
##   62.768 μs (1003 allocations: 117.28 KiB)

## julia> @btime xs_ys13akw(vs) setup=(vs=[randn(1000) for i in 1:3]);
##   65.905 μs (1003 allocations: 117.28 KiB)

## make a matrix n x N, then go down 1:n
function xs_ys1(vs)
    A=hcat(vs...)
    Tuple([A[i,:] for i in eachindex(first(vs))])
end

## broadcast push!
function xs_ys2(vs)
    u = first(vs); N = length(vs)
    T = eltype(u); n = length(u)
    v0 = evec(T,n)
    for v in vs
        push!.(v0, v)
    end
    v0
end

## broadcast push!
function xs_ys2a(vs)
    u = first(vs); N = length(vs)
    n = length(u)
    v0 = Tuple(eltype(u)[] for _ in eachindex(u))
    for v in vs
        push!.(v0, v)
    end
    v0
end

## broadcast setindex!
function xs_ys3(vs)
  u = first(vs); N = length(vs)
    T = eltype(u); n = length(u)
    v0 = evec(T,N,n)
    for (i,v) in enumerate(vs)
        setindex!.(v0, v, i)
    end

    v0
end

## 10 times faster ~77mus avoiding passing T
function xs_ys3a(vs)
  u = first(vs); N = length(vs)
    n = length(u)
    v0 = Tuple(Vector{eltype(u)}(undef, N) for _ in eachindex(u))
    for (i,v) in enumerate(vs)
        setindex!.(v0, v, i)
    end

    v0
end

function xs_ys3b(vs)
  u = first(vs); N = length(vs)
    n = length(u)
    v0 = ntuple(_ -> Vector{eltype(u)}(undef, N), n)
    for (i,v) in enumerate(vs)
        setindex!.(v0, v, i)
    end

    v0
end

## loop N n
function xs_ys4(vs)
    u = first(vs); N = length(vs)
    T = eltype(u); n = length(u)
    v0 = evec(T,N,n)

    for i in 1:N
        for j in 1:n
            v0[j][i] = vs[i][j]
        end
    end
    v0
end

## loop N n
function xs_ys4a(vs)
    u = first(vs); N = length(vs)
    T = eltype(u); n = length(u)
    v0 = evec(T,N,n)

    for (i,v) in enumerate(vs)
        for j in 1:n
            v0[j][i] = v[j]
        end
    end
    v0
end

## fast 67mus
function xs_ys4b(vs)
    u = first(vs); N = length(vs)
    n = length(u)
    v0 = Tuple(Vector{eltype(u)}(undef, N) for _ in eachindex(u))

    for (i,v) in enumerate(vs)
        for j in 1:n
            v0[j][i] = v[j]
        end
    end
    v0
end

## loop n N
function xs_ys5(vs)
    u = first(vs); N = length(vs)
    T = eltype(u); n = length(u)
    v0 = evec(T,N,n)

    for j in 1:n
        for i in 1:N
            v0[j][i] = vs[i][j]
        end
    end
    v0
end

function xs_ys6(vs)
    u = first(vs); N = length(vs)
    T = eltype(u); n = length(u)
    A = Matrix{T}(undef, (N,n))
    for (i,v) in enumerate(vs)
        A[i,:] = v
    end
    Tuple(A[:,i] for i in 1:n)
end

function xs_ys7(vs)
    u = first(vs); N = length(vs)
    T = eltype(u); n = length(u)
    A = Matrix{T}(undef, (n, N))
    for (i,v) in enumerate(vs)
        A[:,i] = v
    end
    Tuple(A[i, :] for i in 1:n)
end

# faster but doesn't wotk with plot recipes
# and may be slower once realized
function xs_ys8(vs)
    N = length(vs)
    u = first(vs); T = eltype(u); n = length(u)
    Tuple((vs[j][i] for j in 1:N) for i in 1:n)
end

function xs_ys9(vs)
    N = length(vs)
    u = first(vs); T = eltype(u); n = length(u)
    Tuple(collect(vs[j][i] for j in 1:N) for i in 1:n)
end

function xs_ys10(vs)
    N = length(vs)
    u = first(vs); T = eltype(u); n = length(u)
    v0 = evec(T,N, n)
    for j in 1:n
        v0[j][:] .= (v[j] for v in vs)
    end
    v0
end

# mauro3 https://github.com/JuliaDiffEq/ODE.jl/issues/80
_pluck(y,i) =  eltype(first(y))[el[i] for el in y]
xs_ys11(vs) = Tuple(_pluck(vs, i) for i in eachindex(first(vs)))

# slower
xs_ys11a(vs) = ntuple(i->_pluck(vs, i), length(first(vs)))

# one liner
xs_ys11b(vs) = Tuple(eltype(first(vs))[el[i] for el in vs]  for i in eachindex(first(vs)))

function xs_ys11c(vs)
    u = first(vs)
    Tuple(eltype(u)[el[i] for el in vs] for i in eachindex(u))
end
xs_ys11d(vs) = (u=first(vs);  Tuple(eltype(u)[el[i] for el in vs] for i in eachindex(u)))
xs_ys11e(vs) = (u=first(vs);  ntuple(i->eltype(u)[v[i] for v in vs], length(u)))
xs_ys11f(vs) = (u=first(vs);n::Int=length(u);T::DataType=eltype(u);ntuple(i->eltype(u)[v[i] for v in vs], n))
xs_ys11g(vs::Vector{Vector{T}}) where {T} = (u=first(vs);n::Int=length(u);ntuple(i->T[v[i] for v in vs], n))

@inline _pluck(T, y, i) =  T[el[i] for el in y]
function xs_ys11b(vs)
    T = eltype(first(vs))
    Tuple(_pluck(T, vs, i) for i in eachindex(first(vs)))
end

function xs_ys12(vs)
    N = length(vs)
    u = first(vs); T = eltype(u); n = length(u)
    Tuple(T[el[i] for el in vs] for i in eachindex(first(vs)))
    end

function xs_ys12a(vs)
    N = length(vs)
    u = first(vs); T = eltype(u); n = length(u)
    ntuple( i -> T[el[i] for el in vs], n)
end

function xs_ys11h(vs)
    u = first(vs)
    T = eltype(u)
    Tuple(T[el[i] for el in vs] for i in eachindex(u))
    end

function xs_ys11i(vs)
    u = first(vs)
    Tuple(eltype(u)[el[i] for el in vs] for i in eachindex(u))
end

function _xs_ys12(vs, u::Vector{T}) where {T}
    Tuple(T[el[i] for el in vs] for i in eachindex(u))
end

xs_ys13(vs, u::Vector{T}=first(vs)) where {T} = Tuple(T[el[i] for el in vs] for i in eachindex(u))

xs_ys13a(vs, u::Vector{T}=first(vs), n::Val{N}=Val(length(u))) where {T,N} = ntuple(i -> T[el[i] for el in vs], n)

## cleaned up
function xs_ys13a(vs, u::Vector{T}=first(vs), n::Val{N}=Val(length(u))) where {T,N}
   plucki = i -> T[el[i] for el in vs]
   ntuple(plucki, n)
end

function xs_ys13akw(vs; u::Vector{T}=first(vs), n::Val{N}=Val(length(u))) where {T,N}
   plucki = i -> T[el[i] for el in vs]
   ntuple(plucki, n)
end

function xs_ys13b(vs, u::Vector{T}=first(vs), n::Val{N}=Val(length(u))) where {T,N}
    Tuple(T[el[i] for el in vs] for i in eachindex(u))
end

xs_ys14(vs) = Tuple(eltype(vs[1])[vs[i][j] for i in 1:length(vs)] for j in 1:length(vs[1]))

xs_ys14a(vs) = Tuple([vs[i][j] for i in 1:length(vs)] for j in 1:length(first(vs)))


text:
export PresenceUpdate,
    TypingStart,
    UserUpdate

"""
Sent when a [`User`](@ref)'s [`Presence`](@ref) is updated.
"""
struct PresenceUpdate <: AbstractEvent
    presence::Presence
end
@boilerplate PresenceUpdate :docs
PresenceUpdate(d::Dict{String, Any}) = PresenceUpdate(Presence(d))

"""
Sent when a [`User`](@ref) begins typing.
"""
struct TypingStart <: AbstractEvent
    channel_id::Snowflake
    guild_id::Union{Snowflake, Missing}
    user_id::Snowflake
    timestamp::DateTime
end
@boilerplate TypingStart :dict :docs

"""
Sent when a [`User`](@ref)'s details are updated.
"""
struct UserUpdate <: AbstractEvent
    user::User
end
@boilerplate UserUpdate :docs
UserUpdate(d::Dict{String, Any}) = UserUpdate(User(d))

ext: jl


text:
include("required.jl")

## Get data
@load joinpath("data", "jld2", "comparison-results.jld2") raw sdm
results = CSV.read(joinpath("data", "proc", "comparison-results.csv"), DataFrame)

# Check correlation
cor(results.richness_raw, results.richness_sdm)
cor(results.lcbd_raw, results.lcbd_sdm)

# Check distribution
histogram(results.richness_raw; title="richness_raw", legend=:none)
histogram(results.richness_sdm; title="richness_sdm", legend=:none)
histogram(results.lcbd_raw; title="lcbd_raw", legend=:none)
histogram(results.lcbd_sdm; title="lcbd_sdm", legend=:none)

mean(results.richness_raw)
std(results.richness_raw)
variation(results.richness_raw)
var(results.richness_raw) / mean(results.richness_raw)

mean(results.richness_sdm)
std(results.richness_sdm)
variation(results.richness_sdm)

## Test linear regression in Julia
using GLM

# LM
lm_richness = lm(@formula(richness_sdm ~ richness_raw), results)
lm_lcbd = lm(@formula(lcbd_sdm ~ lcbd_raw), results)

# Test utility functions
coeftable(lm_richness)
coef(lm_richness)
deviance(lm_richness)
dof_residual(lm_richness)
r2(lm_richness)
stderror(lm_richness)
vcov(lm_richness)
predict(lm_richness)

function glance(model::StatsModels.TableRegressionModel{<:LinearModel})
    return DataFrame(;
        r_squared=r2(model),
        adj_r_squared=adjr2(model),
        # sigma(model),
        sigma=missing,
        # ftest(model),
        statistic=missing,
        # pvalue(model),
        pvalue=missing,
        df=dof(model),
        loglik=loglikelihood(model),
        AIC=aic(model),
        BIC=bic(model),
        deviance=deviance(model),
        df_residual=dof_residual(model),
        nobs=nobs(model),
    )
end
glance(lm_richness)
glance(lm_lcbd)

# Test GLMs
# Poisson
glm_richness = glm(@formula(richness_sdm ~ richness_raw), results, Poisson())
deviance(glm_richness) / dof_residual(glm_richness) # Overdispersion

function glance(model::StatsModels.TableRegressionModel{<:GeneralizedLinearModel})
    return glancedf = DataFrame(;
        # null_deviance = nulldeviance(model) # not working
        null_deviance=missing,
        df_null=dof_residual(model),
        loglik=loglikelihood(model),
        AIC=aic(model),
        BIC=bic(model),
        deviance=deviance(model),
        df_residual=nobs(model) - dof(model),
        nobs=nobs(model),
    )
end
glance(glm_richness)

# Negative binomial
negb_richness = negbin(@formula(richness_sdm ~ richness_raw), results, LogLink())
glance(negb_richness)

# LCBD Gamma
glm_lcbd = glm(@formula(lcbd_sdm ~ lcbd_raw), results, Gamma())
glance(glm_lcbd)

# residuals(glm_richness)

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
mutable struct MetaGraph{T <: Integer,U <: Real} <: AbstractMetaGraph{T}
    graph::SimpleGraph{T}
    vprops::Dict{T,PropDict}
    eprops::Dict{SimpleEdge{T},PropDict}
    gprops::PropDict
    weightfield::Symbol
    defaultweight::U
    metaindex::MetaDict
    indices::Set{Symbol}
end

function MetaGraph(x, weightfield::Symbol, defaultweight::U) where U <: Real
    T = eltype(x)
    g = SimpleGraph(x)
    vprops = Dict{T,PropDict}()
    eprops = Dict{SimpleEdge{T},PropDict}()
    gprops = PropDict()
    metaindex = MetaDict()
    idxs = Set{Symbol}()
    MetaGraph(g, vprops, eprops, gprops, weightfield, defaultweight, metaindex, idxs)
end

MetaGraph() = MetaGraph(SimpleGraph())
MetaGraph{T,U}() where T <: Integer where U <: Real = MetaGraph(SimpleGraph{T}(), one(U))
MetaGraph{T,U}(x::Integer) where T <: Integer where U <: Real = MetaGraph(T(x), :weight, U(1.0))
MetaGraph(x) = MetaGraph(x, :weight, 1.0)
MetaGraph(x, weightfield::Symbol) = MetaGraph(x, weightfield, 1.0)
MetaGraph(x, defaultweight::Real) = MetaGraph(x, :weight, defaultweight)

# converts MetaGraph{Int, Float64} to MetaGraph{UInt8, Float32}
function MetaGraph{T,U}(g::MetaGraph) where T <: Integer where U <: Real
    newg = SimpleGraph{T}(g.graph)
    return MetaGraph(newg, U(g.defaultweight))
end

function MetaGraph{T,U}(g::SimpleGraph, weightfield::Symbol=:weight, defaultweight::Real=1.0) where T <: Integer where U <: Real
    newg = SimpleGraph{T}(g)
    return MetaGraph(newg, weightfield, U(defaultweight))
end

function MetaGraph{T,U}(g::SimpleGraph, defaultweight::Real) where T <: Integer where U <: Real
    newg = SimpleGraph{T}(g)
    return MetaGraph(newg, :weight, U(defaultweight))
end

function MetaGraph(g::MetaDiGraph{T,U}) where T <: Integer where U <: Real
    return MetaGraph(Graph(g.graph), deepcopy(g.vprops), deepcopy(g.eprops), deepcopy(g.gprops), g.weightfield, g.defaultweight, deepcopy(g.metaindex), deepcopy(g.indices))
end

SimpleGraph(g::MetaGraph) = g.graph

is_directed(::Type{MetaGraph}) = false
is_directed(::Type{MetaGraph{T,U}}) where T where U = false
is_directed(g::MetaGraph) = false

weighttype(g::MetaGraph{T,U}) where T where U = U
function props(g::MetaGraph, _e::SimpleEdge)
    e = LightGraphs.is_ordered(_e) ? _e : reverse(_e)
    get(g.eprops, e, PropDict())
end

function set_props!(g::MetaGraph, _e::SimpleEdge, d::Dict)
    e = LightGraphs.is_ordered(_e) ? _e : reverse(_e)
    if has_edge(g, e)
        if !_hasdict(g, e)
            g.eprops[e] = d
        else
            merge!(g.eprops[e], d)
        end
        return true
    end
    return false
end

zero(g::MetaGraph{T,U}) where T where U = MetaGraph{T,U}(SimpleGraph{T}())


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
include("../common.jl")

gap_version = v"400.1191.001"
gap_lib_version = v"400.1191.000"
name = "orb"
upstream_version = v"4.8.4" # when you increment this, reset offset to v"0.0.0"
offset = v"0.0.0" # increment this when rebuilding with unchanged upstream_version, e.g. gap_version changes
version = offset_version(upstream_version, offset)

# Collection of sources required to build libsingular-julia
sources = [
    ArchiveSource("https://github.com/gap-packages/$(name)/releases/download/v$(upstream_version)/$(name)-$(upstream_version).tar.gz",
                  "5d44d52cf7ae7ab6ad1208d6ed563f3dc73df10e31a71024a1f9b7ba7dd71427"),
]

# Bash recipe for building across all platforms
script = raw"""
cd orb*
./configure ${prefix}/share/gap
make -j${nproc}

# copy the loadable module
mkdir -p ${prefix}/lib/gap
cp bin/*/*.so ${prefix}/lib/gap/

install_license LICENSE
"""

name = gap_pkg_name(name)
platforms, dependencies = setup_gap_package(gap_version, gap_lib_version)

# The products that we will ensure are always built
products = [
    FileProduct("lib/gap/orb.so", :orb),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", preferred_gcc_version=v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Autogenerated wrapper script for SLEEF_jll for x86_64-w64-mingw32
export libsleef, libsleefquad

using CompilerSupportLibraries_jll
using MPFR_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("SLEEF")
JLLWrappers.@declare_library_product(libsleef, "libsleef.dll")
JLLWrappers.@declare_library_product(libsleefquad, "libsleefquad.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, MPFR_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libsleef,
        "bin\\libsleef.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsleefquad,
        "bin\\libsleefquad.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
include("basis/basis.jl")
include("linear.jl")

export LinearFunction

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
import Base.length
mutable struct StateGrid{D}
    cutPoints::Vector{Vector{Float64}}
    indicies::Vector{Int}

    function StateGrid{D}(cutPoints...) where D
        newCutPoints = Vector{Vector{Float64}}(undef, length(cutPoints))
        for i = 1:D
            if length(Set(cutPoints[i])) != length(cutPoints[i])
                error(@sprintf("Duplicates cutpoints are not allowed (duplicates observed in dimension %d)",i))
            end
            if !issorted(cutPoints[i])
                error("Cut points must be sorted")
            end
            newCutPoints[i] = cutPoints[i]
        end
        indicies = cumprod([1; [length(points)+1 for points in newCutPoints[1:end-1]]...])
        return new(newCutPoints, indicies)
    end
end
StateGrid(cutPoints...) = StateGrid{Base.length(cutPoints)}(cutPoints...)
Base.length(grid::StateGrid{D}) where D = D
Base.size(grid::StateGrid{D}) where D = NTuple{D, Int}(length(points)+1 for points in grid.cutPoints)

zeros_like(grid::StateGrid{D}) where D = zeros(Int, NTuple{D, Int}(length(points)+1 for points in grid.cutPoints))

function access(grid::StateGrid{D}, access_cnt::Array{Int,D}, s::S, pomdp::POMDP{S}) where {S,D}
    s = convert(SVector{D,Float64}, s)
    index = 1
    for d in 1:D
        cutPoints = grid.cutPoints[d]
        # Binary search for the apt grid
        start_ind = 1
        end_ind = length(cutPoints) + 1
        mid_ind = div(start_ind+end_ind, 2)
        while start_ind < end_ind
            cutPoint = cutPoints[mid_ind]
            if s[d] < cutPoint
                end_ind = mid_ind
            else
                start_ind = mid_ind + 1
            end
            mid_ind = div(start_ind+end_ind, 2)
        end
        index += (mid_ind-1) * grid.indicies[d]
    end
    access_cnt[index] += 1
    return access_cnt[index] == 1 ? true : false
end


text:
module TestLAE

using QuantEcon
using Base.Test
using FactCheck
using Distributions

# copied from the lae lecture
s = 0.2
δ = 0.1
a_σ = 0.4       # A = exp(B) where B ~ N(0, a_sigma)
α = 0.4         # We set f(k) = k**alpha
ϕ = LogNormal(a_σ)

function p(x, y)
    d = s * x.^α

    pdf_arg = clamp((y .- (1-δ) .* x) ./ d, eps(), Inf)
    return pdf(ϕ, pdf_arg) ./ d
end

# other data
n_a, n_b, n_y = 50, (5, 5), 20
a = rand(n_a) + 0.01
b = rand(n_b...) + 0.01

y = linspace(0, 10, 20)

lae_a = LAE(p, a)
lae_b = LAE(p, b)

laes = [lae_a, lae_b]

facts("Testing lae.jl") do
    # test stuff
    for l in laes
        # make sure X is made a column vector by constructor
        @fact size(l.X, 2) --> 1

        # make sure X is indeed a column vector (2 dims)
        @fact ndims(l.X) --> 2

        # Make sure we get 1d estimate out
        @fact size(lae_est(l, y)) --> (n_y, )
    end
end  # facts
end  # module


text:
setCableLengthFeet(rhd::FPGA,port,lengthInFeet::Float64)=setCableLengthMeters(rhd,port,.3048*lengthInFeet)

function setCableLengthMeters(rhd::FPGA,port, lengthInMeters::Float64)

    tStep=1.0 / (2800.0 * rhd.sampleRate)

    distance = 2.0 * lengthInMeters

    timeDelay = (distance / cableVelocity) + xilinxLvdsOutputDelay + rhd2000Delay + xilinxLvdsInputDelay + misoSettleTime

    delay = convert(Int32,floor(((timeDelay / tStep) + 1.0) +0.5))

    if delay < 1
        delay=1
    end

    setCableDelay(rhd,port, delay)
    nothing
end

function estimateCableLengthMeters(fpga,delay)

    tStep=1.0 / (2800.0 * fpga.sampleRate)
    cableVelocity = 0.555 * speedOfLight
    distance = ((delay - 1.0) * tStep - (xilinxLvdsOutputDelay + rhd2000Delay + xilinxLvdsInputDelay + misoSettleTime))
    if distance <0.0
        distance = 0.0
    end
    distance
end

function approxCableLengthFeet(fpga,delay)
    estimateCableLengthMeters(fpga,delay) * 3.2808
end

function setCableDelay(rhd::FPGA,port, delay)

    #error checking goes here

    if delay<0
        delay=0
    elseif delay>15
        delay=15
    end

    #here i should update the bit shift int and cableDelay vector of ints appropriately. I have no idea what the cableDelay vector does

    if port=="PortA"
        bitShift=0;
    elseif port=="PortB"
        bitShift=4
    elseif port=="PortC"
        bitShift=8
    elseif port=="PortD"
        bitShift=12
    end

    bitShift=convert(Int32, bitShift)

    SetWireInValue(rhd,WireInMisoDelay, delay << bitShift, 0x0000000f << bitShift)
    UpdateWireIns(rhd)

    nothing
end

function check_delay_output(fpga::FPGA,port,output)

    data_stream_inds=find(fpga.dataStreamEnabled.==1)
    if port=="PortA"
        outinds=find((data_stream_inds.==1).|(data_stream_inds.==2).|(data_stream_inds.==9).|(data_stream_inds.==10))
    elseif port=="PortB"
        outinds=find((data_stream_inds.==3).|(data_stream_inds.==4).|(data_stream_inds.==11).|(data_stream_inds.==12))
    elseif port=="PortC"
        outinds=find((data_stream_inds.==5).|(data_stream_inds.==6).|(data_stream_inds.==13).|(data_stream_inds.==14))
    elseif port=="PortD"
        outinds=find((data_stream_inds.==7).|(data_stream_inds.==8).|(data_stream_inds.==15).|(data_stream_inds.==16))
    end

    hits=0
    for j in outinds
        hits+=output[33,j]==UInt16("I"[1])
        hits+=output[34,j]==UInt16("N"[1])
        hits+=output[35,j]==UInt16("T"[1])
        hits+=output[36,j]==UInt16("A"[1])
        hits+=output[37,j]==UInt16("N"[1])
        hits+=output[25,j]==UInt16("R"[1])
        hits+=output[26,j]==UInt16("H"[1])
        hits+=output[27,j]==UInt16("D"[1])
    end

    hits==length(outinds)*8
end

function determine_delay(fpga::FPGA,port)

    setMaxTimeStep(fpga,SAMPLES_PER_DATA_BLOCK)
    setContinuousRunMode(fpga,false)
    selectAuxCommandBank(fpga,port,"AuxCmd3",0)

    output_delay=falses(16)

    for delay=0:15

        setCableDelay(fpga,port,delay)

	flushBoard(fpga)

        runBoard(fpga)

        while isRunning(fpga)
        end

	if fpga.usb3
	    SetWireInValue(fpga,WireInResetRun, 1 << 16, 1 << 16)
	    UpdateWireIns(fpga)
	    ReadFromBlockPipeOut(fpga,PipeOutData,2*convert(Clong,fpga.numWords),fpga.usbBuffer)
	    SetWireInValue(fpga,WireInResetRun, 0 << 16, 1 << 16);
	    UpdateWireIns(fpga)
	else
	    ReadFromPipeOut(fpga,PipeOutData, 2*convert(Clong, fpga.numWords * 1), fpga.usbBuffer)
	end
        index=1
        output=zeros(UInt16,60,fpga.numDataStreams)
        for t=1:60

            index+=12

            for j=1:3
                for i=1:fpga.numDataStreams
                    if j==3
                        output[t,i]=convertUsbWordu(fpga.usbBuffer,index)
                    end
                    index+=2
                end
            end

            #Amplifier

	    for i=1:32
	        for j=1:fpga.numDataStreams
		    index+=2
	        end
	    end

	    #skip 36 filler word
        if fpga.numDataStreams>0
	       index += 2
       end

	    #ADCs
            for i=1:8
                index+=2
            end

	    #TTL
	    index += 4
        end

        if check_delay_output(fpga,port,output)
            output_delay[delay+1]=true
        end
    end

    if all(.!output_delay)
        println("No delay setting produces optimum results")
	setCableLengthFeet(fpga,port, 6.0)
    else
        setCableDelay(fpga,port,find(output_delay.==true)[1]-1)
        println("Optimum delay on ", port, " is ", find(output_delay.==true)[1]-1)
        println("Approx ", approxCableLengthFeet(fpga,find(output_delay.==true)[1]-1), " feet")
    end

    nothing
end


text:
# ir_dump.jl
# 2019-03-05 Jeff Fessler

using Printf

"""
`ir_dump(x::DataType)`

show all the filelds of a structure more nicely than dump() does
"""
function ir_dump(y::Any)
	x = typeof(y)
	print(x)
	fields = fieldnames(x)
	fieldtypes = x.types
	for idx in 1:length(fields)
		println()
		fd = fields[idx]
		print(stdout, " ", fd, "::")
		ft = fieldtypes[idx]
		print(ft)
		if ft <: Number
			print(stdout, " ", getfield(y, fd))
		end
		if ft == String
			print(stdout, " '", getfield(y, fd), "'")
		end
	end
	println()
	nothing
end

"""
`ir_dump(:test)`
"""
function ir_dump(test::Symbol)
	@assert test == :test
	x = (a=1, b=2)
	ir_dump(x)
	true
end


text:
@testset "Fit2B" begin include("test_fit2b.jl") end
@testset "Fit3B" begin include("test_fit3b.jl") end


text:
"Your optional docstring here"
function distance(a, b)
    distance = 0
    if sizeof(a) == sizeof(b)
        for (i,j) in zip(a,b)
            if i != j
                distance += 1
            end
        end
        return distance
    else
        throw(ArgumentError("input strands are not of the same size"))
    end

end

max_stars_repo_name: AbstractGeek/exercism-expeditions

max_issues_repo_name: AbstractGeek/exercism-expeditions

max_forks_repo_name: AbstractGeek/exercism-expeditions


text:
A = [4 5 31 1 2 5 3252 23 1 5236 1]

#Done
function insertionSort!(x::Array)
    for j in 2:length(x)
        key = x[j]
        i = j-1
        while i>0 && x[i]< key
            x[i+1] = x[i]
            i = i-1
        end
        x[i+1] = key
    end
end

#Done
function mergeSort!(x::Array)
    mergeSortStep!(x,1, length(x))
end

function mergeSortStep!(x::Array, p::Int, r::Int)
    if p < r
        q = Int(floor((p+r)/2))
        mergeSortStep!(x,p,q)
        mergeSortStep!(x,q+1,r)
        mergeParts!(x,p,q,r)
    end
end

function mergeParts!(x::Array, p::Int, q::Int, r::Int)
    L = copy(x[p:q])
    R = copy(x[q+1:r])
    push!(L, typemin(Int64))
    push!(R, typemin(Int64))
    i = j = 1
    for k = p:r
        if L[i] >= R[j]
            x[k] = L[i];
            i = i+1
        else
            x[k] = R[j]
            j = j+1
        end
    end
end

#Done
function selectionSort!(x)
    n = length(x)
    for i in 1:n-1
        minI = i
        for j in i+1:n
            if x[j] > x[minI]
                minI = j
            end
        end
        if i!= minI
            x[i], x[minI] = x[minI], x[i]
        end
    end
end

#Done
function quickSort!(x)
    quickSortStep!(x,1,length(x))
end

function quickSortStep!(x,p,r)
    if p<r
        q = partition!(x,p,r)
        quickSortStep!(x,p,q-1)
        quickSortStep!(x,q+1,r)
    end
end

function partition!(x,p,r)
    pivot = x[r]
    i = p - 1
    for j in p:r-1
        if x[j] >= pivot
            i = i+1
            x[i],x[j] = x[j],x[i]
        end
    end
    x[i+1], x[r] = x[r], x[i+1]
    return i+1
end

selectionSort!(A)
print(A)


text:
# This file is part of FlagSets.jl project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

using FlagSets
using Test, Serialization

using FlagSets: deprecated_syntax_message, undef_var_error_hint
using FlagSets: invalid_bit_error, invalid_flagspec_error, invalid_type_error
using FlagSets: overflow_error, not_unique_error
using FlagSets: BitMask

macro macrocall(ex)
    @assert Meta.isexpr(ex, :macrocall)
    ex.head = :call
    for i = 2:length(ex.args)
        ex.args[i] = QuoteNode(ex.args[i])
    end
    insert!(ex.args, 3, __module__)
    return esc(ex)
end

if VERSION >= v"1.6"
    @testset "FlagSets of symbols" verbose=true begin
        include("symbol_flagsets.jl")
    end

    @testset "FlagSets of other types" verbose=true begin
        include("flagset_type.jl")
    end
else
    @testset "FlagSets of symbols" begin
        include("symbol_flagsets.jl")
    end

    @testset "FlagSets of other types" begin
        include("flagset_type.jl")
    end
end


text:
"""
A structure that is assumed to hold a homogenous transform.
"""
struct basic_dh{T}
    mat::SMatrix{4,4,T,16}
    function basic_dh(t::SVector{3,T}) where {T}
        mat = SMatrix{4,4,T,16}( one(T), zero(T), zero(T), zero(T),
                                zero(T),  one(T), zero(T), zero(T),
                                zero(T), zero(T),  one(T), zero(T),
                                   t[1],    t[2],    t[3],  one(T))
        return new{T}(mat)
    end
    function basic_dh(R::Rotation{3,T}) where {T}
        R = R[:]
        mat = SMatrix{4,4,T,16}(   R[1],    R[2],    R[3], zero(T),
                                   R[4],    R[5],    R[6], zero(T),
                                   R[7],    R[8],    R[9], zero(T),
                                zero(T), zero(T), zero(T),  one(T))
        return new{T}(mat)
    end
    function basic_dh(R::Rotation{3,T}, t::SVector{3,T}) where {T}
        R = R[:]
        mat = SMatrix{4,4,T,16}(   R[1],    R[2],    R[3], zero(T),
                                   R[4],    R[5],    R[6], zero(T),
                                   R[7],    R[8],    R[9], zero(T),
                                   t[1],    t[2],    t[3],  one(T))
        return new{T}(mat)
    end
    function basic_dh(R::SMatrix{3,3,T,9}, t::SVector{3,T}) where {T}
        R = R[:]
        mat = SMatrix{4,4,T,16}(   R[1],    R[2],    R[3], zero(T),
                                   R[4],    R[5],    R[6], zero(T),
                                   R[7],    R[8],    R[9], zero(T),
                                   t[1],    t[2],    t[3],  one(T))
        return new{T}(mat)
    end
    function basic_dh(R::SMatrix{3,3,T,9}) where {T}
        R = R[:]
        mat = SMatrix{4,4,T,16}(   R[1],    R[2],    R[3], zero(T),
                                   R[4],    R[5],    R[6], zero(T),
                                   R[7],    R[8],    R[9], zero(T),
                                zero(T), zero(T), zero(T), one(T))
        return new{T}(mat)
    end
    # function basic_dh(;scale::SVector{3,T}) where {T}
    #     return basic_dh(SMatrix{3,3,T,9}(scale[1], zero(T),  zero(T),
    #                                      zero(T),  scale[2], zero(T),
    #                                      zero(T),  zero(T),  scale[3]))
    # end
    function basic_dh(d::Diagonal{T,SVector{3,T}}) where {T}
        return basic_dh(SMatrix{3,3,T,9}(d))
    end
    basic_dh(s::T) where {T} = basic_dh(s * one(SMatrix{3,3,T,9}))
    basic_dh(mat::SMatrix{4,4,T,16}) where {T} = new{T}(mat)
end

"""
Extracts the rotational and translational part of the transformation matrix.
"""
function dh_R_t(a::basic_dh{T}) where {T}
    a_top = getTop(a.mat)
    t = SVector{3,T}(a_top[10], a_top[11], a_top[12])
    R = SMatrix{3,3,T,9}(a_top[1], a_top[2], a_top[3], a_top[4], a_top[5], a_top[6], a_top[7], a_top[8], a_top[9])
    return R, t
end

function povray_12(a::basic_dh)
    # http://www.povray.org/documentation/view/3.6.0/49/
    # qx = Val00 * px + Val10 * py + Val20 * pz + Val30
    # qy = Val01 * px + Val11 * py + Val21 * pz + Val31
    # qz = Val02 * px + Val12 * py + Val22 * pz + Val32
    #
    # matrix <Val00, Val01, Val02, Val10, Val11, Val12, Val20, Val21, Val22, Val30, Val31, Val32>
    return getTop(a.mat)[:]
end

Base.:*(a::basic_dh, b::basic_dh) = basic_dh(a.mat * b.mat)

function Base.inv(a::basic_dh)
    R, t = dh_R_t(a)
    return basic_dh(R', -R' * t)
end

Base.one(::Type{basic_dh{T}}) where {T} = basic_dh(one(SMatrix{4,4,T,16}))

dh_vector_mul(dh::basic_dh{T1}, v::SVector{3,T2}) where {T1,T2} = unPad(dh.mat * onePad(v))


text:
#!/usr/bin/env julia6
using ReadGlobal

function VectorDirec_part!(w₁::A,w₂::A,w₃::A,out::A) where {A<:Vector{<:Number}}
    Threads.@threads for i in 1:length(w₁)
        @inbounds out[i] = w₃[i]/sqrt((w₁[i]^2 + w₂[i]^2 + w₃[i]^2))
    end
    return 0
end

function VectorDirec(w1::String,w2::String,w3::String,out::String)
    doinchunks(VectorDirec_part!,input=(w1,w2,w3),output=(out,))
end

VectorDirec(ARGS...)


text:
const vfname = "puppersim.mp4"      # Video file name

##################################################### functions
function finish_recording(s::mjSim)
    # Primarily see avio.jl reference in render function below. Also of potential interest:
    # https://github.com/JuliaIO/VideoIO.jl/tree/master/examples
    # https://discourse.julialang.org/t/creating-a-video-from-a-stack-of-images/646/7

    @static if use_VideoIO
        println("Saving video to: $vfname")
        props = [:priv_data => ("crf"=>"22","preset"=>"medium")]
        @time encodedvideopath = VideoIO.encodevideo(vfname, s.imgstack, framerate=30, AVCodecContextProperties=props, silent=false)
        s.imgstack = []
    else
        println("Closing $vfname")
        close(s.record)
    end

    println("Done writing video!")
    s.record = nothing
end

function record_video(s::mjSim)
    if s.record === nothing
        println("Recording")
        @static if use_VideoIO
            s.record = 0
        else
            println("Saving video to $vfname")
            w, h = GLFW.GetFramebufferSize(s.window)

            # -y overwrite output files
            # -f force format
            @ffmpeg_env s.record = open(`ffmpeg -y
                            -f rawvideo -pixel_format rgb24
                            -video_size $(w)x$(h) -framerate $(s.refreshrate)
                            -i pipe:0
                            -preset fast -threads 0
                            -vf "vflip" $vfname`, "w")
        end
    else
        finish_recording(s)
    end
end

function record_video_frame(s::mjSim, width, height)
    @static if use_VideoIO
    if s.record <= max_video_frames
        s.record += 1

        # Image dims must be a multiple of two
        width  = div(width,  2) * 2 # ensure that width is even
        height = div(height, 2) * 2 # ensure that height is even
        buflen = width * height * sizeof(TPixel)

        # If user has resized the window, we may need to allocate a new buffer
        if length(s.vidbuf) != buflen
            s.vidbuf = Vector{UInt8}(undef, buflen)
        end

        # Get the pixels from MuJoCo
        viewrect = mjrRect(0, 0, width, height)
        mjr_readPixels(s.vidbuf, C_NULL, viewrect, s.con)

        # Reference: @testset "Encoding video across all supported colortypes" block in file avio.jl:
        # (https://github.com/JuliaIO/VideoIO.jl/blob/master/test/avio.jl)

        # Reinterpret the video buffer as pixels
        pixels = reinterpret(TPixel, s.vidbuf)

        # Shape the buffer into an image array
        image = reshape(pixels, width, height)

        # Allocate an uninitialized frame on the image stack
        push!(s.imgstack, Array{TPixel,2}(undef, height, width))

        # Permute image array from column major to row major and write the
        # result to the uninitialized memory at the top of the image stack
        permutedims!(s.imgstack[end], image, (2,1))

        # Flip the image in place on the image stack in the vertical direction
        vflip!(s.imgstack[end])
    else    # s.record <= max_video_frames
        finish_recording(s)
    end     # s.record <= max_video_frames
    else    # @static if use_VideoIO
        buflen = width * height * sizeof(TPixel)

        # If user has resized the window, we may need to allocate a new buffer
        if length(s.vidbuf) != buflen
            s.vidbuf = Vector{UInt8}(undef, buflen)
        end

        # Get the pixels from MuJoCo
        viewrect = mjrRect(0, 0, width, height)
        mjr_readPixels(s.vidbuf, C_NULL, viewrect, s.con)
        write(s.record, s.vidbuf)
    end # @static if use_VideoIO
end


text:
include("xml_parameters.jl")

seed_time = @elapsed rand(100, 100) * rand(100, 100)
seed = Int(round(seed_time * 1e10))
srand(seed)

L = 8
W = L
betas = 1./collect(0.1:0.1:10)
betas = [5, 10, 15, 20]
betas = [20]
len = length(betas)

for beta in betas
   for Delta in 0.1:0.1:7.0
       mkpath("jobs/L_$(L)_W_$(W)/L_$(L)_W_$(W)_beta_$(beta)/L_$(L)_W_$(W)_beta_$(beta)_Delta_$(Delta)")
       cd("jobs/L_$(L)_W_$(W)/L_$(L)_W_$(W)_beta_$(beta)/L_$(L)_W_$(W)_beta_$(beta)_Delta_$(Delta)")

       for h in -13.0:0.1:13.0
            prefix = "L_$(L)_W_$(W)_beta_$(beta)_Delta_$(Delta)_h_$(h)"
            try
                mkdir(prefix)
    	    catch end
            cd(prefix)

            p = Dict{Any, Any}()
            p["L"] = L
            p["W"] = W
            p["Delta"] = Delta
            p["EPS"] = 0.1
            p["h"] = h
            p["BETA"] = [beta]
            p["THERMALIZATION"] = 2^13
            p["SWEEPS"] = 2^14
            p["SEED"] = convert(Int, round(rand() * 1e4))

            parameters2xml(p, prefix)
            cd("..")
        end
    end
end

max_stars_repo_name: pebroecker/CythonHardcoreBosons

max_issues_repo_name: pebroecker/CythonHardcoreBosons

max_forks_repo_name: pebroecker/CythonHardcoreBosons


text:
module BenchSimulation

using BenchmarkTools
using LatticeBoltzmann

import LatticeBoltzmann: collide!, stream!, apply_boundary_conditions!, next!, simulate

suite = BenchmarkGroup()

function initialize_model(q = D2Q9(), τ = 1.0, scale = 2)
    cs = 1 / q.speed_of_sound_squared
    benchmark_problem = LatticeBoltzmann.CouetteFlow(cs * (τ - 0.5), scale)
    return LatticeBoltzmann.LatticeBoltzmannModel(
        benchmark_problem,
        q,
        process_method = LatticeBoltzmann.ProcessingMethod(benchmark_problem, false, 0),
    )
end
for q in LatticeBoltzmann.Quadratures
    suite[string(q)] = BenchmarkGroup([string(q), "simulation"])

    model = initialize_model(q)

    suite[string(q)]["collision"] =
        @benchmarkable(collide!(model, time = 0.0), setup = (model = initialize_model($q)))
    suite[string(q)]["stream"] =
        @benchmarkable(stream!(model), setup = (model = initialize_model($q)))
    suite[string(q)]["boundary conditions"] = @benchmarkable(
        apply_boundary_conditions!(model, time = 0.0),
        setup = (model = initialize_model($q))
    )
    suite[string(q)]["processing"] =
        @benchmarkable(next!(model, 0), setup = (model = initialize_model($q)))
    suite[string(q)]["simulate"] =
        @benchmarkable(simulate(model, 0:10), setup = (model = initialize_model($q)))
end

end  # module
BenchSimulation.suite


text:
export PRISMA, PrismaParams, PrismaStepsize

function PRISMA(x, # starting point
				Lf, # Lipshitz constant of f
				grad_f,
				prox_g,
				prox_h,
				obj,
				params::OptParams=PrismaParams(PrismaStepsize(), 10))

	# initialize
	ssr      = params.stepsizerule
	if ssr.initial_stepsize == Inf
		R = sqrt(obj(x)) # estimate of distance to solution
		ssr.initial_stepsize = 2*R/Lf
	end
	y        = copy(x)
	xkpp     = copy(x)
	betakpp  = step(ssr)
	Lkpp     = Lf + 1/betakpp
	thetak   = 1
	thetakpp = 1 # XXX see orabona's code

	if params.verbose > 0
		@printf("%10s%12s\n", "iter", "obj")
		@printf("%10d%12.4e\n", 0, obj(xkpp))
	end

	# iterate
	for k=1:params.maxiter
		betak, betakpp    = betakpp, step(ssr)
		Lk, Lkpp          = Lkpp, Lf + 1/betakpp
		thetak, thetakpp  = thetakpp, 2/(1+sqrt(1+4Lkpp/thetak^2/Lk))
		xk, xkpp          = xkpp, prox_h((1-1/Lk/betak)*y - grad_f(y)/Lk + prox_g(y, betak)/Lk/betak, 1/Lk)
		y                 = xkpp + thetakpp*(1/thetak - 1)*(xkpp - xk)
		if params.verbose > 0
			@printf("%10d%12.4e\n", k, obj(xkpp))
		end
		if stop(params, xk, xkpp)
			break
		end
	end

	return xkpp
end

type PrismaStepsize<:StepSizeRule
	initial_stepsize::AbstractFloat
	iteration_counter::Int
end
function step(s::PrismaStepsize, args...; kwargs...)
	s.iteration_counter += 1
	return s.initial_stepsize/(s.iteration_counter)
end
PrismaStepsize() = PrismaStepsize(Inf)
PrismaStepsize(initial_stepsize) = PrismaStepsize(initial_stepsize, 0)

type PrismaParams<:OptParams
	stepsizerule::StepSizeRule
	maxiter::Int
	verbose::Int
	reltol::Float64
end
PrismaParams(;stepsize::StepSizeRule=PrismaStepsize(),
	          maxiter::Int=100,
	          verbose::Int=1,
	          reltol::Float64=1e-5) = PrismaParams(stepsize,maxiter,verbose,reltol)
stop(s::PrismaParams, xk, xkpp) = vecnorm(xk - xkpp)/vecnorm(xk) < s.reltol ? true : false

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
module TestReturns

using Base.Test
using TimeData
using Econometrics
using DataFrames
using Dates
using DataArrays

###########
## ranks ##
###########

x = [0.4, 0.8, 0.2, 0.1]
@test Econometrics.ranks(x) == [0.6, 0.8, 0.4, 0.2]

xx = [x x]
@test Econometrics.ranks(xx) == [Econometrics.ranks(x)
                                 Econometrics.ranks(x)]

## DataArrays
##-----------

da = @data([3, 4, NA, 2.3, 1 ])
expOut = @data([0.6, 0.8, NA, 0.4, 0.2])
@test isequal(ranks(da), daOut)

## Timenum
##--------

xIn = testcase(Timenum, 4)
da1 = @data([1:5]/6)
da2 = @data([0.2, 0.4, NA, 0.6, 0.8])
expOutDf = DataFrame(prices1 = da1, prices2 = da2)

actOut = ranks(xIn)
@test isequal(actOut.vals, expOutDf)
@test isequal(names(actOut), names(xIn))
@test idx(actOut) == idx(xIn)

## Timematr
##---------

xIn = testcase(Timematr, 1)
da1 = @data([0.2, 0.6, 0.4, 0.8])
da2 = @data([0.4, 0.6, 0.2, 0.8])
expOutDf = DataFrame(prices1 = da1, prices2 = da2)

actOut = ranks(xIn)
@test isequal(actOut.vals, expOutDf)
@test isequal(names(actOut), names(xIn))
@test idx(actOut) == idx(xIn)

end


text:
# v0.12 deprecations
@deprecate Dropout(p, dims) Dropout(p; dims=dims)
@deprecate InstanceNorm(λ, β, γ, μ, σ², ϵ, momentum) InstanceNorm(λ, β, γ, μ, σ², ϵ, momentum, nothing)
@deprecate BatchNorm(λ, β, γ, μ, σ², ϵ, momentum) BatchNorm(λ, β, γ, μ, σ², ϵ, momentum, nothing)
@deprecate GroupNorm(G, λ, β, γ, μ, σ², ϵ, momentum) GroupNorm(G, λ, β, γ, μ, σ², ϵ, momentum, nothing)
@deprecate outdims(f, inputsize) outputsize(f, inputsize)
@deprecate Conv(; weight,  bias, activation=identity, kws...) Conv(weight, bias, activation; kws...)
@deprecate ConvTranspose(; weight, bias, activation=identity, kws...) ConvTranspose(weight, bias, activation; kws...)
@deprecate DepthwiseConv(; weight, bias, activation=identity, kws...) DepthwiseConv(weight, bias, activation; kws...)


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

mutable struct TreeJointPrismatic <: Modia3D.AbstractPrismatic
   _internal::ModiaMath.ComponentInternal  # Data common to assembly component classes
   obj1::Object3D                        # Translation of obj1 to obj2 along axis of obj1
   obj2::Object3D
   canCollide::Bool                        # = true, if collision detection is enabled between the two Object3Ds connected by the Prismatic joint (default = false)
   isDriven::Bool                          # = true, if setDistance!(..) can be called on the joint

   # Definitions special to Prismatic joints
   axis::Int                       #  axis (= 1,2,3 or -1,-2,-3) of translation
   eAxis::ModiaMath.Vector3D   # Unit vector in direction of axis

   s::ModiaMath.RealScalar
   v::ModiaMath.RealScalar
   a::ModiaMath.RealScalar
   f::ModiaMath.RealScalar
   residue::ModiaMath.RealScalar

   P::ModiaMath.RealScalar       # Total power flowing in obj1, obj2, axis computed at communication points (must be zero; = computed with interpolated values of x,derx)

   flange::PrismaticFlange

   TreeJointPrismatic(_internal, obj1, obj2, canCollide, isDriven, axis, eAxis) = new(_internal, obj1, obj2, canCollide, isDriven, axis, eAxis)
end

mutable struct CutJointPrismatic <: Modia3D.AbstractPrismatic
   _internal::ModiaMath.ComponentInternal  # Data common to assembly component classes
   obj1::Object3D                                # Translation of obj1 to obj2 along axis of obj1
   obj2::Object3D

   # Definitions special to Prismatic joints
   visited::Bool                                   # = true, if already visited when traversing the spanning tree

   # Definitions special to Prismatic joints
   axis::Int                          # axis (= 1,2,3, -1,-2,-3) of translation
   eAxis::ModiaMath.Vector3D      # Unit vector in direction of axis
   e_lambda1::ModiaMath.Vector3D  # Unit vector in direction of lambda1
   e_lambda2::ModiaMath.Vector3D  # Unit vector in direction of lambda2

   lambda1::ModiaMath.RealScalar
   lambda2::ModiaMath.RealScalar
   mue1::ModiaMath.RealScalar
   mue2::ModiaMath.RealScalar
   residue1::ModiaMath.RealScalar
   residue2::ModiaMath.RealScalar
   residue1d::ModiaMath.RealScalar
   residue2d::ModiaMath.RealScalar

   CutJointPrismatic(_internal, obj1, obj2, visited, axis, eAxis, e_lambda1, e_lambda2) = new(_internal, obj1, obj2, visited, axis, eAxis, e_lambda1, e_lambda2)
end

function Base.show(io::IO, joint::TreeJointPrismatic)
   print(io,"Prismatic(", ModiaMath.fullName(joint.obj1),
                   ", ", ModiaMath.fullName(joint.obj2),
                      ", axis = ", joint.axis,
                      ", s_start = ", joint.s.start, " m)")
end

#=
function JSON.show_json(io::JSON.StructuralContext,
                         s::JSON.CommonSerialization, joint::TreeJointPrismatic)
   JSON.begin_object(io)
      JSON.show_pair(io, s, "obj_a"  , joint.obj1.name)
      JSON.show_pair(io, s, "obj_b"  , joint.obj2.name)
      JSON.show_pair(io, s, "driven"   , joint.driven)
      JSON.show_pair(io, s, "cutJoint" , joint.cutJoint)
      JSON.show_pair(io, s, "axis"     , joint.axis)
      JSON.show_pair(io, s, "s_start", joint.s.start)
   JSON.end_object(io)
end
=#

get_eAxis(axis::Int) = axis== 1 ? ModiaMath.Vector3D( 1.0,  0.0,  0.0) :
                       axis== 2 ? ModiaMath.Vector3D( 0.0,  1.0,  0.0) :
                       axis== 3 ? ModiaMath.Vector3D( 0.0,  0.0,  1.0) : error("Modia3D.Prismatic: axis = ", axis, " but must be 1, 2 or 3.")

"""
    joint = Modia3D.Prismatic(obj1::Object3D, obj2::Object3D;
                              axis=1, s_start=0, v_start=0, canCollide=false)

Return a `joint` object that constrains the movement of `obj2::`[`Object3D`](@ref)
with respect to `obj1::`[`Object3D`](@ref) along
coordinate axis `axis` (`axis = 1,2,3,-1,-2,-3`). The initial position/velocity of `obj2` with respect
to `obj1` along `axis` is `s_start` [m] and `v_start` [m/s], respectively.
If `canCollide=false`, no collision detection will occur between `obj1` and `obj2`
(and `Object3D`s that are directly or indirectly rigidly fixed to `obj1` or `obj2`).

If a `Prismatic` joint *closes a kinematic loop*, then the already present objects must be consistent
to the `Prismatic` joint that is the frames of `obj1` and `obj2` must be *parallel* to each other and
movement of `obj1` along its axis `axis` with `s_start` results in `obj2`. If `s_start=NaN`,
its value is computed in this case.

# Examples
```julia
using Modia3D
import ModiaMath

@assembly FallingBall(;h=1.0) begin
   world  = Object3D()
   sphere = Object3D( Modia3D.Solid(Modia3D.SolidSphere(0.1), "Aluminium") )

   # Constrain sphere movement (initial placement at position [0,h,0])
   prismatic = Modia3D.Prismatic(world, sphere, axis=2, s_start=h)
end

simulationModel = Modia3D.SimulationModel( FallingBall(h=1.5), stopTime=1.0 )
result          = ModiaMath.simulate!(simulationModel)
ModiaMath.plot(result, ("prismatic.s", "prismatic.v"))
```
"""
function Prismatic(obj1::Object3D, obj2::Object3D;
                   axis::Int = 1,
                   s_start::Number = 0.0,
                   v_start::Number = 0.0,
                   canCollide::Bool = false)
   eAxis = get_eAxis(axis)
   rs_start = Float64(s_start)
   rv_start = Float64(v_start)
   (obj_a,obj_b,cutJoint) = attach(obj1, obj2)
   @assert(1 <= abs(axis) <= 3)
   if obj_b === obj2
      axis2 = axis
   else
      axis2 = -axis
      eAxis = -eAxis
   end

   if cutJoint  # Currently only planar cut-joint, without checking that it is planar
      println("... Prismatic joint connecting ", ModiaMath.fullName(obj1), " with ", ModiaMath.fullName(obj2), " is a cut-joint")

      # Check that obj1 and obj2 are parallel to each other
      R_rel = obj_b.R_abs*obj_a.R_abs
      diff1 = norm(R_rel - ModiaMath.NullRotation)
      if diff1 > Basics.neps
         error("\nError from Modia3D.Prismatic(", ModiaMath.fullName(obj1), ", ", ModiaMath.fullName(obj2), ", ...):\n",
               "The coordinate systems of the provided Object3Ds are not parallel to each other.\n",
               "( norm(R_rel - nullRotation()) = ", string(diff1), " > ", Basics.neps, ")")
      end
      s_start2 =

      if isnan(rs_start)
         # Compute distance along axis
         r_rel = obj_b.r_abs - obj_a.r_abs
         diff2 = cross(eAxis, r_rel)
         if diff2 > Basics.neps
            error("\nError from Modia3D.Prismatic(", ModiaMath.fullName(obj1), ", ", ModiaMath.fullName(obj2), ", ...):\n",
                  "Moving from obj1 along its axis axis does not arrive at obj2.\n",
                  "( sin( angle(obj2.r_abs - obj1.r_abs, axis(obj1) )))  = ", string(diff2), " > ", Basics.neps, ")")
         end

      else
         r2 = obj_a.r_abs + obj_a.R_abs'*(eAxis*s_start)
         diff3 = obj_b.r_abs - r2

         if diff3 > Basics.neps
            error("\nError from Modia3D.Prismatic(", ModiaMath.fullName(obj1), ", ", ModiaMath.fullName(obj2), ", ...):\n",
                 "The origins of the provided Object3Ds are not along axis/s_start.\n",
                 "( norm( <obj2 position> - <obj1.position+s_start along axis>) = ", string(diff33), " > ", Basics.neps,")")
         end
      end

      # Define variables of cut-joint
      joint    =  CutJointPrismatic(ModiaMath.ComponentInternal(),obj_a, obj_b, false, axis, e_lambda1, e_lambda2)
      lambda1   = ModiaMath.RealScalar("lambda1"  , joint, unit="N", info="Constraint force in direction1"   , numericType=ModiaMath.LAMBDA, analysis=ModiaMath.QuasiStaticAndDynamicAnalysis)
      lambda2   = ModiaMath.RealScalar("lambda2"  , joint, unit="N", info="Constraint force in direction2"   , numericType=ModiaMath.LAMBDA, analysis=ModiaMath.QuasiStaticAndDynamicAnalysis)
      mue1      = ModiaMath.RealScalar("mue1"     , joint,           info="Stabilizing mue in direction1"    , numericType=ModiaMath.MUE   , analysis=ModiaMath.OnlyDynamicAnalysis)
      mue2      = ModiaMath.RealScalar("mue2"     , joint,           info="Stabilizing mue in direction2"    , numericType=ModiaMath.MUE   , analysis=ModiaMath.OnlyDynamicAnalysis)
      residue1  = ModiaMath.RealScalar("residue1" , joint,           info="Position constraint in direction1", numericType=ModiaMath.FC)
      residue2  = ModiaMath.RealScalar("residue2" , joint,           info="Position constraint in direction2", numericType=ModiaMath.FC)
      residue1d = ModiaMath.RealScalar("residue1d", joint,           info="Velocity constraint in direction1", numericType=ModiaMath.FC    , analysis=ModiaMath.OnlyDynamicAnalysis)
      residue2d = ModiaMath.RealScalar("residue2d", joint,           info="Velocity constraint in direction2", numericType=ModiaMath.FC    , analysis=ModiaMath.OnlyDynamicAnalysis)

      # direction1/direction2
      if axis == 1
         e_lambda1 = ModiaMath.Vector3D(0,1,0)
         e_lambda2 = ModiaMath.Vector3D(0,0,1)
      elseif axis == 2
         e_lambda1 = ModiaMath.Vector3D(0,0,1)
         e_lambda2 = ModiaMath.Vector3D(1,0,0)
      else
         e_lambda1 = ModiaMath.Vector3D(1,0,0)
         e_lambda2 = ModiaMath.Vector3D(0,1,0)
      end

      # Add joint to the frames
      push!(obj_a.twoObject3Dobject, joint)
      push!(obj_b.twoObject3Dobject, joint)
      obj_a.hasCutJoint = true
      obj_b.hasCutJoint = true
   else
      # Joint in the spanning tree
      if isnan(rs_start)
         error("\nerror from Modia3D.Prismatic(", ModiaMath.fullName(obj1), ", ", ModiaMath.fullName(obj2), ", s_start=NaN):\n",
              "Argument s_start = NaN (Not-a-Number). However, this is not allowed if the joint is in the spanning tree")
      end
      joint   = TreeJointPrismatic(ModiaMath.ComponentInternal(), obj_a, obj_b, canCollide, false, axis, eAxis)
      s       = ModiaMath.RealScalar(:s      , joint, unit="m"    , start=rs_start, fixed=true, info="Relative distance along axis         "    , numericType=ModiaMath.XD_EXP)
      v       = ModiaMath.RealScalar(:v      , joint, unit="m/s"  , start=rv_start, fixed=true, info="=der(s): Relative velocity along axis"    , numericType=ModiaMath.XD_IMP    , integral=s, analysis=ModiaMath.OnlyDynamicAnalysis)
      a       = ModiaMath.RealScalar(:a      , joint, unit="m/s^2", start=0.0     ,             info="=der(v): Relative acceleration along axis", numericType=ModiaMath.DER_XD_IMP, integral=v, analysis=ModiaMath.OnlyDynamicAnalysis)
      f       = ModiaMath.RealScalar(:f      , joint, unit="N"    , start=0.0     ,             info="Driving force"                            , numericType=ModiaMath.WR                    , analysis=ModiaMath.QuasiStaticAndDynamicAnalysis)
      residue = ModiaMath.RealScalar(:residue, joint,               start=0.0     ,             info="Prismatic residue"                        , numericType=ModiaMath.FD_IMP                , analysis=ModiaMath.OnlyDynamicAnalysis)

      P       = ModiaMath.RealScalar(:P      , joint, unit="J", info="Total power computed with interpolated x,derx, values (should be zero)", numericType=ModiaMath.WC, analysis=ModiaMath.OnlyDynamicAnalysis)
      joint.flange = PrismaticFlange()

      obj_a.hasChildJoint = true
      obj_b.joint = joint
      obj_b.r_rel = eAxis*rs_start
      obj_b.r_abs = obj_a.r_abs + obj_a.R_abs'*obj_b.r_rel
      obj_b.R_rel = ModiaMath.NullRotation
      obj_b.R_abs = obj_a.R_abs
   end
   return joint
end

function driveJoint!(joint::TreeJointPrismatic)
   joint.isDriven = true
   joint.s.numericType = ModiaMath.WR
   joint.v.numericType = ModiaMath.WR
   joint.a.numericType = ModiaMath.WR
end

function setDistance!(joint::TreeJointPrismatic, s::Float64)
   @assert( joint.isDriven )
   joint.s.value = s
   joint.obj2.r_rel = joint.eAxis*s
   return joint
end

function computeKinematics!(joint::TreeJointPrismatic, obj::Object3D, analysis::ModiaMath.AnalysisType, time::Float64)::NOTHING
   parent::Object3D = obj.parent

   obj.r_rel = joint.eAxis*joint.s.value
   obj.r_abs = parent.r_abs + parent.R_abs'*obj.r_rel
   obj.R_abs = parent.R_abs

   if analysis == ModiaMath.DynamicAnalysis
      dynamics::Object3Ddynamics       = obj.dynamics
      parentDynamics::Object3Ddynamics = parent.dynamics

      v_rel       = joint.eAxis*joint.v.value
      a_rel       = joint.eAxis*joint.a.value
      dynamics.v0 = parentDynamics.v0 + parent.R_abs'*(v_rel + cross(parentDynamics.w, obj.r_rel))
      dynamics.a0 = parentDynamics.a0 + parent.R_abs'*(a_rel + cross(parentDynamics.z, obj.r_rel) +
                    cross(parentDynamics.w, v_rel + cross(parentDynamics.w, obj.r_rel)) )
      dynamics.w  = parentDynamics.w
      dynamics.z  = parentDynamics.z
   end
   return nothing
end

function computeForceTorqueAndResidue!(joint::TreeJointPrismatic, obj::Object3D, analysis::ModiaMath.AnalysisType, time::Float64)::NOTHING
   parent::Object3D                 = obj.parent
   dynamics::Object3Ddynamics       = obj.dynamics
   parentDynamics::Object3Ddynamics = parent.dynamics
   parentDynamics.f  += dynamics.f
   parentDynamics.t  += dynamics.t + cross(obj.r_rel, dynamics.f)
   dynamics.f         = -dynamics.f
   dynamics.t         = -dynamics.t

   joint.residue.value = -joint.f.value + dot(joint.eAxis,dynamics.f)

   # For checking: compute total power
   joint.P.value = dot(parentDynamics.f, joint.obj1.R_abs*parentDynamics.v0) +
                   dot(dynamics.f      , joint.obj2.R_abs*dynamics.v0) +
                   dot(parentDynamics.t, parentDynamics.w) +
                   dot(dynamics.t      , dynamics.w) -
                   joint.f.value*joint.v.value
   return nothing
end

function computePositionResidues!(joint::CutJointPrismatic, time::Float64)
   r_rel = joint.obj1.R_abs*(joint.obj2.r_abs - joint.obj1.r_abs)
   joint.residue_x.value = dot(joint.e_lambda1, r_rel)
   joint.residue_y.value = dot(joint.e_lambda2, r_rel)
   return nothing
end

function computeVelocityResidues!(joint::CutJointPrismatic, time::Float64)
   v_rel = joint.obj1.R_abs*(joint.obj2.dynamics.v0 - joint.obj1.dynamics.v0)
   joint.residue_vx.value = dot(joint.e_lambda1, v_rel)
   joint.residue_vy.value = dot(joint.e_lambda2, v_rel)
   return nothing
end

function computeCutForcesAndToques!(joint::CutJointPrismatic, time::Float64)
   dynamics1.f += joint.e_lambda1*lambda1 + joint.e_lambda2*lambda2
   dynamics2.f += dynamics.f
end

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

function assignPrismaticFlange(flange::PFlange, prismFlange::PrismaticFlange, joint::Modia3D.AbstractPrismatic, assembly::Modia3D.Composition.AssemblyInternal)
  names = fieldnames(typeof(flange))
  for val in names # loops over s,v,a,f
    flange_val    = getfield(flange, val)
    prismFlange_val = getfield(prismFlange,val)
    joint_val     = getfield(joint,val)
    if causalitiesAreOk(flange_val, joint_val)
        setFlangeVariable!(flange_val, prismFlange_val)
        setJointVariable!(joint_val, prismFlange_val)
        addVariableToArray(flange_val, joint_val, assembly)
    else
      return false
    end
  end
  return true
end

function connect(torqueObj::Modia3D.AbstractForceAdaptor, joint::Modia3D.AbstractPrismatic)
  actualAssembly = torqueObj._internal.within._internal

  flangeTorque  = torqueObj.flange
  flangeJoint   = joint.flange

  if !assignPrismaticFlange(flangeTorque, flangeJoint, joint, actualAssembly) # !assignRevoluteFlangeAndSignal(signal, flangeSig, joint, flangeJoint, assembly)
    error("The causalities are different. Please, check all flanges which are connected with joint = ", joint ,".")
  end
  if !in(torqueObj.forceElement, actualAssembly.uniqueForceTorques)
    push!(actualAssembly.uniqueForceTorques, torqueObj.forceElement)
  end
  if !isSInput(flangeJoint) && isFInput(flangeJoint)
    joint.isDriven = false
  elseif isSInput(flangeJoint) && !isFInput(flangeJoint)
    driveJoint!(joint)
  end
end
connect(joint::Modia3D.AbstractPrismatic, torqueObj::Modia3D.AbstractForceAdaptor) = connect(torqueObj::Modia3D.AbstractForceAdaptor, joint::Modia3D.AbstractPrismatic)

function connect(signal::Modia3D.AbstractSignalAdaptor, joint::Modia3D.AbstractPrismatic)
    actualAssembly = signal._internal.within._internal

    flangeSig   = signal.flange
    flangeJoint = joint.flange

    if !assignPrismaticFlange(flangeSig, flangeJoint, joint, actualAssembly) # !assignRevoluteFlangeAndSignal(signal, flangeSig, joint, flangeJoint, assembly)
      error("The causalities are different. Please, check all flanges which are connected with joint = ", joint ,".")
    end
    if !in(signal.signal, actualAssembly.uniqueSignals)
      push!(actualAssembly.uniqueSignals, signal.signal)
    end
    if !isSInput(flangeJoint) && isFInput(flangeJoint)
      joint.isDriven = false
    elseif isSInput(flangeJoint) && !isFInput(flangeJoint)
      driveJoint!(joint)
    end
end
connect(joint::Modia3D.AbstractPrismatic, signal::Modia3D.AbstractSignalAdaptor) = connect(signal::Modia3D.AbstractSignalAdaptor, joint::Modia3D.AbstractPrismatic)


text:
function δ(X, discrete = true)
    mn, mx = extrema(Iterators.flatten(X))
    return mx - mn + discrete
end

function δ(X, Y, discrete = true)
    mnx, mxx = extrema(Iterators.flatten(X))
    mny, mxy = extrema(Iterators.flatten(Y))
    return max(mxx, mxy) - min(mnx, mny) + discrete
end

to_domains(domain_sizes::Vector{Int}) = map(ds -> domain(0:ds), domain_sizes)

function to_domains(X, ds::Int = δ(X))
    d = domain(0:ds-1)
    return fill(d, length(first(X)))
end

function to_domains(X, d::D) where {D <: DiscreteDomain}
    n::Int = length(first(X)) / domain_size(d)
    return fill(d, n)
end


text:
mutable struct GaussianProcessRegressor
    _kernel::Kernel
    _beta::Float64
    # cache the result of prior training
    _cov::AbstractArray{Float64,2}
    _precision::AbstractArray{Float64,2}
    _X::AbstractArray{Float64,2}
    _t::AbstractArray{Float64,1}
    function GaussianProcessRegressor(kernel::Kernel, beta = 1.0)
        new(kernel, beta, zeros(1, 1), zeros(1, 1), zeros(1, 1), zeros(1))
    end
end

function log_likelihood(gp::GaussianProcessRegressor)
    return -0.5 * log(abs(det(gp._cov))) -
           0.5 * (transpose(gp._t) * gp._precision * gp._t) -
           0.5 * length(gp._t) * log(2.0 * pi)
end

"""
    fitting(gp::GaussianProcessRegressor, X::AbstractArray{Float64, 2}, t::AbstractArray{Float64, 1}, iter_max=0, learning_rate=0.1)

For training data X = [x_1,.... x_N] of size (n_feature, n_samples) and t = [t_1,... t_N] of size(n_samples)
"""
function fitting(
    gp::GaussianProcessRegressor,
    X::AbstractArray{Float64,2},
    t::AbstractArray{Float64,1},
    iter_max = 0,
    learning_rate = 0.1,
)
    n_samples = size(X)[2]
    n_samples == size(t)[1] || error("size does not match")

    gp._X = copy(X)
    gp._t = copy(t)
    E = Matrix(I, n_samples, n_samples) * 1.0
    K = kernel(gp._kernel, X, X)
    gp._cov = K + E / gp._beta
    gp._precision = inv(gp._cov)

    log_likes = [-Inf]
    # learn hyper parameter
    precision = gp._precision
    for i = 1:iter_max
        n_params = gp._kernel._n_params
        grads = derivatives(gp._kernel, X, X)
        updates = [
            -tr(precision * grads[n, :, :]) +
            transpose(t) * precision * grads[n, :, :] * precision * t for n = 1:n_params
        ]
        for j = 1:iter_max
            update_parameters(gp._kernel, learning_rate * updates)
            K = kernel(gp._kernel, X, X)
            gp._cov = K + E / gp._beta
            gp._precision = inv(gp._cov)
            log_like = log_likelihood(gp)
            if log_like > log_likes[end]
                push!(log_likes, log_like)
                break
            else
                update_parameters(gp._kernel, -learning_rate * updates)
                learning_rate *= 0.9
            end
        end
    end
    popfirst!(log_likes)
    return log_likes
end

"""
    predict(gp::GaussianProcessRegressor, X::AbstractArray{Float64, 2}, with_error=false)

For input X = [x_1,... x_N] of size (n_feature, n_samples), return mus = [m(x_1),... m(x_N)] of size (n_feature, n_samples) and sigmas = [sigma(x_1),... sigma(x_N)] of size(n_feature, n_feature, n_samples) if with_error == true. Otherwise return mus.
"""
function predict(
    gp::GaussianProcessRegressor,
    X::AbstractArray{Float64,2},
    with_error = false,
)
    n_feature = size(X)[1]
    n_samples = size(X)[2]

    mus = zeros(1, n_samples)
    sigmas = zeros(1, 1, n_samples)

    for n = 1:n_samples
        x_n = X[:, n]
        x_n = reshape(x_n, size(x_n)[1], 1)
        k_n = kernel(gp._kernel, gp._X, x_n)
        size(k_n) == (size(gp._X)[2], 1) ||
            error("size does not match: size(k_n) = $(size(k_n))")
        c = kernel(gp._kernel, x_n, x_n)
        size(c) == (1, 1) || error("size does not match")
        c = c .+ 1.0 / gp._beta
        mus[:, n] = transpose(k_n) * gp._precision * gp._t
        if with_error
            k_n = reshape(k_n, length(k_n))
            sigmas[:, :, n] = c .- transpose(k_n) * gp._precision * k_n
        end
    end

    if with_error
        return mus, sigmas
    else
        return mus
    end
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
include("DataTypes.jl")

WGS84_wkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"

"""
    GeoArray(A::AbstractArray{T,3} where T <: Union{Real,Union{Missing,Real}})

Construct a GeoArray from any Array. A default `AffineMap` and `CRS` will be generated.

# Examples
```julia-repl
julia> GeoArray(rand(10,10,1))
10x10x1 Array{Float64, 3} with AffineMap([1.0 0.0; 0.0 1.0], [0.0, 0.0]) and undefined CRS
```
"""
GeoArray(A::AbstractArray{T, 3}) where T<:Union{Real, Union{Missing, Real}} =
    GeoArray(A, geotransform_to_affine(SVector(0.,1.,0.,0.,0.,1.)), "")

GeoArray(A::AbstractArray{T, 3}, f::AffineMap) where T<:Union{Real, Union{Missing, Real}} =
    GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), ""))

GeoArray(A::AbstractArray{T, 3}, f::AffineMap, crs::String) where T<:Union{Real, Union{Missing, Real}} = begin
    if crs == ""; crs = WGS84_wkt; end
    GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), crs))
end

function GeoArray(A::AbstractArray{T, 3}, x::AbstractRange, y::AbstractRange, args...) where T<:Union{Real, Union{Missing, Real}}
    size(A)[1:2] != (length(x), length(y)) &&
        error("Size of `GeoArray` $(size(A)) does not match size of (x,y): $((length(x),length(y))). Note that this function takes *center coordinates*.")
    f = unitrange_to_affine(x, y)
    GeoArray(A, f, args...)
end

function GeoArray(A::AbstractArray{T, 3}, b::bbox; proj = 4326) where T<:Union{Real, Union{Missing, Real}}
    ga = GeoArray(A)
    st_bbox!(ga, b)
    epsg!(ga, proj)
    ga
end

GeoArray(A::AbstractArray{T, 2} where T<:Union{Real, Union{Missing, Real}}, args...) =
    GeoArray(reshape(A, size(A)..., 1), args...)

Base.size(ga::GeoArray) = size(ga.A)
Base.IndexStyle(::Type{T}) where {T<:GeoArray} = IndexLinear()
Base.iterate(ga::GeoArray) = iterate(ga.A)
Base.length(ga::GeoArray) = length(ga.A)
Base.parent(ga::GeoArray) = ga.A
# Base.map(f, ga::GeoArray) = GeoArray(map(f, ga.A), ga.f, ga.crs)
# Base.convert(::Type{Array{T, 3}}, A::GeoArray{T}) where {T} = convert(Array{T,3}, ga.A)
Base.eltype(::Type{GeoArray{T}}) where {T} = T

Base.show(io::IO, ::MIME"text/plain", ga::GeoArray) = show(io, ga)
function Base.show(io::IO, ga::GeoArray)
    crs = GeoFormatTypes.val(ga.crs)
    wkt = length(crs) == 0 ? "undefined CRS" : "CRS $crs"
    print(io, "$(join(size(ga), "x")) $(typeof(ga.A)) with $(ga.f) and $(wkt)")
end

include("st_[.jl")
include("st_Ops.jl")


text:
h5file  = "temp.h5"

@testset "TrainingSet" begin

    nd = 2
    np = 10
    nh = 5
    ns = 11

    μ = Parameter(:μ, 0.0, 1.0, 3)
    ν = Parameter(:ν, 1.0, 1.0, 1)
    σ = Parameter(:σ, 0.0, 4.0, 2)

    ps = ParameterSpace(μ, ν, σ)

    X = rand(nd, np, ns, length(ps))
    V = rand(nd, np, ns, length(ps))
    A = rand(nd, np, ns, length(ps))
    Φ = rand(nd, nh, ns, length(ps))
    W = rand(ns, length(ps))
    K = rand(ns, length(ps))
    M = rand(ns, length(ps))

    ts1 = TrainingSet(ps, X, V, A, Φ, W, K, M)

    @test ts1.paramspace == ps

    @test ts1.X == X
    @test ts1.V == V
    @test ts1.A == A
    @test ts1.Φ == Φ
    @test ts1.W == W
    @test ts1.K == K
    @test ts1.M == M

    h5save(ts1, h5file; mode="w")
    @test isfile(h5file)

    ts2 = h5load(TrainingSet, h5file)
    rm(h5file)
    @test ts1 == ts2

end


text:
function complete_cases(x::DataArray)
    xna=sum(!isna(x), 2)
    ind=Int[]
    n=size(x,1)
    ncols=size(x,2)
    vec=1:n
    for i=1:n
        if xna[i]==ncols
            push!(ind, i)
        end
    end
    return vec[ind]
end

function deleterows!(x::DataArray)
    keep_inds::Array{Int}=complete_cases(x)
    output=DataArray(Any, length(keep_inds), size(x,2))
    for i=1:size(x,2)
        for j=1:length(keep_inds)
            output[j,i]=x[keep_inds[j],i]
        end
    end
    return output
end

#Finding duplicated values
function duplicated{S <: Real}(x::Vector{S})
    xsort=sort(x)
    dup=false
    for i=1:(length(xsort)-1)
        if xsort[i]==xsort[i+1]
            dup=true
            break
        end
    end
    return dup
end
#duplicated(x::Vector{Int})=duplicated(convert(Vector{Float64}, x))
#duplicated(x::DataVector)=duplicated(removeNA(x))

function seq(Start::Number, End::Number, increment::Number)
    delta=End-Start
    if delta.<=0
        error("`start` must be larger than `end`")
    else
        n::Int=floor(delta/increment)+1
        output=zeros(n)
        for i=1:n
            output[i]=(i-1).*increment+Start
        end
    end
    return output
end

function trimcibt_loop(x, n, nboot, tr, side, randid, test)
    xbar=tmean(x, tr=tr)
    temp=zeros(n)
    tval=zeros(nboot)
    if side
        test=abs(test)
        pval=0.0
        for i=1:(nboot*n)
            if (i%n)!=0
                temp[i%n]=x[randid[i]]-xbar
            else
                temp[n]=x[randid[i]]-xbar
                tval[div(i, n)]=abs(tmean(temp, tr=tr)./trimse(temp, tr=tr))
                pval += tval[div(i, n)]>=test? 1/nboot :0.0
            end
        end
        return sort!(tval), pval
    else
        for i=1:(nboot*n)
            if (i%n)!=0
                temp[i%n]=x[randid[i]]-xbar
            else
                temp[n]=x[randid[i]]-xbar
                tval[div(i, n)]=tmean(temp, tr=tr)./trimse(temp, tr=tr)
            end
        end
        return sort!(tval)
    end
end

function akerd_loop(x, pts, hval, alam)
    npts  = length(pts)
    dhat  = zeros(npts)
    sqrt5 = sqrt(5)
    n     = length(x)
    temp  = zeros(n)
    for i = 1:npts
        epan_alam_hval = 0.0
        for j = 1:n
            temp[j] = (pts[i]-x[j])/(hval*alam[j])
            if abs(temp[j]) < sqrt5
                epan_alam_hval += (0.75*(1-0.2*temp[j]*temp[j])/sqrt5)/alam[j]
            end
        end
        dhat[i]=epan_alam_hval/npts
    end
    return dhat./hval
end

#Faster Pearson's r for
function cor!(x, y)
    n=length(x)
    xbar=mean(x); ybar=mean(y)
    top=bot1=bot2=0.0
    for i=1:n
        top+=(x[i]-xbar).*(y[i]-ybar)
        bot1+=(x[i]-xbar)*(x[i]-xbar)
        bot2+=(y[i]-ybar)*(y[i]-ybar)
    end
    top./sqrt(bot1.*bot2)
end

function regts1(vstar, yhat, res, mflag, x, tr)
    const n=length(res)
    ystar=zeros(n)
    [ystar[i]=yhat+res[i].*vstar[i] for i=1:n]
    rval=zeros(n)
    ystar_bar=tmean(ystar, tr=tr)
    rval=zeros(n)
    [rval[i]= sum(ystar[mflag[:,i]]-ystar_bar) for i=1:n]
    return rval
end

function indt_mflag(x)
    n=size(x, 1)
    np=size(x, 2)
    output=falses(n, n)
    for i=1:n
        for j=1:n
            total=0
            for k=1:np
                total+=x[i,k]<=x[j,k]?1:0
            end
            if total==np
                output[i,j]=true
            end
        end
    end
    return output
end

function bootindirect{S <: Real, T <: Real, W <: Real}(x::Vector{S}, y::Vector{T}, m::Vector{W}, nboot::Integer)
    n = length(x)
    tempx  = hcat(rep(1.0, n), zeros(n))
    tempy  = zeros(n)
    tempm  = zeros(n)
    tempmx = hcat(rep(1.0, n), zeros(Real, n), zeros(Real, n))
    bvec   = zeros(nboot)
    randid = rand(1:n, n*nboot)
    for i = 1:nboot
        for j = 1:n
            tempx[j,2] = tempmx[j,2] = x[randid[(i-1)*n + j]]
            tempy[j]   = y[randid[(i-1)*n + j]]
            tempm[j]   = tempmx[j,3] = m[randid[(i-1)*n + j]]
        end
        bvec[i] = (inv(tempx'tempx)'*(tempx'tempm))[2]*(inv(tempmx'tempmx)'*(tempmx'tempy))[3]
    end
    bvec
end

function mad!(v::AbstractArray)
    v=copy(v)
    center = median!(v, checknan=false)
    for i in 1:length(v)
        v[i] = abs(v[i]-center)
    end
    1.4826 * median!(v, checknan=false)
end

function t1waycore(args...)
    x, grp, tr, method = args[1], args[2], args[3], args[4]
    g = unique(grp)
    J = length(g)
    h, w, xbar, n = zeros(J), zeros(J), zeros(J), zeros(J)

    for j = 1:J
        n[j] = sum(grp.==g[j])
        temp = x[grp.==g[j]]
        h[j] = n[j] - 2*floor(tr*n[j])
        w[j] = h[j]*(h[j] - 1)/((n[j] -1)*winvar(temp, tr=tr))
        xbar[j] = tmean(temp, tr=tr)
    end
    u = sum(w)
    xtil = sum(w .* xbar)/u
    A = sum(w .* (xbar - xtil).*(xbar - xtil))/(J - 1)
    B = 2*(J - 2)*sum((1 - w./u).*(1 - w./u)./(h - 1))/(J*J - 1)
    statistic = A/(B + 1)
    nu1 = J - 1
    nu2 = 1/(3 * B/2/(J-2))
    p   = 1 - Rmath.pf(statistic, nu1, nu2)
    if method
        METHOD="Heteroscedastic one-way ANOVA for trimmed means\nusing a generalization of Welch's method.\n"
    else
        METHOD=nothing
    end
    output    = testOutput()
    output.n  = n
    output.df = [nu1, nu2]
    output.statistic = statistic
    output.p  = p
    output.method = METHOD
    output
end


text:
module DashDraggable
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.1.1a"

include("''_dashboarditem.jl")
include("''_dashboarditemresponsive.jl")
include("''_draggabledashboard.jl")
include("''_draggabledashboardresponsive.jl")
include("''_gridlayout.jl")
include("''_responsivegridlayout.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dash_draggable",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "dash_draggable.min.js",
    external_url = nothing,
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_draggable.min.js.map",
    external_url = nothing,
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end

max_stars_repo_name: lonlazer/dash-draggable

max_issues_repo_name: lonlazer/dash-draggable

max_forks_repo_name: lonlazer/dash-draggable


text:
function image_widget(fn)
  img = Gtk.GtkImageLeaf(fn)
  vbox = Gtk.GtkBoxLeaf(:v)
  push!(vbox, Gtk.GtkLabelLeaf(fn))
  push!(vbox, img)
  show(img)
  vbox
end

function replace_refimg(tmpfn, reffn)
  cmd = `cp $tmpfn $reffn`
  run(cmd)
  info("Replaced reference image with: $cmd")
end

"Show a Gtk popup with both images and a confirmation whether we should replace the new image with the old one"
function replace_refimg_dialog(tmpfn, reffn_old; reffn_new = reffn_old)

  # add the images
  imgbox = Gtk.GtkBoxLeaf(:h)
  push!(imgbox, image_widget(tmpfn))
  push!(imgbox, image_widget(reffn_old))

  win = Gtk.GtkWindowLeaf("Should we make this the new reference image?")
  push!(win, Gtk.GtkFrameLeaf(imgbox))

  showall(win)

  # now ask the question
  if Gtk.ask_dialog("Should we make this the new reference image?", "No", "Yes")
    replace_refimg(tmpfn, reffn_new)
  end

  destroy(win)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule OpenJpeg_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("OpenJpeg")
JLLWrappers.@generate_main_file("OpenJpeg", UUID("643b3616-a352-519d-856d-80112ee9badc"))
end  # module OpenJpeg_jll


text:
function dataProvider()
    L=4;T=10
    t=0:0.5:10000;x=cos(t);y=x[[2:end;1]]
    nData=1000;inputBase=zeros(1,nData,T);outputBase=zeros(1,nData,T);state=1
    for data=1:nData,t=1:T
      inputBase[:,data,t]=x[state];outputBase[:,data,t]=y[state]; state+=1
    end
  return inputBase,outputBase
end


text:
# Load Julia packages (libraries) needed  for the snippets in chapter 0

using StatisticalRethinking
using CmdStan
#gr(size=(600,600));

# CmdStan uses a tmp directory to store the output of cmdstan

ProjDir = rel_path("..", "scripts", "04")
cd(ProjDir)

# ### Preliminary snippets

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

# Use only adults

df2 = filter(row -> row[:age] >= 18, df);

# Center weight and store as weight_c

mean_weight = mean(df2[:weight])
df2 = hcat(df2, df2[:weight] .- mean_weight)
rename!(df2, :x1 => :weight_c); # Rename our col :x1 => :weight_c
first(df2, 5)

# Define the Stan language model

weightsmodel = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] height; // Predictor
 vector[N] weight; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 height ~ normal(alpha + weight * beta , sigma);
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

stanmodel = Stanmodel(name="weights", monitors = ["alpha", "beta", "sigma"],model=weightsmodel,
  output_format=:mcmcchains);

# Input data for cmdstan.

heightsdata = Dict("N" => length(df2[:height]), "height" => df2[:height], "weight" => df2[:weight_c]);

# Sample using cmdstan

rc, chn, cnames = stan(stanmodel, heightsdata, ProjDir, diagnostics=false,
  summary=false, CmdStanDir=CMDSTAN_HOME);

# ### Snippet 4.47

# Show first 5 draws of correlated parameter values in chain 1

chn.value[1:5,:,1]

# ### Snippets 4.48 & 4.49

# Plot estimates using the N = [10, 50, 150, 352] observations

nvals = [10, 50, 150, 352];

# Create the 4 nvals plots

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 4)
for i in 1:length(nvals)
  N = nvals[i]
  heightsdataN = [
    Dict("N" => N, "height" => df2[1:N, :height], "weight" => df2[1:N, :weight_c])
  ]
  rc, chnN, cnames = stan(stanmodel, heightsdataN, ProjDir, diagnostics=false,
    summary=false, CmdStanDir=CMDSTAN_HOME)

  rws, vars, chns = size(chnN)
  xi = -15.0:0.1:15.0
  alpha_vals = convert(Vector{Float64}, reshape(chnN.value[:, 1, :], (rws*chns)))
  beta_vals = convert(Vector{Float64}, reshape(chnN.value[:, 2, :], (rws*chns)))

  p[i] = scatter(df2[1:N, :weight_c], df2[1:N, :height], leg=false, xlab="weight_c")
  for j in 1:N
    yi = alpha_vals[j] .+ beta_vals[j]*xi
    plot!(p[i], xi, yi, title="N = $N")
  end
end
plot(p..., layout=(2, 2))

# ### Snippet 4.50

# Get dimensions of chains

rws, vars, chns = size(chn)
mu_at_50 = link(50:10:50, chn, [1, 2], mean_weight);
density(mu_at_50)

# ### Snippet 4.54

# Show posterior density for 6 mu_bar values

mu = link(25:10:75, chn, [1, 2], mean_weight);

q = Vector{Plots.Plot{Plots.GRBackend}}(undef, size(mu, 1))
for i in 1:size(mu, 1)
  q[i] = density(mu[i], ylim=(0.0, 1.5),
    leg=false, title="mu_bar = $(round(mean(mu[i]), digits=1))")
end
plot(q..., layout=(2, 3), ticks=(3))

# End of `clip_48_54s.jl`


text:
using DrWatson
@quickactivate "EPExperiment"

using EPExperiments

using FileIO
using JLD2
using StatsPlots
using AnnealedIS
using StatsFuns: logsumexp

function AnnealedIS.effective_sample_size(log_weights)
    denominator = logsumexp(2 * log_weights)
    numerator = 2 * logsumexp(log_weights)
    return exp(numerator - denominator)
end

function post_hoc_analysis(
    data_filename,
    is_ground_truth_name,
    taanis_name;
    mcmc_baseline_name=nothing,
    make_joint_plot=false
)
    is_ground_truth_dir = datadir("is_ground_truth", is_ground_truth_name)
    taanis_dir = datadir("inferences", taanis_name)
    input_data = load(datadir("sims", data_filename))

    true_gamma = input_data["params"].γ

    is_results = load(joinpath(is_ground_truth_dir, "results.jld2"))
    taanis_results = load(joinpath(taanis_dir, "results.jld2"))

    # experiment_config = taanis_results["experiment_config"]
    # @unpack cost_fn_name = experiment_config
    cost_fn_name = :cost_fn_sir
    @show cost_fn_name
    cost_fn = COST_FUNCTIONS[cost_fn_name]

    true_expectation_value = is_results["expectation_estimate"]

    num_xlogspaced = 100
    ep, cp = convergence_plot(
        "sir",
        taanis_results["ode_tabi"][:Z1_positive_info],
        taanis_results["ode_tabi"][:Z2_info],
        (i, b) -> cost_fn(b, true_gamma),
        true_expectation_value;
        thinning=1,
        num_xlogspaced=num_xlogspaced
    )
    savefig(ep, joinpath(taanis_dir, "test_error_plot_log.png"))
    #savefig(cp, joinpath(taanis_dir, "convergence_plot.png"))
    # ep, cp = convergence_plot(
    #     "sir",
    #     taanis_results["ode_tabi"][:Z1_positive_info],
    #     taanis_results["ode_tabi"][:Z2_info],
    #     (i, b) -> cost_fn(b, true_gamma),
    #     true_expectation_value;
    #     thinning=1,
    #     num_xlogspaced=num_xlogspaced,
    #     mean_and_std=false
    # )
    # savefig(ep, joinpath(taanis_dir, "error_plot_means.png"))

    Z1_lws = Array(taanis_results["ode_tabi"][:Z1_positive_info][:log_weight])
    Z1_betas = Array(taanis_results["ode_tabi"][:Z1_positive_info][:β])

    Z2_lws = Array(taanis_results["ode_tabi"][:Z2_info][:log_weight])
    Z2_betas = Array(taanis_results["ode_tabi"][:Z2_info][:β])

    # Z1_ess = effective_sample_size(Z1_lws)
    # Z2_ess = effective_sample_size(Z2_lws)

    Z2_lws_retargeted = Z2_lws .+ log.(cost_fn.(Z2_betas, true_gamma))
    Z1_lws_retargeted = Z1_lws .- log.(cost_fn.(Z1_betas, true_gamma))

    # Z2_ess_retargeted = effective_sample_size(Z2_lws_retargeted)
    # Z1_ess_retargeted = effective_sample_size(Z1_lws_retargeted)

    ess_plot = make_ess_plot(
        Z1_lws,
        Z2_lws,
        Z2_lws_retargeted;
        yaxisscale=:log10,
        take_min=false,
        legend=false,
        num_xlogspaced=num_xlogspaced
    )
    savefig(ess_plot, joinpath(taanis_dir, "sir_ess_not_min.pdf"))

    if !isnothing(mcmc_baseline_name)
        mcmc_results = load(datadir(
            "sir_nuts_baseline",
            mcmc_baseline_name,
            "ess.jld2"
        ))
        mcmc_estimates = mcmc_results["exp_estimates"]
        num_intermediates = mcmc_results["num_intermediates"]
        fake_ess = mcmc_results["fake_ess"]
        @show mcmc_estimates[end]
        @show true_expectation_value
        @show relative_squared_error(
            true_expectation_value, mcmc_estimates[end]
        )
        mcmc_rse = relative_squared_error.(
            true_expectation_value, mcmc_estimates
        )

        @show size(mcmc_estimates)
        @show mcmc_rse[end,:]
        @show mcmc_estimates[end,:]

        mds, qs = get_median_and_quantiles(mcmc_rse)

        thinning = 100
        num_samples = size(Z2_lws, 1)
        # TODO: Remove magic numbers
        # ixs = Array(0:4500:size(mcmc_rse, 1))
        # ixs[1] = 90
        # @show size(mcmc_rse[ixs])
        # plot_ixs = Array(0:50:num_samples)
        # plot_ixs[1] = 1
        # @show num_samples
        # @show size(plot_ixs)
        if !isnothing(num_xlogspaced)
            log_ixs = unique(round.(Int, geomspace(1, num_samples, num_xlogspaced)))
            plot_ixs = log_ixs
            @show 90*log_ixs
            mds = mds[90*log_ixs]
            qs = qs[:,90*log_ixs]
        else
            ixs = Array(0:4500:size(mcmc_rse, 1))
            ixs[1] = 90
            @show size(mcmc_rse[ixs])
            plot_ixs = Array(0:50:num_samples)
            plot_ixs[1] = 1
            mds = mds[ixs]
            qs = qs[ixs]
        end
        plot!(
            ep,
            plot_ixs,
            mds,
            ribbon=[(mds.-qs[1,:]), (qs[2,:]-mds)],
            label="MCMC",
            legend=:bottomleft
        )
        savefig(ep, joinpath(taanis_dir, "sir_error_plot_mcmc_log.pdf"))

        # plot!(
        #     ess_plot,
        #     1:num_samples,
        #     fake_ess,
        #     label="MCMC"
        # )
        # savefig(ess_plot, joinpath(taanis_dir, "sir_ess_test_mcmc.png"))
    end

    if make_joint_plot
        @unpack i_0_true, β_true = input_data["params"]
        joint_plot = plot_sir_samples(
            taanis_results["ode_tabi"][:Z2_info],
            color=palette(:default)[2]
        )
        savefig(joint_plot, joinpath(taanis_dir, "anis_Z2_joint_samples.png"))

        joint_plot = plot_sir_samples(
            taanis_results["ode_tabi"][:Z1_positive_info],
            color=palette(:default)[1]
        )
        savefig(joint_plot, joinpath(taanis_dir, "Z1_joint_samples.png"))
    end

    return
    num_samples = length(Z2_lws)
    # TODO: Is the logevidence field in the Chains object still correct if we merged several chains?
    Z2_weight_normalisation = taanis_results["ode_tabi"][:Z2_info].logevidence + log(num_samples)
    Z1_weight_normalisation = taanis_results["ode_tabi"][:Z1_positive_info].logevidence + log(num_samples)
    savefig(
        plot_is_samples(Z2_betas, Z2_lws .- Z2_weight_normalisation),
        joinpath(taanis_dir, "Z2_betas.png")
    )
    savefig(
        plot_is_samples(Z1_betas, Z1_lws .- Z1_weight_normalisation),
        joinpath(taanis_dir, "Z1_betas.png")
    )
end

post_hoc_analysis(
    "neg_binomial_data.jld2",
    GROUND_TRUTH_FOLDER, # TODO: Give name of folder with ground truth.
    RESULTS_FOLDER; # TODO: Give name of folder which contains (TA)AnIS results.
    # mcmc_baseline_name="combined_baselines_5replications",
    make_joint_plot=true
)


text:
using Base.Test
include("UTF16.jl")
using UTF16

u8 = "𝕥𝟶f𠂊"
u16 = utf16(u8)
@test length(u16.data) == 7
@test length(u16) == 4
@test utf8(u16) == u8
@test collect(u8) == collect(u16)


text:
import MathProgBase

export Problem, Solution, minimize, maximize, satisfy, add_constraint!, add_constraints!
export Float64OrNothing

const Float64OrNothing = Union{Float64, Nothing}

# TODO: Cleanup
mutable struct Solution{T<:Number}
    primal::Array{T, 1}
    dual::Array{T, 1}
    status::Symbol
    optval::T
    has_dual::Bool
end

Solution(x::Array{T, 1}, status::Symbol, optval::T) where {T} =
    Solution(x, T[], status, optval, false)
Solution(x::Array{T, 1}, y::Array{T, 1}, status::Symbol, optval::T) where {T} =
    Solution(x, y, status, optval, true)

mutable struct Problem
    head::Symbol
    objective::AbstractExpr
    constraints::Array{Constraint}
    status::Symbol
    optval::Float64OrNothing
    model::Union{MathProgBase.AbstractConicModel, Nothing}
    solution::Solution

    function Problem(head::Symbol, objective::AbstractExpr,
                     model::Union{MathProgBase.AbstractConicModel, Nothing},
                     constraints::Array=Constraint[])
        if sign(objective)== Convex.ComplexSign()
            error("Objective can not be a complex expression")
        else
            return new(head, objective, constraints, Symbol("not yet solved"), nothing, model)
        end
    end
end

# constructor if model is not specified
function Problem(head::Symbol, objective::AbstractExpr, constraints::Array=Constraint[],
                 solver::Union{MathProgBase.AbstractMathProgSolver, Nothing}=nothing)
    model = solver !== nothing ? MathProgBase.ConicModel(solver) : solver
    Problem(head, objective, model, constraints)
end

# If the problem constructed is of the form Ax=b where A is m x n
# returns:
# index: n
# constr_size: m
# var_to_ranges a dictionary mapping from variable id to (start_index, end_index)
# where start_index and end_index are the start and end indexes of the variable in A
function find_variable_ranges(constraints)
    index = 0
    constr_size = 0
    var_to_ranges = Dict{UInt64, Tuple{Int, Int}}()
    for constraint in constraints
        for i = 1:length(constraint.objs)
            for (id, val) in constraint.objs[i]
                if !haskey(var_to_ranges, id) && id != objectid(:constant)
                    var = id_to_variables[id]
                    if var.sign == ComplexSign()
                        var_to_ranges[id] = (index + 1, index + 2*length(var))
                        index += 2*length(var)
                    else
                        var_to_ranges[id] = (index + 1, index + length(var))
                        index += length(var)
                    end
                end
            end
            constr_size += constraint.sizes[i]
        end
    end
    return index, constr_size, var_to_ranges
end

function vexity(p::Problem)
    bad_vex = [ConcaveVexity, NotDcp]

    obj_vex = vexity(p.objective)
    if p.head == :maximize
        obj_vex = -obj_vex
    end
    typeof(obj_vex) in bad_vex && @warn "Problem not DCP compliant: objective is not DCP"

    constr_vex = ConstVexity()
    for i in 1:length(p.constraints)
        vex = vexity(p.constraints[i])
        typeof(vex) in bad_vex && @warn "Problem not DCP compliant: constraint $i is not DCP"
        constr_vex += vex
    end
    problem_vex = obj_vex + constr_vex
    # this check is redundant
    # typeof(problem_vex) in bad_vex && warn("Problem not DCP compliant")
    return problem_vex
end

function conic_form!(p::Problem, unique_conic_forms::UniqueConicForms=UniqueConicForms())
    objective_var = Variable()
    objective = conic_form!(objective_var, unique_conic_forms)
    conic_form!(p.objective - objective_var == 0, unique_conic_forms)
    for constraint in p.constraints
        conic_form!(constraint, unique_conic_forms)
    end
    return objective, objective_var.id_hash
end

function conic_problem(p::Problem)
    if length(p.objective) != 1
        error("Objective must be a scalar")
    end

    # conic problems have the form
    # minimize c'*x
    # st       b - Ax \in cones
    # our job is to take the conic forms of the objective and constraints
    # and convert them into vectors b and c and a matrix A
    # one chunk of rows in b and in A corresponds to each constraint,
    # and one chunk of columns in b and A corresponds to each variable,
    # with the size of the chunk determined by the size of the constraint or of the variable

    # A map to hold unique constraints. Each constraint is keyed by a symbol
    # of which atom generated the constraints, and a integer hash of the child
    # expressions used by the atom
    unique_conic_forms = UniqueConicForms()
    objective, objective_var_id = conic_form!(p, unique_conic_forms)
    constraints = unique_conic_forms.constr_list
    # var_to_ranges maps from variable id to the (start_index, stop_index) pairs of the columns of A corresponding to that variable
    # var_size is the sum of the lengths of all variables in the problem
    # constr_size is the sum of the lengths of all constraints in the problem
    var_size, constr_size, var_to_ranges = find_variable_ranges(constraints)
    c = spzeros(var_size, 1)
    objective_range = var_to_ranges[objective_var_id]
    c[objective_range[1]:objective_range[2]] .= 1

    # slot in all of the coefficients in the conic forms into A and b
    A = spzeros(constr_size, var_size)
    b = spzeros(constr_size, 1)
    cones = Tuple{Symbol, UnitRange{Int}}[]
    constr_index = 0
    for constraint in constraints
        total_constraint_size = 0
        for i = 1:length(constraint.objs)
            sz = constraint.sizes[i]
            for (id, val) in constraint.objs[i]
                if id == objectid(:constant)
                    for l in 1:sz
                        b[constr_index + l] = val[1][l] == 0 ? val[2][l] : val[1][l]
                    end
                    #b[constr_index + sz + 1 : constr_index + 2*sz] = val[2]
                else
                    var_range = var_to_ranges[id]
                    if id_to_variables[id].sign == ComplexSign()
                        A[constr_index + 1 : constr_index + sz, var_range[1] : var_range[1] + length(id_to_variables[id])-1] = -val[1]
                        A[constr_index + 1 : constr_index + sz, var_range[1] + length(id_to_variables[id]) : var_range[2]] = -val[2]
                    else
                        A[constr_index + 1 : constr_index + sz, var_range[1] : var_range[2]] = -val[1]
                    end
                end
            end
            constr_index += sz
            total_constraint_size += sz
        end
        push!(cones, (constraint.cone, constr_index - total_constraint_size + 1 : constr_index))
    end

    # find integral and boolean variables
    vartypes = fill(:Cont, length(c))
    for var_id in keys(var_to_ranges)
        variable = id_to_variables[var_id]
        if :Int in variable.sets
            startidx, endidx = var_to_ranges[var_id]
            for idx in startidx:endidx
                vartypes[idx] = :Int
            end
        end
        if :Bin in variable.sets
            startidx, endidx = var_to_ranges[var_id]
            for idx in startidx:endidx
                vartypes[idx] = :Bin
            end
        end
    end

    if p.head == :maximize
        c = -c
    end

    return c, A, b, cones, var_to_ranges, vartypes, constraints
end

Problem(head::Symbol, objective::AbstractExpr, constraints::Constraint...) =
    Problem(head, objective, [constraints...])

# Allow users to simply type minimize
minimize(objective::AbstractExpr, constraints::Constraint...) =
    Problem(:minimize, objective, collect(constraints))
minimize(objective::AbstractExpr, constraints::Array{<:Constraint}=Constraint[]) =
    Problem(:minimize, objective, constraints)
minimize(objective::Value, constraints::Constraint...) =
    minimize(convert(AbstractExpr, objective), collect(constraints))
minimize(objective::Value, constraints::Array{<:Constraint}=Constraint[]) =
    minimize(convert(AbstractExpr, objective), constraints)

# Allow users to simply type maximize
maximize(objective::AbstractExpr, constraints::Constraint...) =
    Problem(:maximize, objective, collect(constraints))
maximize(objective::AbstractExpr, constraints::Array{<:Constraint}=Constraint[]) =
    Problem(:maximize, objective, constraints)
maximize(objective::Value, constraints::Constraint...) =
    maximize(convert(AbstractExpr, objective), collect(constraints))
maximize(objective::Value, constraints::Array{<:Constraint}=Constraint[]) =
    maximize(convert(AbstractExpr, objective), constraints)

# Allow users to simply type satisfy (if there is no objective)
satisfy(constraints::Constraint...) = Problem(:minimize, Constant(0), [constraints...])
satisfy(constraints::Array{<:Constraint}=Constraint[]) =
    Problem(:minimize, Constant(0), constraints)
satisfy(constraint::Constraint) = satisfy([constraint])

# +(constraints, constraints) is defined in constraints.jl
add_constraints!(p::Problem, constraints::Array{<:Constraint}) =
    +(p.constraints, constraints)
add_constraints!(p::Problem, constraint::Constraint) = add_constraints!(p, [constraint])
add_constraint! = add_constraints!

# caches conic form of x when x is the solution to the optimization problem p
function cache_conic_form!(conic_forms::UniqueConicForms, x::AbstractExpr, p::Problem)
    objective = conic_form!(p.objective, conic_forms)
    for c in p.constraints
        conic_form!(c, conic_forms)
    end
    cache_conic_form!(conic_forms, x, objective)
end


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

# Simple coupling interface of four elements

using FEMBase, MortarContact2D, Test
using MortarContact2D: get_slave_dofs, get_master_dofs,
                       get_mortar_matrix_D, get_mortar_matrix_M,
                       get_mortar_matrix_P

# Slave side coordinates:
Xs = Dict(
    1 => [0.0, 1.0],
    2 => [1.0, 1.0],
    3 => [2.0, 1.0])

# Master side coordinates:
Xm = Dict(
    4 => [0.0, 1.0],
    5 => [5/4, 1.0],
    6 => [2.0, 1.0])

# Define elements, update geometry:
X = merge(Xm , Xs)
es1 = Element(Seg2, [1, 2])
es2 = Element(Seg2, [2, 3])
em1 = Element(Seg2, [4, 5])
em2 = Element(Seg2, [5, 6])
elements = [es1, es2, em1, em2]
update!(elements, "geometry", X)

# Define new problem `Mortar2D`, coupling elements using mortar methods
problem = Problem(Mortar2D, "test interface", 1, "u")
@test_throws ErrorException add_elements!(problem, [es1, es2])
add_slave_elements!(problem, [es1, es2])
add_master_elements!(problem, [em1, em2])

# Assemble problem
assemble!(problem, 0.0)
s = get_slave_dofs(problem)
m = get_master_dofs(problem)
@test s == [1, 2, 3]
@test m == [4, 5, 6]

P = get_mortar_matrix_P(problem)
D = get_mortar_matrix_D(problem)
M = get_mortar_matrix_M(problem)

# From my thesis
D_expected = [1/3 1/6 0; 1/6 2/3 1/6; 0 1/6 1/3]
M_expected = [11/30 2/15 0; 41/160 13/20 3/32; 1/480 13/60 9/32]
P_expected = 1/320*[329 -24 15; 46 304 -30; -21 56 285]

# Results:

@test isapprox(D, D_expected)
@test isapprox(M, M_expected)
@test isapprox(P, P_expected)


text:
module FisherMetrics

using Manifolds
using ManifoldsBase
using Distributions
using Optim
using OrdinaryDiffEq
using BlockDiagonals

import Manifolds: local_metric
import ManifoldsBase: inner, log, log!

# Fisher information matrix for common distributions
include("fisher_information.jl")

# Fisher metric
include("fisher_metrics.jl")

export fim, FisherMetricManifold, default_basis

end

max_stars_repo_name: ap-browning/FisherMetrics

max_issues_repo_name: ap-browning/FisherMetrics

max_forks_repo_name: ap-browning/FisherMetrics


text:
## a type to indicate zero vector

immutable ZeroVector{T}
    len::Int
end

ZeroVector{T}(::Type{T}, n::Int) = ZeroVector{T}(n)

eltype{T}(v::ZeroVector{T}) = T
length(v::ZeroVector) = v.len
full(v::ZeroVector) = zeros(T, v.len)

convert{T}(::Type{Vector{T}}, v::ZeroVector{T}) = full(v)

+ (x::DenseArray, v::ZeroVector) = x
- (x::DenseArray, v::ZeroVector) = x
.+ (x::DenseArray, v::ZeroVector) = x
.- (x::DenseArray, v::ZeroVector) = x

##### Utility functions

type NoArgCheck end

isunitvec{T}(v::AbstractVector{T}) = (vecnorm(v) - 1.0) < 1.0e-12

function allfinite{T<:Real}(x::Array{T})
    for i = 1 : length(x)
        if !(isfinite(x[i]))
            return false
        end
    end
    return true
end

function allzeros{T<:Real}(x::Array{T})
    for i = 1 : length(x)
        if !(x == zero(T))
            return false
        end
    end
    return true
end

allzeros(x::ZeroVector) = true

function isprobvec(p::Vector{Float64})
    s = 0.
    for i = 1:length(p)
        pi = p[i]
        s += pi
        if pi < 0
            return false
        end
    end
    return abs(s - 1.0) <= 1.0e-12
end

function pnormalize!{T<:FloatingPoint}(v::AbstractVector{T})
    s = 0.
    n = length(v)
    for i = 1:n
        @inbounds s += v[i]
    end
    for i = 1:n
        @inbounds v[i] /= s
    end
    v
end

add!(x::DenseArray, y::DenseVector) = broadcast!(+, x, x, y)
add!(x::DenseVecOrMat, y::ZeroVector) = x

function multiply!(x::AbstractArray, c::Number)
    for i = 1:length(x)
        @inbounds x[i] *= c
    end
    x
end

function exp!(x::AbstractArray)
    for i = 1:length(x)
        @inbounds x[i] = exp(x[i])
    end
    x
end

# for checking the input range of quantile functions
# comparison with NaN is always false, so no explicit check is required
macro checkquantile(p,ex)
    :(zero($p) <= $p <= one($p) ? $ex : NaN)
end
macro checkinvlogcdf(lp,ex)
    :($lp <= zero($lp) ? $ex : NaN)
end

# because X == X' keeps failing due to floating point nonsense
function isApproxSymmmetric(a::Matrix{Float64})
    tmp = true
    for j in 2:size(a, 1)
        for i in 1:(j - 1)
            tmp &= abs(a[i, j] - a[j, i]) < 1e-8
        end
    end
    return tmp
end

# because isposdef keeps giving the wrong answer for samples
# from Wishart and InverseWisharts
hasCholesky(a::Matrix{Float64}) = isa(trycholfact(a), Cholesky)

function trycholfact(a::Matrix{Float64})
    try cholfact(a)
    catch e
        return e
    end
end


text:
"""
`SymmetricDoubleWell` - The classic symmetric doublewell potential.

### Fields
* `x` - Position x
"""
function SymmetricDoubleWell(x)
    return (x^2 -1)^2;
end

"""
`AsymmetricDoubleWell` - An asymmetric double well potential requires |δ|< 1 for
multiple minima.

### Fields
* `x` - Position x
### Optional Fields
* `δ = 0.5` - Asymmetry parameter
"""
function AsymmetricDoubleWell(x; δ = 0.5)
    return x^4 - 1.5 * x^2 - δ * x;
end

"""
`FatSkinnyDoubleWell10` - The fat-skinny double well of degree 10. This introduces
entropic effects in dimension 1.

### Fields
* `x` - Position x
"""
function FatSkinnyDoubleWell10(x)
    return ((8-5*x)^8 * (2+5*x)^2)/(2^26);
end

max_stars_repo_name: gideonsimpson/TestLandscapes

max_issues_repo_name: gideonsimpson/TestLandscapes

max_forks_repo_name: gideonsimpson/TestLandscapes