text:
"""
Create an expression a method for `Static` input to `fun` if the result
is among the `targets`
"""
function genstaticmethod1(fun, arg1, targets)
    r = try
        @eval $fun($arg1)
    catch
        return nothing
    end
    ts = filter(isequal(r), targets)
    if length(ts) == 1
        return :( $fun(::$(typeof(static(arg1)))) = $(static(first(ts))) )
    end
    return nothing
end

"""
Create an expression a method for two `Static` inputs to `fun` if the result
is among the `targets`
"""
function genstaticmethod2(fun, arg1, arg2, targets)
    r = try
        @eval $fun($arg1, $arg2)
    catch
        return nothing
    end
    ts = filter(isequal(r), targets)
    if length(ts) == 1
        return :( $fun(::$(typeof(static(arg1))),::$(typeof(static(arg2)))) = $(static(first(ts))) )
    end
    return nothing
end

"""
Create an array of all the methods for `Static` 1-argument inputs to `fun`
where the result is either among the possible inputs, `args`, or among
the additional `targets`.
"""
function genstaticmethods1(funs, args, targets)
    m = Vector{Expr}()
    tu = union(args, targets)
    for f in funs
        for a in args
            e = genstaticmethod1(f, a, tu)
            if e isa Expr
                push!(m, e)
            end
        end
    end
    return m
end

"""
Create an array of all the methods for `Static` 2-argument inputs to `fun`
where the result is either among the possible inputs, `args`, or among
the additional `targets`.
"""
function genstaticmethods2(funs, args, targets)
    m = Vector{Expr}()
    tu = union(args, targets)
    for f in funs
        for a1 in args
            for a2 in args
                e = genstaticmethod2(f, a1, a2, tu)
                if e isa Expr
                    push!(m, e)
                end
            end
        end
    end
    return m
end

"""
@staticnumbers numbers 1argfuns 2argfuns

This macro creates methods that return `Static` numbers when
functions are called with only `Static` arguments.

The inputs should be a list of literal numbers that will be tested as
inputs and outputs of all functions.

Optionally a fourth argument can give a list of numbers that will only
be considered as results, but not as inputs.

Example:
```
@staticnumbers (0, 1) (sin, cos) (+, -)
```
will create all the following method definitions:
```
sin(::StaticInteger{0}) = StaticInteger{0}()
cos(::StaticInteger{0}) = StaticInteger{1}()
+(::StaticInteger{0}, ::StaticInteger{0}) = StaticInteger{0}()
+(::StaticInteger{0}, ::StaticInteger{1}) = StaticInteger{1}()
+(::StaticInteger{1}, ::StaticInteger{0}) = StaticInteger{1}()
-(::StaticInteger{0}, ::StaticInteger{0}) = StaticInteger{0}()
-(::StaticInteger{1}, ::StaticInteger{0}) = StaticInteger{0}()
-(::StaticInteger{1}, ::StaticInteger{1}) = StaticInteger{1}()
```
(Note: The macro will run in the local scope. Functions from `Base`
must be imported before they can be extended.)
"""
macro staticnumbers(args::Expr, funs1::Expr, funs2::Expr, targets::Expr=:(()))
    args.head == :tuple || error("Expected a Tuple of numbers")
    funs1.head == :tuple || error("Expected a Tuple of 1-arg functions")
    funs2.head == :tuple || error("Expected a Tuple of 2-arg functions")
    targets.head == :tuple || error("Expected a Tuple of target numbers")
    return esc(Expr(:block,vcat(
       genstaticmethods1(funs1.args, args.args, targets.args),
       genstaticmethods2(funs2.args, args.args, targets.args)
       )...))
end


text:
using VisualCircuits, Compose

for logic in [:NOT, :AND, :NAND, :OR, :NOR, :XOR, :XNOR]
    g = Gate(1, 1, logic)
    c = compose(
        context(units = UnitBox(0.25, 0.25, 1.6, 1.5)),
        render(g),
        # Draw input & output lines (connected to x = 0; x = 2)
        connect(g[:in1], 0),
        connect(g[:in2], 0),
        connect(g[:out], 2),
        stroke("black"), linewidth(1mm)
    )
    draw(PDF(string(logic) * ".pdf", 40mm, 40mm), c)
end

ext: jl


text:
@testset "Parsimonious flux balance analysis with StandardModel" begin
    model = test_toyModel()

    d = parsimonious_flux_balance_analysis_dict(
        model,
        Tulip.Optimizer;
        modifications = [
            change_constraint("EX_m1(e)", lb = -10.0),
            change_optimizer_attribute("IPM_IterationsLimit", 500),
        ],
        qp_modifications = [
            change_optimizer(OSQP.Optimizer),
            change_optimizer_attribute("polish", true),
            silence,
        ],
    )

    # The used optimizer doesn't really converge to the same answer everytime
    # here, we therefore tolerate a wide range of results.
    @test isapprox(d["biomass1"], 10.0, atol = QP_TEST_TOLERANCE)
end


text:
n = 10
using Test

pairwise_funcs = [
    StatsBase.cov,
    StatsBase.cor,
    StatsBase.countne,
    StatsBase.counteq,
    StatsBase.corkendall,
    StatsBase.corspearman,
    StatsBase.maxad,
    StatsBase.meanad,
    StatsBase.msd,
    StatsBase.psnr,
    StatsBase.rmsd,
    StatsBase.sqL2dist,
    StatsBase.crosscor,
    StatsBase.crosscov
]

uvals = example_uvals

@testset "Pairwise statistics on datasets" begin
    @testset "$(pairwise_funcs[i])" for i = 1:length(pairwise_funcs)
        f = pairwise_funcs[i]
        if f == StatsBase.psnr
            maxv = 100
            @test f(uvals, uvals, maxv, n) isa Vector{T} where T <: Real

        elseif f ∈ [StatsBase.crosscor, StatsBase.crosscov]
            @test f(uvals, uvals, n) isa Vector{Vector{T}} where T <: Real
            @test f(uvals, uvals, 1:5, n) isa Vector{Vector{T}} where T <: Real
        else
            @test f(uvals, uvals, n) isa Vector{T} where T <: Real
        end
    end
end

UV = UncertainValueDataset(example_uvals)

@testset "Pairwise statistics on datasets" begin
    @testset "$(pairwise_funcs[i])" for i = 1:length(pairwise_funcs)
        f = pairwise_funcs[i]
        if f == StatsBase.psnr
            maxv = 100
            @test f(UV, UV, maxv, n) isa Vector{T} where T <: Real

        elseif f ∈ [StatsBase.crosscor, StatsBase.crosscov]
            @test f(UV, UV, n) isa Vector{Vector{T}} where T <: Real
            @test f(UV, UV, 1:5, n) isa Vector{Vector{T}} where T <: Real
        else
            @test f(UV, UV, n) isa Vector{T} where T <: Real
        end
    end
end

UI = UncertainIndexDataset(example_uvals)

@testset "Pairwise statistics on datasets" begin
    @testset "$(pairwise_funcs[i])" for i = 1:length(pairwise_funcs)
        f = pairwise_funcs[i]
        if f == StatsBase.psnr
            maxv = 100
            @test f(UI, UI, maxv, n) isa Vector{T} where T <: Real

        elseif f ∈ [StatsBase.crosscor, StatsBase.crosscov]
            @test f(UI, UI, n) isa Vector{Vector{T}} where T <: Real
            @test f(UI, UI, 1:5, n) isa Vector{Vector{T}} where T <: Real
        else
            @test f(UI, UI, n) isa Vector{T} where T <: Real
        end
    end
end

# Functions that under the hood use functions with strictly positive domains
# special_pairwise_funcs = [
#     StatsBase.gkldiv,
#     StatsBase.kldivergence,
# ]

# @testset "$(pairwise_funcs[i])" for i = 1:length(pairwise_funcs)
#     f = pairwise_funcs[i]
#     @testset for (i, uval) in enumerate(example_uvals)
#         @test f(uval, uval, n) isa T where T <: Real
#     end
# end;


text:
using Optim, Test
@testset "#711" begin
	# make sure it doesn't try to promote df
	dof = 7

	fun(x) = 0.0; x0 = fill(0.1, dof)
	df = TwiceDifferentiable(fun, x0)

	lx = fill(-1.2, dof); ux = fill(+1.2, dof)
	dfc = TwiceDifferentiableConstraints(lx, ux)

	res = optimize(df, dfc, x0, IPNewton(); autodiff=:forward)
	res = optimize(df, dfc, x0, IPNewton())
end

@testset "#600" begin
	function exponential(x)
		return exp((2.0 - x[1])^2) + exp((3.0 - x[2])^2)
	end

	function exponential_gradient!(storage, x)
		storage[1] = -2.0 * (2.0 - x[1]) * exp((2.0 - x[1])^2)
		storage[2] = -2.0 * (3.0 - x[2]) * exp((3.0 - x[2])^2)
		storage
	end
	function exponential_hessian!(storage, x)
		Optim.NLSolversBase.ForwardDiff.hessian!(storage, exponential, x)
	end

	function exponential_gradient(x)
		storage = similar(x)
		storage[1] = -2.0 * (2.0 - x[1]) * exp((2.0 - x[1])^2)
		storage[2] = -2.0 * (3.0 - x[2]) * exp((3.0 - x[2])^2)
		storage
	end

	initial_x = [0.0, 0.0]
	optimize(exponential, exponential_gradient!, initial_x, BFGS())
	lb = fill(-0.1, 2)
	ub = fill(1.1, 2)
	od = OnceDifferentiable(exponential, initial_x)
	optimize(od, lb, ub, initial_x, IPNewton())
	optimize(od, lb, ub, initial_x, IPNewton(), Optim.Options())
	optimize(exponential, lb, ub, initial_x, IPNewton())
	optimize(exponential, lb, ub, initial_x, IPNewton(), Optim.Options())
	optimize(exponential, exponential_gradient!, lb, ub, initial_x, IPNewton())
	optimize(exponential, exponential_gradient!, lb, ub, initial_x, IPNewton(), Optim.Options())
	optimize(exponential, exponential_gradient!, exponential_hessian!, lb, ub, initial_x, IPNewton())
	optimize(exponential, exponential_gradient!, exponential_hessian!, lb, ub, initial_x, IPNewton(), Optim.Options())
	optimize(TwiceDifferentiable(od, initial_x), lb, ub, initial_x)
	optimize(TwiceDifferentiable(od, initial_x), lb, ub, initial_x, Optim.Options())
end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    varshape(f::Function)::AbstractValueShape

Get the value shape of the input/argument of an unary function `f`. `f`
should support call syntax

    f(x::T)

with `valshape(x) == varshape(f)` as well all

    f(x::AbstractVector{<:Real})

with `length(eachindex(x)) == vardof(f)` (see [`vardof`](@ref)).
"""
function varshape end
export varshape

"""
    vardof(f::Function)::Integer

Get the number of degrees of freedom of the input/argument of f.

