text:
using Knet
using Knet: Net, params, axpy!
if !isdefined(:MNIST)
    include(Pkg.dir("Knet/examples/mnist.jl"))
    mnist = minibatch(MNIST.xtrn, MNIST.ytrn, 100)
end

@knet function model02(x; winit=Gaussian(0,.1), hidden=32, f1=:relu, o...)
    h    = wf(x; out=hidden, f=f1, winit=winit)
    return wf(h; out=10, f=:soft, winit=winit)
end

# We'll try to create the two weight matrices from a single array

function initmodel{T}(w::CudaArray{T}; X=784, H=32, Y=10, B=100)
    @assert length(w) == X*H+H*Y
    f = compile(:model02; hidden=H)
    forw(f, zeros(T, X, B))
    f.reg[2].out0 = CudaArray(w.ptr, (H,X), w.dev)
    f.reg[5].out0 = CudaArray(w.ptr+H*X*sizeof(T), (Y,H), w.dev)
    return f
end

# Given a minibatch (x,y), eval model and return softloss
# Note that each minibatch gives a different loss function

evalmodel(f, x, y)=softloss(forw(f,x),y)

# We can find the right epsilon by finding the max radius at which 50%
# of the steps are positive.  For mnist, this radius seems to be about
# 0.01-0.05 starting at 0.

function deneme1(r; N=100)
    global w = fill!(CudaArray(Float32,794*32),0)
    global f = initmodel(w)
    (x,y) = first(mnist)
    @show f0 = softloss(forw(f,x),y)
    w0 = copy(w)
    global dw = copy(w)
    better=0
    for i=1:N
        scale!(r, randn!(dw))
        axpy!(1, dw, w)
        f1 = softloss(forw(f,x),y)
        f1 < f0 && (better += 1)
        copy!(w, w0)
    end
    better/N
end

# Simple algorithm:
# Repeat last step at increased size if success.
# Try random step at decreased size if failure.
#
# SGD gets to 0.36 in 1 epoch.
# This gets to around 0.60 sloss in 1000 epochs.
# using step=0.01 grow=1.001
# There is a long period where model is stuck around 0.70
#
# It turns out the grow factor is very important
# grow  sloss@64epoch
# 1.001	0.91
# 1.01	0.61
# 1.1	0.96 (though very high variance)
#
# (epoch,sloss,stepsize) using step=0.01 grow=1.01
# (epochs,step,grow) = (1024,0.01,1.01)
# (epochs,step,grow) = (1024,0.01,1.01)
# (1,1.7985317f0,0.0051531769028983964)
# (2,1.3769107f0,0.002341992948164361)
# (4,1.0358456f0,0.0011349164150302124)
# (8,0.8762574f0,0.00146262243400234)
# (16,0.7483307f0,0.0008924781269490016)
# (32,0.6375183f0,0.0008380797458788005)
# (64,0.5952995f0,0.0007391564809047644)
# (128,0.5699461f0,0.0003445230915309831)
# (256,0.59447426f0,0.0005493959444114676)
# (512,0.57121503f0,0.0002871048146276168)
# (1024,0.6157344f0,0.0008438423214864373)

function dhc01(; epochs=64, step=0.01, grow=1.01)
    @show (epochs,step,grow)
    global w = fill!(CudaArray(Float32,794*32),0)
    global f = initmodel(w)
    global dw = randn!(similar(w))
    global w0 = copy(w)
    nextn = 1
    for epoch in 1:epochs
        sloss = iter = ssize = 0
        for (x,y) in mnist
            f0 = softloss(forw(f,x),y)
            f1 = 0
            while true
                axpy!(step, dw, w)
                f1 = softloss(forw(f,x),y)
                if f1 < f0
                    sloss += f1; ssize += step; iter += 1
                    copy!(w0, w)
                    step *= grow
                    break
                else
                    copy!(w, w0)
                    step > 1e-6 && (step /= grow)
                    randn!(dw)
                end
            end
        end
        if epoch == nextn || epoch == epochs
            println((epoch,sloss/iter,ssize/iter))
            nextn*=2
        end
    end
end

# DHC1 is basically Procedure 3.1 from AITR1569.
# Let's try the original DHC algorithm, Procedure 4.1.

function dhc02(; epochs=64, step=0.01, grow=1.01, maxiter=5)
    @show (epochs,step,grow)
    global w = fill!(CudaArray(Float32,794*32),0)
    global w0 = copy(w)
    global u = fill!(similar(w), 0)
    global v = randn!(similar(w))
    global f = initmodel(w)
    nextn = 1
    for epoch in 1:epochs
        sloss = ssize = steps = 0
        for (x,y) in mnist
            f0 = softloss(forw(f,x),y)
            axpy!(step, v, w)
            f1 = softloss(forw(f,x),y)
            iter = 0
            while f1 > f0 && iter < maxiter
                iter += 1
                copy!(w, w0)
                randn!(v)
                axpy!(step, v, w)
                f1 = softloss(forw(f,x),y)
            end
            sloss += f1
            ssize += step
            steps += 1
            if f1 > f0
                step /= grow
            elseif iter == 0
                copy!(w0, w)
                axpy!(step, v, u)
                step *= grow
            else
                copy!(w0, w)    # save w+v
                axpy!(1, u, w)  # try w+v+u
                f2 = softloss(forw(f,x),y)
                if f2 > f1      # u+v did not work
                    copy!(w, w0)
                    fill!(u, 0) # reset memory
                    axpy!(step, v, u) # u is scaled v is not
                    step *= grow
                else
                    copy!(w0, w)
                    axpy!(step, v, u)
                    step *= grow
                end
            end
        end
        if epoch == nextn || epoch == epochs
            println((epoch,sloss/steps,ssize/steps))
            nextn*=2
        end
    end
end

# function train02(; epochs=1000, lr=0.5, adam=false, batchsize=100, seed=0, target=nothing, o...)
#     seed > 0 && setseed(seed)
#     global f = compile(:model02; o...)
#     global d = minibatch(MNIST.xtrn, MNIST.ytrn, batchsize)

#     # algorithm parameters
#     nhistory = 10                # remember last history steps
#     ihistory = 0
#     stepsize = 0                # dx = stepsize*gradient + (minstep + rndstep*stepsize)*random
#     minstep = 1e-4
#     rndstep = 0.1

#     for epoch=1:epochs
#         zloss = sloss = iter = 0
#         for (x,ygold) in d
#             ypred = forw(f, x)
#             sl = softloss(ypred, ygold)
#             zl = zeroone(ypred, ygold)
#             iter += 1; zloss += zl; sloss += sl

#         end
#         println((epoch, sloss/iter, zloss/iter))
#         zloss == 0 && break
#     end
#     return f
# end

# weights(f::Net)=map(r->r.out, params(f))
# gradients(f::Net)=map(r->r.dif, params(f))

# function vnorm(w)
#     sqnorm = 0
#     for u in w; sqnorm += vecnorm(u)^2; end
#     sqrt(sqnorm)
# end

# function vcos(u,v)
#     uv = 0
#     for i=1:length(u); uv += CUBLAS.dot(u[i],v[i]); end
#     uv / (vnorm(u)*vnorm(v))
# end

# function test02(f, data, loss)
#     sumloss = numloss = 0
#     for (x,ygold) in data
#         ypred = forw(f, x)
#         sumloss += loss(ypred, ygold)
#         numloss += 1
#     end
#     sumloss / numloss
# end

# f1 = train02(seed=1)
# f2 = train02(seed=1, target=f1)

# Can get to zloss=0 with hidden=24 using adam=true lr=0.001

#     setp(f, lr=lr, adam=adam)
            # # back(f, ygold, softloss)

            # if target != nothing
            #     # we start at w0 and want to reach w2
            #     # currently we are at w1 and about to take the step dw
            #     # compare dw and w2-w1
            #     # compare w2-w1 and w2-w0
            #     w1 = weights(f)
            #     dw = gradients(f)
            #     if w0 == nothing
            #         w0 = deepcopy(w1)
            #         w2 = weights(target)
            #         w2_w0 = deepcopy(w2)
            #         w2_w1 = deepcopy(w2)
            #         for i=1:length(w2); axpy!(-1,w0[i],w2_w0[i]); end
            #     end
            #     for i=1:length(w2); copy!(w2_w1[i],w2[i]); axpy!(-1,w1[i],w2_w1[i]); end
            #     @printf("w2-w0=%g,1 w2-w1=%g,%g dw=%g,%g\n",
            #             vnorm(w2_w0), vnorm(w2_w1), vcos(w2_w1, w2_w0),
            #             lr*vnorm(dw), -vcos(dw, w2_w1))
            # end

            # # update!(f)

# This doesn't work because initforw copies the arrays.
# @knet function model02(x; ftype=Float32, xinput=784, hidden=32, output=10,
#                        wvec=fill!(CudaArray(ftype,xinput*hidden+hidden*output),0))
#     h    = wf(x; f=:relu, winit=CudaArray(wvec.ptr, (hidden,xinput), wvec.dev))
#     return wf(h; f=:soft, winit=CudaArray(wvec.ptr+hidden*xinput*sizeof(ftype), (output,hidden), wvec.dev))
# end


text:
function scale_down_camera(camera, FACTOR)
    camera_modified = GL.CameraIntrinsics(
        width=round(Int,camera.width/FACTOR), height=round(Int,camera.height/FACTOR),
        fx=camera.fx/FACTOR, fy=camera.fy/FACTOR, cx=camera.cx/FACTOR, cy=camera.cy/FACTOR,
        near=camera.near, far=camera.far)
    camera_modified
end

export scale_down_camera

ext: jl


text:
struct AveragedPotential{M,T} <: AdiabaticModel
    models::M
    tmp_derivative::Matrix{T}
end

AveragedPotential(models, r) = AveragedPotential(models, zero(r))

NQCModels.ndofs(model::AveragedPotential) = NQCModels.ndofs(model.models[1])

function NQCModels.potential(model::AveragedPotential, r::AbstractMatrix)
    V = zero(eltype(r))
    for m in model.models
        V += NQCModels.potential(m, r)
    end
    return V / length(model.models)
end

function NQCModels.derivative!(model::AveragedPotential, D::AbstractMatrix, r::AbstractMatrix)
    fill!(D, zero(eltype(D)))
    for m in model.models
        NQCModels.derivative!(m, model.tmp_derivative, r)
        D .+= model.tmp_derivative
    end
    D ./= length(model.models)
    return D
end


text:
using SHA
using Random

block_string = "test5"
block_string = randstring(10)
#println(bytes2hex(sha256("test")))

max_nonce = 2 ^ 32 # 4 billion
max_nonce = 2 ^ 16

is_hash_found = 0
for nonce in 0:max_nonce
    hash = bytes2hex(sha2_256(string(nonce) * block_string))
    #println(hash)
    if startswith(hash, "0000")
        global is_hash_found = 1
        println(hash)
        break
    end
end

if is_hash_found == 0
	println("no proper hash found, please try again later")
end

max_stars_repo_name: marchunter/litcoin

max_issues_repo_name: marchunter/litcoin

max_forks_repo_name: marchunter/litcoin


text:
module DielectricFunctions

export piecewise_constant_ε,
piecewise_constant_F

"""
    piecewise_constant_ε(x,y,params)
"""
function piecewise_constant_ε(x,y,params::Dict)
    n2 = get!(params,:n₂,0)
    n1 = get!(params,:n₁,1)
    return complex(n1,n2)^2
end

"""
    piecewise_constant_F(x,y,params)
"""
piecewise_constant_F(x,y,params) = get!(params,:F,0.0)

end # module


text:
# -*- coding: utf-8 -*-
# # Polynomial Ring of Transformations
#
# One goal of SemanticModels is to lift model augmentation into the world of category theory and algebra.
# This is analogous to a long tradition in the *symbolification* of mathematics, for example geometry is the study of shapes and their relations, which gets lifted to group theory by the algebraists who want to to make everything symbolic reasoning.
#
# He we will represent a model class *monomial regression* with a data structure `MonomialRegression` and define a set group of transformations that can act on models from this class.

using Polynomials
using SemanticModels
using SemanticModels.ModelTools
using SemanticModels.ModelTools.Transformations
using SemanticModels.ModelTools.MonomialRegressionModels

# ## Monomial Regression
#
# Here we take our class of models to be monomial regression and write a script that implements an experiment with this model.
#
# 1. Define a function $f(a, x) = ax^p$ for some integer $p$
# 2. Generate data from the distribution $x, f(a, x)+Normal(0,\sigma)$
# 3. Given the sample data, use least squares to solve for $a$.
#
# Note that the value of $p$ when we generate the data does not need to match the value of $p$ when we solve the least squares problem.

# Here is an implementation of our model.

expr = quote
    module Regression
    using Random
    mid(a,b) = (a+b)/2
    function f(x, a)
        y = first(a)*x^0
        return y
    end

    function optimize(l, interval)
        left = interval[1]
        right = interval[2]
        if left > right - 1e-8
            return (interval)
        end
        midp = (left+right)/2
        if l(mid(left, midp)) > l(mid(midp, right))
            return optimize(l, (midp, right))
        else
            return optimize(l, (left, midp))
        end
    end

    function sample(g::Function, n::Int)
        x = randn(Float64, n)
        target = g.(x) .+ randn(Float64, n)./8
        return x, target
    end

    #setup

    Random.seed!(42)
    a = 1/2
    n = 10
    g(x) = a*(1 + x^1 + x^2)
    x, target = sample(g, n)
    loss(a) = sum((f.(x, a) .- target).^2)
    a₀ = [-1, 1]
    # @show loss.([-1,-1/2,-1/4, 0, 1/4,1/3,1/2,2/3, 1])

    #solving
    ahat = optimize(loss, a₀)

    end
end;

