text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using DocOpt
using Images

function blur(img, iterations)
    w, h = size(img)
    for i = 1:iterations
      img[3:w-2,3:h-2] =
           img[3-2:w-4,3-2:h-4] * 0.003  + img[3-1:w-3,3-2:h-4] * 0.0133 + img[3:w-2,3-2:h-4] * 0.0219 + img[3+1:w-1,3-2:h-4] * 0.0133 + img[3+2:w,3-2:h-4] * 0.0030 +
           img[3-2:w-4,3-1:h-3] * 0.0133 + img[3-1:w-3,3-1:h-3] * 0.0596 + img[3:w-2,3-1:h-3] * 0.0983 + img[3+1:w-1,3-1:h-3] * 0.0596 + img[3+2:w,3-1:h-3] * 0.0133 +
           img[3-2:w-4,3+0:h-2] * 0.0219 + img[3-1:w-3,3+0:h-2] * 0.0983 + img[3:w-2,3+0:h-2] * 0.1621 + img[3+1:w-1,3+0:h-2] * 0.0983 + img[3+2:w,3+0:h-2] * 0.0219 +
           img[3-2:w-4,3+1:h-1] * 0.0133 + img[3-1:w-3,3+1:h-1] * 0.0596 + img[3:w-2,3+1:h-1] * 0.0983 + img[3+1:w-1,3+1:h-1] * 0.0596 + img[3+2:w,3+1:h-1] * 0.0133 +
           img[3-2:w-4,3+2:h-0] * 0.003  + img[3-1:w-3,3+2:h-0] * 0.0133 + img[3:w-2,3+2:h-0] * 0.0219 + img[3+1:w-1,3+2:h-0] * 0.0133 + img[3+2:w,3+2:h-0] * 0.0030
    end
    return img
end

function main()
    doc = """gaussian-blur.jl

Gaussian blur image processing.

Usage:
  gaussian-blur.jl -h | --help
  gaussian-blur.jl [--img-file=<img-file>] [--iterations=<iterations>]

Options:
  -h --help                  Show this screen.
  --img-file=<img-file>      Specify a path to an input image file [default: ../../example.jpg].
  --iterations=<iterations>  Specify a number of iterations [default: 100].
"""
    arguments = docopt(doc)

    dir = dirname(@__FILE__)
    img_file = joinpath(dir, arguments["--img-file"])
    iterations = parse(Int, arguments["--iterations"])

    (fname, ext) = splitext(img_file)
    out_file = string(fname, "-blur", ".jpg")

    println("input file = ", img_file)
    println("iterations = ", iterations)
    println("output file = ", out_file)

    function gaussian_blur(input_fname, output_fname, iterations)
        local img = convert(Array{Float32, 2}, load(input_fname))
        tic()
        blur(img, 0)
        println("SELFPRIMED ", toq())

        tic()
        img = blur(img, iterations)
        println("SELFTIMED ", toq())
        save(output_fname, ufixed8sc(convert(Image, img)))
    end

    gaussian_blur(img_file, out_file, iterations)

end

main()


text: Dict{String, Any}("people" => Any[Dict{String, Any}("first_name" => Dict{String, Any}("value" => "Bruce", "type" => "string"), "last_name" => Dict{String, Any}("value" => "Springsteen", "type" => "string")), Dict{String, Any}("first_name" => Dict{String, Any}("value" => "Eric", "type" => "string"), "last_name" => Dict{String, Any}("value" => "Clapton", "type" => "string")), Dict{String, Any}("first_name" => Dict{String, Any}("value" => "Bob", "type" => "string"), "last_name" => Dict{String, Any}("value" => "Seger", "type" => "string"))])

ext: jl


text:
# using Distances

const start = time()

const DEBUG = length(ARGS) >=1 && ARGS[1] == "1"

function distance(solution::Array{Int32})
    cost = 0
    for i in 1:length(solution)-1
        cost += d[solution[i],solution[i+1]]
    end
    return cost
end

tic()
tic()
# file lines iterator object
# input = eachline(STDIN)

# state = start(input)
# read data size N
const n = parse(Int32, readline(STDIN))
const k = max(0.01*n, 25.0)

# 2xN node's coordinates matrix
data = Array{Tuple{Float64, Float64}}(n)

# read file data to Data matrix
for i in 1:n
    v = readline(STDIN)
    (i_v, x_v, y_v) = split(v)
    data[parse(Int32, i_v)] = (parse(Float64, x_v), parse(Float64, y_v))
end
# calculate nodes distances as NxN matrix

maxtime = parse(Float64, readline(STDIN))

# d = Array{Float64}(n, n)
# pairwise!(d, Euclidean(1e-12), data)

d = [hypot(x2-x1, y2-y1) for (x2, y2) in data, (x1, y1) in data]

currentsol = Array{Int32}(n+1)
currentsol[1] = currentsol[n+1] = Int32(1)

for i in 1:n
    d[i,i] = Inf
end

d1 = copy(d)
j = 1
for i in 2:length(currentsol)-1
    currentsol[i] = findfirst(d1[j,:].==min(d1[j,:]...))
    j = currentsol[i]
    d1[j,1] = Inf
    for visited_node in 1:i
        d1[j, currentsol[visited_node]] = Inf
    end
end

currentcost = distance(currentsol)

optimumsol = copy(currentsol)
optimumcost = currentcost

####### SETTINGS #######
# how many times to iterate (minimum)
const minItrs     = 10000000

# fraction of melting point for starting temperature
const meltPointF  = 0.7

# fraction of melting point for ending temperature
const targetTempF = 0.01

# fraction of stagnant minItrs allowed before reheating
const stagItrsF   = 0.1
########################

function iterate(t::Float64, xsol::Array{Int32}, xcost::Float64)
	ysol = copy(xsol)
    ycost = xcost

	t1 = rand(2:n-2)
	t2 = rand(t1+2:n)

    ycost = ycost -
    		d[ysol[t1-1], ysol[t1]] -
    		d[ysol[t1], ysol[t1+1]] +
    		d[ysol[t1-1], ysol[t2]] +
    		d[ysol[t2], ysol[t1+1]] -
    		d[ysol[t2-1], ysol[t2]] -
    		d[ysol[t2], ysol[t2+1]] +
    		d[ysol[t2-1], ysol[t1]] +
    		d[ysol[t1], ysol[t2+1]]

	ysol[t1], ysol[t2] = ysol[t2], ysol[t1]

	if ycost <= xcost
		return ysol, ycost
	elseif rand() < exp((xcost-ycost)/t)
		return ysol, ycost
	end
	return xsol, xcost
end

function iterateP(p::Float64, xsol::Array{Int32}, xcost::Float64)
	ysol = copy(xsol)
    ycost = xcost

	t1 = rand(2:n-2)
    t2Max = round(min(t1+k, n))
	t2 = rand(Int64(t1+2):Int64(t2Max))

    ycost = ycost -
    		d[ysol[t1-1], ysol[t1]] -
    		d[ysol[t1], ysol[t1+1]] +
    		d[ysol[t1-1], ysol[t2]] +
    		d[ysol[t2], ysol[t1+1]] -
    		d[ysol[t2-1], ysol[t2]] -
    		d[ysol[t2], ysol[t2+1]] +
    		d[ysol[t2-1], ysol[t1]] +
    		d[ysol[t1], ysol[t2+1]]

	ysol[t1], ysol[t2] = ysol[t2], ysol[t1]

	if ycost <= xcost
		return ysol, ycost
	elseif rand() < p
		return ysol, ycost
	end
	return xsol, xcost
end

const initialcost = currentcost
if DEBUG
    println(string("initial cost: ", initialcost))
    println(string("min iterations: ", minItrs))
    println(string("initial time: ", toq(), " seconds"))
end

testsol = copy(currentsol)
testcost = currentcost
minT = 10.0e10
maxT = 0

for i in 1:max(0.01*minItrs, 2.0)
    testsol, testcost = iterateP(0.001, testsol, testcost)

    minT = min(minT, testcost)
    maxT = max(maxT, testcost)
end

const meltPoint = (maxT - minT) * 10.0^(-1.0*log10(n))

const t0 = meltPoint * meltPointF

if DEBUG
    @printf("T0 = %.2f\n", t0)
end
t = t0

const tDecay = targetTempF ^ (1.0/minItrs)

itr = 1
optItr = 1
stagItrs = 0

function printresult()
    if DEBUG
        @printf("T = %.2f\n", t)
    	@printf("Optimum: %.5f\n", optimumcost)
    	@printf("%.2f%% improvement\n", 100.0-100*optimumcost/initialcost)
    	@printf("after %.2f%% of iterations\n", 100*optItr/itr)
    	@printf("done %.2f%% of minimum iterations\n", 100*itr/minItrs)
    	toc()
    else
        println(STDOUT, optimumcost)
        for city in optimumsol
            print(STDERR, city)
            print(STDERR, " ")
        end
        println(STDERR)
    end
end

atexit(printresult)

while t > t0*targetTempF

    if time() - start > maxtime - 1.0
        exit()
    end

	currentsol, currentcost = iterate(t, currentsol, currentcost)
	if currentcost < optimumcost
		optimumsol = copy(currentsol)
        optimumcost = currentcost
		optItr = itr
	else
		stagItrs += 1
	end

	if stagItrs == Int(stagItrsF*minItrs) && itr <= minItrs
		stagItrs = 0
        if DEBUG
            println("reheating")
        end
		t *= 1.0 + 0.7*(minItrs-itr)/minItrs
	end

	t *= tDecay

	if DEBUG && itr%(minItrs/5) == 0
		@printf("T = %.2f\n", t)
		@printf("cost: %.5f\n", currentcost)
		@printf("optimum: %.5f\n", optimumcost)
	end
    itr += 1
end


text:
#
# fmin-funcs.jl -
#
# Define functions for testing univariate optimization.
#
module FminTestFunctions

# Counter of function evaluations.
const cnt = Ref{Int}(0)
incr_nevals(n::Integer=1) = (cnt[] += n; nothing)
reset_nevals() =  (cnt[] = 0; nothing)
get_nevals() = cnt[]

"""
    obj = TestFunc{T}(f; id="...", a=..., b=.., x0=..., f0=NaN, ga=a, gb=b)

yields an object that wraps objective function `f` to be tested for univariate
minimization.  Keywords must be used to specify the identifier `id`, the bounds
`a` and `b` of the search interval, the (approximate) correct solution `x0`,
and corresponding minimum function value `f0`.  Keywords `ga` and `gb` are to
specify other endpoints for the search interval in global optimization (or
plotting).  All keywords but `f0`, `ga` and `gb` are mandatory.  Parameter `T`
is the floating-point type for computations, `Float64` by default.

The wrapper is intended to bundle the settings for testing the minimization of
the function and to improve type stability when the wrapper is called as
function to evaluate the objective function: `obj(x)` yields `f(convert(T,x))`.
This improved type-stability yields no regressions, on the contrary, it may reduce
allocations and speed-up tests in some cases (e.g., the AMPGO problem 18).

Also when the objective function is called via the wrapper, the global variable
storing the number of calls is automatically incremented.

"""
struct TestFunc{T<:AbstractFloat,F}
    f::F
    id::String
    a::T
    b::T
    ga::T
    gb::T
    x0::T
    f0::T
    function TestFunc{T}(f::F;
                         id::AbstractString,
                         a::Real, ga::Real=a,
                         b::Real, gb::Real=b,
                         x0::Real, f0::Real=NaN) where {T,F}
        new{T,F}(f, id, a, b, ga, gb, x0,
                 isnan(f0) ? f(T(x0)) : f0)
    end