Equivalent to `totalndof(varshape(f))` (see [`varshape`](@ref)).
"""
function vardof end
export vardof

vardof(f::Function) = totalndof(varshape(f))

struct FuncWithVarShape{F<:Function,VS<:AbstractValueShape} <: Function
    f::F
    varshape::VS
end

varshape(fws::FuncWithVarShape) = fws.varshape

Base.@propagate_inbounds (fws::FuncWithVarShape)(x::Any) = fws.f(x)

Base.@propagate_inbounds function (fws::FuncWithVarShape)(x::AbstractVector{<:Real})
    # ToDo: Check performance:
    x_shaped = stripscalar(varshape(fws)(x))
    fws.f(x_shaped)
end

import Base.>>
@inline >>(varshape::AbstractValueShape, f::Function) = FuncWithVarShape(f, varshape)


text:
# Adapted from https://raw.githubusercontent.com/JuliaData/Tables.jl/main/find_integrations.jl
####
#### automatically generate a list of integrations
####

###
### Usage
###
### 1. ensure a development version of GeoInterface.jl (`pkg> add GeoInterface`)
### 2. make sure the General registry is up to date (`pkg> up`)
### 3. run this script, which uses the first depot from DEPOT_PATH

DEPOT = first(DEPOT_PATH)
REGISTRIES = joinpath(DEPOT, "registries")
@info DEPOT
# find each package w/ a direct dependency on GeoInterface.jl
general = joinpath(DEPOT, "General")
mkpath(general)
# run(`tar -xzf General.tar.gz -C $general`)
pkgs = cd(REGISTRIES) do
    dirname.(readlines(`grep -rl ^GeoInterface General/. --include=Deps.toml`))
end

pkgnames = sort([splitpath(x)[end] for x in pkgs])

function parseurl(file)
    repo = readlines(file)[end]
    return strip(split(repo, " = ")[end], '"')
end

urls = [parseurl(joinpath(REGISTRIES, "General", string(first(nm)), nm, "Package.toml"))
        for nm in pkgnames]

open(joinpath(dirname(@__DIR__), "GeoInterface.jl", "INTEGRATIONS.md"), "w+") do io
    println(io, "# Packages\nPackages currently integrating with GeoInterface.jl:")
    for (nm, url) in zip(pkgnames, urls)
        println(io, "* [$nm]($url)")
    end
end


text:
@nograd readline, Base.gc_num, Base.time_ns, Base.print, Base.println, Base.show,
  Core.show, Core.print, Core.println, string, repr, Threads.nthreads, Threads.threadid

# Gradient of AD stacks

grad_mut(::AbstractVector) = []

@adjoint! function _push!(a::Vector, x)
  _push!(a, x), function (y)
    dstk = grad_mut(__context__, a)
    return (nothing, pop!(dstk))
  end
end

@adjoint! function pop!(stk::Stack)
  pop!(stk), function (Δ)
    dstk = grad_mut(__context__, stk.data)
    push!(dstk, Δ)
    return
  end
end

# Dictionaries

grad_mut(d::AbstractDict) = Dict()
grad_mut(d::IdDict) = IdDict()

# TODO perhaps look up mutable gradients in `pullback`
function accum(a::AbstractDict, b::AbstractDict)
  @assert a === b
  return a
end

@adjoint function getindex(d::AbstractDict, k)
  d[k], function (Δ)
    grad = grad_mut(__context__, d)
    grad[k] = accum(get(grad, k, nothing), Δ)
    return (grad, nothing)
  end
end

@adjoint! function setindex!(d::AbstractDict, v, k)
  setindex!(d, v, k), function (_)
    Δ = get(grad_mut(__context__, d), k, nothing)
    delete!(grad_mut(__context__, d), k)
    (nothing, Δ, nothing)
  end
end

@nograd haskey

# Channels

@nograd Channel, schedule

grad_mut(ch::Channel) = Channel(ch.sz_max)

@adjoint! function put!(ch::Channel, x)
  put!(ch, x), function (ȳ)
    x̄ = take!(grad_mut(__context__, ch))
    (nothing, accum(x̄, ȳ), nothing)
  end
end

@adjoint! function take!(ch::Channel)
  take!(ch), function (x̄)
    put!(grad_mut(__context__, ch), x̄)
    return
  end
end

@adjoint! function Task(f)
  t = Task(f)
  t.code = function ()
    y, back = _pullback(__context__, f)
    cache(__context__)[t] = Task(back)
    return y
  end
  t, _ -> fetch(cache(__context__)[t])
end

function runadjoint(cx, t, ȳ = nothing)
  t̄ = cache(cx)[t]
  f = t̄.code
  t̄.code = () -> f(ȳ)
  @static if VERSION > v"1.3-"
    t̄.sticky = t.sticky
  end
  schedule(t̄)
end

@adjoint! function wait(t::Task)
  wait(t), _ -> (runadjoint(__context__, t); nothing)
end

@adjoint! function fetch(t::Task)
  fetch(t), ȳ -> (runadjoint(__context__, t, ȳ); nothing)
end

@adjoint! function Base.sync_end(refs)
  Base.sync_end(refs), _ -> foreach(t -> runadjoint(__context__, t), refs)
end

# Make @sync work
# Safe as long as other kinds of mutation are disallowed
@adjoint push!(refs::Vector{Any}, t::Task) = push!(refs, t), _ -> nothing


text:
using Revise, MDToolbox, Plots, JLD2, Printf, DelimitedFiles, BenchmarkTools

include("../src/mathTool.jl")
include("../src/plotTool.jl")
include("../src/BioAFMResult.jl")
include("../src/BioAFM_Q.jl")
include("../src/BioAFM_QD.jl")
include("../src/BioAFM_QD_FFT.jl")
include("../src/BioAFMize_QD_FFT.jl")
include("../src/afm/afmize.jl")

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
function jacobi_preconditioned_gmres(A::SparseMatrixCSC, x0::Vector{Float64}, b::Vector{Float64}, Nz, pchunks; w=2/3, Na=10, Ns=100, m=20, err_thresh=1e-8)

    Nk = div(length(x0), Nz)

    J = JacobiSmoother(A, b::Vector{Float64}, w, pchunks, Nk, Nz, Na)
    x = x0[:]
    x[:] = J*x
    for t=1:Ns
        x[:] = J * x
    end

    err = norm(A*x - b) / norm(b)

    if err > err_thresh
        zslices = collect(Iterators.partition(1:Nk*Nz, Nz))
        c = pmap(zslices) do inds
            Aloc = A[inds, inds]
            bloc = b[inds]
            gmres_solve(Aloc, x[inds], bloc, m; err_thresh=err_thresh)
        end
        return vcat(c...)
    else
        return x
    end

end

#
#
# struct LinearLeftHandSide
#
#     A :: SparseMatrixCSC
#     M :: SparseMatrixCSC
#     D :: SparseMatrixCSC
#
# end
# LinearLeftHandSide(A::SparseMatrixCSC, w::Float64) = LinearLeftHandSide(A, compute_jacobi_matrices(A, w)...)
#
#
#
# function jacobi_preconditioned_gmres(L::LinearLeftHandSide, x0::Vector{Float64}, b::Vector{Float64}; n=100, m=10, err_thresh=1e-10)
#
#     # precodition using jacobi_relaxation()
#     x = jacobi_relaxation(n, L.A, L.M, L.D, x0, b)
#
#     # finish with gmres
#     x = gmres_solve(L.A, x, b, m; err_thresh=err_thresh)
#
#     return x
# end
#
#
# function parallel_jacobi_preconditioned_gmres(L::LinearLeftHandSide, x0::Vector{Float64}, b::Vector{Float64}, pchunks::Int64, Nz::Int64; nsynch=100, nasynch=5, m=10, err_thresh=1e-10)
#
#     # jacobi preconditioning
#     x = distribute(x0; dist=pchunks*Nz)
#     bdist = distribute(b; dist=pchunks*Nz)
#     for _=1:nsynch
#         x = asynch_jacobi_relaxation(nasynch, L.M, L.D, x, bdist)
#     end
#
#     # do gmres in parallel
#     x = DArray((size(x)[1],), workers()[1:Nz], Nz) do inds
#
#         xs = Array(x0[inds[1],1])
#         bs = b[inds[1]]
#         gmres_solve(L.A, xs, bs, m, err_thresh=err_thresh)
#     end
#
#     return Array(x)
# end


text:
__precompile__()

module PhyloTrees

  # Methods expanded
  import
    Base.show,
    Base.push!,
    Base.append!,
    Base.getindex

  # Functions provided
  export
    # Trees
    Tree,

    # Utilities
    addnode!,
    addnodes!,
    addbranch!,
    branch!,
    setsource!,
    settarget!,
    setlength!,
    setheight!,
    changesource!,
    changetarget!,
    changelength!,
    changeheight!,
    deletenode!,
    deletebranch!,
    indegree,
    outdegree,
    isroot,
    isleaf,
    isinternal,
    findroots,
    findroot,
    findleaves,
    findinternals,
    findexternals,
    findnonroots,
    findnonleaves,
    areconnected,
    nodepath,
    branchpath,
    parentnode,
    childnodes,
    descendantcount,
    descendantnodes,
    ancestorcount,
    ancestornodes,
    nodetype,

    # Distance
    distance,

    # Traversal
    postorder

  # Package files
  include("structure.jl")
  include("show.jl")
  include("construction.jl")
  include("utilities.jl")
  include("distance.jl")
  include("traversal.jl")
end


text:
function getUMD()
    umdUrl = "http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Momentum_Factor_daily_TXT.zip"
    (data, dataNames, varnames) = readFamaFrenchRaw(umdUrl)

    newVarnames = ["UMD"]
    finalData = TimeSeries.rename(data[1], newVarnames)

    fileName = joinpath(Pkg.dir("EconDatasets"), "data", "UMD.csv")
    writetimearray(finalData, fileName)

end


text:
include("../common.jl")
build_harfbuzz(ARGS, "HarfBuzz_ICU")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function gencenters(side, λ)
    n = pois_rand(λ * side^2)
    return reduce(hcat, (rand(1:side, 2) for i in 1:n))
end

function gendisks(side, R, λ)
    spheres = zeros(Int8, (side + 2R + 1, side + 2R + 1))
    sphere  = zeros(Int8, (2R + 1, 2R + 1))
    centers = gencenters(side, λ)
    for i in -R:R
        for j in -R:R
            dist = i^2 + j^2
            if dist < R^2
                sphere[j+R+1, i+R+1] = 1
            end
        end
    end
    for center in (centers[:,i] for i in 1:size(centers,2))
        x = center[1]
        y = center[2]
        spheres[x:x + 2R, y:y + 2R] .|= sphere
    end
    return spheres[R+1:end-R-1, R+1:end-R-1]
end

function s2_theory(r, R, λ)
    tmp = r/(2R)
    tmp2 = (r > 2R) ? 2 : 2/π*(π + tmp*sqrt(1 - tmp^2) - acos(tmp))
    η = λ * π * R^2
    return exp(-η*abs(tmp2))
end

function l2_theory(r, R, λ)
    η = λ * π * R^2
    p = exp(-η)
    return p^(1 + 2r/(π*R))
end

function ss_theory(r, R, λ)
    if r < 2R
        A = 4R^2 - r^2
        B = acos(r/(2R))
        part = exp(-λ*(r*sqrt(A)/2 + 2π*R^2 - 2*B*R^2))
        return part * ((4A*B^2 - 8π*A*B + 4*π^2*A)*R^2*λ^2*r + 4*sqrt(A)*R^2*λ) / (A*r)
    else
        return (2π*λ*R)^2*exp(-2π*λ*R^2)
    end
end

function sv_theory(r, R, λ)
     if r < 2R
        A = 4R^2 - r^2
        B = acos(r/(2R))
        return 2R*λ*(π-B)exp(-sqrt(A)*r*λ/2 + 2B*R^2*λ - 2π*R^2*λ)
    else
        return 2π*R*λ*exp(-2π*λ*R^2)
    end
end

function pore_size_theory(r, R, λ)
    s(x) = 2π * x
    v(x) = π * x^2
    η  = v(R) * λ
    p  = exp(-η)
    s1 = s(r + R)
    v1 = v(r + R)
    return λ * s1 * exp(-λ * v1) / p
end

function chord_length_theory(r, R, λ)
    η = λ * π * R^2
    p = exp(-η)
    return 2η/(π*R)*p^(2r/(π*R))
end

mean_chord_length(R, λ) = 1/(2λ*R)

@testset "S2 on random overlapping disks generated by Poisson process" begin
    # Area = S^2, radius of a disk = R
    # Poisson parameter = λ
    S = 7000; R = 40; λ = 5e-5;

    disks = gendisks(S, R, λ)
    calc = Directional.s2(disks, 0; periodic = true) |> mean
    theory = [s2_theory(r-1, R, λ) for r in 0:length(calc) - 1]

    err = relerr.(calc, theory)
    @test maximum(err) < 0.15
end

@testset "L2 on random overlapping disks generated by Poisson process" begin
    # Area = S^2, radius of a disk = R
    # Poisson parameter = λ
    S = 7000; R = 40; λ = 5e-5; N = 700

    disks = gendisks(S, R, λ)
    calc = Directional.l2(disks, 0; len = N, periodic = true) |> mean .|> log
    theory = [(log ∘ l2_theory)(r-1, R, λ) for r in 0:N - 1]

    err = relerr.(calc, theory)
    @test maximum(err) < 0.15
end

@testset "SS on random overlapping disks generated by Poisson process" begin
    # Area = S^2, radius of a disk = R
    # Poisson parameter = λ
    S = 7000; R = 30; λ = 4e-4;

    disks = gendisks(S, R, λ)
    calc = Directional.surfsurf(disks, 0; periodic = true) |> mean
    theory = [ss_theory(r, R, λ) for r in 0:length(calc) - 1]

    err = relerr.(calc, theory)
    # Surface-surface function with r < 2R has a discontinuity and two
    # points where it goes to +Inf. Do not test it.
    @test maximum(err[2R + 10:end]) < 0.15
end

@testset "SV on random overlapping disks generated by Poisson process" begin
    # Area = S^2, radius of a disk = R
    # Poisson parameter = λ
    S = 7000; R = 30; λ = 4e-4;

    disks = gendisks(S, R, λ)
    calc = Directional.surfvoid(disks, 0; periodic = true) |> mean
    theory = [sv_theory(r, R, λ) for r in 0:length(calc) - 1]

    err = relerr.(calc, theory)
    @test maximum(err) < 0.2
end

@testset "Pore size on random overlapping disks generated by Poisson process" begin
    # Area = S^2, radius of a disk = R, Poisson parameter = λ
    S = 7000; R = 50; λ = 5e-5
    disks = gendisks(S, R, λ)

    calc = pore_size(disks; nbins = 20)
    edges = calc.edges[1]
    s = step(edges)
    theory = [integrate(x -> pore_size_theory(x, R, λ), n:0.05:n+s)
              for n in 0:s:s*(length(edges) - 2)]

    # Compare cummulative distributions instead of probability
    # densities because density is almost zero for big radii.
    calc_cdf = scan(calc.weights)
    theory_cdf = scan(theory)

    err = relerr.(calc_cdf, theory_cdf)
    @test maximum(err) < 0.1
end

@testset "Chord length on random overlapping disks generated by Poisson process" begin
    # Area = S^2, radius of a disk = R, Poisson parameter = λ
    S = 7000; R = 60; λ = 1e-4
    disks = gendisks(S, R, λ)

    calc, mc = Directional.chord_length(disks, 0; nbins = 35)
    edges = calc.edges[1]
    s = step(edges)
    theory = [integrate(x -> chord_length_theory(x, R, λ), n:0.05:n+s)
              for n in 0:s:s*(length(edges) - 2)]

    # Compare cummulative distributions instead of probability
    # densities because density is almost zero for big lengths.
    calc_cdf = scan(calc.weights)
    theory_cdf = scan(theory)

    err = relerr.(calc_cdf, theory_cdf)
    @test maximum(err) < 0.1
    @test relerr(mc, mean_chord_length(R, λ)) < 0.1
end


text:
using Colors
import Compose: compose, context, polygon
using Lazy
using Gadfly

import Escher: @d

# Home page.

pkgname(name="Escher") =
    hbox(title(3, name), hskip(1em), Escher.tex("\\beta")) |> fontcolor("#999")

*(f::Function, g::Function) = x -> f(g(x))

sidenote(note, icon="info", iconcolor="#aaa") =
    hbox(Escher.icon(icon) |>
            size(100Escher.px, 100Escher.px) |>
            fontcolor(iconcolor),
        hskip(4em),
        note |>
           fontsize(0.9em) * fontcolor("#666")) |>
        pad(1em) |>
        bordercolor("#e1e1e1") |>
        Escher.borderstyle([top, bottom], solid) |>
        Escher.borderwidth([top, bottom], 1Escher.px)

function sierpinski(n)
    if n == 0
        compose(context(), polygon([(1,1), (0,1), (1/2, 0)]))
    else
        t = sierpinski(n - 1)
        compose(context(),
                (context(1/4,   0, 1/2, 1/2), t),
                (context(  0, 1/2, 1/2, 1/2), t),
                (context(1/2, 1/2, 1/2, 1/2), t))
    end
end

nᵗ = Signal(5) # The angle at any given time
connected_slider = subscribe(nᵗ, slider(0:7, value=5))

reactive_eg = map(nᵗ) do n
    vbox(
        connected_slider,
        sierpinski(n)
    )
end

tex_eg = """f(x) = \\int_{-\\infty}^\\infty
    \\hat f(\\xi)\,e^{2 \\pi i \\xi x}
    \\,d\\xi
"""
md_example = md"""**Things to do:**

- Create *universe*
- Make a *pie*
"""

part1 = md"""

$(vskip(1em))

$(

vbox(
hline(color=colorant"#e1e1e1"),
vskip(2em),
Escher.fontsize(2em,
     "With Escher you can build beautiful Web UIs entirely in Julia.") |>
        textalign(centertext) |>
        fontweight(200) |>
        lineheight(1.2em) |>
        maxwidth(15em) |>
        x -> hbox(flex(), x, flex()),
vskip(2em),
hline(color=colorant"#e1e1e1"),
)
)
$(vskip(3em))