# ## Transforming the model
#
# Valid transformations for this model include increasing and decreasing the power of the monomial by 1. And you can apply those transformations repeatedly.
#
# Mathematically we can go from $f(a,x) = ax^p$ to $f(a,x) = ax^{p+1}$ or $f(a,x) = ax^{p-1}$. If we allow repeated application of these two generator transformations as a free monoid, we find that the set of transformations we get is isomorphic to the *group* $(Z, +, 0)$ which is the integers with addition. In Julia, we have decided to represent elements from this transformation group with a struct `Pow{Int} <: Transformation`. The composition operation for these transformations can be found in the source file `groups.jl`.
#
# Composition of `Pow{Int}` maps naturally to addition of the power to transform with. You can see how transformations are applied below.

m = model(MonomialRegressionModel, deepcopy(expr))
@show bodyblock(m.f)

# With this machinery in place we can think about the action of $Z$ on a `MonomialRegression` model to create a family of models with different powers of monomoial. Each monomial model finds the best fit to a particular data set, and by leveraging this group action, we can systematically explore the space of all models in the class.

eval(m.expr);
â = Regression.ahat

function (p::Poly)(m::MonomialRegressionModel)
    if m.expr.head == :module
        # we need to add the using statement to get polyval defined
        pushfirst!(m.expr.args[3].args, :(using Polynomials))
    end
    bodyblock(m.f)[end-2] = Parsers.replace(bodyblock(m.f)[end-2], :(x^0), :(polyval($p,x)))
end
Poly([1,0,1])(m)
m

eval(m.expr);
â = Regression.ahat

# ## Chosing the right polynomial
# Given a fixed dataset, we want to loop over all the models in a class and fit the best coefficient, then we will decide which polynomial allows for the best fit to this data.

m = model(MonomialRegressionModel, deepcopy(expr))
# println(m)
# sol = eval(m.expr)
# @show sol.ahat
# @show sol.loss(sol.ahat[1])
results = Any[]
for i in 0:1
    for j in 0:1
        for k in 0:2
            m = model(MonomialRegressionModel, deepcopy(expr))
            (i*Poly([1]) + j*Poly([0, 1]) + k*Poly([0,0,1]))(m)
            @show m.f.args[2].args[2]
            sol = eval(m.expr)
            p = bodyblock(m.f)
            ahat = @show sol.ahat[1]
            lhat = @show sol.loss(sol.ahat[1])
            push!(results, (i, ahat, lhat, p))
        end
    end
end

using Printf
fmt(i::Int) = i
fmt(f::Real) = Printf.@sprintf "% .4f" round(f, sigdigits=5)
fmt(f::Any) = f[end-2].args[2].args[3].args[2]
println("\nResults:\n\n\ni\tâ         ℓ⋆ \t\tpoly\n----------------------------------------------")
for r in results
    println(join(fmt.(r), "     "))
end
best = sort(results, by=x->x[end-1])[1][2:end]
println("Model order $(fmt(best[end])) is the best with a=$(fmt(best[1])) and loss $(fmt(best[2]))")

# ## Conclusions
#
# We can see that this sweep over polynomials, recovers the right model. We can see some fundamental algebra at work here.
#
#
# This example shows that you can represent the transformations on a model as a group and the action of the group on a model is the way to represent transformations on models. Insights from the algebraic structures homormorphic/isomorphic to the algebra of transformations provide insights on the class of models.

div(Poly([1,2,3])*Poly([0,1]), Poly([0,1]))

divrem((Poly([1,2,3])*Poly([0,1,2])), Poly([1,1,2,1]))

gcd((Poly([1,2,3])*Poly([1,3,4])), Poly([1,2,3]))


text:
#=
Aodwt.jl
2019-02-23 Jeff Fessler, University of Michigan
=#

export Aodwt

using Plots
using LinearMapsAA: LinearMapAA
using Wavelets: dwt, idwt, wavelet, WT
#using MIRT: jim

"""
    A, levels, mfun = Aodwt(dims ; level::Int=3, wt=wavelet(WT.haar))

Create orthogonal discrete wavelet transform (ODWT) `LinearMapAA`

in
- `dims::Dims` tuple of dimensions

option
- `level::Int` # of levels; default 3
- `wt` wavelet transform type (see `Wavelets` package); default Haar

out
- `A` a `LinearMapAA` object
- `scales` array of size `dims` showing the scale of each coefficient
which is useful when imposing scale-dependent regularization
- `mfun` convenience function for A*X when X is a Matrix or Array (not vector)

2019-02-23 Jeff Fessler, University of Michigan
"""
function Aodwt(dims::Dims ; level::Int=3, wt=wavelet(WT.haar))

	Afun = (level) -> LinearMapAA(
		x -> vec(dwt(reshape(x, dims), wt, level)),
		y -> vec(idwt(reshape(y, dims), wt, level)),
		(prod(dims), prod(dims)), (wt=wt, level=level))

	A = Afun(level)

	mfun = (A, x) -> reshape(A * vec(x), dims)

	scales = zeros(dims)
	for il=1:level
		Al = Afun(il)
		tmp = mfun(Al, ones(dims)) .== 0
		scales += il * (tmp .& (scales .== 0))
	end

	return A, scales, mfun
end

function Aodwt()
	@doc Aodwt
end

"""
    Aodwt_show( ; M::Int=32, N::Int=64)
show scales
"""
function Aodwt_show( ; dims::Dims = (32, 64))
	W, scales, mfun = Aodwt(dims)
	jim(scales)
end

"""
    Aodwt(:test)
self test

`Aodwt(:show)`
visualize
"""
function Aodwt(test::Symbol)
	if test === :show
		return Aodwt_show()
	end

	test != :test && throw(ArgumentError("test $test"))
	Aodwt()

	W,_,_ = Aodwt((8,16), level=2)
	@test Matrix(W)' == Matrix(W') # check adjoint
	isinteractive() && (@show W.wt)
	isinteractive() && (@show propertynames(W.wt))
	@test W.level == 2
	@test W.wt.name == "haar"

	Aodwt(:show)
	true
end


text:
struct PingFirst{P<:Policy} <: Policy
    p::P
end

function POMDPs.action(p::PingFirst, b::AbstractParticleBelief)
    if first(particles(b)).aware
        return action(p.p, b)
    else
        return PING
    end
end


text:
export SoftPlus

"""
`SoftPlus([domainType=Float64::Type,] dim_in::Tuple)`

Creates the softplus non-linear operator with input dimensions `dim_in`.
```math
\\sigma(\\mathbf{x}) = \\log (1 + e^{x} )
```

"""
struct SoftPlus{T,N} <: NonLinearOperator
	dim::NTuple{N,Int}
end

function SoftPlus(DomainType::Type, DomainDim::NTuple{N,Int}) where {N}
	SoftPlus{DomainType,N}(DomainDim)
end

SoftPlus(DomainDim::NTuple{N,Int}) where {N} = SoftPlus{Float64,N}(DomainDim)

function mul!(y::AbstractArray{T,N}, L::SoftPlus{T,N}, x::AbstractArray{T,N}) where {T,N}
	y .= log.(1 .+exp.(x))
end

function mul!(y::AbstractArray,
              J::AdjointOperator{Jacobian{A,TT}},
              b::AbstractArray) where {T, N, A <: SoftPlus{T,N}, TT <: AbstractArray{T,N} }
    L = J.A
	y .= 1 ./(1 .+exp.(-L.x)).*b
end

fun_name(L::SoftPlus) = "σ"

size(L::SoftPlus) = (L.dim, L.dim)

domainType(L::SoftPlus{T,N}) where {T,N} = T
codomainType(L::SoftPlus{T,N}) where {T,N} = T


text:
@mlj_model mutable struct APIx0 <: Model
    f0::Int
end
@mlj_model mutable struct APIx0b <: Model
    f0::Int
end

mutable struct APIx1 <: Model end

@testset "selectrows(model, data...)" begin
    X = (x1 = [2, 4, 6],)
    y = [10.0, 20.0, 30.0]
    @test selectrows(APIx0(), 2:3, X, y) == ((x1 = [4, 6],), [20.0, 30.0])
end

@testset "fit-x" begin
    m0 = APIx0(f0=1)
    m1 = APIx0b(f0=3)
    # no weight support: explicit fallback
    M.fit(m::APIx0, v::Int, X, y) = (5, nothing, nothing)
    M.fit(m::APIx0, v::Int, X, y, w) = (5, nothing, nothing)
    @test fit(m0, 1, randn(2), randn(2), 5) == (5, nothing, nothing)
    # with weight support: use
    M.fit(m::APIx0b, v::Int, X, y, w) = (7, nothing, nothing)
    @test fit(m1, 1, randn(2), randn(2), 5) == (7, nothing, nothing)
    # default fitted params
    @test M.fitted_params(m1, 7) == (fitresult=7,)
    # default iteration_parameter
    @test M.training_losses(m0, nothing) === nothing
    #update fallback = fit
    @test update(m0, 1, 5, nothing, randn(2), 5) == (5, nothing, nothing)
    # training losses:
    f, c, r = MLJModelInterface.fit(m0, 1, rand(2), rand(2))
    @test M.training_losses(m0, r) === nothing
end

struct DummyUnivariateFinite end

mutable struct UnivariateFiniteFitter <: Probabilistic end

@testset "models fitting a distribution to data" begin
    MMI = MLJModelInterface

    function MMI.fit(model::UnivariateFiniteFitter, verbosity::Int, X, y)
        fitresult = DummyUnivariateFinite()
        report = nothing
        cache = nothing

        verbosity > 0 && @info "Fitted a $fitresult"

        return fitresult, cache, report
    end

    function MMI.predict(model::UnivariateFiniteFitter, fitresult, X)
        return fill(fitresult, length(X))
    end

    MMI.input_scitype(::Type{<:UnivariateFiniteFitter}) = Nothing

    MMI.target_scitype(::Type{<:UnivariateFiniteFitter}) = AbstractVector{<:Finite}

    y = categorical(collect("aabbccaa"))
    X = nothing
    model = UnivariateFiniteFitter()
    fitresult, cache, report = MMI.fit(model, 1, X, y)

    @test cache === nothing
    @test report === nothing

    ytest = y[1:3]
    yhat = predict(model, fitresult, fill(nothing, 3))
    @test yhat == fill(DummyUnivariateFinite(), 3)

end


text:
@testset "util" begin

@testset "argmax, argmin" begin
    v = Float64[11, 13, 14, 10, 12]

    # Test the behavior for a nonempty index vector.
    indv = [1,2,5]
    @test @inferred(argmax(identity, v)) == 3
    @test @inferred(argmax(identity, v, indv)) == 2
    @test @inferred(argmin(identity, v)) == 4
    @test @inferred(argmin(identity, v, indv)) == 1

    # Test the behavior for an empty index vector.
    indv = Int[]
    @test @inferred(argmax(identity, v, indv)) == 0
    @test @inferred(argmin(identity, v, indv)) == 0

end  # @testset "argmax, argmin"

end  # @testset "util"


text:
#=
lower bound on relative entropy of entanglement (PPT relaxation)
adapted from
https://github.com/hfawzi/cvxquad/blob/master/examples/rel_entr_entanglement.m
=#

struct RelEntrEntanglementJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    na::Int
    nb::Int
end

function build(inst::RelEntrEntanglementJuMP{T}) where {T <: Float64}
    (na, nb) = (inst.na, inst.nb)
    side = na * nb
    Rho = randn(T, side, side)
    Rho = Rho * Rho'
    Rho = Symmetric(Rho / tr(Rho))
    vec_dim = Cones.svec_length(side)
    rho_vec = zeros(T, vec_dim)
    Cones.smat_to_svec!(rho_vec, Rho, sqrt(T(2)))

    model = JuMP.Model()
    JuMP.@variable(model, tau_vec[1:vec_dim])
    Tau = zeros(JuMP.AffExpr, side, side)
    Cones.svec_to_smat!(Tau, one(T) * tau_vec, sqrt(T(2)))
    JuMP.@constraint(model, tr(Tau) == 1)

    JuMP.@variable(model, y)
    JuMP.@objective(model, Min, y / log(T(2)))
    JuMP.@constraint(model, vcat(y, tau_vec, rho_vec) in
        Hypatia.EpiTrRelEntropyTriCone{T}(1 + 2 * vec_dim))
    pt = partial_transpose(Symmetric(Tau), 2, [na, nb])
    JuMP.@SDconstraint(model, Symmetric(pt) >= 0)

    return model
end

# partial transpose of Q over system i given subsystem dimensions subs
function partial_transpose(
    Q::Symmetric,
    i::Int,
    subs::Vector{Int},
    )
    @assert 1 <= i <= length(subs)
    @assert size(Q, 1) == prod(subs)
    n = length(subs)
    s = n + 1 - i
    p = collect(1:2n)
    p[s] = n + s
    p[n + s] = s
    rev = reverse(subs)
    permQ = permutedims(reshape(Q, (rev..., rev...)), p)
    side = prod(subs)
    return reshape(permQ, (side, side))
end


text:
function barzanCurveFit(lockType, trainX, trainY, lowConf, upConf, minStep, maxSteps)
	totalSteps = 0;
	conf = 0;
	if !(size(lowConf, 1) == 1 && size(upConf, 1) == 1 && size(minStep, 1) == 1 && size(maxSteps, 1) == 1)
		error("Config should be a row-vector");
	end
	if !(size(trainY,2)==1)
		error("We only accept 1-D as target values");
	end

	lowVal = implicitUseLockModel(lockType,lowConf, trainX);
	upVal = implicitUseLockModel(lockType,upConf, trainX);

	for col=1:size(lowConf,2)
		steps = 0;
		while abs(upConf[col]-lowConf[col]) > minStep[col] && steps <= maxSteps[col]
			steps = steps+1;
			totalSteps = totalSteps+1;

			lowErr = mean(abs(lowVal-trainY));
			if isnan(lowErr)
				lowConf[col] = (lowConf[col] + upConf[col])/2;
				lowVal = implicitUseLockModel(lockType,lowConf, trainX);
				continue;
			end

			upErr = mean(abs(upVal-trainY));
			if isnan(upErr)
				upConf[col] = (lowConf[col] + upConf[col])/2;
				upVal = implicitUseLockModel(lockType,upConf, trainX);
				continue;
			end

			if lowErr>upErr
				lowConf[col] = (lowConf[col] + upConf[col])/2;
				lowVal = implicitUseLockModel(lockType,lowConf, trainX);
			else
				upConf[col] = (lowConf[col] + upConf[col])/2;
				upVal = implicitUseLockModel(lockType,upConf, trainX);
			end
			println(string("steps=",steps,", low=",valueToString(lowConf)", up=",valueToString(upConf),", lowErr=",lowErr,", upErr=",upErr));
		end

		if (steps > maxSteps[col])
			println(string("Reached the maximum number of iterations: ", steps));
		else
			println(string("Converged after ", steps, " iterations for this column"));
		end

		conf = (upConf + lowConf)/2;
	end

	println(string("Overall number of steps: ", totalSteps, " iterations"));
	return conf
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
module Monotonicers