end
TestFunc(args...; kwds...) = TestFunc{Float64}(args...; kwds...)
TestFunc{T}(f::TestFunc{T}) where {T<:AbstractFloat} = f
TestFunc{T}(f::TestFunc) where {T<:AbstractFloat} =
    TestFunc{T}(func(f); id=ident(f), a=local_lower(f), b=local_upper(f),
                ga=global_lower(f), gb=global_upper(f), x0=xsol(f), f0=fsol(f))

@inline (A::TestFunc{T})(x::Real) where {T<:AbstractFloat} = A(convert(T, x))
@inline (A::TestFunc{T})(x::T) where {T<:AbstractFloat} = (incr_nevals(); func(A)(x))
func(A::TestFunc) = A.f
ident(A::TestFunc) = A.id
local_lower(A::TestFunc) = A.a
local_upper(A::TestFunc) = A.b
global_lower(A::TestFunc) = A.ga
global_upper(A::TestFunc) = A.gb
xsol(A::TestFunc) = A.x0
fsol(A::TestFunc) = A.f0

# Problems 4 AMPGO (http://infinity77.net/global_optimization/test_functions_1d.html).
ampgo_2 = TestFunc(x -> sin(x) + sin(10x/3);
                   id="AMPGO_2",
                   a = 4.3,  b = 6.2,  # for local optimization
                   ga = 2.5, gb = 7.5, # for global optimization
                   x0 = 5.1457352882823, f0 = -1.8995993491521)
ampgo_3 = TestFunc(x -> (begin;
                         s = zero(x);
                         # FIXME: The formula (in AMPGO web page) is for
                         #        `k ∈ 1:6` but the figure and the given
                         #        minimum if for `k ∈ 1:5`.
                         for k in 1:5; s -= k*sin((k + 1)*x + k); end
                         return s;
                         end);
                   id = "AMPGO_3",
                   a = -7.3, b = -6.1, # for local optimization
                   ga = -10, gb = 10,  # for global optimization
                   x0 = -6.7745761455604, f0 = -12.0312494421671)
ampgo_4 = TestFunc(x -> -(16x^2 - 24x + 5)*exp(-x);
                   id = "AMPGO_4", a = 1.9, b = 3.9,
                   x0 = 2.8680339892289, f0 = -3.8504507088002)
ampgo_5 = TestFunc(x -> (3x - oftype(x, 1.4))*sin(18x);
                   id = "AMPGO_5",
                   a = 0.8, b = 1.1,   # for local optimization
                   ga = 0.0, gb = 1.2, # for global optimization
                   x0 = 0.9660858038332, f0 = -1.4890725386896)
ampgo_6 = TestFunc(x -> -(x + sin(x))*exp(-x^2);
                   id = "AMPGO_6",
                   a = -0.6, b = 10,  # for local optimization
                   ga = -10, gb = 10, # for global optimization
                   x0 = 0.6795786601089, f0 = -0.8242393984761)
ampgo_7 = TestFunc(x -> sin(x) + sin(10x/3) + log(x) - oftype(x,0.84)*x + 3;
                   id = "AMPGO_7",
                   a = 4.2, b = 6.1,   # for local optimization
                   ga = 2.7, gb = 7.5, # for global optimization
                   x0 = 5.1997783686004, f0 = -1.6013075464944)
ampgo_8 = TestFunc(x -> (begin;
                         s = zero(x);
                         # FIXME: The formula (in AMPGO web page) is for
                         #        `k ∈ 1:6` but the figure and the given
                         #        minimum if for `k ∈ 1:5`.
                         for k in 1:5; s -= k*cos((k + 1)*x + k); end
                         return s;
                         end);
                   id = "AMPGO_8",
                   a = -7.7, b = -6.15, # for local optimization
                   ga = -10, gb = 10,  # for global optimization
                   x0 = -7.0835064075500, f0 = -14.5080079271950)
ampgo_9 = TestFunc(x -> sin(x) + sin(2x/3);
                   id = "AMPGO_9",
                   a = 13.6, b = 20.4,  # for local optimization
                   ga = 3.1, gb = 20.4, # for global optimization
                   x0 = 17.0391989476448, f0 = -1.9059611187158)
ampgo_10 = TestFunc(x -> -x*sin(x);
                    id = "AMPGO_10",
                    a = 5, b = 10,   # for local optimization
                    ga = 0, gb = 10, # for global optimization
                    x0 = 7.9786657125325, f0 = -7.9167273715878)
ampgo_11 = TestFunc(x -> 2cos(x) + cos(2x);
                    id = "AMPGO_11",
                    a = 0, b = 3,       # for local optimization
                    ga = -π/2, gb = 2π, # for global optimization
                    x0 = 2.0943950957161, f0 = -1.5)
ampgo_12 = TestFunc(x -> sin(x)^3 + cos(x)^3;
                    id = "AMPGO_12",
                    a = 1.6, b = 3.9, # for local optimization
                    ga = 0, gb = 2π, # for global optimization
                    x0 = π, f0 = -1)
ampgo_13 = TestFunc(x -> -x^oftype(x,2/3) - (1 - x^2)^oftype(x,1/3);
                    id = "AMPGO_13",
                    a = 0.001, b = 0.99,
                    x0 = 1/sqrt(2), f0 = -1.5874010519682)
ampgo_14 = TestFunc(x -> -exp(-x)*sin(2π*x);
                    id = "AMPGO_14",
                    a = 0, b = 0.7,
                    ga = 0, gb = 4,
                    x0 = 0.2248803847323, f0 = -0.7886853874087)
ampgo_15 = TestFunc(x -> (x^2 - 5x + 6)/(x^2 + 1);
                    id = "AMPGO_15",
                    a = -0.3, b = 5,
                    ga = -5, gb = 5,
                    x0 = 2.4142135734136, f0 = -0.0355339059327)
ampgo_18 = TestFunc(x -> x ≤ 3 ? (x - 2)^2 : 2log(x - 2) + 1;
                    id = "AMPGO_18", a = 0, b = 6, x0 = 2, f0 = 0)
ampgo_20 = TestFunc(x -> (sin(x) - x)*exp(-x^2);
                    id = "AMPGO_20",
                    a = -1.1, b = 6.2,
                    ga = -10, gb = 10,
                    x0 = 1.1951366409628, f0 = -0.0634905289364)
ampgo_21 = TestFunc(x -> (sin(x) + cos(2x))*x;
                    id = "AMPGO_21",
                    a = 3.1, b = 6.5,
                    ga = 0, gb = 10,
                    x0 = 4.7954086865801, f0 = -9.5083504406331)
ampgo_22 = TestFunc(x -> exp(-3x) - sin(x)^3;
                    id = "AMPGO_22",
                    a = 11.1, b = 17.2,
                    ga = 0, gb = 20,
                    x0 = 14.1371669411515, f0 = -1.0)

# Test functions from GSL (GNU Scientific Library).
gsl_fmin_1 = TestFunc(x -> x^4 - 1;
                      id = "GSL_FMIN_1", a = -3, b = 17, x0 = 0, f0 = 0)
gsl_fmin_2 = TestFunc(x -> sqrt(abs(x));
                      id = "GSL_FMIN_2", a = -2.0, b = 1.5, x0 = 0, f0 = 0)
gsl_fmin_3 = TestFunc(x -> x < 1 ? float(one(x)) : -exp(-x);
                      id = "GSL_FMIN_3", a = -2.0, b = 4.0, x0 = 1.0)
gsl_fmin_4 = TestFunc(x -> x - 30/(1 + 100_000*(x - oftype(x,0.8))^2);
                      id = "GSL_FMIN_4", a = 0.72, b = 2.0,
                      x0 = 0.7999998333408, f0 = -29.2000000833333)

end # module


text:
function molecule(::Molecule"methane")
  return """
C	0.0000000	0.0000000	0.0000000
H	0.6268910	0.6268910	0.6268910
H	-0.6268910	-0.6268910	0.6268910
H	-0.6268910	0.6268910	-0.6268910
H	0.6268910	-0.6268910	-0.6268910
  """
end


text:
using Test
using OceanTurbulenceParameterEstimation
using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize
using Oceananigans.TurbulenceClosures: ConvectiveAdjustmentVerticalDiffusivity

#####
##### Parameters
#####

Nz = 16
Lz = 128
Qᵇ = 1e-8
Qᵘ = -1e-5
Δt = 20.0
f₀ = 1e-4
N² = 1e-5

stop_time = 6Δt
save_interval = 2Δt
N_ensemble = 3

experiment_name = "convective_adjustment_test"
data_path = experiment_name * ".jld2"

# "True" parameters to be estimated by calibration
convective_κz = 1.0
convective_νz = 0.9
background_κz = 1e-4
background_νz = 1e-5

@testset "SyntheticObservations tests" begin
    #####
    ##### Generate synthetic observations
    #####

    grid = RectilinearGrid(size=Nz, z=(-Lz, 0), topology=(Flat, Flat, Bounded))
    closure = ConvectiveAdjustmentVerticalDiffusivity(; convective_κz, background_κz, convective_νz, background_νz)
    coriolis = FPlane(f=f₀)

    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
    b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), bottom = GradientBoundaryCondition(N²))

    model = HydrostaticFreeSurfaceModel(grid = grid,
                                        tracers = :b,
                                        buoyancy = BuoyancyTracer(),
                                        boundary_conditions = (; u=u_bcs, b=b_bcs),
                                        coriolis = coriolis,
                                        closure = closure)

    set!(model, b = (x, y, z) -> N² * z)

    simulation = Simulation(model; Δt, stop_time)

    init_with_parameters(file, model) = file["parameters"] = (; Qᵇ, Qᵘ, Δt, N², tracers=keys(model.tracers))

    simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                          schedule = TimeInterval(save_interval),
                                                          prefix = experiment_name,
                                                          field_slicer = nothing,
                                                          force = true,
                                                          init = init_with_parameters)

    run!(simulation)

    #####
    ##### Load truth data as observations
    #####

    data_path = experiment_name * ".jld2"

    b_observations = SyntheticObservations(data_path, field_names=:b)
    ub_observations = SyntheticObservations(data_path, field_names=(:u, :b))
    uvb_observations = SyntheticObservations(data_path, field_names=(:u, :v, :b))

    @test keys(b_observations.field_time_serieses) == tuple(:b)
    @test keys(ub_observations.field_time_serieses) == tuple(:u, :b)
    @test keys(uvb_observations.field_time_serieses) == tuple(:u, :v, :b)