$(title(2, "What's inside") |> Escher.fontsize(1.75em))

$(
vbox(md"**A web server for 2016.** Escher's built-in web server allows you to create interactive UIs with very little code. It takes care of messaging between Julia and the browser under-the-hood. It can also hot-load code: you can see your UI evolve as you save your changes to it.",

md"**A rich functional library of UI components.** the built-in library functions support Markdown, Signal widgets, TeX-style Layouts, Styling, TeX, Code, Behaviors, Tabs, Menus, Slideshows, Plots (via [Gadfly](http://gadfly.org), [Vega](http://johnmyleswhite.github.io/Vega.jl/), [VegaLite](https://github.com/fredo-dedup/VegaLite.jl)) and Vector Graphics (via [Compose](http://composejl.org)) -- everything a Julia programmer would need to effectively visualize data or to create user-facing GUIs. The API comprehensively covers features from HTML and CSS, and also provides advanced features. Its user merely needs to know how to write code in Julia."
) |> pad([left, right], 2em)
)

$(vskip(2em))

# Installation

In a Julia REPL, run:

```julia
Pkg.add("Escher")
```

You might want to link escher executable to `/usr/local/bin` or somewhere in your `PATH`:

```sh
ln -s ~/.julia/v0.4/Escher/bin/escher /usr/local/bin/
```

# Usage

From a directory in which you want to serve Escher UI files, run:

```
<Escher-package-path/bin>/escher --serve
```

This will bring up a web server on port 5555. The `examples/` directory in `Pkg.dir("Escher")` contains a few examples. After running the escher server from this directory, you can visit `http://localhost:5555/<example-file.jl>` to see the output of `<example-file.jl>`. Note that examples containing plots may take a while to load the first time you visit them.

See `escher --help` for other options to the exectuable.

Alternatively, you can start the server from a Julia REPL:

```julia
julia> using Escher
julia> include(Pkg.dir("Escher", "src", "cli", "serve.jl"))
julia> cd(Pkg.dir("Escher", "examples")) # or any other directory
julia> escher_serve()
```

This might be what you need if you installed Escher using a dmg file on OSX.
# An Overview

Escher's API consitently employs some patterns. Understanding these rules gives you a great foundation to understanding Escher's comprehensive API.

## Pattern 1: UIs are immutable values