using Dates
using DataFrames
using Random
using Statistics

export fit!,transform!,ismonotonic,dailyflips
export Monotonicer

export monotonicerrun

using TSML.TSMLTypes
import TSML.TSMLTypes.fit! # to overload
import TSML.TSMLTypes.transform! # to overload
using TSML.Utils

# Transforms instances with nominal features into one-hot form
# and coerces the instance matrix to be of element type Float64.
mutable struct Monotonicer <: Transformer
  model
  args

  function Monotonicer(args=Dict())
    default_args = Dict(
    )
    new(nothing, mergedict(default_args, args))
  end
end

function fit!(st::Monotonicer, features::T, labels::Vector=[]) where {T<:Union{Vector,Matrix,DataFrame}}
  typeof(features) <: DataFrame || error("Monotonicer.fit!: data should be a dataframe: Date,Val ")
  ncol(features) == 2 || error("dataframe must have 2 columns: Date, Val")
  st.model = st.args
end

function transform!(st::Monotonicer, features::T) where {T<:Union{Vector,Matrix,DataFrame}}
  features != [] || return DataFrame()
  typeof(features) <: DataFrame || error("Monotonicer.fit!: data should be a dataframe: Date,Val ")
  ncol(features) == 2 || error("dataframe must have 2 columns: Date, Val")
  sum(names(features) .== (:Date,:Value))  == 2 || error("wrong column names")
  mfeatures=features
  # double check monotonic
  # based on flips and daily flips
  if ismonotonic(features[:Value])
    mfeatures = antimonotonize(features)
  end
  # daily mono reset everynight at most 2x
  ndailyflips = dailyflips(mfeatures)
  if ndailyflips < 0.5
    antimono = antimonotonize(mfeatures)
  elseif 0.5 <= ndailyflips <= 2.0
    antimono = antimonotonizedaily(mfeatures)
  else
    antimono = mfeatures
  end
  return antimono
end

function antimonotonizedaily(dat::DataFrame)
  # get diffs, negative values represent reset
  # replace negative values by succeeding neighbor
  df = deepcopy(dat)
  df[:Value] = [ -1.0; diff(dat[:Value]) ]
  valuelength = length(df[:Value])
  # find locations of negatives
  negatives = findall(x->x<0,df[:Value])
  # if last element is negative, zero it
  if negatives[end] == valuelength
    df[:Value][negative[end]] = 0.0
    pop!(negatives)
  end
  # replace negatives with immediate neighbor
  for val in reverse(negatives)
    df[:Value][val] = df[:Value][val+1]
  end
  return df
end

function antimonotonize(data::Union{Vector,DataFrame,Matrix})
  dat = deepcopy(data)
  typeof(dat) <: DataFrame || error("input must be a dataframe")
  ncol(dat) == 2 || error("input must have two columns")
  sum(names(dat) .== (:Date,:Value))  == 2 || error("wrong column names")
  vals = dat[:Value]
  fvals = diff(vals)
  dvals = [fvals[1];fvals]
  dat[:Value] = dvals
  return dat
end

function dailyflips(dat::DataFrame)
  dates = dat[:Date]
  values = dat[:Value]
  datemin = minimum(dates)
  datemax = maximum(dates)
  ndays = round(datemax-datemin,Dates.Day(1))
  nflips::Float64 = sum((skipmissing(values) |> collect |> diff) .< 0 )
  avg::Float64 = nflips/ndays.value
end

function ismonotonic(dat::Vector{T}) where T
  #check if monotonic increasing
  cdata=dat |> skipmissing |> collect
  maxflips = minimum([10;div(length(cdata),2)]) # flip cutoff to identify monotonic
  nflips = length(cdata) - sum(cdata .== accumulate(max,cdata))
  nflips < maxflips ? true : false
end

end


text:
## Julia-level access of parts (arguments) of Mxpr's by name.
#
# First, recall that Mxpr expressions have the form h(a,b,c,d...)
# Parts can be accessed from the Julia level by eg  h[1] --> a
# They can also be accesed by eg margs(h(a,b,c,d))[1]
# The arguments (parts) have varying semantic interpretation depending on the head 'h'.
# We can get cheap, consistent, semantic access to the parts using multiple dispatch.
# For example if p is an instance of Mxpr{:Power}, we can write base(p) to get the base.
# This works for getting, but not for setting yet.

# There may be a performance penalty, in some cases, for getting the
# args field of the Mxpr each time named access is used, rather
# getting the args once and then indexing directly into them.

base(p::Mxpr{:Power}) = p[1]
expt(p::Mxpr{:Power}) = p[2]
exponent(p::Mxpr{:Power}) = p[2]

expr(p::Mxpr{:Part}) = p[1]
indices(p::Mxpr{:Part}) = (a = margs(p); view(a,2:length(a)))

# We used these a bit
terms(s::Mxpr{:Plus}) = margs(s)
factors(p::Mxpr{:Times}) = margs(p)

lhs(r::RulesT) = r[1]
rhs(r::RulesT) = r[2]


text:
using DataDeps

register(DataDep(
  "image-segmentation",
  "http://archive.ics.uci.edu/ml/datasets/Image+Segmentation",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/image/segmentation.data",
  "ed5b5aee8081bbe875357149f73f50d7085e22cc0897f5b52eb7b52fea969a5b",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:20,
        skipto=6,
        delim=",",
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# Autogenerated wrapper script for Xorg_xcb_util_keysyms_jll for aarch64-linux-gnu
export libxcb_keysyms

using Xorg_xcb_util_jll
JLLWrappers.@generate_wrapper_header("Xorg_xcb_util_keysyms")
JLLWrappers.@declare_library_product(libxcb_keysyms, "libxcb-keysyms.so.1")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_xcb_util_jll)
    JLLWrappers.@init_library_product(
        libxcb_keysyms,
        "lib/libxcb-keysyms.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "Monte Carlo Tests" begin
    @test size(mc_supersonic(10)) == (10, 7)
end;


text:
import DPFEHM
import GaussianRandomFields
import Optim
import PyPlot
import Random
import Zygote

Random.seed!(0)

function getobsnodes(coords, obslocs)
	obsnodes = Array{Int}(undef, length(obslocs))
	for i = 1:length(obslocs)
		obsnodes[i] = findmin(map(j->sum((obslocs[i] .- coords[:, j]) .^ 2), 1:size(coords, 2)))[2]
	end
	return obsnodes
end

mins = [0, 0]#meters
maxs = [100, 100]#meters
allns = [[101, 101], [201, 201], [401, 401], [801, 801]]
allns = [[allns[end]]; allns]
num_eigenvectors = 200
x_true = randn(num_eigenvectors)
x0 = zeros(num_eigenvectors)
sqrtnumobs = 16
obslocs_x = range(mins[1], maxs[1]; length=sqrtnumobs + 2)[2:end - 1]
obslocs_y = range(mins[2], maxs[2]; length=sqrtnumobs + 2)[2:end - 1]
obslocs = collect(Iterators.product(obslocs_x, obslocs_y))[:]
observations = Array{Float64}(undef, length(obslocs))
for (i, ns) in enumerate(allns)
	Random.seed!(0)
	global x_true
	global x0
	global observations
	coords, neighbors, areasoverlengths, _ = DPFEHM.regulargrid2d(mins, maxs, ns, 1.0)
	Qs = zeros(size(coords, 2))
	boundaryhead(x, y) = 5 * (x - maxs[1]) / (mins[1] - maxs[1])
	dirichletnodes = Int[]
	dirichleths = zeros(size(coords, 2))
	for i = 1:size(coords, 2)
		if coords[1, i] == mins[1] || coords[1, i] == maxs[1] || coords[2, i] == mins[2] || coords[2, i] == maxs[2]
			push!(dirichletnodes, i)
			dirichleths[i] = boundaryhead(coords[1:2, i]...)
		end
	end

	lambda = 50.0#meters -- correlation length of log-conductivity
	sigma = 1.0#standard deviation of log-conductivity
	mu = 0.0#mean of log-permeability can be arbitrary because this is steady-state and there are no fluid sources

	cov = GaussianRandomFields.CovarianceFunction(2, GaussianRandomFields.Matern(lambda, 1; σ=sigma))
	x_pts = range(mins[1], maxs[1]; length=ns[1])
	y_pts = range(mins[2], maxs[2]; length=ns[2])
	@time grf = GaussianRandomFields.GaussianRandomField(cov, GaussianRandomFields.KarhunenLoeve(num_eigenvectors), x_pts, y_pts)
	@time logKs = GaussianRandomFields.sample(grf)
	parameterization = copy(grf.data.eigenfunc)
	sigmas = copy(grf.data.eigenval)

	if i == 1#only plot this the first time around
		#plot a realization
		fig, ax = PyPlot.subplots()
		ax.imshow(logKs)
		ax.title.set_text("Random Conductivity Field")
		display(fig)
		println()
		PyPlot.close(fig)

		#plot the first few eigenvectors and some random ones
		fig, axs = PyPlot.subplots(2, 4, figsize=(16, 8))
		rp = Random.randperm(size(parameterization, 2))
		for (i, ax) in enumerate(axs)
			if i <= length(axs) / 2
				ax.imshow(reshape(parameterization[:, i], ns...))
				ax.title.set_text("Eigenvector $i")
			else
				ax.imshow(reshape(parameterization[:, rp[i]], ns...))
				ax.title.set_text("Eigenvector $(rp[i])")
			end
		end
		display(fig)
		println()
		PyPlot.close(fig)
	end

	#logKs2Ks_neighbors(Ks) = map(neighbor->exp(0.5 * (Ks[neighbor[1]] + Ks[neighbor[2]])), neighbors)
	#logKs2Ks_neighbors(Ks) = exp.([0.5 * (Ks[i1] + Ks[i2]) for (i1, i2) in neighbors])
	logKs2Ks_neighbors(Ks) = exp.(0.5 * (Ks[map(p->p[1], neighbors)] .+ Ks[map(p->p[2], neighbors)]))#Zygote differentiates this efficiently but the definitions above are ineffecient with Zygote
	function solveforh(logKs, dirichleths)
		@assert length(logKs) == length(Qs)
		if maximum(logKs) - minimum(logKs) > 25
			return fill(NaN, length(Qs))#this is needed to prevent the solver from blowing up if the line search takes us somewhere crazy
		else
			Ks_neighbors = logKs2Ks_neighbors(logKs)
			return DPFEHM.groundwater_steadystate(Ks_neighbors, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs)
		end
	end
	solveforheigs(x) = solveforh(x2logKs(x), dirichleths)
	x2logKs(x) = reshape(parameterization * (sigmas .* x), ns...)
	logKs_true = x2logKs(x_true)
	h_true = solveforheigs(x_true)

	numobs = 250
	obsnodes = getobsnodes(coords, obslocs)
	obssigma = 1e-3
	observations .= h_true[obsnodes]#set up the observations
	f(logKs) = sum((solveforh(logKs, dirichleths)[obsnodes] - observations) .^ 2 ./ obssigma ^ 2)
	h = solveforh(logKs, dirichleths)
	h = solveforh(logKs, dirichleths)
	print("forward solve time")
	@time h = solveforh(logKs_true, dirichleths)
	zg = Zygote.gradient(f, logKs)[1]
	zg = Zygote.gradient(f, logKs)[1]#work out the precompilation
	print("gradient time")
	@time zg = Zygote.gradient(f, logKs)[1]

	#now reformulate the function in terms of the cofficients of the eigenvectors with some regularization in there
	feigs(x) = f(x2logKs(x)) + sum(x .^ 2)
	print("eigs forward solve time")
	@time feigs(x_true)
	zgeigs = Zygote.gradient(feigs, x_true)[1]
	print("eigs gradient time")
	@time zgeigs = Zygote.gradient(feigs, x_true)[1]

	if i == 1#first time around plot a bunch of stuff and set up the observations
		#plot the solution, the difference between the solution and the solution for a uniform medium, and the logKs
		fig, axs = PyPlot.subplots(1, 4, figsize=(32, 8))
		ims = axs[1].imshow(reshape(h, ns...))
		axs[1].title.set_text("Head for Random Conductivity")
		fig.colorbar(ims, ax=axs[1])
		ims = axs[2].imshow(reshape(h .- map(i->boundaryhead(coords[:, i]...), 1:size(coords, 2)), ns...))
		axs[2].title.set_text("Head Fluctuation for Random Conductivity")
		fig.colorbar(ims, ax=axs[2])
		ims = axs[3].imshow(logKs)
		axs[3].title.set_text("Random Conductivity Field")
		fig.colorbar(ims, ax=axs[3])
		ims = axs[4].imshow(zg)
		axs[4].title.set_text("Gradient of Loss w.r.t. Conductivity")
		fig.colorbar(ims, ax=axs[4])
		fig.tight_layout()
		display(fig)
		println()
		PyPlot.close(fig)

		#plot the solution, the difference between the solution and the solution for a uniform medium, and the logKs
		fig, axs = PyPlot.subplots(1, 4, figsize=(32, 8))
		ims = axs[1].imshow(reshape(h_true, ns...))
		axs[1].title.set_text("Head for Random Conductivity")
		fig.colorbar(ims, ax=axs[1])
		ims = axs[2].imshow(reshape(h_true .- map(i->boundaryhead(coords[:, i]...), 1:size(coords, 2)), ns...))
		axs[2].title.set_text("Head Fluctuation for Random Conductivity")
		fig.colorbar(ims, ax=axs[2])
		ims = axs[3].imshow(x2logKs(x_true))
		axs[3].title.set_text("True Conductivity Field")
		fig.colorbar(ims, ax=axs[3])
		ims = axs[4].plot(zgeigs)
		axs[4].title.set_text("Gradient of Loss w.r.t. Eigenvector Coefficients")
		fig.tight_layout()
		display(fig)
		println()
		PyPlot.close(fig)

	else#otherwise, do the optimization
		options = Optim.Options(iterations=200, extended_trace=false, store_trace=true, show_trace=ifelse(i == length(allns), true, false), x_tol=1e-6, time_limit=60 * 60 * 2)
		print("Optimization time")
		@time opt = Optim.optimize(feigs, x->Zygote.gradient(feigs, x)[1], x0, Optim.LBFGS(), options; inplace=false)
		#@time opt = Optim.optimize(feigs, x->Zygote.gradient(feigs, x)[1], x0, Optim.GradientDescent(), options; inplace=false)
		x_est = opt.minimizer
		logKs_est = x2logKs(x_est)

		#plot results from the inverse analysis
		fig, axs = PyPlot.subplots(2, 3, figsize=(12, 8))
		ims = axs[1].imshow(logKs_true)
		axs[1].title.set_text("True Conductivity")
		fig.colorbar(ims, ax=axs[1])
		ims = axs[2].imshow(logKs_est)
		axs[2].title.set_text("Estimated Conductivity")
		fig.colorbar(ims, ax=axs[2])
		ims = axs[3].imshow(reshape(solveforh(logKs_true, dirichleths), ns...))
		axs[3].title.set_text("True Head")
		#ims = axs[3].imshow(x2logKs(x0))
		#axs[3].title.set_text("Previous Conductivity Estimate")
		fig.colorbar(ims, ax=axs[3])
		#ims = axs[4].imshow(logKs_est - x2logKs(x0))
		#axs[4].title.set_text("Conductivity Update")
		ims = axs[4].imshow(reshape(solveforh(logKs_est, dirichleths), ns...))
		axs[4].title.set_text("Estimated Head")
		fig.colorbar(ims, ax=axs[4])
		axs[5].plot([0, 5], [0, 5], "k", alpha=0.5)
		axs[5].plot(observations, solveforh(logKs_est, dirichleths)[obsnodes], "k.")
		axs[5].set_xlabel("Observed Head")
		axs[5].set_ylabel("Predicted Head")
		axs[6].semilogy(map(i->opt.trace[i].value, 1:length(opt.trace)))
		axs[6].set_xlabel("Iteration")
		axs[6].set_ylabel("Loss")
		fig.tight_layout()
		fig.savefig("figs/level_$(i - 1).pdf")
		display(fig)
		println()
		PyPlot.close(fig)

		x0 = x_est#update x0 for the next optimization
	end
end


text:
using Documenter, Bio

makedocs()
deploydocs(
    deps = Deps.pip("mkdocs", "pygments", "mkdocs-biojulia"),
    repo = "github.com/BioJulia/Bio.jl.git",
    julia = "0.5",
    osname = "linux",
)


text:
# layout-related functions

@compat typealias ContViews{T,N} Union{ContiguousView{T,N},UnsafeContiguousView{T,N}}
@compat typealias NonContViews{T,N,M} Union{StridedView{T,N,M}, UnsafeStridedView{T,N,M}}

@compat typealias ContiguousArray{T,N} Union{Array{T,N}, ContiguousView{T,N}}
typealias ContiguousVector{T} ContiguousArray{T,1}
typealias ContiguousMatrix{T} ContiguousArray{T,2}

if VERSION >= v"0.4.0-dev+1784"
    import Base: linearindexing, LinearFast, LinearSlow

    linearindexing(a::ContViews) = LinearFast()
    linearindexing(a::NonContViews) = LinearSlow()
    linearindexing{T}(a::NonContViews{T,1}) = LinearFast()
end

## strides method

strides{T}(a::ContViews{T,1}) = (1,)
strides{T}(a::ContViews{T,2}) = (1, a.shp[1])
strides{T}(a::ContViews{T,3}) = ((d1, d2, d3) = size(a); (1, d1, d1 * d2))

strides{T}(a::ContViews{T,4}) = ((d1, d2, d3, d4) = size(a); (1, d1, d1 * d2, d1 * d2 * d3))

function strides{T}(a::ContViews{T,5})
    (d1, d2, d3, d4, d5) = a.shp
    s3 = d1 * d2
    s4 = s3 * d3
    s5 = s4 * d4
    (1, d1, s3, s4, s5)
end

strides(a::NonContViews) = a.strides

## stride method

stride{T}(a::ContiguousView{T,1}, d::Integer) =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 : length(a))::Int

stride{T}(a::ContiguousView{T,2}, d::Integer) =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] : length(a))::Int