end


text:
abstract type LaneFollowingDriver <: DriverModel{LaneFollowingAccel} end
track_longitudinal!(model::LaneFollowingDriver, v_ego::Float64, v_oth::Float64, headway::Float64) = model # do nothing by default

function observe!(model::LaneFollowingDriver, scene::Scene{Entity{S, D, I}}, roadway::Roadway, egoid::I) where {S, D, I}

    ego = get_by_id(scene, egoid)

    fore = find_neighbor(scene, roadway, ego, targetpoint_ego=VehicleTargetPointFront(), targetpoint_neighbor=VehicleTargetPointRear())

    v_ego = vel(ego)
    v_oth = NaN
    headway = NaN

    if fore.ind != nothing
        v_oth = vel(scene[fore.ind].state)
        headway = fore.Δs
    end

    track_longitudinal!(model, v_ego, v_oth, headway)

    return model
end

mutable struct StaticLaneFollowingDriver <: LaneFollowingDriver
    a::LaneFollowingAccel
end
StaticLaneFollowingDriver() = StaticLaneFollowingDriver(LaneFollowingAccel(0.0))
StaticLaneFollowingDriver(a::Float64) = StaticLaneFollowingDriver(LaneFollowingAccel(a))
Base.rand(rng::AbstractRNG, model::StaticLaneFollowingDriver) = model.a
Distributions.pdf(model::StaticLaneFollowingDriver, a::LaneFollowingAccel) = isapprox(a.a, model.a.a) ? Inf : 0.0
Distributions.logpdf(model::StaticLaneFollowingDriver, a::LaneFollowingAccel) = isapprox(a.a, model.a.a) ? Inf : -Inf


text:
const APL = AbstractPolynomialLike

function _expectation(μ::Measure, p::APL, f)
    i = 1
    s = 0
    for t in terms(p)
        while i <= length(μ.x) && monomial(t) != μ.x[i]
            i += 1
        end
        if i > length(μ.x)
            error("The polynomial $p has a nonzero term $t with monomial $(t.x) for which the expectation is not known in $μ")
        end
        s += f(μ.a[i], coefficient(t))
        i += 1
    end
    s
end

"""
    MultivariateMoments.expectation(μ::AbstractMeasureLike, p::AbstractPolynomialLike)
    MultivariateMoments.expectation(p::AbstractPolynomialLike, μ::AbstractMeasureLike)

Computes the expectation ``\\mathbb{E}_{\\mu}[p]``.
"""
function expectation end

expectation(μ::Measure, p::APL) = _expectation(μ, p, (*))
expectation(p::APL, μ::Measure) = _expectation(μ, p, (a, b) -> b * a) # a and b may be noncommutative

"""
    dot(μ::AbstractMeasureLike, p::AbstractPolynomialLike)
    dot(p::AbstractPolynomialLike, μ::AbstractMeasureLike)

See [`expectation`](@ref)
"""
LinearAlgebra.dot(μ::AbstractMeasureLike, p::APL) = expectation(μ, p)
LinearAlgebra.dot(p::APL, μ::AbstractMeasureLike) = expectation(p, μ)


text:
using Test

@testset "test one" begin
    @test true
    @test true
    @test true
end


text:
using Hydra
using Documenter

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

makedocs(;
    modules=[Hydra],
    authors="Arthur Newbury",
    repo="https://github.com/EvoArt/Hydra.jl/blob/{commit}{path}#{line}",
    sitename="Hydra.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://EvoArt.github.io/Hydra.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/EvoArt/Hydra.jl",
)

max_stars_repo_name: EvoArt/Hydra

max_issues_repo_name: EvoArt/Hydra

max_forks_repo_name: EvoArt/Hydra


text:
module ConditionalDists

using LinearAlgebra
using Random
using Distributions
using DistributionsAD
using ChainRulesCore
using Functors
using Requires

export condition

export ConditionalDistribution
export ConditionalMvNormal
export ConditionalCategorical
export SplitLayer

include("cond_dist.jl")

include("batch_mvnormal.jl")
include("cond_mvnormal.jl")
include("categorical.jl")
include("utils.jl")

function __init__()
    @require Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" begin
        using .Flux: Dense
        function SplitLayer(in::Int, outs::Vector{Int}, acts::Vector)
            SplitLayer(Tuple(Dense(in,o,a) for (o,a) in zip(outs,acts)))
        end

        function SplitLayer(in::Int, outs::Vector{Int}, act=identity)
            acts = [act for _ in 1:length(outs)]
            SplitLayer(in, outs, acts)
        end
    end
end

end # module


text:
function perm_samp(data::PopData)
    popcounts = countmap(data.metadata.sampleinfo.population)
    pops, counts = keys(popcounts), values(popcounts)
    gdf = groupby(data.genodata, :name)
    perm_idx = shuffle(Xoroshiro128Star(), keys(gdf))

    #return gdf

end


text:
# Autogenerated wrapper script for LibGD_jll for x86_64-apple-darwin
export gd2copypal, gd2togif, gd2topng, gdcmpgif, gdparttopng, gdtopng, giftogd2, libgd, pngtogd, pngtogd2, webpng