A UI in Escher is simply an immutable Julia value of the abstract type `Tile`. A Tile is rendered into a browser [DOM tree](http://en.wikipedia.org/wiki/Document_Object_Model) to actually display the UI.

Escher has all kinds of functions that generate Tiles.

**Example 1.**

`plaintext("Hello, World!")` is a Tile that contains the plain text 'Hello, World!'.

To actually see this in action, create a file called `hello.jl` in an directory where the Escher server is running. Save the following code in it:

```julia
function main(window)
    plaintext("Hello, World!")
end
```

Now if you visit `http://localhost:5555/hello.jl` you should see the text "Hello, World!" on the top left corner of the screen.

$(vskip(1em))
$(sidenote(md"The function `main` *must* take a window argument, and may or may not use it. The `window` object contains some information about the current browser window. Specifically, `window.assets` is an input signal which can be used to load HTML dependencies on-the-fly. `window.alive` is a boolean signal that tells you if the window is still open. `window.dimension` is a 2-tuple representing the current size of the window in pixels."))
$(vskip(1em))

**Example 2.**

The `md""` string macro can generate markdown tiles from a markdown string.
```julia

function main(window)
    md\"\"\"
**Things to do:**

- Create *universe*
- Make a *pie*
\"\"\"
end
```

*Output:*

$(md_example |> pad([left], 4em))

**Example 3.**

The `tex` function creates a TeX tile.
```julia
function main(window)
    push!(window.assets, "tex")

    tex("cos(x)")
end
```
*Output:*

$(Escher.tex(tex_eg, block=true))

**Example 3.**

[Gadfly](http://gadflyjl.org) plots are essentially immutable values too. Escher type-casts Gadfly plots to tiles.

```julia
using Gadfly

function main(window)
    plot(z=(x,y) -> x*exp(-(x-int(x))^2-y^2),
         x=linspace(-8,8,150), y=linspace(-2,2,150), Geom.contour)
end
```
*Output:*

$(
plot(z=(x,y) -> x*exp(-(x-int(x))^2-y^2),
     x=linspace(-8,8,150), y=linspace(-2,2,150), Geom.contour)

)

**Example 4.**

[Compose](http://composejl.org) graphics work the same way.

```julia
using Compose
using Colors

function sierpinski(n)
    if n == 0
        compose(context(), polygon([(1,1), (0,1), (1/2, 0)]))
    else
        t = sierpinski(n - 1)
        compose(context(),
                (context(1/4,   0, 1/2, 1/2), t),
                (context(  0, 1/2, 1/2, 1/2), t),
                (context(1/2, 1/2, 1/2, 1/2), t))
    end
end

main(window) = compose(sierpinski(6))

```
*Output:*

$(vskip(2em))
$(drawing(3Compose.inch, sqrt(3)/2*3Compose.inch, compose(sierpinski(6))) |> pad([left], 8em))

## Pattern 2: functions that modify a tile return a new tile

Library functions that take tiles as input do not modify the input, they return new tiles that contain the modifications intended. This is a result of Pattern 1, since tiles cannot be modified in-place.

**Example 1.**

This example puts a padding of 5mm around a TeX formula, fills the tile with gray color, and changes the font color. We start off by modifying a tex tile and get a series of modified tiles which build up to the final result. Each tile is immutable and can be used later. The `vbox` function stacks many tiles vertically.

```julia
function main(window)
    push!(window.assets, "tex")

    txt = tex("T = 2\\pi\\sqrt{L\\over g}")
    txt1 = fontcolor("#499", txt)
    txt2 = pad(5mm, txt1)
    txt3 = fillcolor("#eeb", txt2)

    vbox(txt, txt1, txt2, txt3)
end
```
*Output:*

$(begin
    txt = Escher.tex("T = 2\\pi\\sqrt{L\\over g}")
    txt1 = fontcolor("#499", txt)
    txt2 = pad(5Escher.mm, txt1)
    txt3 = fillcolor("#eeb", txt2)
    vbox(txt, txt1, txt2, txt3)
end)

You can of course chain these function calls if you just want the end result.

```julia
function main(window)
    txt = tex("T = 2\\pi\\sqrt{L\\over g}")
    fillcolor("#eeb", fontcolor("#499", pad(5mm, txt)))
end
```

$(fillcolor("#eeb", fontcolor("#499", Escher.pad(5Escher.mm, Escher.tex("T = 2\\pi\\sqrt{L\\over g}")))))

## Pattern 3: Escher functions have curried methods

Omitting the last tile argument to escher functions returns a 1-argument function that takes a tile.

For example, `pad(10mm)` returns an anonymous function of 1 argument which must be a tile, and that returns a new tile with the specified 10mm of padding.

Therefore, `pad(10mm, txt)` is equivalent to `pad(10mm)(txt)` or `txt |> pad(10mm)`. This is helpful when you want to apply, for example, the same padding to a all the tiles in a vector. e.g. `map(pad(10mm), [tile1, tile2])` will return a vector of two tiles with 10mm padding each.

Moreover, using the curried version with the `|>` infix operator makes for code that reads better.

**Example 1.**

The following is equivalent to the last example.

```julia
function main(window)
    tex("T = 2\\pi\\sqrt{L\\over g}") |>
        fontcolor("#499") |>
        pad(5mm) |>
        fillcolor("#eeb")
end
```

You can mentally read this as: to the plaintext Hello, World, apply font color red and then apply a padding of `10mm`.

**Example 2.**

Since the curried methods return 1-argument lambdas, they can be readily used in a call to `map` or any function that expects a 1-argument function as an argument. Here we are going to use `map` and a curried method of `pad` to pad a number of tiles.

```julia
function main(window)
    padded = map(pad([left, right], 1em), ["A", "B", "C", "D"])
    colors = ["#837", "#859", "#892", "#875"]
    tiles = map(fillcolor, colors, padded)

    hbox(intersperse(hskip(1em), tiles))
end
```

*Output:*

$(begin
    tiles = map(fillcolor, ["#837", "#859", "#892", "#875"],
        map(pad([left, right], 1em), ["A", "B", "C", "D"]))
    hbox(intersperse(hskip(1em), tiles))
end)

$(vskip(1em))

*Explanation:*

`map(pad([left, right], 1em), ["A", "B", "C", "D"])` returns a vector of 4 tiles that each have a padding of `1em` to their left and right. Note that `pad([left, right], 1em)` returns a 1-argument function, making this possible.
`map(fillcolor, colors, padded)` takes corresponding elements from the `colors` vector and `padded` vector and combines them using `fillcolor`, resulting in 4 colored tiles.  `intersperse(hskip(1em), tiles))` takes the 4-element vector `tiles` and returns a 7-element vector where `hskip(1em)` is *interspered* between each pair of adjacent tiles in `tiles`. `hskip(1em)` is simply a `1em` empty space along the horizontal axis. Finally, `hbox` puts the result of `intersperse` in a horizontal layout.

## Pattern 4: layout functions combine many Tiles into one

**Example 1.**

To stack tiles vertically, use `vbox`. To stack horizontally, use `hbox`.

```julia
# Four tiles from the previous example.
a,b,c,d = map(fillcolor, ["#837", "#859", "#892", "#875"],
    map(pad([left, right], 1em), ["A", "B", "C", "D"]))

function main(window)
    x = vbox(
        a, b, c, d
    )
    y = vbox(
        d, c, b, a
    )

    hbox(x, y)
end
```

$(vskip(1em))
Should result in an arrangement like this:

$(begin
    a,b,c,d = map(fillcolor, ["#837", "#859", "#892", "#875"],
        map(pad([left, right], 1em), ["A", "B", "C", "D"]))

    x = vbox(
        a, b, c, d
    )
    y = vbox(
        d, c, b, a
    )

    hbox(x, y)
end)

`x` and `y` are vertial arrangements of 4 tiles each, these arrangements are themselves put in a `hbox` to place `x` next to `y`.

**Example 2.**

`The `tabs` function combines tiles into a set of tabs.
```julia
tabs([
    hbox(icon("face"), hskip(1em), "Tab 1"),
    hbox(icon("explore"), hskip(1em), "Tab 2"),
    hbox(icon("extension"), hskip(1em), "Tab 3"),
])
```
*Output:*
$(
tabs([
    hbox(icon("face"), hskip(1em), "Tab 1"),
    hbox(icon("explore"), hskip(1em), "Tab 2"),
    hbox(icon("extension"), hskip(1em), "Tab 3"),
])
)
$(vskip(1em))

**Example 3.**

The `pages` function combines tiles into a set of *pages* - only a single page is visible at a time. Pages can be combined with tabs to allow switching between pages.

```julia
tabbar = tabs([
    hbox(icon("face"), hskip(1em), "Tab 1"),
    hbox(icon("explore"), hskip(1em), "Tab 2"),
    hbox(icon("extension"), hskip(1em), "Tab 3"),
])

tabcontent = pages([
    sierpinski(5),
    plot([sin, cos], 0, 25),
    title(3, "web component all the things"),
])

t, p = wire(tabbar, tabcontent, :tab_channel, :selected)
       # ^^^ returns a pair of "connected" tab set and pages
vbox(t, p)

```

$(begin

tabbar = tabs([
    hbox(icon("face"), hskip(1em), "Tab 1"),
    hbox(icon("explore"), hskip(1em), "Tab 2"),
    hbox(icon("extension"), hskip(1em), "Tab 3"),
])

tabcontent = pages([
    drawing(3Compose.inch, 3*sqrt(3)/2*Compose.inch, sierpinski(5)),
    drawing(5Compose.inch, 3Compose.inch, plot([sin, cos], 0, 25)),
    title(3, "web component all the things"),
])
t, p = wire(tabbar, tabcontent, :tab_channel, :selected)
vbox(
    t, p |> pad(1em)
) |> size(40em, 22em)
end)

Other higher-order layout functions to try are: `menu`, `submenu`, `slideshow`.

## Pattern 5: An interactive UI is a Signal of UIs

[Reactive.jl](http://julialang.org/Reactive.jl) package allows "reactive programming" in Julia. Reactive programming is a style of event-driven programming with signals of data. A signal is a value that can change over time.  [Reactive.jl's documentation](http://julialang.org/Reactive.jl) provides an overview of the signal framework. At this point it is highly recommended that you read it.

There are two facets to this rule:

* Getting the input from tiles
* Creating a signal of UI using these signals

Firstly, some Tiles (particularly those that are subtypes of `Behavior` which is in turn a subtype of `Tile`) can write to Reactive's `Signal` signals. Widgets such as sliders, buttons, dropdown menus are subtypes of `Behavior`. The function `subscribe` lets you pipe updates from a behavior into a signal.

**Example 1.**

Let's create a slider and subscribe its current value to a Reactive signal.

```julia

iterations = Signal(0) # The angle at any given time
connected_slider = subscribe(slider(0:7), iterations)
```

*Output:*

$(slider(0:7))

`connected_slider` renders as a slider and updates the signal `iterations` when the slider's knob is moved by the user.

Secondly, you can use these input signals to create a signal of Escher UIs.

**Example 2.**

Let's now use the `iterations` signal to show an interactive Sierpinski's triangle. Note that you cannot meaningfully interact with the graphic if you are reading this on Github pages. Starting escher server from the `docs/` directory should get you interactivity.

```julia

using Compose

function main(window)
    # Load HTML dependencies related to the slider
    push!(window.assets, "widgets")

    iterations = Signal(5) # The number of iterations to show
    connected_slider = subscribe(slider(0:7, value=5), iterations)

    map(iterations) do n
        vbox(
            connected_slider,
            sierpinski(n)
        )
    end
end
```

"""

part2 = md"""

**Example 3.**

This is Minesweeper in 70 SLOC! About 89 with comments.

```julia
using Colors
using Lazy

#### Model ####

@defonce immutable Board{lost}
    uncovered::AbstractMatrix
    mines::AbstractMatrix
end

newboard(m, n, minefraction=0.2) =
    Board{false}(fill(-1, (m, n)), rand(m, n) .< minefraction)

function mines_around(board, i, j)
    m, n = size(board.mines)

    a = max(1, i-1)
    b = min(i+1, m)
    c = max(1, j-1)
    d = min(j+1, n)

    sum(board.mines[a:b, c:d])
end

### Update ###

next(board::Board{true}, move) = board
function next(board, move)
    i, j = move
    if board.mines[i, j]
        return Board{true}(board.uncovered, board.mines) # Game over
    else
        uncovered = copy(board.uncovered)
        if uncovered[i, j] == -1
            uncovered[i, j] = mines_around(board, i, j)
        end
        return Board{false}(uncovered, board.mines)
    end
end

movesᵗ = Signal((0, 0))
initial_boardᵗ = Signal(Board, newboard(10, 10))
boardᵗ = flatten(
    map(initial_boardᵗ) do b
        foldp(next, b, movesᵗ; typ=Board)
    end
)

### View ###

colors = ["#fff", colormap("reds", 7)]

box(content, color) =
    inset(Escher.middle,
        fillcolor(color, size(4em, 4em, empty)),
        Escher.fontsize(2em, content)) |> paper(1) |> pad(0.2em)

number(x) = box(x == -1 ? "" : string(x) |> fontweight(800), colors[x+2])
mine = box(icon("report"), "#e58")
tile(board::Board{true}, i, j) =
    board.mines[i, j] ? mine :
        number(board.uncovered[i, j])

tile(board, i, j) =
     constant((i, j), clickable(number(board.uncovered[i, j]))) >>> movesᵗ

gameover = vbox(
        title(2, "Game Over!") |> pad(1em),
        addinterpreter(_ -> newboard(10, 10), broadcast(button("Start again"))) >>> initial_boardᵗ
    ) |> pad(1em) |> fillcolor("white")

function showboard{lost}(board::Board{lost})
    m, n = size(board.mines)
    b = hbox([vbox([tile(board, i, j) for j in 1:m]) for i in 1:n])
    lost ? inset(Escher.middle, b, gameover) : b
end

function main(window)
    push!(window.assets, "widgets")

    map(boardᵗ) do board
        vbox(
           vskip(2em),
           title(3, "minesweeper"),
           vskip(2em),
           showboard(board),
        ) |> packacross(center)
    end
end
```

*Output:*

Below is a screen grab of a game that was just lost.

(Run `examples/minesweeper.jl` via escher server to play the game)

![](pkg/Escher/img/minesweeper.png)

# Documentation

The documentation is a work in progress. Here are the other bits available for now:

### API Reference

* [Layout API](layout-api.html)
* [Higher Order Layouts API](layout2-api.html)
* [Embellishment API](embellishment-api.html)
* [Typography API](typography-api.html)
* [Content API](content-api.html)
* [Widgets API](widgets-api.html)
* [Behavior API](behavior-api.html)
* [Signal API](signal-api.html)
* [Slideshow API](slideshow-api.html)
* [TeX API](tex-api.html)
* [Util API](util-api.html)

### WIP: User Guides

* [Layout Guide](layout-guide.html)
* [Typography](typography.html)
* [Reactive programming Guide](reactive.html)

Any help with documenting Escher will be appreciated. Take a look at [this issue](https://github.com/shashi/Escher.jl/issues/26) to see where you can contribute.

"""

include("helpers/page.jl")

function main(window)
    push!(window.assets, "widgets")
    push!(window.assets, "tex")
    push!(window.assets, "codemirror")
    push!(window.assets, "layout2")

    map(reactive_eg) do example
        vbox(
            pkgname(),
            vskip(1em),
            part1,
            example,
            part2,
        ) |> docpage
    end
end


text:
const cosmic_fl = joinpath(homedir(), "haplox/Bone/data/cosmic/CosmicCompleteExport.tsv")

function numsamplesofcancer(cancer)
    sampleids =  Set{ASCIIString}()
    open(cosmic_fl) do file
        header = readline(file)
        header  = split(header, "\t")
        while !eof(file)
            line = readline(file)
            row  = split(line, "\t")
            if row[8] != cancer
                continue
            end
            sample_id = convert(ASCIIString,row[6])
            push!(sampleids, sample_id)
        end
    end
    #@show sampleids
    length(sampleids)
end

function samplerecordsofcancer(cancer, threshold)
    sampleid_numrecords =  Dict{ASCIIString,Int64}()
    open(cosmic_fl) do file
        header = readline(file)
        header  = split(header, "\t")
        while !eof(file)
            line = readline(file)
            row  = split(line, "\t")
            if row[8] != cancer
                continue
            end
            sample_id = convert(ASCIIString,row[6])
            if in(sample_id, keys(sampleid_numrecords))
                sampleid_numrecords[sample_id] += 1
            else
                sampleid_numrecords[sample_id] = 1
            end
        end
    end
    #@show sampleids
    length(filter(x->x[2]>threshold, sort(collect(sampleid_numrecords),by=x->x[2],rev=true)))
end

@show samplerecordsofcancer("small_intestine", 4)


text:
matmul_num_threads = Ref(1)

"""
    ThreadedSparseCSR.set_num_threads(n::Int)

Sets the number of threads used in sparse csr matrix - vector multiplication.
"""
function set_num_threads(n::Int)

    0 < n <=  Threads.nthreads() || throw(DomainError("The numbers of threads must be > 0 and <= $(Threads.nthreads())."))

    global matmul_num_threads[] = n

end

"""
    ThreadedSparseCSR.get_num_threads()

Gets the number of threads used in sparse csr matrix - vector multiplication.
"""
function get_num_threads()

    return matmul_num_threads[]

end


text:
"""
    deleteCols(x, cols)

Delete columns from a 2x2 matrix. Use Float or Int array for x.
"""
function deleteCols(xIn, cols::Array{Int64, 1})

    # get dimensions of input
    (numRows, numCols) = size(xIn)

    # for calculating new dimensions
    numColDelete = length(cols)

    # uninitialized array with new dimensions
    xOut = Array{Float64}(undef, numRows, (numCols - numColDelete))

    # dynamically build new matrix without columns to delete
    xIndexIn  = 1
    xIndexOut = 1
    for ii in cols
        xColsIn =  xIndexIn : (ii - 1)
        xColsOut =  xIndexOut : xIndexOut + length(xColsIn) - 1

        # setting xOut cols to xIn, xColsOut always <= xColsIn
        xOut[:, xColsOut] = xIn[:, xColsIn]

        xIndexIn = ii + 1
        xIndexOut = xIndexOut + length(xColsIn)
    end

    # need one more iteration to get last section of cols
    if xIndexIn <= numCols
        xColsIn = xIndexIn : numCols
        xColsOut = xIndexOut : xIndexOut + length(xColsIn) - 1

        xOut[:, xColsOut] = xIn[:, xColsIn]
    end

    return xOut
end


text:
import ThresholdStability: perms

@testset "Ordering" begin
    @test perms(1) == [[1.],[-1.]]
    @test perms(2) == [[1.,1.],[1.,-1.],[-1.,1.],[-1.,-1.]]
    @test perms(3) == [[1.,1.,1.],[1.,1.,-1.],[1.,-1.,1.],[-1.,1.,1.],[1.,-1.,-1.],[-1.,1.,-1.],[-1.,-1.,1.],[-1.,-1.,-1.]]
end

function Base.:≈(X::T, Y::T) where T <: LightAutomaton
    if X.G == Y.G && X.Σ == Y.Σ && X.nt == Y.nt && X.next_id == Y.next_id
        return true
    end
    false
end

Σ = [randn(2,2) for _ in 1:4]
G = LightAutomaton(4)
add_transition!(G, 1, 1, 1)
add_transition!(G, 1, 3, 1)
add_transition!(G, 2, 1, 2)
add_transition!(G, 2, 3, 2)
add_transition!(G, 3, 2, 3)
add_transition!(G, 3, 4, 3)
add_transition!(G, 4, 2, 4)
add_transition!(G, 4, 4, 4)

H = OneStateAutomaton(4)

@testset "Automata" begin
    @test G ≈ automaton_constructor(Σ)
    @testset "Events" begin
        for q in 1:4
            for r in 1:4
                if has_transition(G, q, r)
                    @test event(G, q, r) == q
                else
                    @test event(G, q, r) == nothing
                end
            end
        end
        @test event(H, 1, 1) == [1,2,3,4]
        @test_throws AssertionError("q == 1") event(H, 2, 3)
        @test_throws AssertionError("r == 1") event(H, 1, 3)
    end
    for n in 2:7
        Σn = [randn(2,2) for _ in 1:2^n]
        @test automaton_constructor(Σn) ≈ automaton_constructor(Σn, perms(n))
    end
    @test islight(G) == true
    @test islight(OneStateAutomaton(2)) == false
end


text:
mutable struct NeutralTracker
    # structure for storing simulation data
    pop::Population
    generations::Int64
    verbose_track::Bool                         # print messages as data are tracked
    ens_strat_dist::Array{Float64, 2}           # track frequencies of strategies
    ens_strat_running::Array{Float64, 2}        # track running average of frequencies of strategies

    ens_y::Array{Float64, 1}                    # track quantity y = P(s_ia == s_ja | i != j)
    ens_y_running::Array{Float64, 1}            # track running average of y
    ens_z::Array{Float64, 1}                    # track quantity z = <h_i h_l | i != j>
    ens_z_running::Array{Float64, 1}            # track running average of z
    ens_g::Array{Float64, 1}                    # track quantity g = <h_i h_l I(s_ia == s_la) | i != l>
    ens_g_running::Array{Float64, 1}            # track running average of g
    ens_h::Array{Float64, 1}                    # track quantity h = <h_j h_l I(s_ia == s_ja) | i != j != l != i>
    ens_h_running::Array{Float64, 1}            # track running average of h

    ens_zp::Array{Float64, 1}                    # track quantity z' = < |h_i||h_l| | i != j>
    ens_zp_running::Array{Float64, 1}            # track running average of z'
    ens_gp::Array{Float64, 1}                    # track quantity g' = < |h_i||h_l| I(s_ia == s_la) | i != l>
    ens_gp_running::Array{Float64, 1}            # track running average of g'
    ens_hp::Array{Float64, 1}                    # track quantity h' = < |h_i||h_l| I(s_ia == s_ja) | i != j != l != i>
    ens_hp_running::Array{Float64, 1}            # track running average of h'

    ens_sia_sjd::Array{Float64, 1}              # track quantity q = <s_ia s_jd | i, j> # not requiring i != j
    ens_sia_sjd_running::Array{Float64, 1}      # track running average of <s_ia s_jd | i, j> # not requiring i != j

    # other options
    ens_sia_sid::Array{Float64, 1}              # track quantity <s_ia s_id>
    ens_sia_sid_running::Array{Float64, 1}      # track quantity <s_ia s_id>
    # ens_bracket_sia_sjd::Array{Float64, 1}   # track quantity <s_ia s_jd>
    # ens_bracket_sia_sja::Array{Float64, 1}   # track quantity <s_ia s_ja>
    # ens_bracket_sid_sjd::Array{Float64, 1}   # track quantity <s_id s_jd> # should be identical to above
    # ens_bracket_sia_hi_hl::Array{Float64, 1}   # track quantity <s_ia h_i h_l>
    # ens_bracket_sia_sid_hi_hl::Array{Float64, 1}   # track quantity <s_ia s_la h_i h_l>
    # ens_bracket_sia_sla_hi_hl::Array{Float64, 1}   # track quantity <s_ia s_la h_i h_l>
    # ens_bracket_sia_sld_hi_hl::Array{Float64, 1}   # track quantity <s_ia s_ld h_i h_l>
    # ens_bracket_sia_sja_hj_hl::Array{Float64, 1}   # track quantity <s_ia s_ja h_i h_l> # triplet
    # ens_bracket_sia_sja_hj_hl::Array{Float64, 1}   # track quantity <s_ia s_ja h_i h_l> # triplet

    # constructor
    function NeutralTracker( pop::Population, generations::Int64, verbose_track::Bool )
        # initialize n_trackers
        ens_strat_dist    = Float64[0.0 for gens=1:generations, strats=1:length(pop.all_strategies)]
        ens_strat_running = Float64[0.0 for gens=1:generations, strats=1:length(pop.all_strategies)]

        ens_y         = Float64[0.0 for gens=1:generations]
        ens_y_running = Float64[0.0 for gens=1:generations]
        ens_z         = Float64[0.0 for gens=1:generations]
        ens_z_running = Float64[0.0 for gens=1:generations]
        ens_g         = Float64[0.0 for gens=1:generations]
        ens_g_running = Float64[0.0 for gens=1:generations]
        ens_h         = Float64[0.0 for gens=1:generations]
        ens_h_running = Float64[0.0 for gens=1:generations]

        ens_zp         = Float64[0.0 for gens=1:generations]
        ens_zp_running = Float64[0.0 for gens=1:generations]
        ens_gp         = Float64[0.0 for gens=1:generations]
        ens_gp_running = Float64[0.0 for gens=1:generations]
        ens_hp         = Float64[0.0 for gens=1:generations]
        ens_hp_running = Float64[0.0 for gens=1:generations]

        ens_sia_sjd         = Float64[0.0 for gens=1:generations]
        ens_sia_sjd_running = Float64[0.0 for gens=1:generations]
        ens_sia_sid         = Float64[0.0 for gens=1:generations]
        ens_sia_sid_running = Float64[0.0 for gens=1:generations]

        return new(pop, generations, verbose_track,
                ens_strat_dist, ens_strat_running,
                ens_y, ens_y_running, ens_z, ens_z_running,
                ens_g, ens_g_running, ens_h, ens_h_running,
                ens_zp, ens_zp_running, ens_gp, ens_gp_running,
                ens_hp, ens_hp_running,
                ens_sia_sjd, ens_sia_sjd_running,
                ens_sia_sid, ens_sia_sid_running
                )
    end
end

function track_neutral!( n_tracker::NeutralTracker, gen::Int64 )
    # main function to track the population dynamics under neutral selection

    if n_tracker.verbose_track println("generation $gen") end

    track_strat_dist!(n_tracker, gen) # track strategy distribution and running average

    # singles
    for i in 1:n_tracker.pop.sets.N n_tracker.ens_sia_sid[gen] += get_sia_times_sid(n_tracker, i) end
    n_tracker.ens_sia_sid[gen] = n_tracker.ens_sia_sid[gen] / n_tracker.pop.sets.N

    # duplets, i != j
    for i in 1:n_tracker.pop.sets.N # keep it simple
        for j in i+1:n_tracker.pop.sets.N # i != j
            n_tracker.ens_y[gen] += get_sia_eq_sja(n_tracker, i, j)
            n_tracker.ens_z[gen] += get_hi_dot_hj(n_tracker, i, j)
            n_tracker.ens_g[gen] += get_hi_dot_hj_sia_eq_sja(n_tracker, i, j)

            n_tracker.ens_zp[gen] += get_hbari_dot_hbarj(n_tracker, i, j)
            n_tracker.ens_gp[gen] += get_hbari_dot_hbarj_sia_eq_sja(n_tracker, i, j)
        end
    end
    n_tracker.ens_y[gen]  = n_tracker.ens_y[gen] / n_tracker.pop.sets.n_pairs
    n_tracker.ens_z[gen]  = n_tracker.ens_z[gen] / n_tracker.pop.sets.n_pairs
    n_tracker.ens_g[gen]  = n_tracker.ens_g[gen] / n_tracker.pop.sets.n_pairs
    n_tracker.ens_zp[gen] = n_tracker.ens_zp[gen] / n_tracker.pop.sets.n_pairs
    n_tracker.ens_gp[gen] = n_tracker.ens_gp[gen] / n_tracker.pop.sets.n_pairs

    # duplets, i == j included
    for i in 1:n_tracker.pop.sets.N # keep it simple
        for j in i:n_tracker.pop.sets.N # i == j in included
            n_tracker.ens_sia_sjd[gen] += get_sia_times_sjd(n_tracker, i, j)
        end
    end
    n_tracker.ens_sia_sjd[gen] = n_tracker.ens_sia_sjd[gen] / (n_tracker.pop.sets.n_pairs + n_tracker.pop.sets.N)

    # triplets, i != j != l (!= i)
    triplet_counter = 0 # eventually store this in the n_tracker struct somewhere
    for i in 1:n_tracker.pop.sets.N # keep it simple
        for j in i+1:n_tracker.pop.sets.N # i != j
            for l in j+1:n_tracker.pop.sets.N
                n_tracker.ens_h[gen] += get_hi_dot_hl_sia_eq_sja(n_tracker, i, j, l)
                n_tracker.ens_h[gen] += get_hi_dot_hl_sia_eq_sja(n_tracker, j, l, i)
                n_tracker.ens_h[gen] += get_hi_dot_hl_sia_eq_sja(n_tracker, l, i, j)

                n_tracker.ens_hp[gen] += get_hbari_dot_hbarl_sia_eq_sja(n_tracker, i, j, l)
                n_tracker.ens_hp[gen] += get_hbari_dot_hbarl_sia_eq_sja(n_tracker, j, l, i)
                n_tracker.ens_hp[gen] += get_hbari_dot_hbarl_sia_eq_sja(n_tracker, l, i, j)

                triplet_counter += 3
            end
        end
    end
    n_tracker.ens_h[gen]  = n_tracker.ens_h[gen] / triplet_counter
    n_tracker.ens_hp[gen] = n_tracker.ens_hp[gen] / triplet_counter

    # update running average
    if gen == 1
        # at gen = 1, running average = dist at gen = 1 itself
        n_tracker.ens_y_running[gen] = n_tracker.ens_y[gen]
        n_tracker.ens_z_running[gen] = n_tracker.ens_z[gen]
        n_tracker.ens_g_running[gen] = n_tracker.ens_g[gen]
        n_tracker.ens_h_running[gen] = n_tracker.ens_h[gen]
        n_tracker.ens_zp_running[gen] = n_tracker.ens_zp[gen]
        n_tracker.ens_gp_running[gen] = n_tracker.ens_gp[gen]
        n_tracker.ens_hp_running[gen] = n_tracker.ens_hp[gen]
        n_tracker.ens_sia_sjd_running[gen] = n_tracker.ens_sia_sjd[gen]
        n_tracker.ens_sia_sid_running[gen] = n_tracker.ens_sia_sid[gen]
    else
        # at gen > 1, running average = ( (running avg y at gen-1) * (gen-1) + y at gen) / gen
        n_tracker.ens_y_running[gen] = ( n_tracker.ens_y_running[gen-1]*(gen-1) + n_tracker.ens_y[gen] ) / gen
        n_tracker.ens_z_running[gen] = ( n_tracker.ens_z_running[gen-1]*(gen-1) + n_tracker.ens_z[gen] ) / gen
        n_tracker.ens_g_running[gen] = ( n_tracker.ens_g_running[gen-1]*(gen-1) + n_tracker.ens_g[gen] ) / gen
        n_tracker.ens_h_running[gen] = ( n_tracker.ens_h_running[gen-1]*(gen-1) + n_tracker.ens_h[gen] ) / gen
        n_tracker.ens_zp_running[gen] = ( n_tracker.ens_zp_running[gen-1]*(gen-1) + n_tracker.ens_zp[gen] ) / gen
        n_tracker.ens_gp_running[gen] = ( n_tracker.ens_gp_running[gen-1]*(gen-1) + n_tracker.ens_gp[gen] ) / gen
        n_tracker.ens_hp_running[gen] = ( n_tracker.ens_hp_running[gen-1]*(gen-1) + n_tracker.ens_hp[gen] ) / gen
        n_tracker.ens_sia_sjd_running[gen] = ( n_tracker.ens_sia_sjd_running[gen-1]*(gen-1) + n_tracker.ens_sia_sjd[gen] ) / gen
        n_tracker.ens_sia_sid_running[gen] = ( n_tracker.ens_sia_sid_running[gen-1]*(gen-1) + n_tracker.ens_sia_sid[gen] ) / gen
    end
end

function get_sia_times_sid( n_tracker::NeutralTracker, i::Int64 )
    # for <s_ia s_id>
    return Float64( n_tracker.pop.strategies[i,1] * n_tracker.pop.strategies[i,2] )
end

function get_sia_eq_sja( n_tracker::NeutralTracker, i::Int64, j::Int64 )
    # for y = P(s_ia == s_ja | i != j)
    return Float64( n_tracker.pop.strategies[i,1] == n_tracker.pop.strategies[j,1] )
end

function get_sia_times_sjd( n_tracker::NeutralTracker, i::Int64, j::Int64 )
    # for q = <s_ia s_jd | i, j>, including i == j
    return Float64( n_tracker.pop.strategies[i,1] * n_tracker.pop.strategies[j,2] )
end

function get_hi_dot_hj( n_tracker::NeutralTracker, i::Int64, j::Int64 )
    # for z = <h_i h_l | i != j>
    product = 0.
    for k in 1:n_tracker.pop.sets.M
        product += Float64( n_tracker.pop.sets.h[i,k] * n_tracker.pop.sets.h[j,k] )
    end
    return product / n_tracker.pop.sets.M
    # return Float64( n_tracker.pop.sets.h[i] * n_tracker.pop.sets.h[j] )
end

function get_hi_dot_hj_sia_eq_sja( n_tracker::NeutralTracker, i::Int64, j::Int64 )
    # for g = <h_i h_l I(s_ia == s_la) | i != l>
    # ie h_i h_l if s_ia = s_la is true, else 0
    product = 0.
    if n_tracker.pop.strategies[i,1]==n_tracker.pop.strategies[j,1]
        for k in 1:n_tracker.pop.sets.M
            product += Float64( n_tracker.pop.sets.h[i,k] * n_tracker.pop.sets.h[j,k] )
        end
    end
    return product / n_tracker.pop.sets.M
    # return n_tracker.pop.strategies[i,1]==n_tracker.pop.strategies[j,1] ? Float64( n_tracker.pop.sets.h[i] * n_tracker.pop.sets.h[j] ) : 0.0
end

function get_hi_dot_hl_sia_eq_sja( n_tracker::NeutralTracker, i::Int64, j::Int64, l::Int64)
    # for h = <h_j h_l I(s_ia == s_ja) | i != j != l != i>
    # ie h_i h_l if s_ia = s_ja is true, else 0
    product = 0.
    if n_tracker.pop.strategies[i,1]==n_tracker.pop.strategies[j,1]
        for k in 1:n_tracker.pop.sets.M
            product += Float64( n_tracker.pop.sets.h[i,k] * n_tracker.pop.sets.h[l,k] )
        end
    end
    return product / n_tracker.pop.sets.M
    # return n_tracker.pop.strategies[i,1]==n_tracker.pop.strategies[j,1] ? Float64( n_tracker.pop.sets.h[i] * n_tracker.pop.sets.h[l] ) : 0.0
end

function get_hbari_dot_hbarj( n_tracker::NeutralTracker, i::Int64, j::Int64 )
    # for z' = < |h_i||h_l| | i != j>
    product = 0.
    for k in 1:n_tracker.pop.sets.M
        product += Float64( n_tracker.pop.sets.h_bar[i,k] * n_tracker.pop.sets.h_bar[j,k] )
    end
    return product / n_tracker.pop.sets.M
end

function get_hbari_dot_hbarj_sia_eq_sja( n_tracker::NeutralTracker, i::Int64, j::Int64 )
    # for g' = < |h_i||h_l| I(s_ia == s_la) | i != l>
    product = 0.
    if n_tracker.pop.strategies[i,1]==n_tracker.pop.strategies[j,1]
        for k in 1:n_tracker.pop.sets.M
            product += Float64( n_tracker.pop.sets.h_bar[i,k] * n_tracker.pop.sets.h_bar[j,k] )
        end
    end
    return product / n_tracker.pop.sets.M
end

function get_hbari_dot_hbarl_sia_eq_sja( n_tracker::NeutralTracker, i::Int64, j::Int64, l::Int64)
    # for h' = < |h_j||h_l| I(s_ia == s_ja) | i != j != l != i>
    product = 0.
    if n_tracker.pop.strategies[i,1]==n_tracker.pop.strategies[j,1]
        for k in 1:n_tracker.pop.sets.M
            product += Float64( n_tracker.pop.sets.h_bar[i,k] * n_tracker.pop.sets.h_bar[l,k] )
        end
    end
    return product / n_tracker.pop.sets.M
end

function track_strat_dist!( n_tracker::NeutralTracker, gen::Int64 )
    # compute and store frequency of individuals with each strategy
    # strategies_base10           = [2*tracker.pop.strategies[i,1] + n_tracker.pop.strategies[i,2] for i in 1:tracker.pop.sets.N]
    n_tracker.ens_strat_dist[gen,:] = [count(x->x==strat, n_tracker.pop.strategies_base10) / n_tracker.pop.sets.N for strat in n_tracker.pop.all_strategies]

    # track running average
    if gen == 1 # at g = 1, running average = dist at g = 1 itself
        n_tracker.ens_strat_running[gen,:] = n_tracker.ens_strat_dist[gen,:]
    else # at g > 1, running average = (dist at g-1 + dist at g) / g
        n_tracker.ens_strat_running[gen,:] = ( n_tracker.ens_strat_running[gen-1,:]*(gen-1) + n_tracker.ens_strat_dist[gen,:] ) / gen
    end
end


text:
using Test
using Interpolations
using IncrementalInference
using Distributions
using TensorCast
using Distributions

##
@testset "Test HeatmapGridDensity" begin
##

x = -10:0.1:10;
y = -10:0.1:10;

img = zeros(length(x),length(y))

# lambda dispatch in this scope was REALLY slow
mv = MvNormal(Diagonal([1.0;1.0]))
g(x,y) = pdf(mv,[x;y])

for (i,x_) in enumerate(x), (j,y_) in enumerate(y)
  img[i,j] = g(x_,y_)
end

##
println("build a HeatmapGridDensity")
hgd = IIF.HeatmapGridDensity(img, (x,y), nothing, 0.07; N=1000)

println("test packing converters")
# check conversions to packed types
phgd = convert(PackedSamplableBelief, hgd)
hgd_ = convert(SamplableBelief, phgd)

@test isapprox( hgd, hgd_ )

# use in graph

## Check that sampling of the HMD is correct

pts_ = sample(hgd,1000)[1]

@cast pts_x[i] := pts_[i][1]
@cast pts_y[i] := pts_[i][2]

f = fit(MvNormal, hcat(pts_x,pts_y)')

@test isapprox([0;0], f.μ; atol=0.15)
@test isapprox([1 0; 0 1], f.Σ.mat; atol=0.25)

##
end


text:
include("../../wrapper/preprocWrapper.jl")
include("../../../src/preproc.jl")

# using BenchmarkTools
using Plots
using DelimitedFiles

function preprocplot(m, p, n)
    # params: row(A): row(B): col ratio

    maxI = 10
    maxJ = 100
    x = zeros(maxI)
    twrapper = zeros(maxI)
    tnative = zeros(maxI)

    for i = 1:maxI
        x[i] = m*i*2
        time1 = 0.0
        time2 = 0.0
        for j = 1:maxJ
            A1 = randn(m*i*2, n*i*2)
            B1 = randn(p*i*2, n*i*2)
            A2 = deepcopy(A1)
            B2 = deepcopy(B1)
            ans1 = @timed dggsvp3!(A1, B1)
            ans2 = @timed preproc!(A2, B2)
            time1 = time1 + ans1[2]
            time2 = time2 + ans2[2]
        end
            twrapper[i] = time1/maxJ
            tnative[i] = time2/maxJ

        # e = eps(Float64)
        # res_a1 = opnorm(ans1[1][1]'*A_*ans1[1][3] - ans1[1][4]*ans1[1][6], 1)/(max(10*m*i,10*n*i)*opnorm(A_, 1)*e)
        # res_a2 = opnorm(ans2[1].U'*A_*ans2[1].Q - ans2[1].D1*ans2[1].R0, 1)/(max(10*m*i,10*n*i)*opnorm(A_, 1)*e)
        println("m: ", x[i])
        # println("residual err of A by new gsvd: ", res_a1)
        # println("residual err of A by current gsvd in Julia: ", res_a2)
        println("--------------------------------------------------------")
    end

    y = [twrapper tnative];
    gr();
    pl = plot(x, y,
    title = "CPU timing of preproc, m:p:n = $m:$p:$n",
    label = ["LAPACK-wrapper" "Native-Julia"],
    xlabel = "m: number of rows of A",
    ylabel = "elapsed time in seconds",
    yformatter = :scientific,
    markershape = [:circle :xcross],
    markersize = 3.5,
    markeralpha = 0.75,
    markerstrokewidth=0,
    lw = 1.5,
    legend = :topleft)
    # savefig("myplot.png") # Saves the CURRENT_PLOT as a .png
    # savefig(pl, "myplot.pdf") # Saves the plot from p as a .pdf vector graphic
    # return y, pl
    open("preproc-m-p-n=$m-$p-$n.txt", "w") do io
        writedlm(io, y)
    end
    savefig(pl, "preproc-m-p-n=$m-$p-$n.pdf")
    savefig("preproc-m-p-n=$m-$p-$n.png")
end

function preprocbarplot(m, p, n)
    # y = randn(20);
    y = readdlm("preproc-m-p-n=$m-$p-$n.txt")
    maxI = 10
    gr();
    pl = bar(
    string.(m*2:m*2:m*maxI*2),
    y[:,1]./y[:,2],
    xlabel = "m: number of rows of A",
    xticks = :all,
    ylabel = "factors of speedup",
    title = "preproc: t(LAPACK-wrapper)/t(Native-Julia), m:p:n=$m:$n:$p",
    fillcolor = :royalblue1,
    fillalpha=0.25,
    legend = nothing)
    savefig(pl, "speedup_preproc-m-p-n=$m-$p-$n.pdf")
    savefig("speedup_preproc-m-p-n=$m-$p-$n.png")
end


text:
using Econlite
using Test
using ForwardDiff

@testset "Econlite.jl" begin

    # utility.jl
    tol = 1e-7
    l = LogUtility()
    crra = IsoElastic(0.42)
    cara = CARA(1.2)

    for ps in (l, crra, cara)
        c = 1.23
        u = util(c, ps)
        cc = invutil(u, ps)

        @test abs(c - cc) < tol

        du=dutil(c, ps)
        ccc=invdutil(du, ps)

        @test abs(c - ccc) < tol

        # check with forward diff
        du_fd = ForwardDiff.derivative(t -> util(t, ps), c)

        @test abs(du_fd - du) < tol
    end

    # aggregators.jl

end


text:
include("base.jl")
include("agents/agents.jl")
include("q_based_policies/q_based_policies.jl")
include("v_based_policies.jl")
include("tabular_random_policy.jl")
include("random_policy.jl")
include("random_start_policy.jl")


text:
using Test

# Use the following command to clean the notebooks if the tests fail
# jupyter nbconvert --ClearOutputPreprocessor.enabled=True --inplace *.ipynb
#
# judge if notebook is clean
@testset "Clean Notebooks" begin
    for _file in readdir(joinpath(@__DIR__, "../notebooks"))
        if _file[end-4:end] == "ipynb"
            for _line in eachline(joinpath(@__DIR__, "../notebooks", _file))
                # execution count is null?
                if occursin("   \"execution_count\":", _line)
                    @test _line == "   \"execution_count\": null,";
                end
                # metadata is empty?
                if occursin("   \"metadata\":", _line)
                    @test _line == "   \"metadata\": {},"
                end
                # output is empty?
                if occursin("   \"outputs\":", _line)
                    @test _line == "   \"outputs\": [],";
                end
            end
        end
    end
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
#===============================================================================

    Newtonian Method

    Based on code on Ben Moll's webite:
            www.princton.edu/~moll/HACTproject.htm

===============================================================================#

#Need to figure out best way to transfer inputs in julia
include("Huggett_initial.jl")
include("Huggett_terminal.jl")
include("Huggett_subroutine.jl")

# Get values for g, r, A from initial
g0, gg0, r, A = initial_val()

r00 = r
A0=A
AT0 = A'

# Get values for g_st, v_st, r_st from terminal
g_st, v_st, r_st = terminal()

plot(a,v_st)

plot(a, g)
plot!(a,g0, xlims=(-0.15,1))

# Create time and time increment variables
global T = 20
global N = 100
global H =1000
global dt = T/N

f= zeros(N,1)
f1=zeros(N,1)

# Create and initial guess of the interest rate sequence
r0 = r_st*ones(N,1)

# Set up space to save variables
S = zeros(N,1)
global SS = zeros(N,1)
global dS = zeros(N,1)

v = zeros(H,2,N)

data=[]

push!(data, v_st)
push!(data, gg0)
push!(data, Aswitch)

maxit=100
convergence_criterion = 10^(-6)

# Set an intial guess of the interest rate squence
x = r_st.*ones(N,1)

#evaluate market clearing at inital guess
global A_t=[]
global gg=zeros(N+1,2*H,1)
f = huggett_subroutine(x, data)
println("Convergence Criterion at inital guess is: $(findmax(abs.(f))[1])")

# Compute the Jacobian at the inital guess
using SparseArrays
h = diagm(0=>sqrt(eps()).*max.(abs.(x),1)[:])

# Implement the Newtonian scheme
J = zeros(N,N)

# For loop for Newtonian method
for k in 1:N
    println("Coloumn of Jacobian = $(k)")

    x1 = x + h[:,k]
    global A_t=[]
    global gg=zeros(N+1,2*H,1)
    f1=huggett_subroutine(x1,data)
    J[:,k] = (f1-f)/h[k,k]
end

# Double check that the Jacobian has full rank
println("Rank of Jacobian = $(rank(J))")
J0 = J # save this inital Jacobian

# Dampening factor
ξ = 0.1

# main loop
Sdist=[]

for it in 1:maxit
    global J, f, x = J, f, x
    dx = -J\f

    x = x + ξ.*dx
    f = huggett_subroutine(x, data)
    push!(Sdist, findmax(abs.(f))[1])

    J = J + f *dx'./(dx'*dx)

    println("Iteration=$(it), Convergence Criterion=$(Sdist[it])")

    if Sdist[it]<convergence_criterion
        break
    end
end

plot(Sdist[:])

# Plot of r_t over time, note the oscillation at the end.
    # this is a common issue with Newtonian methods
r_t = x
plot(1:N, r_t[:], legend=false, xlabel="Period", ylabel="r")
plot!(1:N, r_st.*ones(N,1), line=:dash)

# More plots from tranistion dynamics code

N1 =4
T1=-N1*dt
time = (1:N)'*dt
time1 = T1 .+ (1:N1)'*dt
time2 = [time1'; time']
r_t2 = [r00*ones(N1,1); r_t]

plot(time2, r_t2, xlims=(T1, 10), ylims=(-0.05, 0.035), legend=:false)
plot!(time2, r_st*ones(N1+N,1), line=:dash,
        xlabel="Year", title="Equilibrium Interest Rate")

amax1 = 0.5
gmax=3
n=2

p1 = plot(a, gg[n,1:H,:],xlims=(amin,amax1),
        ylims=(0, gmax), xlabel="Wealth", ylabel="Densities",legend=false)
plot!(a, gg[n,H+1:end,:], title="\$t=0.1\$")

# set up new index
t=2
n=convert(Int64, t/dt)

p2 = plot(a, gg[n,1:H,:],xlims=(amin,amax1),
        ylims=(0, gmax), xlabel="Wealth", ylabel="Densities",legend=false)
plot!(a, gg[n,H+1:end,:], title="\$t=2\$")

t=5
n=convert(Int64, t/dt)

p3 = plot(a, gg[n,1:H,:],xlims=(amin,amax1),
        ylims=(0, gmax), xlabel="Wealth", ylabel="Densities",legend=false)
plot!(a, gg[n,H+1:end,:], title="\$t=5\$")

t=T
n=convert(Int64, t/dt)

p4 = plot(a, gg[n,1:H,:],xlims=(amin,amax1),
        ylims=(0, gmax), xlabel="Wealth", ylabel="Densities",legend=false)
plot!(a, gg[n,H+1:end,:], title="\$t=\\infty\$")

plot(p1,p2,p3,p4,legend=false,layout = (2,2))


text:
# ---
# title: 1655. Distribute Repeating Integers
# id: problem1655
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Backtracking
# link: <https://leetcode.com/problems/distribute-repeating-integers/description/>
# hidden: true
# ---
#
# You are given an array of `n` integers, `nums`, where there are at most `50`
# unique values in the array. You are also given an array of `m` customer order
# quantities, `quantity`, where `quantity[i]` is the amount of integers the
# `ith` customer ordered. Determine if it is possible to distribute `nums` such
# that:
#
#   * The `ith` customer gets **exactly** `quantity[i]` integers,
#   * The integers the `ith` customer gets are **all equal** , and
#   * Every customer is satisfied.
#
# Return `true` _if it is possible to distribute_`nums` _according to the above
# conditions_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,4], quantity = [2]
#     Output: false
#     Explanation: The 0th customer cannot be given two different integers.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3,3], quantity = [2]
#     Output: true
#     Explanation: The 0th customer is given [3,3]. The integers [1,2] are not used.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,1,2,2], quantity = [2,2]
#     Output: true
#     Explanation: The 0th customer is given [1,1], and the 1st customer is given [2,2].
#
#
# **Example 4:**
#
#
#
#     Input: nums = [1,1,2,3], quantity = [2,2]
#     Output: false
#     Explanation: Although the 0th customer could be given [1,1], the 1st customer cannot be satisfied.
#
# **Example 5:**
#
#
#
#     Input: nums = [1,1,1,1,1], quantity = [2,3]
#     Output: true
#     Explanation: The 0th customer is given [1,1], and the 1st customer is given [1,1,1].
#
#
#
#
# **Constraints:**
#
#   * `n == nums.length`
#   * `1 <= n <= 105`
#   * `1 <= nums[i] <= 1000`
#   * `m == quantity.length`
#   * `1 <= m <= 10`
#   * `1 <= quantity[i] <= 105`
#   * There are at most `50` unique values in `nums`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Dice Simulator

using Random

function normaldice()
	return Random.rand(collect(1:6))
end

function specialdice()
	special = [1,6,2,6,3,4,5,6]
	return Random.rand(special)
end

println("You need to roll Normal Dice or Special Dice (s/n): ")
choice = readline()

if choice == "n"
	println(normaldice())
elseif choice == "s"
	println(specialdice())
else
	println("Wrong Choice !")
end

max_stars_repo_name: cashewcodes/Basic-Julia-Projects

max_issues_repo_name: cashewcodes/Basic-Julia-Projects

max_forks_repo_name: cashewcodes/Basic-Julia-Projects


text:
########################################################################################################
# Generates replicates from two base model files for group comparison testing
########################################################################################################

# Dependencies
using Distributions
using DelimitedFiles

# Args
n_reps_g1 = parse(Int64, ARGS[1])
n_reps_g2 = parse(Int64, ARGS[2])

# Read base models
g1_mod = readdlm("examples/grp_comp_example/mod_files/g1_base_model.txt")
g2_mod = readdlm("examples/grp_comp_example/mod_files/g2_base_model.txt")

## Replicates

# Group 1
for rep = 1:n_reps_g1
    rep_mod = deepcopy(g1_mod)
    for i = 1:size(rep_mod)[1]
        θ = parse.(Float64, split(rep_mod[i,4], ','))
        rep_mod[i,4] = "$(θ[1] + rand(Normal())),$(θ[2] + rand(Normal())),$(θ[3] + rand(Normal()))"
    end
    writedlm("examples/grp_comp_example/mod_files/g1_rep$(rep).txt", rep_mod)
end

# Group 2
for rep = 1:n_reps_g2
    rep_mod = deepcopy(g2_mod)
    for i = 1:size(rep_mod)[1]
        θ = parse.(Float64, split(rep_mod[i,4], ','))
        rep_mod[i,4] = "$(θ[1] + rand(Normal())),$(θ[2] + rand(Normal())),$(θ[3] + rand(Normal()))"
    end
    writedlm("examples/grp_comp_example/mod_files/g2_rep$(rep).txt", rep_mod)
end


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# ╔═╡ 142b62e1-364e-45d6-9bca-7c69b794f8ce
md"# Tutorial 1"

# ╔═╡ 27f5936d-4043-4cb6-b795-033f6f2a0674
md"Crash course in Julia basics:"

# ╔═╡ 3f08e241-e53d-4d65-935f-ddf6a6bfbbdd
md"""
Arithmetic, arrays, tuples, strings, dictionaries, functions,
iteration, random numbers, package loading, plotting
"""

# ╔═╡ 3e3bcc53-7840-4af7-af2a-25235e544a31
md"(40 min)"

# ╔═╡ 96774cd2-130c-4455-8af5-148d95ea2900
md"## Julia is a calculator:"

# ╔═╡ 8e8a0cbc-96b9-490c-a6bf-47194d7e8e12
1 + 2^3

# ╔═╡ a65fb83a-f977-47ce-828a-4cc485149963
sqrt(1 + 2^3) # do `sqrt(ans)` in REPL

# ╔═╡ dc16182e-b6b3-47ee-9e54-3f7b3ca87ecb
sin(pi)

# ╔═╡ bc1912b8-c66a-4ab5-ba1f-6363f43ec697
md"""
Query a function's document string using `?sin` at the REPL, or in a
notebook:
"""

# ╔═╡ 048a1195-1107-46e9-95ea-2955abd45275
@doc sin

# ╔═╡ ef72056a-9321-4290-b1b3-dd47e367ba54
md"""
I've forgotten how the arcsin is called. Is it `asin` or `arcsin`? I
can search all doc-strings containing "sine" to locate the method:
"""

# ╔═╡ af710d1b-f8c2-4219-8d7e-8e449afd1c48
apropos("sine")

# ╔═╡ 3e690f1c-4e80-4016-a949-7f3bd292fe31
md"Okay, I see it's called `asin`."

# ╔═╡ b541c5ac-3c15-4e71-8514-99938a2932db
asin(1 + 3*im)

# ╔═╡ c00b48c5-d35b-409a-a0de-1721c1bc2df2
md"## Arrays"

# ╔═╡ 3694faf5-d9a4-4089-bca8-7eec7950fd82
md"One dimensional vectors:"

# ╔═╡ cad4d998-43d7-421b-9873-1b1430e635cc
v = [3, 5, 7]

# ╔═╡ 5c43b0f6-30f4-4ef7-b43e-d0ebe87da176
length(v)

# ╔═╡ 26ce1939-e595-496a-9009-4cfb2012998f
md"A \"row vector\" is a 1 x n array:"

# ╔═╡ d66140c7-958b-4347-abd3-cb77d7a79683
row = [3 5 7]

# ╔═╡ 7aee99b6-3c6e-4a71-879e-dc128f3db7b3
md"Multiple row vectors separated by semicolons or new lines define matrices:"

# ╔═╡ 96c1adb9-1b29-433d-a367-9aa3c6cf34d0
A = [3 5 7
     2 4 6
     1 3 5]

# ╔═╡ ef6a776c-8c63-4d30-bf32-94657e65284e
size(A)

# ╔═╡ 4bb2001b-3b46-40f3-9afd-65a5b5facac9
length(A)

# ╔═╡ b97d12e6-2346-410b-b6c8-54f5ed90a964
md"Accessing elements (Julia indices start at 1 not 0):"

# ╔═╡ 177a30bb-473a-4e5e-9292-fe822525fc77
A[1, 2]

# ╔═╡ 109564cb-ea6f-4264-ae5e-4b83dcbc9675
A[1, 2] == A[2]

# ╔═╡ 26714555-a4ad-49cb-8a29-32f11452654a
A[:, 2] # 2nd column

# ╔═╡ d00b2507-e891-4ec5-a5f1-151fcbe229a2
md"Changing elements:"

# ╔═╡ 4d1dd37a-194e-4ed2-81bc-2d5603785a09
A[1, 1] = 42

# ╔═╡ 38477997-58dd-4248-9d87-1d1fbb0e3997
inv(A) # inverse

# ╔═╡ e61c7ae2-8204-4a20-b950-fddef2a1fb10
isapprox(inv(A)*v, A\v) # but RHS more efficient

# ╔═╡ d38b7040-bb41-44be-951b-64cb2a36c8e3
md"## \"Variables\" in Julia *point* to objects"

# ╔═╡ 9f1e3326-b738-4b60-b0e6-441461cc8770
md"Corollary: all passing of function arguments is pass by reference"

# ╔═╡ 494f9e7c-e53a-4c65-8cb0-de601961bc02
md"Like Python; Unlike R, C or FORTRAN."

# ╔═╡ 259bf65c-6868-47b0-a87d-950c50fb2955
begin
  w = v

  w
end

# ╔═╡ 6ef57fdd-3188-4b19-8448-4bcb089096cd
v[1] = 42

# ╔═╡ ef2473e2-bb35-4bd7-a012-a84240254fb6
v

# ╔═╡ 76dc6ac8-03fe-4c6e-bbdd-9799f7bb2e60
w

# ╔═╡ 5b3a3587-31f7-4a00-97a7-88e3af4f10ee
md"## Tuples"

# ╔═╡ 82f84505-336f-4845-b372-8c9866e51852
md"Similar to vectors but of fixed length and immutable (cannot be changed)"

# ╔═╡ 5b71f2a1-5567-4e59-8f3d-b82f1e7b6f7a
begin
  t1 = (1, 2.0, "cat")
  typeof(t1)
end

# ╔═╡ 249c3a08-3ca6-463e-ab04-579e5608ae53
t1[3]

# ╔═╡ f85e0bbc-9db0-4a9f-86cf-35420d9e6995
md"Tuples also come in a *named* variety:"

# ╔═╡ b8b91bc0-0a63-4695-a29a-8d0445351914
t2 = (i = 1, x = 2.0, animal="cat")

# ╔═╡ e0e4f0da-69c7-4941-be65-28d0fcca50a8
t2.x

# ╔═╡ 166a1850-cddb-4f72-9a2e-c468345d87d1
md"## Strings and relatives"

# ╔═╡ 44ef9d71-f75b-4140-b5fa-cb79ebf595ef
begin
  a_string = "the cat"
  a_character = 't'
  a_symbol = :t
end

# ╔═╡ ff4b5618-4be7-4111-91c5-524da38aa391
a_string[1] == a_character

# ╔═╡ fb16016a-87d7-4e9f-ad90-b4405b216771
md"""
A `Symbol` is string-like but
[interned](https://en.wikipedia.org/wiki/String_interning). Generally
use `String` for ordinary textual data, but use `Symbol` for
language reflection (metaprogramming). For example:
"""

# ╔═╡ 20ba7245-15ed-417e-895b-997a92b731e0
isdefined(Main, :z)

# ╔═╡ 542ec737-2c88-4293-a525-e9b04a4bd246
begin
  z = 1 + 2im
  isdefined(Main, :z)
end

# ╔═╡ 8b1e2b91-947c-4046-80f0-ba7781e173cf
z.im

# ╔═╡ 3dcab89a-20ca-462a-9cbc-9a2c39793333
fieldnames(typeof(z))

# ╔═╡ 65eb310b-82e8-4a1d-b887-2b16710e5502
md"Symbols are generalized by *expressions*:"

# ╔═╡ 35998a8b-6fda-4a2f-9452-5298a8a4a401
begin
  ex = :(z == 3)
  eval(ex)
end

# ╔═╡ 2db5e586-237b-45aa-b017-96f1602f2cad
md"If this is confusing, forget it for now."

# ╔═╡ 8b1c08d9-f8dd-4ec9-8be2-2a4017c45345
md"## Dictionaries"

# ╔═╡ 4942a061-18dd-41d0-a7ad-80d84f5b0070
d = Dict('a' => "ant", 'z' => "zebra")

# ╔═╡ d1fb9147-4a06-44fa-8378-5aa686f0b407
d['a']

# ╔═╡ 4f4bf27c-b2f9-4c15-9f43-49763e8691a1
begin
  d['b'] = "bat"
  d
end

# ╔═╡ 29d71401-66e5-484e-bb0e-3a45761c733a
keys(d)

# ╔═╡ 4b80c58d-646a-405b-9979-7044b2f2df33
md"The expression 'a' => \"ant\" is itself a stand-alone object:"

# ╔═╡ 44914181-9dd8-49b4-b2a1-fbbde543f620
begin
  pair = 'a' => "ant"

  first(pair)
end

# ╔═╡ 63badfd3-c501-4d7c-910b-f8a9a217eff2
md"## Functions"

# ╔═╡ e089aa5e-8d6e-48a4-aa36-d2e8546da46a
md"Three ways to define a generic function:"

# ╔═╡ c4c0e393-0c62-4e8e-88a1-a2bf91434413
begin
  foo(x) = x^2 # METHOD 1 (inline)
  foo(3)
end

# ╔═╡ 115b6c86-a0f2-4546-a1cb-54e34396a855
md"or"

# ╔═╡ 3c8e1237-ff55-4dd0-8037-0de5806e1a54
3 |> foo

# ╔═╡ 115b6c86-a0f2-4546-9961-a09632c33fb0
md"or"

# ╔═╡ 169c078b-2f81-4b5e-b7cc-46f4ef988c68
3 |> x -> x^2 # METHOD 2 (anonymous)

# ╔═╡ 115b6c86-a0f2-4546-90fb-7f2721f6fcc7
md"or"

# ╔═╡ 1d690f10-1e44-4e9a-af65-ccd25ecb9818
begin
  function foo2(x) # METHOD 3 (verbose)
      y = x
      z = y
      w = z
      return w^2
  end

  foo2(3)
end

# ╔═╡ 7a568a60-5212-4c87-8890-979691212d9b
md"## Basic iteration"

# ╔═╡ 77d49aef-9644-41c2-a6fa-eaf14df5d44b
md"Here are three ways to square the integers from 1 to 10."

# ╔═╡ c3d721a1-ae10-4c0b-8025-5084804a9f6c
md"METHOD 1 (explicit loop):"

# ╔═╡ d9016c54-9656-4684-9e90-09023d201062
begin
  squares = [] # or Int[] if performance matters
  for x in 1:10
      push!(squares, x^2)
  end

  squares
end

# ╔═╡ a707446f-b9f3-4405-b7bb-2f33ef765cc0
md"METHOD 2 (comprehension):"

# ╔═╡ f86cc305-d2b6-4be3-9625-01172c4a0081
[x^2 for x in 1:10]

# ╔═╡ a0c15cdb-6294-47bc-af4f-11c7de9e21dd
md"METHOD 3 (map):"

# ╔═╡ 7d2bb0cb-fe63-443b-8dba-241d9b744683
map(x -> x^2, 1:10)

# ╔═╡ 4455b804-b818-4fa0-a6e5-2f0b4dca5c59
md"METHOD 4 (broadcasting with dot syntax):"

# ╔═╡ 87188f88-6519-4ffd-8550-20498aa03ed0
(1:10) .^ 2

# ╔═╡ d1f615d1-644a-4358-9e7b-5b943cf6b560
md"## Random numbers"

# ╔═╡ da033108-0b83-4ee8-bce6-0a1379cc1259
typeof(2)

# ╔═╡ c494e441-4853-4a6d-9608-b5032c116833
rand() # sample a Float64 uniformly from interval [0, 1]

# ╔═╡ 2063cdaf-dbb9-4587-b473-4aa968e6937a
rand(3, 4) # do that 12 times and put in a 3 x 4 array

# ╔═╡ a2a2be72-f9fe-4e2c-8d9e-644a1b3cc6b6
randn(3, 4) # use normal distribution instead

# ╔═╡ 6deb462e-668e-475d-ac0a-23ded81445da
rand(Int8) # random elment of type Int8

# ╔═╡ 112438a5-ef2f-4c5d-8535-1a088a6a3228
rand(['a', 'b', 'c'], 10) # 10 random elements from a vector

# ╔═╡ b87a9e77-1e26-4a1f-a39f-7893c73eef39
md"Some standard libraries are needed to do more, for example:"

# ╔═╡ cc2371a4-26be-4325-bcca-51a27994a151
using Random

# ╔═╡ 9ac1f1af-bfdc-4499-9b33-09b4b6661170
randstring(30)

# ╔═╡ cf7971b2-b57f-4386-b45d-88d068bb0fa2
using Statistics  # part of standard library

# ╔═╡ 01f39450-9ce3-46c1-92c8-6db3a590d963
begin
  y = rand(30)
  @show mean(y) quantile(y, 0.75);
end

# ╔═╡ bd256b23-4c27-4ca0-abf5-96f457eb2bdf
md"""
(Use the macro @show before stuff you want printed prefixed by
*what* it is that is being printed.)
"""

# ╔═╡ 133e6f93-bb0b-460d-8a60-2e0214c059f5
md"""
For sampling from more general distributions we need
Distributions.jl package which is not part of the standard library.
"""

# ╔═╡ e7c2edc1-91cb-4ef1-a38a-beb5c7157b57
md"## Loading packages"

# ╔═╡ ab20d999-4cfa-4104-81f5-507803ea9ed6
md"If not in the REPL:"

# ╔═╡ c5e1e96d-691f-4366-9b21-6d4eb642d87e
begin
  using Pkg                        # built-in package manager
  Pkg.activate("env", shared=true) # create a new pkg env
end

# ╔═╡ f7215fcd-0f87-4b4a-b98c-3df1f31879bf
md"""
Add some packages to your environment (latest compatible versions
added by default):
"""

# ╔═╡ fdd4c037-11f7-4777-92b7-4ceba56e97ad
begin
  Pkg.add("Distributions")
  Pkg.add("Plots")
end

# ╔═╡ 64a59295-2630-40c2-b122-214de244406c
md"To load the code for use:"

# ╔═╡ b480add5-8651-4720-8a4b-c83694978e38
begin
  using Distributions
  using Plots

  N = 1000
  samples = rand(Normal(), N);   # equivalent to Julia's built-in `randn(d)`
  samples = (samples).^2;        # square element-wise
end

# ╔═╡ f618c762-a936-433f-a8b6-f7ff516dedc4
g = fit(Gamma, samples)

# ╔═╡ 5db01f65-3d0f-4ba4-81e1-74ef03c2e79e
@show mean(g) median(g) pdf(g, 1)

# ╔═╡ 2a96dfa7-cf5a-4e7f-a04e-0b76409c14a7
md"## Plotting"

# ╔═╡ c89ff154-d5a1-401f-b979-3458f2f26667
begin
  f(x) = pdf(g, x)

  xs = 0:0.1:4 # floats from 0 to 4 in steps of 0.1
  ys = f.(xs)  # apply f element-wise to xs

  plot(xs, ys, xrange=(0,4), yrange=(0,0.2))
  histogram!(samples , normalize=true, alpha=0.4)
end

# ╔═╡ ebe6467e-e767-481b-97e3-6677afc6ca9f
savefig("my_first_plot.png")

# ╔═╡ a7f061b8-d1ed-4b1f-b10e-562ce21caa04
md"# Exercises"

# ╔═╡ 19dae74b-24bd-428f-8f78-ab549ef1544e
md"## Exercise 1"

# ╔═╡ 3af8467f-4988-44dc-a8a4-a53f5149481e
md"Write a function named `total` that adds the elements of its vector input."

# ╔═╡ b04e2fdc-8f1a-4964-8704-dbab8e09b4f1
md"## Exercise 2"

# ╔═╡ 2b88d9c5-61a4-4e1f-a02f-27b6c05e4d02
md"""
Generate a 1000 random samples from the standard normal
distribution. Create a second such sample, and add the two samples
point-wise.  Compute the (sample) mean and variance of the combined
samples. In the same
plot, show a frequency-normalized histogram of the combined samples
and a plot of the pdf for normal distribution with zero mean and
variance `2`.
"""

# ╔═╡ b899ecf2-1e43-4783-bea4-3d467d48781c
md"""
You can use `std` to compute the standard deviation and `sqrt` to
compute square roots.
"""

# ╔═╡ 4f939d4a-e802-4b7d-97c5-2979af8a507d
md"## Exercise 3"

# ╔═╡ 63f410ea-c37a-4761-b639-63c76c72c513
md"The following shows that named tuples share some behaviour with dictionaries:"

# ╔═╡ cb2c42bf-21b5-4e04-8f5a-e1ee9eb5c15c
begin
  t = (x = 1, y = "cat", z = 4.5)
  keys(t)
end

# ╔═╡ b8996fa3-c046-4fb1-adc5-e671db8bca5d
t[:y]

# ╔═╡ 6584c100-ffd6-4fa1-86f1-0ea38de21abb
md"""
Write a function called `dict` that converts a named tuple to an
actual dictionary. You can create an empty dictionary using `Dict()`.
"""

# ╔═╡ 135dac9b-0bd9-4e1d-a55c-0cb44ab816d7
md"""
---

*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
"""

# ╔═╡ Cell order:
# ╟─142b62e1-364e-45d6-9bca-7c69b794f8ce
# ╟─27f5936d-4043-4cb6-b795-033f6f2a0674
# ╟─3f08e241-e53d-4d65-935f-ddf6a6bfbbdd
# ╟─3e3bcc53-7840-4af7-af2a-25235e544a31
# ╟─96774cd2-130c-4455-8af5-148d95ea2900
# ╠═8e8a0cbc-96b9-490c-a6bf-47194d7e8e12
# ╠═a65fb83a-f977-47ce-828a-4cc485149963
# ╠═dc16182e-b6b3-47ee-9e54-3f7b3ca87ecb
# ╟─bc1912b8-c66a-4ab5-ba1f-6363f43ec697
# ╠═048a1195-1107-46e9-95ea-2955abd45275
# ╟─ef72056a-9321-4290-b1b3-dd47e367ba54
# ╠═af710d1b-f8c2-4219-8d7e-8e449afd1c48
# ╟─3e690f1c-4e80-4016-a949-7f3bd292fe31
# ╠═b541c5ac-3c15-4e71-8514-99938a2932db
# ╟─c00b48c5-d35b-409a-a0de-1721c1bc2df2
# ╟─3694faf5-d9a4-4089-bca8-7eec7950fd82
# ╠═cad4d998-43d7-421b-9873-1b1430e635cc
# ╠═5c43b0f6-30f4-4ef7-b43e-d0ebe87da176
# ╟─26ce1939-e595-496a-9009-4cfb2012998f
# ╠═d66140c7-958b-4347-abd3-cb77d7a79683
# ╟─7aee99b6-3c6e-4a71-879e-dc128f3db7b3
# ╠═96c1adb9-1b29-433d-a367-9aa3c6cf34d0
# ╠═ef6a776c-8c63-4d30-bf32-94657e65284e
# ╠═4bb2001b-3b46-40f3-9afd-65a5b5facac9
# ╟─b97d12e6-2346-410b-b6c8-54f5ed90a964
# ╠═177a30bb-473a-4e5e-9292-fe822525fc77
# ╠═109564cb-ea6f-4264-ae5e-4b83dcbc9675
# ╠═26714555-a4ad-49cb-8a29-32f11452654a
# ╟─d00b2507-e891-4ec5-a5f1-151fcbe229a2
# ╠═4d1dd37a-194e-4ed2-81bc-2d5603785a09
# ╠═38477997-58dd-4248-9d87-1d1fbb0e3997
# ╠═e61c7ae2-8204-4a20-b950-fddef2a1fb10
# ╟─d38b7040-bb41-44be-951b-64cb2a36c8e3
# ╟─9f1e3326-b738-4b60-b0e6-441461cc8770
# ╟─494f9e7c-e53a-4c65-8cb0-de601961bc02
# ╠═259bf65c-6868-47b0-a87d-950c50fb2955
# ╠═6ef57fdd-3188-4b19-8448-4bcb089096cd
# ╠═ef2473e2-bb35-4bd7-a012-a84240254fb6
# ╠═76dc6ac8-03fe-4c6e-bbdd-9799f7bb2e60
# ╟─5b3a3587-31f7-4a00-97a7-88e3af4f10ee
# ╟─82f84505-336f-4845-b372-8c9866e51852
# ╠═5b71f2a1-5567-4e59-8f3d-b82f1e7b6f7a
# ╠═249c3a08-3ca6-463e-ab04-579e5608ae53
# ╟─f85e0bbc-9db0-4a9f-86cf-35420d9e6995
# ╠═b8b91bc0-0a63-4695-a29a-8d0445351914
# ╠═e0e4f0da-69c7-4941-be65-28d0fcca50a8
# ╟─166a1850-cddb-4f72-9a2e-c468345d87d1
# ╠═44ef9d71-f75b-4140-b5fa-cb79ebf595ef
# ╠═ff4b5618-4be7-4111-91c5-524da38aa391
# ╟─fb16016a-87d7-4e9f-ad90-b4405b216771
# ╠═20ba7245-15ed-417e-895b-997a92b731e0
# ╠═542ec737-2c88-4293-a525-e9b04a4bd246
# ╠═8b1e2b91-947c-4046-80f0-ba7781e173cf
# ╠═3dcab89a-20ca-462a-9cbc-9a2c39793333
# ╟─65eb310b-82e8-4a1d-b887-2b16710e5502
# ╠═35998a8b-6fda-4a2f-9452-5298a8a4a401
# ╟─2db5e586-237b-45aa-b017-96f1602f2cad
# ╟─8b1c08d9-f8dd-4ec9-8be2-2a4017c45345
# ╠═4942a061-18dd-41d0-a7ad-80d84f5b0070
# ╠═d1fb9147-4a06-44fa-8378-5aa686f0b407
# ╠═4f4bf27c-b2f9-4c15-9f43-49763e8691a1
# ╠═29d71401-66e5-484e-bb0e-3a45761c733a
# ╟─4b80c58d-646a-405b-9979-7044b2f2df33
# ╠═44914181-9dd8-49b4-b2a1-fbbde543f620
# ╟─63badfd3-c501-4d7c-910b-f8a9a217eff2
# ╟─e089aa5e-8d6e-48a4-aa36-d2e8546da46a
# ╠═c4c0e393-0c62-4e8e-88a1-a2bf91434413
# ╟─115b6c86-a0f2-4546-a1cb-54e34396a855
# ╠═3c8e1237-ff55-4dd0-8037-0de5806e1a54
# ╟─115b6c86-a0f2-4546-9961-a09632c33fb0
# ╠═169c078b-2f81-4b5e-b7cc-46f4ef988c68
# ╟─115b6c86-a0f2-4546-90fb-7f2721f6fcc7
# ╠═1d690f10-1e44-4e9a-af65-ccd25ecb9818
# ╟─7a568a60-5212-4c87-8890-979691212d9b
# ╟─77d49aef-9644-41c2-a6fa-eaf14df5d44b
# ╟─c3d721a1-ae10-4c0b-8025-5084804a9f6c
# ╠═d9016c54-9656-4684-9e90-09023d201062
# ╟─a707446f-b9f3-4405-b7bb-2f33ef765cc0
# ╠═f86cc305-d2b6-4be3-9625-01172c4a0081
# ╟─a0c15cdb-6294-47bc-af4f-11c7de9e21dd
# ╠═7d2bb0cb-fe63-443b-8dba-241d9b744683
# ╟─4455b804-b818-4fa0-a6e5-2f0b4dca5c59
# ╠═87188f88-6519-4ffd-8550-20498aa03ed0
# ╟─d1f615d1-644a-4358-9e7b-5b943cf6b560
# ╠═da033108-0b83-4ee8-bce6-0a1379cc1259
# ╠═c494e441-4853-4a6d-9608-b5032c116833
# ╠═2063cdaf-dbb9-4587-b473-4aa968e6937a
# ╠═a2a2be72-f9fe-4e2c-8d9e-644a1b3cc6b6
# ╠═6deb462e-668e-475d-ac0a-23ded81445da
# ╠═112438a5-ef2f-4c5d-8535-1a088a6a3228
# ╟─b87a9e77-1e26-4a1f-a39f-7893c73eef39
# ╠═cc2371a4-26be-4325-bcca-51a27994a151
# ╠═9ac1f1af-bfdc-4499-9b33-09b4b6661170
# ╠═cf7971b2-b57f-4386-b45d-88d068bb0fa2
# ╠═01f39450-9ce3-46c1-92c8-6db3a590d963
# ╟─bd256b23-4c27-4ca0-abf5-96f457eb2bdf
# ╟─133e6f93-bb0b-460d-8a60-2e0214c059f5
# ╟─e7c2edc1-91cb-4ef1-a38a-beb5c7157b57
# ╟─ab20d999-4cfa-4104-81f5-507803ea9ed6
# ╠═c5e1e96d-691f-4366-9b21-6d4eb642d87e
# ╟─f7215fcd-0f87-4b4a-b98c-3df1f31879bf
# ╠═fdd4c037-11f7-4777-92b7-4ceba56e97ad
# ╟─64a59295-2630-40c2-b122-214de244406c
# ╠═b480add5-8651-4720-8a4b-c83694978e38
# ╠═f618c762-a936-433f-a8b6-f7ff516dedc4
# ╠═5db01f65-3d0f-4ba4-81e1-74ef03c2e79e
# ╟─2a96dfa7-cf5a-4e7f-a04e-0b76409c14a7
# ╠═c89ff154-d5a1-401f-b979-3458f2f26667
# ╠═ebe6467e-e767-481b-97e3-6677afc6ca9f
# ╟─a7f061b8-d1ed-4b1f-b10e-562ce21caa04
# ╟─19dae74b-24bd-428f-8f78-ab549ef1544e
# ╟─3af8467f-4988-44dc-a8a4-a53f5149481e
# ╟─b04e2fdc-8f1a-4964-8704-dbab8e09b4f1
# ╟─2b88d9c5-61a4-4e1f-a02f-27b6c05e4d02
# ╟─b899ecf2-1e43-4783-bea4-3d467d48781c
# ╟─4f939d4a-e802-4b7d-97c5-2979af8a507d
# ╟─63f410ea-c37a-4761-b639-63c76c72c513
# ╠═cb2c42bf-21b5-4e04-8f5a-e1ee9eb5c15c
# ╠═b8996fa3-c046-4fb1-adc5-e671db8bca5d
# ╟─6584c100-ffd6-4fa1-86f1-0ea38de21abb
# ╟─135dac9b-0bd9-4e1d-a55c-0cb44ab816d7


text:
"""
    initial(::Type{RandomInitialLayout}, N::T) where {T <: EcologicalNetworks.AbstractEcologicalNetwork}

Random disposition of nodes in a circle. This is a good starting point for any
force-directed layout. The circle is scaled so that its radius is twice the
square root of the network richness, which helps most layouts converge faster.
"""
function initial(::Type{RandomInitialLayout}, N::T) where {T <: EcologicalNetworks.AbstractEcologicalNetwork}
  L = Dict([s => NodePosition() for s in species(N)])
  _adj = 2sqrt(richness(N))
  for s in species(N)
    L[s].x *= _adj
    L[s].y *= _adj
  end
  return L
end

"""
    initial(::Type{BipartiteInitialLayout}, N::T) where {T <: EcologicalNetworks.AbstractBipartiteNetwork}

Random disposition of nodes on two levels for bipartite networks.
"""
function initial(::Type{BipartiteInitialLayout}, N::T) where {T <: EcologicalNetworks.AbstractBipartiteNetwork}
  level = NodePosition[]
  for (i, s) in enumerate(species(N))
    this_level = s ∈ species(N; dims=1) ? 1.0 : 0.0
    push!(level, NodePosition(rand(), this_level, 0.0, 0.0))
  end
  return Dict(zip(species(N), level))
end

"""
    initial(::Type{FoodwebInitialLayout}, N::T) where {T <: EcologicalNetworks.AbstractUnipartiteNetwork}

Random disposition of nodes on trophic levels for food webs. Note that the
continuous trophic level is used, but the layout can be modified afterwards to
use another measure of trophic rank.
"""
function initial(::Type{FoodwebInitialLayout}, N::T) where {T <: EcologicalNetworks.AbstractUnipartiteNetwork}
  L = initial(RandomInitialLayout, N)
  tl = trophic_level(N)
  for s in species(N)
    L[s].y = tl[s]
  end
  return L
end

"""
    initial(::Type{CircularInitialLayout}, N::T) where {T <: EcologicalNetworks.AbstractEcologicalNetwork}

Random disposition of nodes on a circle. This is the starting point for
circle-based layouts.
"""
function initial(::Type{CircularInitialLayout}, N::T) where {T <: EcologicalNetworks.AbstractEcologicalNetwork}
  level = NodePosition[]
  n = richness(N)
  for (i, s) in enumerate(species(N))
    θ = 2i * π/n
    x, y = cos(θ), sin(θ)
    push!(level, NodePosition(x, y, i))
  end
  return Dict(zip(species(N), level))
end

"""
    initial(::Type{UnravelledInitialLayout}, N::T) where {T <: EcologicalNetworks.AbstractUnipartiteNetwork}

Unravelled disposition of nodes on trophic levels for food webs, where the x
axis is the omnivory index. Note that the *fractional* trophic level is used,
but the layout can be modified afterwards to use the continuous levels. See the
documentation for `UnravelledLayout` to see how.
"""
function initial(::Type{UnravelledInitialLayout}, N::T) where {T <: EcologicalNetworks.AbstractUnipartiteNetwork}
  layout = Dict([s => NodePosition() for s in species(N)])
  tl = fractional_trophic_level(N)
  oi = omnivory(N)
  for s in species(N)
    layout[s].x = float(oi[s])
    layout[s].y = float(tl[s])
  end
  return layout
end


text:
using Base.Test
using EGR
using Redis

println("TestFindBest")

client = RedisConnection()

# run(`redis-cli keys "Test*"` |> `xargs redis-cli del`)
run(pipeline(`redis-cli keys "Test*"`, `xargs redis-cli del`))

createOracleOutputLevel = 1
numEquivalentPasses = 1
algOutputLevel = 0
maxOutputNum=20
constStepSize(k)=1

sds = [
(n,ntp,dt)->uQuadBeta1(1.0,n,1,false),
(n,ntp,dt)->uExpBeta1(2.0,n,1,false),
(n,ntp,dt)->SGsd( "SG")
]

myWriteFunction(problem, sd, opts, k, gnum, origWant, fromOutputsFunction) = writeFunction(client, problem,  opts, sd,k, gnum, origWant, fromOutputsFunction; outputLevel = 0)

for (gradientOracle, numVars, numTP, csDataType, LossFunctionString, myOutputter, L2reg, thisDataName, thisProblem) in Oracles
	println(" $thisDataName $LossFunctionString L2reg = $L2reg")
	maxG  = Int(round(numEquivalentPasses*numTP))

	expIndices=unique(round(Int, logspace(0,log10(maxG+1),maxOutputNum)))-1

	myOutputOpts =  OutputOpts(myOutputter,expIndices)
	myOpts(stepSizePower) = Opts(zeros(numVars,1);  stepSizePower=stepSizePower, maxG=maxG, outputLevel=algOutputLevel)
	for sd in sds

		thisSD = sd(numVars,numTP,csDataType)

		println("  stepString = $(thisSD.stepString)")

		algForSearch(stepSizePower) =alg(thisProblem(Task(() -> getSequentialFinite(numTP, gradientOracle))), myOpts(stepSizePower),thisSD, myOutputOpts, myWriteFunction)

		returnResultIfExists(stepSizePower) = returnIfExists(client,thisProblem(Task(() -> getSequentialFinite(numTP, gradientOracle))), myOpts(stepSizePower), thisSD,expIndices,myOutputOpts.outputter.numOutputsFromOutputsFunction)

		findBests = [((res)->getF(res,maxG),0.0,"getF"),((res)->getPCC(res,maxG),-1.0,"getPCC"), ((res)->getMCC(res,maxG),-1.0,"getMCC")]

		for findBest in findBests
			println("   findBest = $(findBest[3])")
			(bestPower, bestVal) = findBestStepsizeFactor(algForSearch,findBest...,returnResultIfExists;outputLevel = 0)

			println("   Best stepsize power = $bestPower, best value = $bestVal")
		end
	end
end

println("TestFindBest successful!")
println()