stride{T}(a::ContiguousView{T,3}, d::Integer) =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] :
     d == 3 ? a.shp[1] * a.shp[2] : length(a))::Int

 stride{T}(a::ContiguousView{T,4}, d::Integer) =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] :
     d == 3 ? a.shp[1] * a.shp[2] :
     d == 4 ? a.shp[1] * a.shp[2] * a.shp[3] : length(a))::Int

stride{T}(a::ContiguousView{T,5}, d::Integer) =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] :
     d == 3 ? a.shp[1] * a.shp[2] :
     d == 4 ? a.shp[1] * a.shp[2] * a.shp[3] :
     d == 5 ? a.shp[1] * a.shp[2] * a.shp[3] * a.shp[4] : length(a))::Int

stride{T,N}(a::NonContViews{T,N}, d::Integer) =
    (d > 0 || error("dimension out of range.");
     d <= N ? a.strides[d] : a.strides[N])

 # test contiguousness

iscontiguous(a::ContViews) = true

iscontiguous{T,N}(a::NonContViews{T,N}) = _iscontiguous(size(a), strides(a))

_iscontiguous(shp::NTuple{0,Int}, strides::NTuple{0,Int}) = true
_iscontiguous(shp::NTuple{1,Int}, strides::NTuple{1,Int}) = (strides[1] == 1)
_iscontiguous(shp::NTuple{2,Int}, strides::NTuple{2,Int}) = (strides[1] == 1 && strides[2] == shp[1])
_iscontiguous(shp::NTuple{3,Int}, strides::NTuple{3,Int}) =
    (strides[1] == 1 &&
     strides[2] == shp[1] &&
     strides[3] == shp[1] * shp[2])

function _iscontiguous(shp::NTuple{4,Int}, strides::NTuple{4,Int})
    s2 = shp[1]
    s3 = s2 * shp[2]
    s4 = s3 * shp[3]
    strides[1] == 1 && strides[2] == s2 && strides[3] == s3 && strides[4] == s4
end

function _iscontiguous(shp::NTuple{5,Int}, strides::NTuple{5,Int})
    s2 = shp[1]
    s3 = s2 * shp[2]
    s4 = s3 * shp[3]
    s5 = s4 * shp[4]
    strides[1] == 1 && strides[2] == s2 && strides[3] == s3 && strides[4] == s4 && strides[5] == s5
end


text:
"""
    shiftlocus(locus::Locus, x)

Shift the index of `x` from the current locus to the new locus.

We only shift `Samped`, `Regular` or `Explicit`, `Intervals`.
"""
function shiftlocus(locus::Locus, lookup::LookupArray)
    samp = sampling(lookup)
    samp isa Intervals || error("Cannot shift locus of $(nameof(typeof(samp)))")
    newindex = _shiftindexlocus(locus, lookup)
    newlookup = rebuild(lookup; data=newindex)
    return set(newlookup, locus)
end

# Fallback - no shifting
_shiftindexlocus(locus::Locus, lookup::LookupArray) = index(lookup)
# Sampled
function _shiftindexlocus(locus::Locus, lookup::AbstractSampled)
    _shiftindexlocus(locus, span(lookup), sampling(lookup), lookup)
end
# TODO:
_shiftindexlocus(locus::Locus, span::Irregular, sampling::Sampling, lookup::LookupArray) = index(lookup)
# Sampled Regular
function _shiftindexlocus(destlocus::Center, span::Regular, sampling::Intervals, dim::LookupArray)
    if destlocus === locus(sampling)
        return index(dim)
    else
        offset = _offset(locus(sampling), destlocus)
        shift = ((index(dim) .+ abs(step(span))) .- index(dim)) .* offset
        return index(dim) .+ shift
    end
end
function _shiftindexlocus(destlocus::Locus, span::Regular, sampling::Intervals, lookup::LookupArray)
    index(lookup) .+ (abs(step(span)) * _offset(locus(sampling), destlocus))
end
# Sampled Explicit
_shiftindexlocus(::Start, span::Explicit, sampling::Intervals, lookup::LookupArray) = val(span)[1, :]
_shiftindexlocus(::End, span::Explicit, sampling::Intervals, lookup::LookupArray) = val(span)[2, :]
function _shiftindexlocus(destlocus::Center, span::Explicit, sampling::Intervals, lookup::LookupArray)
    _shiftindexlocus(destlocus, locus(lookup), span, sampling, lookup)
end
_shiftindexlocus(::Center, ::Center, span::Explicit, sampling::Intervals, lookup::LookupArray) = index(lookup)
function _shiftindexlocus(::Center, ::Locus, span::Explicit, sampling::Intervals, lookup::LookupArray)
    # A little complicated so that DateTime works
    (view(val(span), 2, :)  .- view(val(span), 1, :)) ./ 2 .+ view(val(span), 1, :)
end

_offset(::Start, ::Center) = 0.5
_offset(::Start, ::End) = 1
_offset(::Center, ::Start) = -0.5
_offset(::Center, ::End) = 0.5
_offset(::End, ::Start) = -1
_offset(::End, ::Center) = -0.5
_offset(::T, ::T) where T<:Locus = 0

maybeshiftlocus(locus::Locus, l::LookupArray) = _maybeshiftlocus(locus, sampling(l), l)

_maybeshiftlocus(locus::Locus, sampling::Intervals, l::LookupArray) = shiftlocus(locus, l)
_maybeshiftlocus(locus::Locus, sampling::Sampling, l::LookupArray) = l

"""
    basetypeof(x) => Type

Get the "base" type of an object - the minimum required to
define the object without it's fields. By default this is the full
`UnionAll` for the type. But custom `basetypeof` methods can be
defined for types with free type parameters.

In DimensionalData this is primariliy used for comparing `Dimension`s,
where `Dim{:x}` is different from `Dim{:y}`.
"""
@inline basetypeof(x::T) where T = basetypeof(T)
@generated function basetypeof(::Type{T}) where T
    if T isa Union
        T
    else
        getfield(parentmodule(T), nameof(T))
    end
end

# unwrap
# Unwrap Val and Type{Val}, or return unchanged
unwrap(::Val{X}) where X = X
unwrap(::Type{Val{X}}) where X = X
unwrap(x) = x

# orderof
# Detect the order of an abstract array
orderof(A::AbstractUnitRange) = ForwardOrdered()
orderof(A::AbstractRange) = _order(A)
function orderof(A::AbstractArray)
    local sorted, indord
    # This is awful. But we don't know if we can
    # call `issorted` on the contents of `A`.
    # This may be resolved by: https://github.com/JuliaLang/julia/pull/37239
    try
        indord = _order(A)
        sorted = issorted(A; rev=LookupArrays.isrev(indord))
    catch
        sorted = false
    end
    return sorted ? indord : Unordered()
end

_order(A) = first(A) <= last(A) ? ForwardOrdered() : ReverseOrdered()


text:
using BlackBoxOptimizationBenchmarking
using Documenter

DocMeta.setdocmeta!(BlackBoxOptimizationBenchmarking, :DocTestSetup, :(using BlackBoxOptimizationBenchmarking); recursive=true)

makedocs(;
    modules=[BlackBoxOptimizationBenchmarking],
    authors="Jonathan Bieler <jonathan.bieler@alumni.epfl.ch> and contributors",
    repo="https://github.com/jonathanBieler/BlackBoxOptimizationBenchmarking.jl/blob/{commit}{path}#{line}",
    sitename="BlackBoxOptimizationBenchmarking.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jonathanBieler.github.io/BlackBoxOptimizationBenchmarking.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jonathanBieler/BlackBoxOptimizationBenchmarking.jl",
    devbranch = "main",
)


text:
include("L_DVBN.jl")

f = readcsv("data_auto_mpg.csv")

data = Array(Any,392,8)

discrete_index = [2,7,8]
continuous_index = [1,3,4,5,6]

for i = 1 : 392
    for j = 1 : 8
        if j in discrete_index
            data[i,j] = round(Int64,f[i,j])
        else
            data[i,j] = f[i,j]
        end
    end
end

##### Network Structure is known in advance #####
graph = [2,(2,3),(3,5),(5,1),(1,5,3,7),(3,4),(4,6),8]
order_continuous_order = graph_to_reverse_conti_order(graph,continuous_index)
cut_time = 8

e1 = Discretize_All(data,graph,continuous_index,cut_time)
##### Network structure is not known in advance #####

u = 2
cut_time = 5
times = 3

order = [1,2,3,4,5,6,7,8]
e2 = Learn_DVBN(data,continuous_index,order,u,cut_time)

e3 = Learn_Discrete_Bayesian_Net(data,continuous_index,u,cut_time,times)


text:
"""
    matcher(t::AbstractTerm) -> AbstractMatcher

Build the matcher for `t`.
"""
matcher(t) = matcher(t, Set{Variable}())[1]

const Substitution = Dict{Variable,AbstractTerm}

struct Matches{S<:AbstractSubproblem}
    p::Substitution
    s::S
end
Base.IteratorSize(::Type{<:Matches}) = Base.SizeUnknown()

struct Fail end
Base.iterate(::Fail) = nothing
Base.length(::Fail) = 0
const fail = Fail()