using JpegTurbo_jll
using Zlib_jll
using libpng_jll
using Libtiff_jll
using Xorg_libXpm_jll
using Libiconv_jll
using libwebp_jll
JLLWrappers.@generate_wrapper_header("LibGD")
JLLWrappers.@declare_executable_product(gd2copypal)
JLLWrappers.@declare_executable_product(gd2togif)
JLLWrappers.@declare_executable_product(gd2topng)
JLLWrappers.@declare_executable_product(gdcmpgif)
JLLWrappers.@declare_executable_product(gdparttopng)
JLLWrappers.@declare_executable_product(gdtopng)
JLLWrappers.@declare_executable_product(giftogd2)
JLLWrappers.@declare_library_product(libgd, "@rpath/libgd.3.dylib")
JLLWrappers.@declare_executable_product(pngtogd)
JLLWrappers.@declare_executable_product(pngtogd2)
JLLWrappers.@declare_executable_product(webpng)
function __init__()
    JLLWrappers.@generate_init_header(JpegTurbo_jll, Zlib_jll, libpng_jll, Libtiff_jll, Xorg_libXpm_jll, Libiconv_jll, libwebp_jll)
    JLLWrappers.@init_executable_product(
        gd2copypal,
        "bin/gd2copypal",
    )

    JLLWrappers.@init_executable_product(
        gd2togif,
        "bin/gd2togif",
    )

    JLLWrappers.@init_executable_product(
        gd2topng,
        "bin/gd2topng",
    )

    JLLWrappers.@init_executable_product(
        gdcmpgif,
        "bin/gdcmpgif",
    )

    JLLWrappers.@init_executable_product(
        gdparttopng,
        "bin/gdparttopng",
    )

    JLLWrappers.@init_executable_product(
        gdtopng,
        "bin/gdtopng",
    )

    JLLWrappers.@init_executable_product(
        giftogd2,
        "bin/giftogd2",
    )

    JLLWrappers.@init_library_product(
        libgd,
        "lib/libgd.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        pngtogd,
        "bin/pngtogd",
    )

    JLLWrappers.@init_executable_product(
        pngtogd2,
        "bin/pngtogd2",
    )

    JLLWrappers.@init_executable_product(
        webpng,
        "bin/webpng",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module TableReader

export readdlm, readcsv, readtsv

using Dates:
    Date,
    DateTime,
    @dateformat_str
using Unicode:
    isletter
using DataFrames:
    DataFrame
using TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    NoopStream,
    Memory,
    Buffer,
    State,
    Noop,
    fillbuffer,
    makemargin!,
    buffermem
using CodecZlib:
    GzipDecompressorStream
using CodecZstd:
    ZstdDecompressorStream
using CodecXz:
    XzDecompressorStream

const SP = UInt8(' ')
const CR = UInt8('\r')
const LF = UInt8('\n')

include("stringcache.jl")
include("tokenizer.jl")
include("parser.jl")
include("download.jl")

# logarithm of base two
const MINIMUM_CHUNK_BITS = 14  # 16 KiB
const DEFAULT_CHUNK_BITS = 20  #  1 MiB
const MAXIMUM_CHUNK_BITS = 36  # 64 GiB

# Printable characters
const CHARS_PRINT = ' ':'~'

# Whitelist of delimiters and quotation marks
const ALLOWED_DELIMITERS = tuple(['\t'; CHARS_PRINT[.!(isletter.(CHARS_PRINT) .| isdigit.(CHARS_PRINT))]]...)
const ALLOWED_QUOTECHARS = tuple(CHARS_PRINT[.!(isletter.(CHARS_PRINT) .| isdigit.(CHARS_PRINT))]...)

# Special mark that never occurs in valid UTF-8 text.
const NO_QUOTECHAR = 0xff

"""
    readdlm(filename, command, or IO object;
            delim = nothing,
            quot = '"',
            trim = true,
            lzstring = true,
            skip = 0,
            skipblank = true,
            colnames = nothing,
            normalizenames = false,
            hasheader = (colnames === nothing),
            chunkbits = 20  #= 1 MiB =#)

Read a character delimited text file.

[`readcsv`](@ref) and [`readtsv`](@ref) call this function behind. To read a
CSV or TSV file, consider to use these dedicated function instead.

## Data source

The first (and the only positional) argument specifies the source to read data
from there.

If the argument is a string, it is considered as a local file name or the URL
of a remote file. If the name matches with `r"^\\w+://.*"` in regular
expression, it is handled as a URL. For example,
`"https://example.com/path/to/file.csv"` is regarded as a URL and its content
is streamed using the `curl` command.

If the argument is a command object, it is considered as a source whose
standard output is text data to read. For example, `unzip -p path/to/file.zip
somefile.csv` can be used to extract a file from a zipped archive. It is also
possible to pipeline several commands using `pipeline`.

If the arguments is an object of the `IO` type, it is considered as a direct
data source. The content is read using `read` or other similar functions. For
example, passing `IOBuffer(text)` makes it possible to read data from the raw
text object.

The data source is transparently decompressed if the compression format is
detectable. Currently, gzip, zstd, and xz are supported. The format is detected
by the magic bytes of the stream header, and therefore other information such
as file names does not affect the detection.

## Parser parameters

`delim` specifies the field delimiter in a line.  This cannot be the same
character as `quot`.  If `delim` is `nothing`, the parser tries to guess a
delimiter from data.  Currently, the following delimiters are allowed:
$(join(repr.(ALLOWED_DELIMITERS), ", ")).

`quot` specifies the quotation character to enclose a field. This cannot be the
same character as `delim`. If `quot` is `nothing`, no characters are recognized
as a quotation mark. Currently, the following quotation characters are allowed:
$(join(repr.(ALLOWED_QUOTECHARS), ", ")).

`trim` specifies whether the parser trims space (0x20) characters around a field.
If `trim` is true, `delim` and `quot` cannot be a space character.

`lzstring` specifies whether fields with excess leading zeros are treated as
strings.  If `lzstring` is true, fields such as "0003" will be interpreted as
strings instead of integers.

`skip` specifies the number of lines to skip before reading data.  The next
line just after the skipped lines is considered as a header line if the
`colnames` parameter is not specified.

`skipblank` specifies whether the parser ignores blank lines. If `skipblank` is
false, encountering a blank line throws an exception.

`comment` specifies the leading sequence of comment lines. If it is a non-empty
string, text lines that start with the sequence will be skipped as comments.
The default value (empty string) does not skip any lines as comments.

## Column names

`colnames` specifies the column names. If `colnames` is `nothing` (default),
the column names are read from the first line just after skipping lines
specified by `skip` (no lines are skipped by default). Any iterable object is
allowed.

`normalizenames` uses 'safe' names for Julia symbols used in DataFrames.
If `normalizenames` is `false` (default), the column names will be the same
as the source file using basic parsing. If `normalizenames` is `true` then
reserved words and characters will be removed/replaced in the column names.

`hasheader` specified whether the data has a header line or not. The default
value is `colnames === nothing` and thus the parser assumes there is a header
if and only if no column names are specified.

The following table summarizes the behavior of the `colnames` and `hasheader`
parameters.

| `colnames` | `hasheader` | column names |
|:-----------|:------------|:-------------|
| `nothing` | `true`  | taken from the header (default) |
| `nothing` | `false` | automatically generated (X1, X2, ...) |
| specified | `true`  | taken from `colnames` (the header line is skipped) |
| specified | `false` | taken from `colnames` |

If unnamed columns are found in the header, they are renamed to `UNNAMED_{j}`
for ease of access, where `{j}` is replaced by the column number. If the number
of header columns in a file is less than the number of data columns by one, a
column name `UNNAMED_0` will be inserted into the column names as the first
column.  This is useful to read files written by the `write.table` function of
R with `row.names = TRUE`.

## Data types

Integers, floating-point numbers, boolean values, dates, datetimes, missings,
and strings are automatically detected and converted from the text data.  The
following list is the summary of the corresponding data types of Julia and the
text formats described in the regular expression:

- Integer (`Int`): `[-+]?\\d+`
- Float (`Float64`): `[-+]?\\d*\\.?\\d+`, `[-+]?\\d*\\.?\\d+([eE][-+]?\\d+)?`,
                     `[-+]?NaN` or `[-+]?Inf(inity)?` (case-insensitive)
- Bool (`Bool`): `t(rue)?` or `f(alse)?` (case-insensitive)
- Date (`Dates.Date`): `\\d{4}-\\d{2}-\\d{2}`
- Datetime (`Dates.DateTime`): `\\d{4}-\\d{2}-\\d{2}[T ]\\d{2}:\\d{2}:\\d{2}(\\.\\d+)?`
- Missing (`Missing`): empty field or `NA` (case-sensitive)
- String (`String`): otherwise

Integers and floats have some overlap. The parser precedes integers over
floats.  That means, if all values in a column are parsable as integers and
floats, they are parsed as integers instead of floats; otherwise, they are
parsed as floats. Similarly, all the types have higher precedence than strings.

The parser parameter `lzstring` affects interpretation of numbers. If
`lzstring` is true, numbers with excess leading zeros (e.g., "0001", "00.1")
are interpreted as strings. Fields without excess leading zeros (e.g., "0",
"0.1") are interepreted as numbers regardless of this parameter.

## Parsing behavior

The only supported text encoding of a file is UTF-8, which is the default
character encoding scheme of many functions in Julia.  If you need to read text
encoded other than UTF-8, it is required to wrap the data stream with an
encoding conversion tool such as the `iconv` command or StringEncodings.jl.

```julia
# Convert text encoding from Shift JIS (Japanese) to UTF8.
readcsv(`iconv -f sjis -t utf8 somefile.csv`)
```

A text file will be read chunk by chunk to save memory. The chunk size is
specified by the `chunkbits` parameter, which is the base two logarithm of
actual chunk size.  The default value is 20 (i.e., 2^20 bytes = 1 MiB).  The
data type of each column is guessed from the values in the first chunk.  If
`chunkbits` is set to zero, it disables chunking and the data types are guessed
from all rows. The chunk size will be automatically expanded when it is
required to store long lines.

A chunk cannot be larger than 64 GiB and a field cannot be longer than 16 MiB.
These limits are due to the encoding method of tokens used by the tokenizer.
Therefore, you cannot parse data larger than 64 GiB without chunking and fields
longer than 16 MiB. Trying to read such a file will result in error.
"""
function readdlm end

"""
    readcsv(filename, command, or IO object; delim = ',', <keyword arguments>)

Read a CSV (comma-separated values) text file.

This function is the same as [`readdlm`](@ref) but with `delim = ','`.
See `readdlm` for details.
"""
function readcsv end

"""
    readtsv(filename, command, or IO object; delim = '\\t', <keyword arguments>)

Read a TSV (tab-separated values) text file.

This function is the same as [`readdlm`](@ref) but with `delim = '\\t'`.
See `readdlm` for details.
"""
function readtsv end

for (fname, delim) in [(:readdlm, nothing), (:readcsv, ','), (:readtsv, '\t')]
    # prepare keyword arguments
    kwargs = Expr[]
    push!(kwargs, Expr(:kw, :(delim::Union{Char,Nothing}), delim))  # delim::Union{Char,Nothing} = $(delim)
    push!(kwargs, Expr(:kw, :(quot::Union{Char,Nothing}), '"'))  # quot::Union{Char,Nothing} = '"'
    push!(kwargs, Expr(:kw, :(trim::Bool), true))  # trim::Bool = true
    push!(kwargs, Expr(:kw, :(lzstring::Bool), true))  # lzstring::Bool = true
    push!(kwargs, Expr(:kw, :(skip::Integer), 0))  # skip::Integer = 0
    push!(kwargs, Expr(:kw, :(skipblank::Bool), true))  # skipblank::Bool = true
    push!(kwargs, Expr(:kw, :(comment::AbstractString), ""))  # command::AbstractString = ""
    push!(kwargs, Expr(:kw, :(colnames), nothing))  # colnames = nothing
    push!(kwargs, Expr(:kw, :(normalizenames), false))  # normalizenames::Bool = false
    push!(kwargs, Expr(:kw, :(hasheader::Bool), :(colnames === nothing)))  # hasheader::Bool = (colnames === nothing)
    push!(kwargs, Expr(:kw, :(chunkbits::Integer), DEFAULT_CHUNK_BITS))  # chunkbits::Integer = DEFAULT_CHUNK_BITS

    # deprecated parameter
    push!(kwargs, Expr(:kw, :(chunksize::Union{Integer,Nothing}), nothing))

    # generate methods
    @eval begin
        function $(fname)(filename::AbstractString; $(kwargs...))
            if chunksize !== nothing
                @warn "the chunksize parameter is deprecated; use chunkbits instead (note that chunksize = 2^chunkbits)"
                chunkbits = chunksize == 0 ? 0 : trailing_zeros(nextpow(2, chunksize))
            end
            params = LexerParameters(delim, quot, trim, lzstring, skip, skipblank, comment, colnames, normalizenames, hasheader, chunkbits)
            if occursin(r"^\w+://", filename)  # URL-like filename
                curl = find_curl()
                if curl === nothing
                    throw(ArgumentError("the curl command is not available"))
                end
                # read a remote file using curl
                return open(
                    proc -> readdlm_internal(wrapstream(proc, params), params),
                    `$(curl) --silent --location --globoff $(filename)`,
                )
            end
            # read a local file
            return open(file -> readdlm_internal(wrapstream(file, params), params), filename)
        end

        function $(fname)(cmd::Base.AbstractCmd; $(kwargs...))
            if chunksize !== nothing
                @warn "the chunksize parameter is deprecated; use chunkbits instead (note that chunksize = 2^chunkbits)"
                chunkbits = chunksize == 0 ? 0 : trailing_zeros(nextpow(2, chunksize))
            end
            params = LexerParameters(delim, quot, trim, lzstring, skip, skipblank, comment, colnames, normalizenames, hasheader, chunkbits)
            return open(proc -> readdlm_internal(wrapstream(proc, params), params), cmd)
        end

        function $(fname)(file::IO; $(kwargs...))
            if chunksize !== nothing
                @warn "the chunksize parameter is deprecated; use chunkbits instead (note that chunksize = 2^chunkbits)"
                chunkbits = chunksize == 0 ? 0 : trailing_zeros(nextpow(2, chunksize))
            end
            params = LexerParameters(delim, quot, trim, lzstring, skip, skipblank, comment, colnames, normalizenames, hasheader, chunkbits)
            return readdlm_internal(wrapstream(file, params), params)
        end
    end
end

# Wrap a stream by TranscodingStream.
function wrapstream(stream::IO, params::LexerParameters)
    bufsize = 2^max(params.chunkbits, MINIMUM_CHUNK_BITS)
    if !applicable(mark, stream) || !applicable(reset, stream)
        stream = NoopStream(stream, bufsize = bufsize)
    end
    format = checkformat(stream)
    if params.chunkbits != 0
        # with chunking
        if format == :gzip
            return GzipDecompressorStream(stream, bufsize = bufsize)
        elseif format == :zstd
            return ZstdDecompressorStream(stream, bufsize = bufsize)
        elseif format == :xz
            return XzDecompressorStream(stream, bufsize = bufsize)
        elseif stream isa TranscodingStream
            return stream
        else
            return NoopStream(stream, bufsize = bufsize)
        end
    end
    # without chunking
    if format == :gzip
        data = read(GzipDecompressorStream(stream))
    elseif format == :zstd
        data = read(ZstdDecompressorStream(stream))
    elseif format == :xz
        data = read(XzDecompressorStream(stream))
    else
        data = read(stream)
    end
    buffer = Buffer(data)
    return TranscodingStream(Noop(), devnull, State(buffer, buffer))
end

# Check the file format of a stream.
function checkformat(stream::IO)
    mark(stream)
    magic = zeros(UInt8, 6)
    nb = readbytes!(stream, magic, 6)
    reset(stream)
    if nb != 6
        return :unknown
    elseif magic[1:6] == b"\xFD\x37\x7A\x58\x5A\x00"
        return :xz
    elseif magic[1:2] == b"\x1f\x8b"
        return :gzip
    elseif magic[1:4] == b"\x28\xb5\x2f\xfd"
        return :zstd
    else
        return :unknown
    end
end

# The main function of parsing a character delimited file.
# `stream` is asuumed to be an input stream of plain text.
function readdlm_internal(stream::TranscodingStream, params::LexerParameters)
    # Skip lines.
    line = skiplines(stream, params.skip) + 1
    line += skip_unwanted_lines(stream, params)

    # Determine delimiter.
    if params.delim === nothing
        params = LexerParameters(
            Char(guessdelimiter(stream, params)),
            Char(params.quot),
            params.trim,
            params.lzstring,
            params.skip,
            params.skipblank,
            params.comment,
            params.colnames,
            params.normalizenames,
            params.hasheader,
            params.chunkbits,
        )
    end
    @assert params.delim isa UInt8

    # Determine column names
    mem = bufferlines(stream)
    if params.hasheader
        if params.colnames === nothing
            # Scan the header line to get the column names
            n, headertokens = scanheader(mem, params)
            skip(stream, n)
            if length(headertokens) == 1 && length(headertokens[1]) == 0  # zero-length token
                throw(ReadError("found no column names in the header at line $(line)"))
            end
            line += 1
            colnames = Symbol[]
            for (i, token) in enumerate(headertokens)
                start, length = location(token)
                if (kind(token) & QSTRING) != 0
                    name = qstring(mem, start, length, params.quot)
                else
                    name = unsafe_string(mem.ptr + start - 1, length)
                end
                if isempty(name)  # unnamed column
                    name = "UNNAMED_$(i)"
                end
                push!(colnames, Symbol(name))
            end
        else
            # skip the header line
            line += skiplines(stream, 1)
            colnames = params.colnames
        end
    else
        # no header line
        if params.colnames === nothing
            # count the number of columns from data
            n_max_cols = countbytesline(mem, params.delim) + 1
            _, n_cols, _ = scanline!(
                Array{Token}(undef, (n_max_cols, 1)), 1, mem, 0, line, params)
            colnames = [Symbol("X", i) for i in 1:n_cols]
        else
            colnames = params.colnames
        end
    end
    ncols = length(colnames)

    # Scan the next line to get the number of data columns
    line += skip_unwanted_lines(stream, params)
    mem = bufferlines(stream)
    tokens = Array{Token}(undef, (ncols + 1, 1))
    _, i, _ = scanline!(tokens, 1, mem, 0, line, params)
    if i == 1 && location(tokens[1,1])[2] == 0
        # no data
        return DataFrame([[] for _ in 1:length(colnames)], colnames, copycols = false, makeunique = true)
    elseif i == ncols
        # the header and the first row have the same number of columns
    elseif i == ncols + 1
        # the first column is supposed to be unnamed
        ncols += 1
        pushfirst!(colnames, :UNNAMED_0)
    else
        throw(ReadError("unexpected number of columns at line $(line)"))
    end

    # Estimate the number of rows.
    nrows_estimated = countbytes(mem, LF)
    if nrows_estimated == 0
        nrows_estimated = countbytes(mem, CR)
    end
    n_chunk_rows = max(nrows_estimated, 5)  # allocate at least five rows

    # Read data.
    tokens = Array{Token}(undef, (ncols, n_chunk_rows))
    columns = Vector[]
    while !eof(stream)
        # Tokenize data.
        mem = bufferlines(stream)
        pos = 0
        n_new_rows = 0
        while pos < lastindex(mem) && n_new_rows < n_chunk_rows
            pos, i, skip = scanline!(tokens, n_new_rows + 1, mem, pos, line, params)
            if pos == 0
                break
            elseif skip  # blank or comment
                line += 1
                continue
            elseif i != ncols
                throw(ReadError("unexpected number of columns at line $(line)"))
            end
            n_new_rows += 1
            line += 1
        end
        if n_new_rows == 0
            # only blank or comment lines
            Base.skip(stream, pos)
            continue
        end

        # Parse data.
        bitmaps = summarizecolumns(tokens, n_new_rows)
        if isempty(columns)
            # infer data types of columns
            resize!(columns, ncols)
            for i in 1:ncols
                T = datatype(bitmaps[i])
                @debug "Filling $(colnames[i])::$(T) column"
                columns[i] = fillcolumn!(
                    Vector{T}(undef, n_new_rows),
                    n_new_rows, mem, tokens, i, params.quot)
            end
        else
            # check existing columns
            for i in 1:ncols
                col = columns[i]
                S = eltype(col)
                T = datatype(bitmaps[i])
                U = Union{S,T}
                if (S <: Union{Int,Missing} && T <: Union{Float64,Missing}) ||
                   (T <: Union{Int,Missing} && S <: Union{Float64,Missing})
                    U = promote_type(S, T)
                elseif S <: Union{String,Missing} && String <: S
                    # any field is interpretable as string
                    if Missing <: U
                        U = Union{String,Missing}
                    else
                        U = String
                    end
                elseif !(U <: S || U <: T)
                    throw(ReadError(string(
                        "type guessing failed at column $(i) ",
                        "(guessed to be $(S) but found records of $(T)); ",
                        "try larger chunkbits or chunkbits = 0 to disable chunking")))
                end
                n_rows = length(col)
                if U <: S
                    resize!(col, n_rows + n_new_rows)
                else
                    col = copyto!(Vector{U}(undef, n_rows + n_new_rows), 1, col, 1, n_rows)
                end
                @debug "Filling $(colnames[i])::$(U) column"
                columns[i] = fillcolumn!(col, n_new_rows, mem, tokens, i, params.quot)
            end
        end
        Base.skip(stream, pos)
    end

    # Parse strings as date or datetime objects.
    for i in 1:ncols
        col = columns[i]
        if eltype(col) <: Union{String,Missing}
            if is_date_like(col)
                try
                    columns[i] = parse_date(col)
                catch
                    # not a date column
                end
            elseif is_datetime_like(col)
                hasT = is_T_delimited(col)  # check delimited by T or space
                try
                    columns[i] = parse_datetime(col, hasT)
                catch
                    # not a datetime column
                end
            end
        end
    end

    # Normalize column names to remove/convert unfriendly characters
    if params.normalizenames
        colnames = [normalizename(String(name)) for name in colnames]
    end

    return DataFrame(columns, colnames, copycols = false, makeunique = true)
end

# Count the number of `byte` in a memory block.
function countbytes(mem::Memory, byte::UInt8)
    n = 0
    @inbounds @simd for i in 1:length(mem)
        n += mem[i] == byte
    end
    return n
end

# Count the number of `byte` in a line.
function countbytesline(mem::Memory, byte::UInt8)
    n = 0
    for i in 1:lastindex(mem)
        @inbounds x = mem[i]
        if x == CR || x == LF
            break
        elseif x == byte
            n += 1
        end
    end
    return n
end

function datatype(bitmap::UInt8)
    if (bitmap & 0b010000) != 0  # all values are missing
        return Missing
    end
    T = (bitmap & INTEGER) != 0 ? Int     :
        (bitmap & FLOAT)   != 0 ? Float64 :
        (bitmap & BOOL)    != 0 ? Bool    : String
    if (bitmap & 0b100000) != 0  # at least one value is missing
        T = Union{T,Missing}
    end
    return T
end

# Summarize columns using bitmaps.
function summarizecolumns(tokens::Matrix{Token}, nrows::Int)
    # From the least significant bit:
    #   1. INTEGER
    #   2. FLOAT
    #   3. BOOL
    #   4. QSTRING
    #   5. ∀ missing
    #   6. ∃ missing
    ncols = size(tokens, 1)
    bitmaps = Vector{UInt8}(undef, ncols)
    fill!(bitmaps, 0b011111)
    @inbounds for j in 1:nrows, i in 1:ncols
        # Note that the tokens matrix is transposed.
        x = kind(tokens[i,j])
        ismissing = x == 0b1111
        y = bitmaps[i]
        bitmap = 0b000000
        bitmap |= (y & 0b100000) | ifelse(ismissing, 0b100000, 0b000000)
        bitmap |= (y & 0b010000) & ifelse(ismissing, 0b010000, 0b000000)
        bitmap |= (y & 0b001111) & x
        bitmaps[i] = bitmap
    end
    return bitmaps
end

# Skip the number of lines specified by `skip`.
function skiplines(stream::TranscodingStream, skip::Int)
    skipped = 0
    while skipped < skip && !eof(stream)
        mem = bufferlines(stream)
        # find the first newline
        i = 1
        while i ≤ lastindex(mem)
            @inbounds x = mem[i]
            if x == CR
                if i + 1 ≤ lastindex(mem) && mem[i+1] == LF
                    i += 1
                end
                break
            elseif x == LF
                break
            end
            i += 1
        end
        Base.skip(stream, i)
        skipped += 1
    end
    return skipped
end

# Skip blank and/or comment lines if required.
function skip_unwanted_lines(stream::TranscodingStream, params::LexerParameters)
    skipped = 0
    @label loop
    n = 0
    if params.skipblank
        n += skipblanklines(stream, params.trim)
    end
    if !isempty(params.comment)
        n += skipcommentlines(stream, params.comment)
    end
    if n > 0
        skipped += n
        @goto loop
    end
    return skipped
end

# Skip blank lines.
function skipblanklines(stream::TranscodingStream, trim::Bool)
    skipped = 0
    while !eof(stream)
        mem = bufferlines(stream)
        @assert length(mem) > 0
        i = 1
        if trim
            while i ≤ lastindex(mem) && mem[i] == SP
                i += 1
            end
        end
        if i + 1 ≤ lastindex(mem) && mem[i] == CR && mem[i+1] == LF
            skip(stream, i + 1)
        elseif i ≤ lastindex(mem) && (mem[i] == CR || mem[i] == LF)
            skip(stream, i)
        else
            # found a non-blank line
            break
        end
        skipped += 1
    end
    return skipped
end

# Skip comment lines.
function skipcommentlines(stream::TranscodingStream, comment::String)
    @assert !isempty(comment)
    skipped = 0
    while !eof(stream)
        mem = bufferlines(stream)
        i = 1
        while i ≤ lastindex(mem) &&
                i ≤ sizeof(comment) &&
                mem[i] != CR &&
                mem[i] != LF &&
                mem[i] == codeunit(comment, i)
            i += 1
        end
        if i ≤ sizeof(comment)
            # found a non-comment line
            break
        end
        # found a comment; skip the current newline
        while i ≤ lastindex(mem) && mem[i] != CR && mem[i] != LF
            i += 1
        end
        if i + 1 ≤ lastindex(mem) && mem[i] == CR && mem[i+1] == LF
            skip(stream, i + 1)
        elseif i ≤ lastindex(mem) && (mem[i] == CR || mem[i] == LF)
            skip(stream, i)
        end
        skipped += 1
    end
    return skipped
end

# Guess the delimiter byte from data.  The current logic is very simple:
# choosing the most frequent one in the first line from the list of candidates.
function guessdelimiter(stream::TranscodingStream, params::LexerParameters)
    mem = bufferlines(stream)
    n_max = 0
    delim = UInt8(',')
    for d in UInt8.((',', '\t', '|', ';', ':'))
        n = countbytesline(mem, d)
        if n > n_max
            n_max = n
            delim = d
        end
    end
    return delim
end

# Buffer lines and return the memory view and the last newline position.
function bufferlines(stream::TranscodingStream)
    # try to find the last newline (LF, CR, or CR+LF)
    @label SEARCH_NEWLINE
    fillbuffer(stream, eager = true)
    buffer = stream.state.buffer1
    mem = buffermem(buffer)
    lastnl = lastindex(mem)
    while lastnl > 0
        @inbounds x = mem[lastnl]
        if x == LF
            break
        elseif x == CR
            if lastnl == lastindex(mem)  # cannot determine CR or CR+LF
                makemargin!(buffer, 1)
                fillbuffer(stream, eager = true)
                mem = buffermem(buffer)
                lastnl = lastindex(mem)
                while lastnl > 0
                    @inbounds x = mem[lastnl]
                    if x == LF || x == CR
                        break
                    end
                    lastnl -= 1
                end
            end
            break
        end
        lastnl -= 1
    end
    if lastnl == 0  # found no newlines in the buffered chunk
        if TranscodingStreams.marginsize(buffer) == 0
            # No newlines in the current buffer but newlines may be in the next.
            bufsize = length(buffer)
            if bufsize ≥ MAX_TOKEN_START
                throw(ReadError("found a too long line to store in a chunk"))
            end
            newbufsize = min(bufsize * 2, MAX_TOKEN_START)
            TranscodingStreams.makemargin!(buffer, newbufsize - bufsize)
            @goto SEARCH_NEWLINE
        else
            # Found EOF without newlines; terminate the buffered data with LF.
            TranscodingStreams.writebyte!(buffer, LF)
            mem = buffermem(buffer)
            lastnl = lastindex(mem)
        end
    end
    @assert lastnl > 0 && (mem[lastnl] == CR || mem[lastnl] == LF)
    return Memory(mem.ptr, lastnl)
end

# Generated from tools/snoop.jl.
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(TableReader.scanline!), Array{TableReader.Token, 2}, Int64, TranscodingStreams.Memory, Int64, Int64, TableReader.LexerParameters})
    precompile(Tuple{typeof(TableReader.scanheader), TranscodingStreams.Memory, TableReader.LexerParameters})
    precompile(Tuple{typeof(TableReader.bufferlines), TranscodingStreams.TranscodingStream{TranscodingStreams.Noop, Base.IOStream}})
    precompile(Tuple{typeof(TableReader.checkformat), Base.IOStream})
    precompile(Tuple{getfield(TableReader, Symbol("#24#29")), Base.Process})
    precompile(Tuple{typeof(TableReader.skipcommentlines), TranscodingStreams.TranscodingStream{TranscodingStreams.Noop, Base.IOStream}, String})
    precompile(Tuple{typeof(TableReader.skipblanklines), TranscodingStreams.TranscodingStream{TranscodingStreams.Noop, Base.IOStream}, Bool})
    precompile(Tuple{typeof(TableReader.fillcolumn!), Array{Int64, 1}, Int64, TranscodingStreams.Memory, Array{TableReader.Token, 2}, Int64, UInt8})
    precompile(Tuple{typeof(TableReader.checkformat), TranscodingStreams.TranscodingStream{TranscodingStreams.Noop, Base.IOStream}})
    precompile(Tuple{typeof(TableReader.fillcolumn!), Array{Float64, 1}, Int64, TranscodingStreams.Memory, Array{TableReader.Token, 2}, Int64, UInt8})
    precompile(Tuple{typeof(TableReader.readdlm_internal), TranscodingStreams.TranscodingStream{TranscodingStreams.Noop, Base.IOStream}, TableReader.LexerParameters})
    precompile(Tuple{typeof(TableReader.qstring), TranscodingStreams.Memory, Int64, Int64, UInt8})
    precompile(Tuple{typeof(TableReader.summarizecolumns), Array{TableReader.Token, 2}, Int64})
    precompile(Tuple{typeof(TableReader.checkformat), TranscodingStreams.TranscodingStream{TranscodingStreams.Noop, Base.Process}})
    precompile(Tuple{typeof(TableReader.parse_date), Array{String, 1}})
    precompile(Tuple{typeof(TableReader.parse_datetime), Array{String, 1}, Bool})
    precompile(Tuple{typeof(TableReader.wrapstream), Base.Process, TableReader.LexerParameters})
    precompile(Tuple{typeof(TableReader.wrapstream), Base.IOStream, TableReader.LexerParameters})
    precompile(Tuple{typeof(TableReader.checkformat), Base.Process})
    precompile(Tuple{getfield(TableReader, Symbol("#25#30")), Base.IOStream})
    precompile(Tuple{typeof(TableReader.allocate!), TableReader.StringCache, Ptr{UInt8}, Int64})
    precompile(Tuple{typeof(TableReader.fillcolumn!), Array{String, 1}, Int64, TranscodingStreams.Memory, Array{TableReader.Token, 2}, Int64, UInt8})
    precompile(Tuple{typeof(TableReader.is_date_like), Array{String, 1}})
    precompile(Tuple{typeof(TableReader.is_datetime_like), Array{String, 1}})
    precompile(Tuple{typeof(TableReader.skip_unwanted_lines), TranscodingStreams.TranscodingStream{TranscodingStreams.Noop, Base.IOStream}, TableReader.LexerParameters})
    precompile(Tuple{typeof(TableReader.countbytes), TranscodingStreams.Memory, UInt8})
    precompile(Tuple{typeof(TableReader.parse_datetime), Array{Union{Base.Missing, String}, 1}, Bool})
    precompile(Tuple{getfield(TableReader, Symbol("##readcsv#14")), Char, Char, Bool, Bool, Int64, Bool, String, Nothing, Bool, Bool, Int64, Nothing, typeof(TableReader.readcsv), String})
    precompile(Tuple{typeof(TableReader.parse_date), Array{Union{Base.Missing, String}, 1}})
    precompile(Tuple{typeof(TableReader.skiplines), TranscodingStreams.TranscodingStream{TranscodingStreams.Noop, Base.IOStream}, Int64})
    precompile(Tuple{getfield(TableReader, Symbol("##readtsv#23")), Char, Char, Bool, Bool, Int64, Bool, String, Nothing, Bool, Bool, Int64, Nothing, typeof(TableReader.readtsv), String})
    precompile(Tuple{typeof(TableReader.normalizename), String})
    precompile(Tuple{getfield(TableReader, Symbol("#15#20")), Base.Process})
    precompile(Tuple{getfield(TableReader, Symbol("#16#21")), Base.IOStream})
    precompile(Tuple{typeof(TableReader.readcsv), String})
    precompile(Tuple{typeof(TableReader.readtsv), String})
end
_precompile_()

end # module


text:
# gll_library
function zwgljd!(z::Vector{Float64}, w::Vector{Float64}, np::Int64, alpha::Float64, beta::Float64)
    # convert to references
    np_ref = Ref{Int32}(Int32(np))
    alpha_ref = Ref{Float64}(alpha)
    beta_ref = Ref{Float64}(beta)

    ccall((:zwgljd_, "../../fortran/lib/gll_library.so"), Cvoid, (Ptr{Float64}, Ptr{Float64}, Ref{Int32}, Ref{Float64}, Ref{Float64}), z, w, np_ref, alpha_ref, beta_ref)
end

# sem_mesh_mod
function anchor_point_index!(iax::Vector{Int32}, iay::Vector{Int32}, iaz::Vector{Int32})
    iax_32 = Int32.(iax)
    iay_32 = Int32.(iay)
    iaz_32 = Int32.(iaz)
    @assert length(iax) == 27 && length(iay) == 27 && length(iaz) == 27
    ccall((:anchor_point_index_, "../../fortran/lib/sem_mesh_mod.so"), Cvoid, (Ptr{Int32}, Ptr{Int32}, Ptr{Int32}), iax_32, iay_32, iaz_32)
    iax .= Int32.(iax_32)
    iay .= Int32.(iay_32)
    iaz .= Int32.(iaz_32)