"""
    match(pattern::AbstractMatcher, term::AbstractTerm)

Match `term` against `pattern`, producing an iterator containing all matches.
"""
function match(p::AbstractMatcher, t::AbstractTerm)
    σ = Substitution()
    s = match!(σ, p, t)
    s === nothing && return fail
    Matches(σ, s)
end
match(p::AbstractTerm, t::AbstractTerm) = match(matcher(p), t)

struct EmptySubproblem <: AbstractSubproblem end
Base.iterate(m::Matches{EmptySubproblem}) = (m.p, nothing)
Base.iterate(::Matches{EmptySubproblem}, ::Any) = nothing


text:
export rrt_trajectory

struct Node
    depth::Int
    score::Float64
    parent::Int     # index of parent
    pose::Vector{Float64}
end

struct Graph
    open_nodes::Vector{Node}
    end_nodes::Vector{Node}
end
function Graph(x0::Vector{Float64})
    n_init = Node(1, 0.0, 0, x0)
    return Graph([n_init], Node[])
end

# K is the number of vertices
function rrt_trajectory(em::ErgodicManager, tm::TrajectoryManager; verbose::Bool=true, logging::Bool=false, K::Int=100, max_depth::Int=100, step_length::Float64=.1)
    # first create the graph
    G = Graph(tm.x0)

    best_depth = 1  # log the deepest we get (mostly for debugging)

    for k = 1:K
        # sample random point
        rx = rand() * em.domain.lengths[1] + em.domain.mins[1]
        ry = rand() * em.domain.lengths[2] + em.domain.mins[2]

        # find nearest node
        nearest_node = G.open_nodes[1]
        nearest_idx = 1
        idx = 1
        d2min = Inf
        for node in G.open_nodes
            dx = node.pose[1] - rx
            dy = node.pose[2] - ry
            d2 = dx*dx + dy*dy
            if d2 < d2min
                d2min = d2
                nearest_node = node
                nearest_idx = idx
            end
            idx += 1
        end

        # step from nearest node to new point
        dx = (rx - nearest_node.pose[1]) * step_length / sqrt(d2min)
        dy = (ry - nearest_node.pose[2]) * step_length / sqrt(d2min)
        new_pose = [nearest_node.pose[1] + dx, nearest_node.pose[2] + dy]

        # compute ergodic score
        traj = [new_pose, nearest_node.pose]
        current_node = nearest_node
        while current_node.parent > 0
            current_node = G.open_nodes[current_node.parent]
            push!(traj, current_node.pose)
        end
        es = ergodic_score(em, traj)

        # create new node and push into graph
        depth = nearest_node.depth + 1
        best_depth = max(depth, best_depth)
        new_node = Node(depth, es, nearest_idx, new_pose)
        if depth < max_depth
            push!(G.open_nodes, new_node)
        else
            push!(G.end_nodes, new_node)
        end

        if verbose
            println("num nodes = ", length(G.open_nodes))
            println("num closed = ", length(G.end_nodes))
        end
    end

    println("best_depth = ", best_depth)

    # now loop through the end nodes, finding the best one
    best_score = Inf
    best_node = G.open_nodes[1]     # just any node
    for node in G.end_nodes
        if node.score < best_score
            best_score = node.score
            best_node = node
        end
    end

    # make the full trajectory
    traj = [best_node.pose]
    current_node = best_node
    while current_node.parent > 0
        current_node = G.open_nodes[current_node.parent]
        push!(traj, current_node.pose)
    end

    # actually, reverse this
    return reverse(traj)
end


text:
function isaksprime(n::Int64)
    if n < 2
        return false
    end
    for c in polycoefs(n)[2:(end-1)]
        if c%n != 0
            return false
        end
    end
    return true
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module CostUtils

import JuMP
export symbol, normalize

symbol(s::Symbol) = s
symbol(s::JuMP.Containers.DenseAxisArrayKey) = only(s.I)

function normalize(weights)
    total = sum(weights)
    map(w -> w / total, weights)
end

function namedtuple(associative_array)
    nt_keys = Tuple(symbol.(keys(associative_array)))
    NamedTuple{nt_keys}(tuple(values(associative_array)...))
end

end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# This file contains configuration-related functions and types

# Dicts to store user-specified configuration
_t = @static Sys.iswindows() ? "windows" : (@static Sys.isapple() ? "aqua" : "qt")
const default_config() = Dict(:mode      => "normal",
                              :timeout   => Sys.isunix() ? 10 : 20,
                              :debug     => false,
                              :term      => _t,
                              :termopts  => "",
                              :preamble  => "",
                              :saveopts  => "",
                              :showable  => "png")

# Set any of Gaston's configuration variables
function set(;reset::Bool = false,
              term        = config[:term],
              termopts    = config[:termopts],
              mode        = config[:mode],
              debug::Bool = config[:debug],
              timeout     = config[:timeout],
              preamble    = config[:preamble],
              saveopts    = config[:saveopts],
              showable    = config[:showable]
            )
    global config

    if reset
        config = default_config()
        return nothing
    end

    config[:mode] = mode
    config[:timeout] = timeout
    debug isa Bool && (config[:debug] = debug)
    config[:term] = term
    config[:termopts] = termopts
    config[:preamble] = preamble
    config[:saveopts] = saveopts
    config[:showable] = showable

    return nothing
end

### Encode terminal capabilities
# supports multiple windows
const term_window = ["qt", "wxt", "x11", "aqua", "windows"]
# outputs text
const term_text = ["dumb", "sixelgd"]

# Validate coordinates
# TODO: validate all possible cases
function valid_coords(x,y,z,supp)
    invalid = false
    if isempty(z)
        length(x) != length(y) && (invalid = true)
    end

    invalid && throw(DimensionMismatch("input vectors must have the same number of elements."))

    return true

end


text:
"""
     AbstractCausalAnalysis{CT <: Union{CausalityTest, MetaCausalityTest}, DT, RT}

Abstract return type for a causality test. `CT` is the type of the causality test,
`DT` is the data type, and `RT` is the result type.

As a minimum, concrete subtypes must implement:

- **`test`**: The parameters of the test.
- **`data`**: The data on which the test was applied.
- **`result`**: The result of applying the relevant causality statistic using
    the parameters in `test` to `data`.
"""
abstract type AbstractCausalAnalysis{CT <: Union{CausalityTest, MetaCausalityTest}, DT, RT} end


text:
## ExoplanetsSysSim/src/io.jl
## (c) 2015 Eric B. Ford

module SysSimIO
using ExoplanetsSysSim
using HDF5, JLD

if VERSION >= v"0.5-"
  import Compat: UTF8String, ASCIIString
end

export save_sim_param, save_sim_results, load_sim_param, load_distances, load_summary_stats

function save_sim_param(filename::ASCIIString, p::SimParam)
 local file
 try
   file = JLD.jldopen(filename,"w")
   write_sim_param(file,p)
 catch
   println("# Problem writing parameters to jld file: ", filename)
   return false
 finally
   close(file)
 end
 return true
end

function save_sim_results(filename::ASCIIString, p::SimParam; distances::Vector{Float64}=Array{Float64}(0), summary_stats::CatalogSummaryStatistics = CatalogSummaryStatistics() )
 local file
 try
   file = JLD.jldopen(filename,"w")
   write_sim_param(file,p)
   write_sim_summary_stats(file,summary_stats)
   write_sim_distances(file,distances)
 catch
   println("# Problem writing data to jld file: ", filename)
 finally
   close(file)
  end
end

function write_sim_distances(file::JLD.JldFile, d::Vector{Float64} )
  JLD.write(file,"distances",d)
end

function write_sim_summary_stats(file::JLD.JldFile, ss::CatalogSummaryStatistics )
  JLD.write(file,"summary_stats",ss.stat)
end

function write_sim_param(file::JLD.JldFile, p::SimParam)
 sim_param_bool = Dict{ASCIIString,Bool}()
 sim_param_int = Dict{ASCIIString,Integer}()
 sim_param_real = Dict{ASCIIString,Real}()
 sim_param_function = Dict{ASCIIString,ASCIIString}()
 sim_param_string = Dict{ASCIIString,ASCIIString}()
 for k in keys(p.param)
   #println("# k=",k,".  v=",p.param[k])
   if typeof(p.param[k]) <: Bool
      sim_param_bool[k] = p.param[k]
   elseif typeof(p.param[k]) <: Integer
      sim_param_int[k] = p.param[k]
   elseif typeof(p.param[k]) <: Real
      sim_param_real[k] = p.param[k]
   elseif typeof(p.param[k]) <: Function
      sim_param_function[k] = string(p.param[k])
   elseif typeof(p.param[k]) <: ASCIIString
      sim_param_string[k] = p.param[k]
   else
	  warn(string("Can't store value of >",k,"< due to type ", typeof(p.param[k])))
   end
 end
   JLD.write(file,"sim_param_int",sim_param_int)
   JLD.write(file,"sim_param_real",sim_param_real)
   JLD.write(file,"sim_param_function",sim_param_function)
   JLD.write(file,"sim_param_string",sim_param_string)
   JLD.write(file,"sim_param_bool",sim_param_bool)
   JLD.write(file,"sim_param_active",p.active)
end

function load_sim_param(filename::ASCIIString)
  local jld_data
  try
    jld_data = load(filename)
  catch
   println("# Problem reading parameters from jld file: ", filename)
  end
  p = SimParam()
  merge!(p.active, jld_data["sim_param_active"])
  merge!(p.param, jld_data["sim_param_int"])
  merge!(p.param, jld_data["sim_param_string"])
  merge!(p.param, jld_data["sim_param_real"])
  # merge!(p.param, jld_data["sim_param_function"])
  df = jld_data["sim_param_function"]
  for k in keys(df)
    p.param[k]::Function = symbol(df[k])
  end
  merge!(p.param, jld_data["sim_param_bool"])
  return p
end

function load_distances(filename::ASCIIString)
  local jld_data
  try
    jld_data = load(filename)
  catch
   println("# Problem reading distances from jld file: ", filename)
  end
  d::Array{Float64,1} = jld_data["distances"]
  return d
end

function load_summary_stats(filename::ASCIIString)
  local jld_data
  try
    jld_data = load(filename)
  catch
   println("# Problem reading parameters from jld file: ", filename)
  end
  s = CatalogSummaryStatistics()
  merge!(s.stat, jld_data["summary_stats"])
  return s
end

function test_io()
  sim_param = setup_sim_param_demo()
  save_sim_param("test.jld",sim_param)
  spd = load_sim_param("test.jld")
  rm("test.jld")
end

end # module SysSimIO


text:
using Flux, CUDA, Test
using Flux: pullback

@testset "CUDNN BatchNorm" begin
    @testset "4D Input" begin
        x = Float64.(collect(reshape(1:12, 2, 2, 3, 1)))
        m = BatchNorm(3)
        cx = gpu(x)
        cm = gpu(m)

        y, back = pullback((m, x) -> m(x), m, x)
        cy, cback = pullback((m, x) -> m(x), cm, cx)

        @test cpu(cy) ≈ y

        Δ = randn(size(y))
        dm, dx = back(Δ)
        cdm, cdx = cback(gpu(Δ))

        @test dm[].γ ≈ cpu(cdm[].γ)
        @test dm[].β ≈ cpu(cdm[].β)
        @test dx ≈ cpu(cdx)
    end

    @testset "2D Input" begin
        x = Float64.(collect(reshape(1:12, 3, 4)))
        m = BatchNorm(3)
        cx = gpu(x)
        cm = gpu(m)

        y, back = pullback((m, x) -> m(x), m, x)
        cy, cback = pullback((m, x) -> m(x), cm, cx)

        @test cpu(cy) ≈ y

        Δ = randn(size(y))
        dm, dx = back(Δ)
        cdm, cdx = cback(gpu(Δ))

        @test dm[].γ ≈ cpu(cdm[].γ)
        @test dm[].β ≈ cpu(cdm[].β)
        @test dx ≈ cpu(cdx)
    end
end


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

JLLWrappers.@generate_main_file_header("yaml_cpp")
JLLWrappers.@generate_main_file("yaml_cpp", UUID("01fea8cc-7d33-533a-824e-56a766f4ffe8"))
end  # module yaml_cpp_jll


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractSampler

A method for sampling from spatial objects.
"""
abstract type AbstractSampler end

"""
    sample(object, sampler)

Sample elements from `object` with `sampler`.
"""
sample(object::AbstractSpatialObject, sampler::AbstractSampler) = @error "not implemented"

#------------------
# IMPLEMENTATIONS
#------------------
include("sampling/uniform_sampler.jl")
include("sampling/ball_sampler.jl")


text:
using Stopping

import Stopping: reinit!, update!

const Iterate           = Union{Number, Vector, Nothing}
const FloatVoid         = Union{Number, Nothing}
const MatrixType        = Union{Number, AbstractArray, AbstractMatrix, Nothing}
"""
Type: MPCCAtX
Methods: update!, reinit!

MPCCAtX contains the information concerning a nonlinear problem at
the iteration x.
min_{x ∈ ℜⁿ} f(x)
subject to lcon <= c(x) <= ucon, lvar <= x <= uvar, lG <= G(x) _|_ H(x) >= lH.

Basic information is:
 - x : the current candidate for solution to our original problem
 - fx : which is the funciton evaluation at x
 - gx : which is the gradient evaluation at x
 - Hx : which is the hessian representation at x

 - mu : Lagrange multiplier of the bounds constraints

 - cx : evaluation of the constraint function at x
 - Jx : jacobian matrix of the constraint function at x
 - lambda : Lagrange multiplier of the constraints

 - current_time : time
 - evals : number of evaluations of the function

Note: * by default, unknown entries are set to nothing (except evals).
      * All these information (except for x and lambda) are optionnal and need to be update when
        required. The update is done trhough the update! function.
      * x and lambda are mandatory entries. If no constraints lambda = [].
      * The constructor check the size of the entries.