end

# sem_mesh_mod
function xyz2cube_bounded!(xyz_anchor::Array{Float64,2}, xyz::Array{Float64,1}, uvw::Array{Float64,1}, misloc::Float64, flag_inside::Bool)
    misloc_ref = Ref{Float64}(misloc)
    flag_inside_ref = Ref{Bool}(flag_inside)
    ccall((:xyz2cube_bounded_, "../../fortran/lib/sem_mesh_mod.so"), Cvoid, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Float64}, Ref{Bool}), xyz_anchor, xyz, uvw, misloc_ref, flag_inside_ref)
    return misloc_ref.x, flag_inside_ref.x
end

# sem_mesh_mod
function lagrange_poly!(x::Float64, ngll::Int64, xgll::Vector{Float64}, lagrange::Vector{Float64})
    x_ref = Ref{Float64}(x)
    ngll_ref = Ref{Int32}(Int32(ngll))
    ccall((:lagrange_poly_, "../../fortran/lib/sem_mesh_mod.so"), Cvoid, (Ref{Float64}, Ref{Int32}, Ptr{Float64}, Ptr{Float64}), x_ref, ngll_ref, xgll, lagrange)
end

max_stars_repo_name: ziyixi/SeisScripts

max_issues_repo_name: ziyixi/SeisScripts

max_forks_repo_name: ziyixi/SeisScripts


text:
using Distances

const DEBUG = length(ARGS) >=1 && ARGS[1] == "1"

function distance(solution::Array{Int32})
    cost = 0
    for i in 1:length(solution)-1
        cost += d[solution[i],solution[i+1]]
    end
    return cost
end

function tab(i::Integer, j::Integer, iteration::Integer)
    tabulist[i] = iteration
    tabulist[j] = iteration
end

function istabbed(i::Integer, j::Integer, iteration::Integer)
    return (iteration - tabulist[i]) < penalty || (iteration - tabulist[j]) < penalty
end

function isaspiring(i::Integer, j::Integer, current::Float64, base::Float64, itr::Integer)
    return base > current + 0.7(penalty-(itr-min(tabulist[i], tabulist[j])))
end

tic()
tic()
# file lines iterator object
input = eachline(STDIN)

state = start(input)
# read data size N
(n, state) = next(input, state)
n = parse(Int32, n)

# 2xN node's coordinates matrix
data = Array{Float64}(2, n)

# read file data to Data matrix
while !done(input, state)
    (v, state) = next(input, state)
    (i_v, x_v, y_v) = split(v)
    data[:,parse(Int32, i_v)] = [parse(Float64, x_v), parse(Float64, y_v)]
end
# calculate nodes distances as NxN matrix

d = Array{Float64}(n, n)
pairwise!(d, Euclidean(1e-12), data)

currentsol = Array{Int32}(n+1)
currentsol[1] = currentsol[n+1] = Int32(1)

for i in 1:n
    d[i,i] = Inf
end

d1 = copy(d)
k = 1
for i in 2:length(currentsol)-1
    currentsol[i] = findfirst(d1[k,:].==min(d1[k,:]...))
    k = currentsol[i]
    d1[k,1] = Inf
    for visited_node in 1:i
        d1[k, currentsol[visited_node]] = Inf
    end
end

currentcost = distance(currentsol)

bestsol = copy(currentsol)
bestcost = currentcost

bestneighborsol = copy(currentsol)
bestneighborcost = currentcost

N = [(j, k) for j=2:length(bestsol)-2 for k=j+1:length(bestsol)-1 if k-j>1]

neighborscount = length(N)

####### SETTINGS #######
const itermax = round(Int32, 1.0e7*round(sqrt(n))/n^2)

const neighborstocheck = round(Int32, sqrt(neighborscount)/log(100, n))

const penalty = round(Int32, 0.7sqrt(n)) # rounds tabbed
########################

tabulist = fill(-penalty, n)
checks = 0

if DEBUG
    println(string("max iterations:\t", itermax))
    println(string("initial cost:\t", bestcost))
    # for n in bestsol
    #     println(STDERR, n)
    # end
    println(string("initial time:\t", toq(), "s"))
end

for iter in 1:itermax
    neighborsol = copy(bestneighborsol)
    neighborcost = bestneighborcost

    checks = 0
    maxchecks = neighborstocheck
    index1 = 0
    index2 = 0

    for (j, k) in shuffle(N)
        city1 = neighborsol[j]
        city2 = neighborsol[k]

        if checks < maxchecks
            currentcost = neighborcost - d[neighborsol[j-1], neighborsol[j]] - d[neighborsol[j], neighborsol[j+1]] + d[neighborsol[j-1], neighborsol[k]] + d[neighborsol[k], neighborsol[j+1]] - d[neighborsol[k-1], neighborsol[k]] - d[neighborsol[k], neighborsol[k+1]] + d[neighborsol[k-1], neighborsol[j]] + d[neighborsol[j], neighborsol[k+1]]

            if !istabbed(city1, city2, iter) || isaspiring(city1, city2, currentcost, neighborcost, iter)
                if currentcost <= bestneighborcost || index1 == 0
                    index1 = j
                    index2 = k
                    bestneighborcost = currentcost
                    if currentcost < bestcost
                        # maxchecks = min(maxchecks, checks+0.01neighborstocheck)
                        break
                    end
                end
                checks += 1
            end
        end
        tab(city1, city2, iter)
    end

    if index1 != 0 && index2 != 0
        bestneighborsol = copy(neighborsol)
        bestneighborsol[index2], bestneighborsol[index1] = bestneighborsol[index1], bestneighborsol[index2]

        if bestneighborcost <= bestcost
            bestsol = copy(bestneighborsol)
            bestcost = bestneighborcost
        end
    else
        break
    end
end

if DEBUG
    println(string("best cost:\t", bestcost))
    println(string("total time:\t", toq(), "s"))
else
    println(STDOUT, bestcost)
    for n in bestsol
        print(STDERR, n)
        print(STDERR, " ")
    end
    println(STDERR)
end

max_stars_repo_name: AfroMetal/algorithms-tsp-tabusearch

max_issues_repo_name: AfroMetal/algorithms-tsp-tabusearch

max_forks_repo_name: AfroMetal/algorithms-tsp-tabusearch


text:
# Autogenerated wrapper script for HelFEM_jll for x86_64-w64-mingw32-cxx03
export libhelfem