"""
mutable struct 	MPCCAtX{T <: AbstractVector} <: AbstractState

#Unconstrained State
    x            :: T     # current point
    fx           :: Union{eltype(T),Nothing}   # objective function
    gx           :: Union{T,eltype(T),Nothing}     # gradient size: x
    Hx           :: MatrixType  # hessian size: |x| x |x|

#Bounds State
    mu           :: Union{T,eltype(T),Nothing}     # Lagrange multipliers with bounds size of |x|

#Constrained State
    cx           :: Union{T,eltype(T),Nothing}     # vector of constraints lc <= c(x) <= uc
    Jx           :: MatrixType  # jacobian matrix, size: |lambda| x |x|
    lambda       :: T    # Lagrange multipliers

#Complementarity State
    cGx           :: Union{T,eltype(T),Nothing}     # vector of constraints lG <= G(x) <= uG
    JGx           :: MatrixType  # jacobian matrix, size: |ncc| x |x|
    lambdaG       :: Union{T,eltype(T),Nothing}    # Lagrange multipliers
    cHx           :: Union{T,eltype(T),Nothing}     # vector of constraints lH <= H(x) <= uH
    JHx           :: MatrixType  # jacobian matrix, size: |ncc| x |x|
    lambdaH       :: Union{T,eltype(T),Nothing}    # Lagrange multipliers

    d            :: Union{T,eltype(T),Nothing} #search direction
    res          :: Union{T,eltype(T),Nothing} #residual

 #Resources State
    current_time   :: Union{eltype(T),Nothing}
    current_score  :: Union{T,eltype(T),Nothing}
    evals          :: MPCCCounters

 function MPCCAtX(x             :: T,
                  lambda        :: T;
                  fx            :: FloatVoid    = nothing,
                  gx            :: Iterate      = nothing,
                  Hx            :: MatrixType   = nothing,
                  mu            :: Iterate      = nothing,
                  cx            :: Iterate      = nothing,
                  Jx            :: MatrixType   = nothing,
                  lambdaG       :: Iterate      = nothing,
                  lambdaH       :: Iterate      = nothing,
                  cGx           :: Iterate      = nothing,
                  JGx           :: MatrixType   = nothing,
                  cHx           :: Iterate      = nothing,
                  JHx           :: MatrixType   = nothing,
                  d             :: Iterate      = nothing,
                  res           :: Iterate      = nothing,
                  current_time  :: FloatVoid    = nothing,
                  current_score :: Iterate     = nothing,
                  evals         :: MPCCCounters = MPCCCounters()) where T <: AbstractVector

  return new{T}(x, fx, gx, Hx, mu,
                cx, Jx, lambda,
                cGx, JGx, lambdaG,
                cHx, JHx, lambdaH,
                d, res, current_time, current_score, evals)
 end
end

"""
reinit!: function that set all the entries at void except the mandatory x

Warning: if x, lambda or evals are given as a keyword argument they will be
prioritized over the existing x, lambda and the default Counters.
"""
function reinit!(stateatx :: MPCCAtX, x :: Vector, l :: Vector; kwargs...)

 for k ∈ fieldnames(MPCCAtX)
   if !(k ∈ [:x,:lambda,:evals]) setfield!(stateatx, k, nothing) end
 end

 return update!(stateatx; x=x, lambda = l, evals = MPCCCounters(), kwargs...)
end

"""
reinit!: short version of reinit! reusing the x in the state

Warning: if x, lambda or evals are given as a keyword argument they will be
prioritized over the existing x, lambda and the default Counters.
"""
function reinit!(stateatx :: MPCCAtX; kwargs...)
 return reinit!(stateatx, stateatx.x, stateatx.lambda; kwargs...)
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

#############
# test_roti #
#############

function test_roti!(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    angle = rand(elty)
    d_x,d_y = CUSPARSE.roti!(d_x,d_y,cos(angle),sin(angle),'O')
    h_x = to_host(d_x)
    h_y = to_host(d_y)
    z = copy(x)
    w = copy(y)
    y[x.nzind] = cos(angle)*w[z.nzind] - sin(angle)*z.nzval
    nx = SparseVector(m,x.nzind,cos(angle)*z.nzval + sin(angle)*w[z.nzind])
    @test_approx_eq(h_y,y)
    @test_approx_eq(h_x,nx)
end

function test_roti(elty)
    x = sparsevec(rand(1:m,k), rand(elty,k), m)
    y = rand(elty,m)
    d_x = CudaSparseVector(x)
    d_y = CudaArray(y)
    angle = rand(elty)
    d_z,d_w = CUSPARSE.roti(d_x,d_y,cos(angle),sin(angle),'O')
    h_w = to_host(d_w)
    h_z = to_host(d_z)
    z = copy(x)
    w = copy(y)
    w[z.nzind] = cos(angle)*y[x.nzind] - sin(angle)*x.nzval
    nz = SparseVector(m,z.nzind, cos(angle)*x.nzval + sin(angle)*y[x.nzind])
    @test_approx_eq(h_w,w)
    @test_approx_eq(h_z,nz)
end

types = [Float32,Float64]
for elty in types
    tic()
    test_roti!(elty)
    println("roti! took ", toq(), " for ", elty)
    tic()
    test_roti(elty)
    println("roti took ", toq(), " for ", elty)
end


text:
include("INCLUDEME.jl")
using Compat
using Compat.Test

using Yao
using Yao.Zoo
using Yao.Blocks

using QCGANS, GradOptim

function train(qcg::QCGAN{N}, g_learning_rate::Real, d_learning_rate::Real, niter::Int, ninner::Int) where N
    for i in 1:niter
        ggrad, dgrad = gradient(qcg)
        dispatch!(+, qcg.generator, -ggrad.*g_learning_rate)
        for j = 1:ninner
            ggrad, dgrad = gradient(qcg)
            dispatch!(-, qcg.discriminator, -dgrad.*d_learning_rate)
        end
        print("dist = $(tracedist(qcg))")
        println("loss = $(loss(qcg)())")
    end
end

function train(qcg::QCGAN{N}, optim, optim_inner, niter::Int, ninner::Int) where N
    for i in 1:niter
        ggrad, dgrad = gradient(qcg)
        params = parameters(qcg.generator)
        update!(params, ggrad, optim)
        dispatch!(qcg.generator, params)
        for j = 1:ninner
            ggrad, dgrad = gradient(qcg)
            params = parameters(qcg.discriminator)
            update!(params, -dgrad, optim_inner)
            dispatch!(qcg.discriminator, params)
        end
        print("dist = $(tracedist(qcg)), ")
        println("loss = $(loss(qcg)())")
    end
end

function train!(qcg::QCGAN{N}, optim, niter::Int) where N
    for i in 1:niter
        ggrad, dgrad = gradient(qcg)
        params = parameters(qcg.circuit)
        update!(params, vcat(ggrad, -dgrad), optim)
        dispatch!(qcg.circuit, params)
        print("dist = $(tracedist(qcg)), ")
        println("loss = $(loss(qcg)())")
    end
end

function train_alter(qcg::QCGAN{N}, g_learning_rate::Real, d_learning_rate::Real, niter::Int) where N
    for i in 1:niter
        ggrad, dgrad = gradient(qcg)
        lossval = loss(qcg)()
        if lossval == 0
            dispatch!(+, qcg.generator, -ggrad.*g_learning_rate)
        else
            ggrad, dgrad = gradient(qcg)
            dispatch!(-, qcg.discriminator, -dgrad.*d_learning_rate)
        end
        #print("dist = $(tracedist(qcg))")
        println("loss = $(lossval)")
    end
end
# TODO using an ideal discriminator!
disc(b) = b in [0, 3, 5, 10, 12, 15]
function train(qcg::QCGAN{N}, g_learning_rate::Real, niter::Int, ninner::Int) where N
    for i in 1:niter
        ggrad, dgrad = gradient(qcg)
        dispatch!(+, qcg.generator, -ggrad.*g_learning_rate)
        print("dist = $(tracedist(qcg))")
        println("loss = $(loss(qcg)())")
    end
end

using BenchmarkTools
pairs_gen = [1=>2, 2=>3, 3=>1]
pairs_dis = [1=>2, 2=>3, 3=>4, 4=>1]
qcg = qcgan(uniform_state(3), 5, 5, pairs_gen, pairs_dis)
qcg.target.state = copy(qcg.reg0) |> qcg.generator |> state
optim = iRProp(2.0, 0.5);
optim_inner = iRProp(2.0, 0.5);
train(qcg, optim, optim_inner, 200, 20)
train(qcg, 0.2, 0.5, 200, 1)
train_alter(qcg, 0.2, 0., 200)
optim = iRProp(2.0, 0.5)
train!(qcg, optim, 200)

reg = rand_state(3)
@code_warntype apply!(reg, qcg.generator)

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
function showblock!(grid, ::ShowMakie, block::Image{2}, obs)
    ax = cleanaxis(grid[1, 1])
    plotimage!(ax, obs)
end
function showblock!(grid, ::ShowMakie, block::Mask{2}, obs)
    ax = cleanaxis(grid[1, 1])
    plotmask!(ax, obs, block.classes)
end

function showblock!(grid, ::ShowMakie, block::Keypoints{2}, obs)
    ax = cleanaxis(grid[1, 1])
    h = maximum(first.(obs))
    ks = [SVector(x, h-y) for (y, x) in obs]
    Makie.scatter!(ax, ks)
end

function showblock!(grid, ::ShowMakie, block::Bounded{2, <:Keypoints{2}}, obs)
    ax = cleanaxis(grid[1, 1])
    h, w = block.size
    ks = [SVector(x, h-y) for (y, x) in obs]
    Makie.xlims!(ax, 0, w)
    Makie.ylims!(ax, 0, h)
    Makie.scatter!(ax, ks)
end

# ## Helpers

@recipe(PlotImage, image) do scene
    Makie.Attributes(
        alpha = 1,
        interpolate = false,
    )
end

function Makie.plot!(plot::PlotImage)
    im = plot[:image]
    rim = @lift alphacolor.(copy(rotr90($im)), $(plot.attributes[:alpha]))
    Makie.image!(plot, rim; plot.attributes...)
    return plot
end

@recipe(PlotMask, mask, classes) do scene
    Makie.Attributes()
end

function Makie.plot!(plot::PlotMask; kwargs...)
    mask = plot[:mask]
    classes = try
        classes = plot[:classes]
    catch
        classes = @lift unique($mask)
    end
    im = @lift _maskimage($mask, $classes)
    plotimage!(plot, im; alpha = 1, plot.attributes...)
    return plot
end

"""
    cleanaxis(f)

Create a `Makie.Axis` with no interactivity, decorations and aspect distortion.
"""
function cleanaxis(f; kwargs...)
    ax = Makie.Axis(f; kwargs...)
    ax.aspect = Makie.DataAspect()
    ax.xzoomlock = true
    ax.yzoomlock = true
    ax.xrectzoom = false
    ax.yrectzoom = false
    ax.panbutton = nothing
    ax.xpanlock = true
    ax.ypanlock = true
    ax.bottomspinevisible = false
    ax.leftspinevisible = false
    ax.rightspinevisible = false
    ax.topspinevisible = false
    Makie.hidedecorations!(ax)

    return ax
end


text:
"""
`FitnessScheme` defines how fitness vectors/values are
compared, presented and aggregated.
`Fitness` represents the score of one and the same individual
on one or a set of evaluations. A scheme is a specific
way to consider the scores in a coherent way.
Type parameter `F` specifies the type of fitness values.

`FitnessScheme` could also be used as a function that defines the fitness
ordering, i.e. `fs(x, y) == true` iff fitness `x` is better than `y`.
"""
abstract type FitnessScheme{F} end

"""
    fitness_type(fs::FitnessScheme)
    fitness_type(fs_type::Type{FitnessScheme})

Get the type of fitness values for fitness scheme `fs`.
"""
fitness_type(::Type{<:FitnessScheme{F}}) where F = F
fitness_type(fs::FitnessScheme) = fitness_type(typeof(fs))
fitness_eltype(::Type{<:FitnessScheme{F}}) where {F<:Number} = F
fitness_eltype(fs::FitnessScheme) = fitness_eltype(typeof(fs))

# trivial convert() between calculated and archived fitness
Base.convert(::Type{F}, fit::F, fit_scheme::FitnessScheme{F}) where F = fit

# ordering induced by the fitness scheme
# FIXME enable once v0.5 issue #14919 is fixed
# (fs::FitnessScheme{F}){F}(x::F, y::F) = is_better(x, y, fs)

"""
In `RatioFitnessScheme` the fitness values can be ranked on a ratio scale so
pairwise comparisons are not required.
The default scale used is the aggregate of the fitness components.

"""
abstract type RatioFitnessScheme{F} <: FitnessScheme{F} end
# FIXME is it necessary?

"""
`Float64`-valued scalar fitness scheme.
The boolean type parameter `MIN` specifies if smaller fitness values
are better (`true`) or worse (`false`).
"""
struct ScalarFitnessScheme{MIN} <: RatioFitnessScheme{Float64}
end

const MinimizingFitnessScheme = ScalarFitnessScheme{true}()
const MaximizingFitnessScheme = ScalarFitnessScheme{false}()

is_minimizing(::ScalarFitnessScheme{MIN}) where {MIN} = MIN
nafitness(::Type{F}) where {F<:Number} = convert(F, NaN)
@inline nafitness(fs::FitnessScheme) = nafitness(fitness_type(fs))
isnafitness(f::F, ::RatioFitnessScheme{F}) where {F<:Number} = isnan(f)
numobjectives(::RatioFitnessScheme) = 1

"""
Aggregation is just the identity function for scalar fitness.
"""
aggregate(fitness::F, ::RatioFitnessScheme{F}) where {F<:Number} = fitness

is_better(f1::Float64, f2::Float64, scheme::ScalarFitnessScheme{true}) = f1 < f2
is_better(f1::Float64, f2::Float64, scheme::ScalarFitnessScheme{false}) = f1 > f2

"""
Complex-valued fitness.
"""
struct ComplexFitnessScheme <: FitnessScheme{ComplexF64}
# FIXME what is isbetter() for ComplexFitnessScheme?
end

# FIXME do we need it? it might be confused with problem's fitness bounds
worst_fitness(fs::FitnessScheme) = is_minimizing(fs) ? Inf : (-Inf)
best_fitness(fs::FitnessScheme) = -worst_fitness(fs)

hat_compare(a1::Number, a2::Number) =
    (a1 < a2) ? -1 : ((a1 > a2) ? 1 : (isnan(a1) ? (isnan(a2) ? 0 : 1) : (isnan(a2) ? -1 : 0)))

"""
Check whether `f1` or `f2` fitness is better.