using libcxxwrap_julia_jll
using armadillo_jll
using GSL_jll
using OpenBLAS_jll
JLLWrappers.@generate_wrapper_header("HelFEM")
JLLWrappers.@declare_library_product(libhelfem, "libhelfem-cxxwrap.dll")
function __init__()
    JLLWrappers.@generate_init_header(libcxxwrap_julia_jll, armadillo_jll, GSL_jll, OpenBLAS_jll)
    JLLWrappers.@init_library_product(
        libhelfem,
        "bin\\libhelfem-cxxwrap.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#time_steps = 1:24

BASE_DIR = abspath(joinpath(dirname(Base.find_package("WaterSystems")), ".."))
TEST_DIR = joinpath(BASE_DIR, "test")
DATA_DIR = joinpath(TEST_DIR, "data")

epanet_test_file = "epanetfiles/VanZyl.inp"
epanet_test_path = joinpath(DATA_DIR, epanet_test_file)


text:
import RestrictedBoltzmannMachines as RBMs
import ConvolutionalRBMs as ConvRBMs
using Test: @test, @testset, @inferred, @test_throws
using Random: bitrand
using LogExpFunctions: logsumexp
using RestrictedBoltzmannMachines: energy, free_energy, interaction_energy
using RestrictedBoltzmannMachines: inputs_v_to_h, inputs_h_to_v
using RestrictedBoltzmannMachines: sample_v_from_v, sample_h_from_h
using RestrictedBoltzmannMachines: sample_v_from_h, sample_h_from_v
using RestrictedBoltzmannMachines: Binary, visible, hidden, weights
using RestrictedBoltzmannMachines: reshape_maybe
using ConvolutionalRBMs: vsizes, hsizes, parts, output_size, inputs_h_to_v_dense, inputs_h_to_v_conv

@testset "DenseConvRBM" begin
    N = (4,)
    C = (2,)
    J = (3,)
    K = (5,)
    M = (7)
    B = (11,)

    rbm = ConvRBMs.DenseConvRBM(
        Binary(randn(N...)), # visible dense
        Binary(randn(C...)), # visible conv
        Binary(randn(M...)), # hidden
        randn(N..., M...), # weights dense
        randn(C..., J..., M...) # weights conv
    )

    v_dense = bitrand(N..., B...)
    v_conv = bitrand(C..., (J .+ K .- 1)..., B...)
    h = bitrand(M..., K..., B...)
    @test size(@inferred inputs_v_to_h(rbm, v_dense, v_conv)) == size(h)
    @test size(@inferred inputs_h_to_v_dense(rbm, h)) == size(v_dense)
    @test size(@inferred inputs_h_to_v_conv(rbm, h)) == size(v_conv)

    h_reduced = reshape(sum(h; dims=length(M) .+ (1:length(K))), M..., B...)
    E_dense_w = interaction_energy(parts(rbm).dense, v_dense, h_reduced)
    E_conv_w = interaction_energy(parts(rbm).conv, v_conv, h)
    @test @inferred(interaction_energy(rbm, v_dense, v_conv, h)) ≈ E_conv_w + E_dense_w

    E_conv = energy(parts(rbm).conv, v_conv, h)
    E_dense = energy(parts(rbm).dense, v_dense, h_reduced)
    E_h = energy(hidden(rbm), h)
    E_h = reshape_maybe(sum(E_h; dims=1:length(J)), B)
    @test @inferred(energy(rbm, v_dense, v_conv, h)) ≈ E_conv .+ E_dense .- E_h

    @test size(@inferred free_energy(rbm, v_dense, v_conv)) == B
end

@testset "pooling DenseConvRBM" begin
    rbm = ConvRBMs.DenseConvRBM(
        Binary(randn(2)), # visible dense
        Binary(randn(3)), # visible conv
        Binary(randn(2)), # hidden
        randn(2, 2), # weights dense
        randn(3, 3, 2); # weights conv
        pool = true
    )

    @test @inferred(free_energy(rbm, bitrand(2), bitrand(3,4))) isa Number
    @test size(@inferred free_energy(rbm, bitrand(2,5), bitrand(3,4,5))) == (5,)
    @test size(@inferred sample_h_from_v(rbm, bitrand(2), bitrand(3,4))) == (2,2)

    h1 = @inferred sample_h_from_v(rbm, bitrand(2), bitrand(3,4))
    h2 = @inferred sample_h_from_h(rbm, h1)
    @test size(h1) == size(h2) == (2,2)
    # pooling: at most one hidden unit active per group
    @test all(sum(h1; dims=2) .≤ 1)
    @test all(sum(h2; dims=2) .≤ 1)
end

@testset "free energy (not pooled)" begin
    rbm = ConvRBMs.DenseConvRBM(
        Binary(randn(2)), # visible dense
        Binary(randn(3)), # visible conv
        Binary(randn(2)), # hidden
        randn(2, 2), # weights dense
        randn(3, 3, 2); # weights conv
        pool = false
    )
    β = rand()
    v_dense, v_conv = bitrand(2), bitrand(3,4)
    hs = [[h1 h2; h3 h4] for h1=0:1, h2=0:1, h3=0:1, h4=0:1]
    @test -β * free_energy(rbm, v_dense, v_conv; β) ≈ logsumexp(-β * energy(rbm, v_dense, v_conv, h) for h in hs)
end

@testset "free energy (pooled)" begin
    rbm = ConvRBMs.DenseConvRBM(
        Binary(randn(2)), # visible dense
        Binary(randn(3)), # visible conv
        Binary(randn(1)), # hidden
        randn(2, 1), # weights dense
        randn(3, 3, 1); # weights conv
        pool = true
    )
    β = rand()
    v_dense, v_conv = bitrand(2), bitrand(3,4)
    hs = [
        [0 1],
        [1 0],
    ]
    K = prod(output_size(rbm, v_conv))
    M = prod(size(hidden(rbm)))
    @test exp(-β * free_energy(rbm, v_dense, v_conv; β)) ≈ (
        sum(exp(-β * energy(rbm, v_dense, v_conv, h)) for h in hs) +
        K^M * exp(-β * energy(rbm, v_dense, v_conv, [0 0]))
    )
end


text:
function fourth(∇f, x)
    fourth(nothing, ∇f, x)
end
function fourth(∇²f, ∇f, x)
    if !(∇²f == nothing)
        ∇²f = 12x^2
    end
    if !(∇f == nothing)
        ∇f = 4x^3
    end

    fx = x^4
    if ∇f == nothing && ∇²f == nothing
        return fx
    elseif ∇²f == nothing
        return fx, ∇f
    else
        return fx, ∇f, ∇²f
    end
end

minimize(fourth, 4.0, BFGS(DirectApprox()))


text:
@inline function reduce{FSA <: Union{FixedArray, Tuple}}(f, a::FSA)
    length(a) == 1 && return a[1]
    @inbounds begin
        red = f(a[1], a[2])
        for i=3:length(a)
            red = f(red, a[i])
        end
    end
    red
end

@inline function reduce(f::Base.Func{2}, a::Mat)
    length(a) == 1 && return a[1,1]
    @inbounds begin
        red = reduce(f, a.(1)[1])
        for i=2:size(a, 2)
            red = f(red, reduce(f, a.(1)[i]))
        end
    end
    red
end

index_expr{T <: Number}(::Type{T}, i::Int, inds::Int...) = :($(symbol("arg$i")))
index_expr{T <: FixedArray}(::Type{T}, i::Int, inds::Int...) = :($(symbol("arg$i"))[$(inds...)])
inner_expr{N}(args::NTuple{N, DataType}, inds::Int...) = :( F($(ntuple(i -> index_expr(args[i], i, inds...), N)...)) )

# This solves the combinational explosion from FixedVectorNoTuple while staying fast.
function constructor_expr{T <: FixedVector}(::Type{T}, tuple_expr::Expr)
    quote
        $(Expr(:boundscheck, false))
        $(Expr(:meta, :inline))
        rvalue = FSA($(tuple_expr))
        $(Expr(:boundscheck,:pop))
        rvalue
    end
end
constructor_expr{T <: Mat}(::Type{T}, tuple_expr::Expr) = quote
    $(Expr(:boundscheck, false))
    $(Expr(:meta, :inline))
    rvalue = Mat($(tuple_expr))
    $(Expr(:boundscheck,:pop))
    rvalue
end
constructor_expr{T <: FixedVectorNoTuple}(::Type{T}, tuple_expr::Expr) = quote
    $(Expr(:boundscheck, false))
    $(Expr(:meta, :inline))
    rvalue = FSA($(tuple_expr)...)
    $(Expr(:boundscheck,:pop))
    rvalue
end
@generated function map{FSA <: FixedArray}(F, arg1::FSA, arg2::FSA)
    inner = fill_tuples_expr((inds...) -> inner_expr((arg1, arg2), inds...), size(FSA))
    constructor_expr(FSA, inner)
end
@generated function map{FSA <: FixedArray}(F, arg1::FSA, arg2::Number)
    inner = fill_tuples_expr((inds...) -> inner_expr((arg1, arg2), inds...), size(FSA))
    constructor_expr(FSA, inner)
end
@generated function map{FSA <: FixedArray}(F, arg1::Number, arg2::FSA)
    inner = fill_tuples_expr((inds...) -> inner_expr((arg1, arg2), inds...), size(FSA))
    constructor_expr(FSA, inner)
end
@generated function map{FSA <: FixedArray}(F, ::Type{FSA})
    inner = fill_tuples_expr((inds...) -> :(F($(inds...))), size(FSA))
    :( FSA($inner) )
end

@generated function map{FSA <: FixedArray}(F, arg1::FSA)
    inner = fill_tuples_expr((inds...) -> :( F(arg1[$(inds...)]) ), size(FSA))
    constructor_expr(FSA, inner)
end

@generated function map{T}(::Type{T}, arg1::FixedArray)
    eltype(arg1) == T && return :(arg1)
    FSA = similar(arg1, T)
    inner = fill_tuples_expr((inds...) -> :( T(arg1[$(inds...)]) ), size(FSA))
    :( $FSA($(inner)) )
end

@inline map{R,C,T}(F::Type{T}, arg1::Mat{R,C,T}) = arg1
@generated function map{R,C,T}(F::DataType, arg1::Mat{R,C,T})
    inner = fill_tuples_expr((inds...) -> :( F(arg1[$(inds...)]) ), (R, C))
    :( Mat{R, C, F}($(inner)) )
end

@generated function map{FSA <: FixedVectorNoTuple}(F::DataType, arg1::FSA)
    eltype(FSA) == F && return :(arg1)
    inner = ntuple(i-> :(F(arg1[$i])), length(FSA))
    :( similar(FSA, F)($(inner...)) )
end


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

struct TenantsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets the tenants for your account.
Param: api_version::String (required)
Return: TenantListResult
"""
function tenantsList(_api::TenantsApi, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", TenantListResult, "/tenants", ["azure_auth"])
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export tenantsList


text:
# moderate modification and adaption to
# https://github.com/JuliaPy/pyjulia/blob/87c669e2729f9743fe2ab39320ec9b91c9300a96/src/julia/juliainfo.jl
export PyJulia_INFO
import Libdl
import Pkg

function PyJulia_INFO()
    println(VERSION)
    println(VERSION.major)
    println(VERSION.minor)
    println(VERSION.patch)
    # binddir
    println(Base.Sys.BINDIR)
    # libjulia_path
    println(Libdl.dlpath(string("lib", splitext(Base.julia_exename())[1])))
    # sysimage
    println(unsafe_string(Base.JLOptions().image_file))
end


text:
# scripts for converting R-tree into `DataFrame`
using DataFrames

# appends node/element information into the `dfcols`
function append_node_info!(dfcols::NamedTuple, node::Any,
                           parent::Union{SI.Node, Missing},
                           ix::Union{Integer, Missing})
    push!(dfcols.type, string(typeof(node).name.name))
    push!(dfcols.level, isa(node, SI.Node) ? SI.level(node) : missing)
    push!(dfcols.ix, ix)
    push!(dfcols.id, isa(node, SI.Node) ? string(Int(pointer_from_objref(node)), base=16) :
         (SI.idtrait(typeof(node)) !== SI.HasNoID ? repr(SI.id(node)) : missing))
    push!(dfcols.parent, !ismissing(parent) ? string(Int(pointer_from_objref(parent)), base=16) : missing)
    push!(dfcols.area, SI.area(SI.mbr(node)))
    push!(dfcols.xmin, SI.mbr(node).low[1])
    push!(dfcols.xmax, SI.mbr(node).high[1])
    push!(dfcols.ymin, SI.mbr(node).low[2])
    push!(dfcols.ymax, SI.mbr(node).high[2])
    return dfcols
end

# appends information about all the nodes in the `node` subtree
# into the `dfcols`
function append_subtree_info!(dfcols::NamedTuple, node::Any,
                              parent::Union{SI.Node, Missing},
                              ix::Union{Integer, Missing})
    append_node_info!(dfcols, node, parent, ix)
    if node isa SI.Node
        for (i, child) in enumerate(SI.children(node))
            append_subtree_info!(dfcols, child, node, i)
        end
    end
    return node
end

# convert `RTree` into `DataFrame`
function Base.convert(::Type{DataFrame}, tree::RTree)
    dfcols = (type = Vector{String}(), level = Vector{Union{Int, Missing}}(),
              id = Vector{Union{String, Missing}}(),
              ix = Vector{Union{Int, Missing}}(),
              parent = Vector{Union{String, Missing}}(),
              area =Vector{Float64}(),
              xmin = Vector{Float64}(), ymin = Vector{Float64}(),
              xmax = Vector{Float64}(), ymax = Vector{Float64}())
    append_subtree_info!(dfcols, tree.root, missing, missing)
    return DataFrame(dfcols)
end


text:
# using Plots
using Omega
# using Parameters
# include(joinpath(dirname(pathof(Omega)), "viz.jl"))

function prob1()
  x = uniform(-1.0, 1.0)
  y = uniform(-1.0, 1.0)
  (x = x, y = y, c = x ==ₛ y, xlims = (-1, 1), ylims = (-1, 1))
end

function prob2()
  x = uniform(-1.0, 1.0)
  y = uniform(-1.0, 1.0)
  (x = x, y = y, c = x >ₛ y, xlims = (-1, 1), ylims = (-1, 1))
end

function prob3()
  x = uniform(-1.0, 1.0)
  y = uniform(-1.0, 1.0)
  (x = x, y = y, c = abs(x) >ₛ abs(y), xlims = (-1, 1), ylims = (-1, 1))
end

function prob4()
  x = uniform(-1.0, 1.0)
  y = uniform(-1.0, 1.0)
  (x = x, y = y, c = x^2 ==ₛ y^2, xlims = (-1, 1), ylims = (-1, 1))
end

function prob5(k = 3, thresh = 0.9999)
  x = uniform(0.0, 1.0)
  y = uniform(0.0, 1.0)
  c = (sin(2π * x * k) * cos(2π * y * k)) >ₛ thresh
  (x = x, y = y, c = c, xlims = (0, 1), ylims = (0, 1))
end

function prob6(k = 3, thresh = 0.8)
  x = uniform(0, 1.0)
  y = uniform(0, 1.0)
  c = (sin(2π * x * k) * cos(2π * y * k)) >ₛ thresh
  (x = x, y = y, c = c, xlims = (0, 1), ylims = (0, 1))
end

function prob7(k = 3, thresh = 0.8)
  x = uniform(-1.0, 1.0)
  y = uniform(-100.0, 100.0)
  c = (x ==ₛ 0.0) & (y ==ₛ 0.0)
  (x = x, y = y, c = c, xlims = (-1, 1), ylims = (-100, 100))
end

# function allsamples(prob, n; algkwargs...)
#   # ωsamples = rand(defΩ(), logerr(prob.c), n, alg, kwargs...)
#   # xy = randtuple((prob.x, prob.yy))
#   # samples = map(ω -> applynotrackerr(x, ω), ωsamples)

#   # ωsamples = rand((prob.x, prob.y), prob.c, n; alg = alg, kwargs...)
#   # (ωsamples = ωsamples, samples = samples)
#   samples = rand((prob.x, prob.y), prob.c, n; algkwargs...)
#   (samples = samples, ωsamples = nothing)
# end

# function scatterxy(samples;
#                    label = nothing, legend = nothing, xlims = (-1, 1),
#                    ylims = (-1, 1), kwargs...)
#   xs, ys = ntranspose(samples)
#   scatter(xs, ys, label = label, legend = legend, xlims = xlims, ylims = ylims)
# end

# function vizall(probs, algs, n)
#   plots = []
#   for prob in probs
#     @unpack x, y, c, xlims, ylims = prob
#     push!(plots, ωcontour(err(c); label = nothing, legend = nothing,
#                                   colorbar = nothing))
#     for alg in algs
#       try
#         @unpack samples, ωsamples = allsamples(prob, n; alg...)
#         push!(plots, scatterxy(samples, xlims = xlims, ylims = ylims))
#       catch e
#         println("Failed")
#         display(e)
#         # rethrow(e)
#         push!(plots, plot())
#       end
#     end
#   end
#   plots
#   # plt = plot(plots..., layout = (length(probs), length(algs) + 1))
# end

# probs = [prob1(), prob2(), prob3(), prob4(), prob5(), prob6()]
# algs = [(alg = SSMH,),
#         (alg = NUTS,),
#         (alg = Replica, nreplicas = 4, inneralg = SSMH)
#         (alg = Replica, nreplicas = 4, inneralg = NUTS)
#         ]

# plots = vizall(probs, algs, 1000)

# function makeplots()
#   plt = plot(plots..., layout = (length(probs), length(algs) + 1),
#              markersize=0.01, tickfontsize=4)
# end

# # n = 8
# # temps = Omega.Inference.logtemps(n)
# # temps = [1e-9, 10000]
# # withkernel(kseα(temps[end])) do
# #   @unpack x, y, c = prob6()
# #   # p1 = ωcontour(err(c); label = nothing, legend = nothing);
# #   p2 = ωcontour(logerr(c); label = nothing, legend = nothing)
# #   # samples = rand((x, y), c, 100000; alg = Replica, nreplicas = 8, temps = temps)
# #   # p3 = scatterxy(samples, xlims = (0, 1), ylims = (0, 1))
# #   # plot(p1, p3, markersize = 0.01, aspectratio = 1)
# # end

x = normal(0, 1)
y = normal(0, 1)
# rand((x, y), x ==ₛ y, alg = Replica, inneralg = HMCFAST)
rand((x, y), x ==ₛ y, 1000, alg = HMCFAST)