Returns
  * `-1` if `f1` is better than `f2`
  * `1` if `f2` is better than `f1`
  * `0` if `f1` and `f2` are equal.
"""
function hat_compare(f1, f2, s::RatioFitnessScheme)
    if is_minimizing(s)
        hat_compare(aggregate(f1, s), aggregate(f2, s))
    else
        hat_compare(aggregate(f2, s), aggregate(f1, s))
    end
end

is_better(f1, f2, scheme::FitnessScheme) = hat_compare(f1, f2, scheme) == -1
is_worse(f1, f2, scheme::FitnessScheme) = hat_compare(f1, f2, scheme) == 1
same_fitness(f1, f2, scheme::FitnessScheme) = hat_compare(f1, f2, scheme) == 0

struct HatCompare{FS<:FitnessScheme}
    fs::FS

    HatCompare(fs::FS) where {FS<:FitnessScheme} = new{FS}(fs)
end

(hc::HatCompare)(x::F, y::F) where {F} = hat_compare(x, y, hc.fs)


text:
push!(LOAD_PATH,"../src/")
using Documenter, ExtendableSparse

makedocs(sitename="ExtendableSparse.jl",
         modules = [ExtendableSparse],
         doctest = true,
         clean = true,
         authors = "J. Fuhrmann",
         repo="https://github.com/j-fu/ExtendableSparse.jl",
         pages=[
             "Home"=>"index.md",
             "example.md",
             "changes.md",
             "api.md"
         ])

deploydocs(repo = "github.com/j-fu/ExtendableSparse.jl.git")


text:
# (1) パッケージをインストールする（初回のみ）
## ] + addで

# (2) ライブラリの読み出し
## dplyr的に操作したいのでDataFramesMetaを使う
using DataFrames
using DataFramesMeta
using HypothesisTests
using RCall
@rlibrary readr

# (3) データの読み込み
## CSV.read()だとネット上のファイルは読み込めないのでRCallを使う。遅い。
email_data =
  rcopy(R"readr::read_csv('http://www.minethatdata.com/Kevin_Hillstrom_MineThatData_E-MailAnalytics_DataMiningChallenge_2008.03.20.csv')")

# (4) データの準備
## 女性向けメールが配信されたデータを削除したデータを作成
@linq male_df =
  email_data |> where(:segment .!= "Womens E-Mail") |> # 女性向けメールが配信されたデータを削除
  transform(treatment = ifelse.(:segment .== "Mens E-Mail", 1, 0)) #介入を表すtreatment変数を追加

# (5) 集計による比較
## byでgroupbyとtransformを同時にできる
@linq summary_by_segment =
  male_df |> by(
    :treatment,
    conversion_rate = mean(:conversion), # グループごとのconversionの平均
    spend_mean = mean(:spend), # グループごとのspendの平均
    count = length(:conversion),
  ) # グループごとのデータ数

# (6) t検定を行う
## (a)男性向けメールが配信されたグループの購買データを得る
@linq mens_mail = male_df |> where(:treatment .== 1)

## (b)メールが配信されなかったグループの購買データを得る
@linq no_mail = male_df |> where(:treatment .== 0)

## (a)(b)の平均の差に対して有意差検定を実行する
## HypothesisTestsの等分散t検定を使う
rct_ttest = EqualVarianceTTest(mens_mail[:, :spend], no_mail[:, :spend])

# (7) セレクションバイアスのあるデータの作成
## seedを固定する
## 書籍とseedを合わせるためRCallを使う
R"set.seed(1)"

## 条件に反応するサンプルの量を半分にする
obs_rate_c = 0.5
obs_rate_t = 0.5

## バイアスのあるデータの作成
@linq biased_data =
  male_df |> transform(
    obs_rate_c = ifelse.(
      (:history .> 300) .| (:recency .< 6) .| (:channel .== "Multichannel"),
      obs_rate_c,
      1,
    ),
    obs_rate_t = ifelse.(
      (:history .> 300) .| (:recency .< 6) .| (:channel .== "Multichannel"),
      1,
      obs_rate_t,
    ),
    random_number = rcopy(R"runif($(nrow(male_df)))"), # 乱数を合わせるためRに作らせる
  ) |> where(
    ((:treatment .== 0) .& (:random_number .< :obs_rate_c)) .|
    ((:treatment .== 1) .& (:random_number .< :obs_rate_t)),
  )

# (8) セレクションバイアスのあるデータで平均を比較
## group_byとsummairseを使って集計(Biased)
## byだけでできる
@linq summary_by_segment_biased =
  biased_data |> by(
    :treatment,
    conversion_rate = mean(:conversion),
    spend_mean = mean(:spend),
    count = length(:conversion),
  )

# (9) Rの関数であるt.testを使ってt検定を行う(Biased)
## (a)男性向けメールが配信されたグループの購買データを得る
@linq mens_mail_biased = biased_data |> where(:treatment .== 1)

## (b)メールが配信されなかったグループの購買データを得る
@linq no_mail_biased = biased_data |> where(:treatment .== 0)

## (a)(b)の平均の差に対して有意差検定を実行
rct_ttest_biased =
  EqualVarianceTTest(mens_mail_biased[:, :spend], no_mail_biased[:, :spend])


text:
module ExtrapTests

using Base.Test, DualNumbers
using Interpolations

f(x) = sin((x-3)*2pi/9 - 1)
xmax = 10
A = Float64[f(x) for x in 1:xmax]

itpg = interpolate(A, BSpline(Linear()), OnGrid())

etpg = extrapolate(itpg, Flat())

@test etpg[-3] == etpg[-4.5] == etpg[0.9] == etpg[1.0] == A[1]
@test etpg[10.1] == etpg[11] == etpg[148.298452] == A[end]

etpf = @inferred(extrapolate(itpg, NaN))

@test @inferred(size(etpf)) == (xmax,)
@test isnan(@inferred(getindex(etpf, -2.5)))
@test isnan(etpf[0.999])
@test_approx_eq @inferred(getindex(etpf, 1)) f(1)
@test_approx_eq etpf[10] f(10)
@test isnan(@inferred(getindex(etpf,10.001)))

@test etpf[2.5,1] == etpf[2.5]   # for show method
@test_throws BoundsError etpf[2.5,2]
@test_throws ErrorException etpf[2.5,2,1]  # this will probably become a BoundsError someday

@test isa(@inferred(getindex(etpf, dual(-2.5,1))), Dual)

etpl = extrapolate(itpg, Linear())
k_lo = A[2] - A[1]
x_lo = -3.2
@test_approx_eq etpl[x_lo] A[1]+k_lo*(x_lo-1)
k_hi = A[end] - A[end-1]
x_hi = xmax + 5.7
@test_approx_eq etpl[x_hi] A[end]+k_hi*(x_hi-xmax)

xmax, ymax = 8,8
g(x, y) = (x^2 + 3x - 8) * (-2y^2 + y + 1)

itp2g = interpolate(Float64[g(x,y) for x in 1:xmax, y in 1:ymax], (BSpline(Quadratic(Free())), BSpline(Linear())), OnGrid())
etp2g = extrapolate(itp2g, (Linear(), Flat()))

@test_approx_eq @inferred(getindex(etp2g, -.5, 4)) itp2g[1,4]-1.5*epsilon(etp2g[dual(1,1),4])
@test_approx_eq @inferred(getindex(etp2g, 5, 100)) itp2g[5,ymax]

etp2ud = extrapolate(itp2g, ((Linear(), Flat()), Flat()))
@test_approx_eq @inferred(getindex(etp2ud, -.5, 4)) itp2g[1,4] - 1.5*epsilon(etp2g[dual(1,1),4])
@test @inferred(getindex(etp2ud, 5, -4)) == etp2ud[5,1]
@test @inferred(getindex(etp2ud, 100, 4)) == etp2ud[8,4]
@test @inferred(getindex(etp2ud, -.5, 100)) == itp2g[1,8] - 1.5 * epsilon(etp2g[dual(1,1),8])

etp2ll = extrapolate(itp2g, Linear())
@test_approx_eq @inferred(getindex(etp2ll, -.5, 100)) itp2g[1,8]-1.5*epsilon(etp2ll[dual(1,1),8]) + (100-8) * epsilon(etp2ll[1,dual(8,1)])

# Allow element types that don't support conversion to Int (#87):
etp87g = extrapolate(interpolate([1.0im, 2.0im, 3.0im], BSpline(Linear()), OnGrid()), 0.0im)
@test @inferred(getindex(etp87g, 1)) == 1.0im
@test @inferred(getindex(etp87g, 1.5)) == 1.5im
@test @inferred(getindex(etp87g, 0.75)) == 0.0im
@test @inferred(getindex(etp87g, 3.25)) == 0.0im

etp87c = extrapolate(interpolate([1.0im, 2.0im, 3.0im], BSpline(Linear()), OnCell()), 0.0im)
@test @inferred(getindex(etp87c, 1)) == 1.0im
@test @inferred(getindex(etp87c, 1.5)) == 1.5im
@test @inferred(getindex(etp87c, 0.75)) == 0.75im
@test @inferred(getindex(etp87c, 3.25)) == 3.25im
@test @inferred(getindex(etp87g, 0)) == 0.0im
@test @inferred(getindex(etp87g, 3.7)) == 0.0im

# Make sure it works with Gridded too
etp100g = extrapolate(interpolate(([10;20],),[100;110], Gridded(Linear())), Flat())
@test @inferred(getindex(etp100g, 5)) == 100
@test @inferred(getindex(etp100g, 15)) == 105
@test @inferred(getindex(etp100g, 25)) == 110
end

include("type-stability.jl")


text:
abstract type Geometry end

struct Sphere <: Geometry
    radius
end

function Sphere(; radius)
    Sphere(radius)
end

struct Point <: Geometry end

struct Cylinder <: Geometry
    radius
    height
    open_ends::Bool
end

function Cylinder(; radius, height, open_ends = true)
    return Cylinder(radius, height, open_ends)
end

height(x::Geometry) = isdefined(x, :height) ? getfield(x, :height) : nothing
radius(x::Geometry) = getfield(x, :radius)
radius(::Point) = 0.0

# TODO: Remove unit stripping when we have proper unit checking implemented
area(x::Sphere) = ustrip(Float64, cm^2, 4*π*radius(x)^2)
area(x::Cylinder) = ustrip(Float64, cm^2, 2*π*radius(x)*(height(x) + (x.open_ends ? 0µm : radius(x))))
area(::Point) = 1.0

#=
@enum StimulusTrigger ContinuousTrigger EdgeTriggered

abstract type Stimulus end

struct CurrentStimulus
    waveform
    offset::Current
    trigger::StimulusTrigger
    delay::TimeF64
    function CurrentStimulus(waveform; offset::Current = 0nA,
                             trigger::StimulusTrigger = ContinuousTrigger,
                             delay::Time = 0ms)
        return new(waveform, offset, trigger, delay)
    end
end

struct VoltageStimulus end
=#

# IfElse.ifelse(t > 100.0, IfElse.ifelse(t <= 250.0, amplitude, 0.0), 0.0)

abstract type AbstractCompartmentSystem <: AbstractTimeDependentSystem end

struct CompartmentSystem <: AbstractCompartmentSystem
    ivs::Num # usually just t
    ## Intrinsic properties
    voltage::Num # symbol that represents membrane voltage
    capacitance::Num # specific membrane capacitance
    geometry::Geometry # compartment geometry metadata (shape, dimensions, etc)
    ## Dynamics
    chans::Set{AbstractConductanceSystem} # conductance systems
    channel_reversals::Set{Num}
    synapses::Set{AbstractConductanceSystem} # synaptic conductance systems
    synaptic_reversals::Set{Num}
    stimuli::Vector{Equation}
    extensions::Vector{ODESystem}
    defaults::Dict
    name::Symbol
end

function CompartmentSystem(
    Vₘ::Num,
    channels,
    reversals;
    capacitance = 1µF/cm^2,
    geometry::Geometry = Point(),
    extensions::Vector{ODESystem} = ODESystem[],
    stimuli::Vector{Equation} = Equation[],
    name::Symbol = Base.gensym("Compartment")
)
    @parameters cₘ = ustrip(Float64, mF/cm^2, capacitance)
    foreach(x -> isreversal(x) || throw("Invalid Equilibrium Potential"), reversals)
    foreach(x -> iscurrent(x.lhs) || throw("Invalid current stimulus"), stimuli)
    return CompartmentSystem(t, Vₘ, cₘ, geometry, Set(channels), Set(reversals), Set(),
                             Set(), stimuli, extensions, Dict(), name)
end

# AbstractSystem interface extensions
get_geometry(x::AbstractCompartmentSystem) = getfield(x, :geometry)
area(x::AbstractCompartmentSystem) = only(@parameters aₘ = area(get_geometry(x)))
capacitance(x::AbstractCompartmentSystem) = getfield(x, :capacitance)
get_output(x::AbstractCompartmentSystem) = getfield(x, :voltage)

get_extensions(x::AbstractCompartmentSystem) = getfield(x, :extensions)
get_channels(x::AbstractCompartmentSystem) = getfield(x, :chans)
get_synapses(x::AbstractCompartmentSystem) = getfield(x, :synapses)
get_stimuli(x::AbstractCompartmentSystem) = getfield(x, :stimuli)
# TODO: define top-level and recursive getters for inputs
function get_inputs(x::AbstractCompartmentSystem) end
function inputs(x::AbstractCompartmentSystem) end

function get_reversals(x::AbstractCompartmentSystem)
    return getfield(x, :channel_reversals), getfield(x, :synaptic_reversals)
end

function build_toplevel(comp_sys)
    dvs = Set{Num}()
    ps = Set{Num}()
    eqs = Equation[]
    defs = Dict()
    currents = Set{Num}()
    build_toplevel!(dvs, ps, eqs, defs, currents, comp_sys)
    return eqs, dvs, ps, defs, currents
end

function build_toplevel!(dvs, ps, eqs, defs, currents, comp_sys::CompartmentSystem)
    aₘ = area(comp_sys)
    cₘ = capacitance(comp_sys)
    Vₘ = get_output(comp_sys)
    push!(dvs, Vₘ)
    push!(ps, aₘ)
    push!(ps, cₘ)
    chan_revs, syn_revs = get_reversals(comp_sys)

    # Check for possible dynamic reversals
    rev_eq_vars = []
    for Erev in union(chan_revs, syn_revs)
        if isparameter(Erev)
            push!(ps, Erev)
        else
            push!(dvs, Erev)
            get_variables!(rev_eq_vars, getdefault(Erev))
            push!(eqs, Erev ~ getdefault(Erev))
        end
    end

    foreach(x -> isparameter(x) && push!(ps, x), rev_eq_vars)

    # Gather channel current equations
    for chan in get_channels(comp_sys)
        ion = permeability(chan)
        Erev = only(filter(x -> isequal(getion(x), ion), chan_revs))
        I = IonCurrent(ion, name = nameof(chan))
        g = renamespace(chan, get_output(chan))
        # TODO: checks for whether we need area scaling
        push!(eqs, I ~ g*aₘ*(Vₘ - Erev))
        push!(dvs, I)
        push!(currents, I)
        merge!(defs, defaults(chan))
    end

    # Gather synaptic current equations
    for synapse in get_synapses(comp_sys)
        ion = permeability(synapse)
        Esyn = only(filter(x -> isequal(getion(x), ion), syn_revs))
        I = IonCurrent(ion, name = nameof(chan))
        g = renamespace(synapse, get_output(synapse))
        push!(eqs, I ~ g*(Vₘ - Esyn))
        push!(dvs, I)
        push!(currents, I)
        merge!(defs, defaults(chan))
    end

    # Gather extension equations
    for extension in get_extensions(comp_sys)
        union!(eqs, equations(extension))
        union!(ps, parameters(extension))
        union!(dvs, states(extension))
        merge!(defs, defaults(extension))
    end

    for stimulus in get_stimuli(comp_sys)
        I = only(get_variables(stimulus.lhs))
        hasdefault(I) || push!(defs, I => stimulus.rhs)
        if isparameter(I)
            push!(ps, I)
        else
            push!(eqs, stimulus)
            push!(dvs, I)
            push!(currents, -I)
        end
    end
    # voltage equation
    push!(eqs, D(get_output(comp_sys)) ~ -sum(currents)/(cₘ*aₘ))
end

# collect _top level_ eqs including from extension + currents + reversals + Vₘ
function get_eqs(x::AbstractCompartmentSystem)
    build_toplevel(x)[1]
end

function get_states(x::AbstractCompartmentSystem)
    collect(build_toplevel(x)[2])
end

MTK.has_ps(x::CompartmentSystem) = true

function get_ps(x::AbstractCompartmentSystem)
    # collect _top level_ parameters from extension + currents + capacitance + area + reversals
    collect(build_toplevel(x)[3])
end

function defaults(x::AbstractCompartmentSystem)
    build_toplevel(x)[4]
end

function get_systems(x::AbstractCompartmentSystem)
    # collect channels + synapses + input systems
    collect(union(getfield(x, :chans), getfield(x, :synapses)))
end

function Base.convert(::Type{ODESystem}, compartment::CompartmentSystem)
    required_states = Set{Num}()
    eqs, dvs, ps, defs, currents = build_toplevel(compartment)
    syss = convert.(ODESystem, collect(get_systems(compartment)))

    # Resolve in/out: "connect" / auto forward cell states to channels
    for x in union(get_channels(compartment), get_synapses(compartment))
        for inp in get_inputs(x)
            # gathering inputs
            # FIXME: add metadata indicating the level we expect to resolve at
            # e.g. Vₘ used by a synapse is resolved at the network level (an external source)
            # therefore it should remain unresolved when building the _host compartment_
            push!(required_states, inp)
            # all inputs are outright connected to the compartment top-level
            subinp = getproperty(x, tosymbol(inp, escape=false))
            push!(eqs, inp ~ subinp)
        end
    end

    int_modified = Set{Num}()
    foreach(x -> get_variables!(int_modified, x.lhs),  eqs)

    # TODO: filter required states from stimuli and extensions
    union!(required_states, setdiff(dvs, int_modified))
    # Resolve unavailable states
    for s in required_states
        # resolvedby(s) !== compartment && continue
        # Handled based on metadata of each state (for now just one case)
        if iscurrent(s) && isaggregate(s)
            push!(eqs, s ~ sum(filter(x -> getion(x) == getion(s), currents)))
            push!(dvs, s)
        end
        # ... other switch cases for required state handlers ...
    end
    return ODESystem(eqs, t, dvs, ps; systems = syss, defaults = defs, name = nameof(compartment))
end


text:
import Base.Dates: Hour

@test Mocking.joinbinding(:Foo) == :(Foo)
@test Mocking.joinbinding(:Foo, :Bar) == :(Foo.Bar)
@test Mocking.joinbinding(:Foo, :Bar, :Baz) == :(Foo.Bar.Baz)

@test Mocking.splitbinding(:(Foo)) == [:Foo]
@test Mocking.splitbinding(:(Foo.Bar)) == [:Foo, :Bar]
@test Mocking.splitbinding(:(Foo.Bar.Baz)) == [:Foo, :Bar, :Baz]

int_expr = is(Int, Int32) ? :(Core.Int32) : :(Core.Int64)
@test Mocking.binding_expr(Int) == int_expr  # typealias. TODO: Change to Core.Int? Shouldn't actually matter
@test Mocking.binding_expr(Int64) == :(Core.Int64)  # concrete type
@test Mocking.binding_expr(Integer) == :(Core.Integer)  # abstract type
@test Mocking.binding_expr(Hour) == :(Base.Dates.Hour)  # unexported type
@test Mocking.binding_expr(rand) == :(Base.Random.rand)  # function
@test Mocking.binding_expr(AbstractArray{Int64}) == :(Core.AbstractArray)  # Core.AbstractArray{Int64}?
# @test Mocking.binding_expr(AbstractArray{T}) == :(Core.AbstractArray{T})

trans = Dict(:Int => Int, :Int64 => Int64, :Integer => Integer)
@test Mocking.adjust_bindings(trans) == Dict(
    :Int => int_expr,
    :Int64 => :(Core.Int64),
    :Integer => :(Core.Integer),
)

expr = :(f(a, b::Int64, c=3, d::Integer=4; e=5, f::Int=6))
@test Mocking.extract_bindings(expr.args[2:end]) == Set([:Int64, :Integer, :Int])
@test Mocking.call_parameters(expr) == [Expr(:parameters, Expr(:kw, :e, :e), Expr(:kw, :f, :f)), :a, :b, :c, :d]

expr = :(f(a..., b::Integer...))
@test Mocking.extract_bindings(expr.args[2:end]) == Set([:Integer])
@test Mocking.call_parameters(expr) == Any[Expr(:..., :a), Expr(:..., :b)]

expr = :(f(h::Hour=Hour(rand(1:24))))
@test Mocking.extract_bindings(expr.args[2:end]) == Set([:Hour, :rand])
@test Mocking.call_parameters(expr) == Any[:h]

expr = :(f(a::AbstractArray{Int64}))
@test Mocking.extract_bindings(expr.args[2:end]) == Set([:(AbstractArray{Int64})])
@test Mocking.call_parameters(expr) == Any[:a]

expr = :(f{T}(a::AbstractArray{T}))
@test Mocking.extract_bindings(expr.args[2:end]) == Set([:(AbstractArray{T})])
@test Mocking.call_parameters(expr) == Any[:a]


text:
function removeTailZeroes!(arr)
    while length(arr)>1 && iszero(arr[end])
        pop!(arr)
    end
    nothing
end

function null2indices(arr)
    res = Int[]
    isempty(arr) && (return res)

    @inbounds for (i, x) in enumerate(arr)
        append!(res, [i-1 for _ in 1:x])
    end
    return res
end

function inf2indices(arr)
    res = Int[]
    isempty(arr) && return res

    @inbounds for (i, x) in enumerate(arr)
        append!(res, [i for _ in 1:x])
    end
    return res
end

function uniquecomplex(arr::Array{T, 1}, tol=0) where {T<:Number}

    iszero(tol) && (tol = maximum(eps.(abs.(arr))))
    res = T[]
    while true
        for (i, c) in enumerate(arr)
            dist = abs.(arr .- c)
            idx = findall(x -> x > tol, dist)
            push!(res, c)
            arr = arr[idx]
            break
        end
        isempty(arr) && break
    end
    return res
end


text:
#### ContinuousIntegration

export ContinuousIntegration

"""
    ContinuousIntegration

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct ContinuousIntegration
    "Azure CI"
    azure::Union{Nothing, Azure}=nothing
    "Travis CI"
    travis::Union{Nothing, Travis}=nothing
    "AppVeyor CI"
    appveyor::Union{Nothing, AppVeyor}=nothing
    "GitHub CI"
    github_ci::Union{Nothing, GitHubCI}=nothing
    "GitHub Docker CI"
    github_docker_ci::Union{Nothing, GitHubDockerCI}=nothing
end

function ci_checks(ci::ContinuousIntegration)
  checks = []
  ci.azure === nothing || push!(checks, ci.azure.checks)
  ci.travis === nothing || push!(checks, ci.travis.checks)
  ci.appveyor === nothing || push!(checks, ci.appveyor.checks)
  ci.github_ci === nothing || push!(checks, ci.github_ci.checks)
  ci.github_docker_ci === nothing || push!(checks, ci.github_docker_ci.checks)
  return checks
end

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates


text:
using BeaData
using Base.Test

b = Bea()
nipatable = get_nipa_table(b, 5, "Q", 2014, 2014)


text:
@testset "Arithmetic" begin
    x = rand(Complex{Double64})
    y = rand(Double64)

    @test x + y isa Complex{Double64}
    @test y + x isa Complex{Double64}
    @test x - y isa Complex{Double64}
    @test y - x isa Complex{Double64}
    @test x * y isa Complex{Double64}
    @test y * x isa Complex{Double64}
    @test x / y isa Complex{Double64}
    @test y / x isa Complex{Double64}
    @test x ^ y isa Complex{Double64}
    @test y ^ x isa Complex{Double64}
    @test x ^ x isa Complex{Double64}
    @test y ^ y isa Double64

    x, y = rand(2)
    @test add2(x, y) == (Double64(x) + Double64(y))
    @test sub2(x, y) == (Double64(x) - Double64(y))
    @test mul2(x, y) == (Double64(x) * Double64(y))
    @test div2(x, y) == (Double64(x) / Double64(y))

    @test x ⊕ y == (Double64(x) + Double64(y))
    @test x ⊖ y == (Double64(x) - Double64(y))
    @test x ⊗ y == (Double64(x) * Double64(y))
    @test x ⊘ y == (Double64(x) / Double64(y))
end

@testset "div,rem,.." begin
    a = 17.0
    b =  5.0
    da = Double64(a)
    db = Double64(b)

    @test div(da, db) == div(a, b)
    @test cld(da, db) == cld(a, b)
    @test fld(da, db) == fld(a, b)
    @test rem(da, db) == rem(a, b)
    @test mod(da, db) == mod(a, b)
    @test divrem(da, db) == divrem(a, b)
    @test fldmod(da, db) == fldmod(a, b)
end

@testset "Trig functions" begin
    x = rand(Double64)
    sinx, cosx = sincos(x)
    @test sinx ≈ sin(x)
    @test cosx ≈ cos(x)
end

@testset "promotion" begin
    @test d64"123" - d64"123" * d64"0.1" * d64"10.0" == 0.0
    @test d64"123" - 123 * d64"0.1" * d64"10.0" == 0.0
end


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

mutable struct CreateOrUpdateVirtualNetworkRuleParameters <: SwaggerModel
    properties::Any # spec type: Union{ Nothing, CreateOrUpdateVirtualNetworkRuleProperties } # spec name: properties

    function CreateOrUpdateVirtualNetworkRuleParameters(;properties=nothing)
        o = new()
        validate_property(CreateOrUpdateVirtualNetworkRuleParameters, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type CreateOrUpdateVirtualNetworkRuleParameters

const _property_map_CreateOrUpdateVirtualNetworkRuleParameters = Dict{Symbol,Symbol}(Symbol("properties")=>Symbol("properties"))
const _property_types_CreateOrUpdateVirtualNetworkRuleParameters = Dict{Symbol,String}(Symbol("properties")=>"CreateOrUpdateVirtualNetworkRuleProperties")
Base.propertynames(::Type{ CreateOrUpdateVirtualNetworkRuleParameters }) = collect(keys(_property_map_CreateOrUpdateVirtualNetworkRuleParameters))
Swagger.property_type(::Type{ CreateOrUpdateVirtualNetworkRuleParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_CreateOrUpdateVirtualNetworkRuleParameters[name]))}
Swagger.field_name(::Type{ CreateOrUpdateVirtualNetworkRuleParameters }, property_name::Symbol) =  _property_map_CreateOrUpdateVirtualNetworkRuleParameters[property_name]

function check_required(o::CreateOrUpdateVirtualNetworkRuleParameters)
    (getproperty(o, Symbol("properties")) === nothing) && (return false)
    true
end

function validate_property(::Type{ CreateOrUpdateVirtualNetworkRuleParameters }, name::Symbol, val)
end


text:
# This file was generated, do not modify it. # hide
function MLJFlux.build(model::MyNetworkBuilder, rng, n_in, n_out)
    init = Flux.glorot_uniform(rng)
    layer1 = Flux.Dense(n_in, model.n1, init=init)
    layer2 = Flux.Dense(model.n1, model.n2, init=init)
    layer3 = Flux.Dense(model.n2, n_out, init=init)
    return Flux.Chain(layer1, layer2, layer3)
end

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials