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

using Markdown
using InteractiveUtils

# ╔═╡ 2bf8e4c2-339b-11eb-1d49-7f1b210a650a
using DelimitedFiles

# ╔═╡ 3e0c149a-3350-11eb-0b17-7bea5d1a28d4
begin
	DAY = 2
	md"# 2020 Day $DAY"
end

# ╔═╡ 18080fe2-339b-11eb-34f0-f393638506fc
md"## Part 1"

# ╔═╡ 4e084a4e-339b-11eb-3ffc-41118215f226
md"*A password is is valid if the character occurs between the lower and higher number of times.*"

# ╔═╡ 24e1613c-339b-11eb-2976-59884053d1ee
input = readdlm("../data/day_$DAY.txt", '\n');

# ╔═╡ 0dcf6e94-346d-11eb-10f4-bb7c5248522b
test_input = split("""1-3 a: abcde
1-3 b: cdefg
2-9 c: ccccccccc""","\n")

# ╔═╡ 94a3659c-346d-11eb-1854-971ac36a57c7
row = test_input[1]

# ╔═╡ 7f6eb578-346d-11eb-050d-8396d5ae6c09
function parse_row(row)
	rest, password = split(row,": ")
	l, rest = split(rest,"-")
	u, char = split(rest," ")
	return parse(Int, l), parse(Int, u), char[1], password
end

# ╔═╡ 9876725e-346d-11eb-3b79-b3f9973d57ca
parse_row(row)

# ╔═╡ 028a6038-346e-11eb-057c-c9b6a9e9a12e
function is_valid_password_part_1(row)
	l, u, char, password = parse_row(row)
	n_chars = sum([char == c for c in password])
	if l ≤ n_chars ≤ u
		return true
	else
		return false
	end
end

# ╔═╡ 2fa28998-346f-11eb-0780-bf2d8ac8f835
part_1(input) = sum(is_valid_password_part_1.(input))

# ╔═╡ 5e2f24aa-346e-11eb-0429-610e74921e6f
@assert part_1(test_input) == 2

# ╔═╡ 438f54ca-346f-11eb-1103-ff18132314b4
part_1(input)

# ╔═╡ 5f797e22-339d-11eb-1ff3-dda90f946838
md"## Part 2"

# ╔═╡ e3177eae-339e-11eb-3f9f-1f7230976832
md"*Now the numbers are indexes and exactly one must match the required character*"

# ╔═╡ 85fa9e14-346f-11eb-3369-7f4fda1d3da2
function is_valid_password_part_2(row)
	l, u, char, password = parse_row(row)
	return (password[l] == char) ⊻ (password[u] == char)
end

# ╔═╡ bebc8de0-346f-11eb-3977-8745fed75644
part_2(input) = sum(is_valid_password_part_2.(input))

# ╔═╡ f1e5075c-346f-11eb-197a-ad88701d918f
@assert part_2(test_input) == 1

# ╔═╡ c6a335e8-346f-11eb-1b2f-ad88443694b6
part_2(input)

# ╔═╡ Cell order:
# ╠═3e0c149a-3350-11eb-0b17-7bea5d1a28d4
# ╠═2bf8e4c2-339b-11eb-1d49-7f1b210a650a
# ╟─18080fe2-339b-11eb-34f0-f393638506fc
# ╟─4e084a4e-339b-11eb-3ffc-41118215f226
# ╠═24e1613c-339b-11eb-2976-59884053d1ee
# ╠═0dcf6e94-346d-11eb-10f4-bb7c5248522b
# ╠═94a3659c-346d-11eb-1854-971ac36a57c7
# ╠═9876725e-346d-11eb-3b79-b3f9973d57ca
# ╠═7f6eb578-346d-11eb-050d-8396d5ae6c09
# ╠═028a6038-346e-11eb-057c-c9b6a9e9a12e
# ╠═2fa28998-346f-11eb-0780-bf2d8ac8f835
# ╠═5e2f24aa-346e-11eb-0429-610e74921e6f
# ╠═438f54ca-346f-11eb-1103-ff18132314b4
# ╟─5f797e22-339d-11eb-1ff3-dda90f946838
# ╠═e3177eae-339e-11eb-3f9f-1f7230976832
# ╠═85fa9e14-346f-11eb-3369-7f4fda1d3da2
# ╠═bebc8de0-346f-11eb-3977-8745fed75644
# ╠═f1e5075c-346f-11eb-197a-ad88701d918f
# ╠═c6a335e8-346f-11eb-1b2f-ad88443694b6

max_stars_repo_name: RensDimmendaal/advent-of-code-julia

max_issues_repo_name: RensDimmendaal/advent-of-code-julia

max_forks_repo_name: RensDimmendaal/advent-of-code-julia


text:
module ImageCoreTests

using ImageCore
using Test, ReferenceTests
using Aqua, Documenter # for meta quality checks

@testset "Project meta quality checks" begin
    # Not checking compat section for test-only dependencies
    Aqua.test_ambiguities(ImageCore)
    Aqua.test_all(ImageCore;
                  ambiguities=false,
                  project_extras=true,
                  deps_compat=true,
                  stale_deps=true,
                  project_toml_formatting=true,
                  unbound_args=false, # FIXME: it fails when this is true
    )
    DocMeta.setdocmeta!(ImageCore, :DocTestSetup, :(using ImageCore); recursive=true)
end

# ReferenceTests uses ImageInTerminal as a default image rendering backend, we need to
# temporarily disable it when we do doctest.
# TODO: ReferenceTests doesn't yet support this switch. That's why we need ImageInTerminal dependency.
using ImageInTerminal
ImageInTerminal.disable_encoding()
doctest(ImageCore, manual = false)
ImageInTerminal.enable_encoding()

include("colorchannels.jl")
include("views.jl")
include("convert_reinterpret.jl")
include("traits.jl")
include("map.jl")
include("functions.jl")
include("show.jl")

# To ensure our deprecations work and don't break code
include("deprecations.jl")

# run these last
isCI = haskey(ENV, "CI") || get(ENV, "JULIA_PKGEVAL", false)
if Base.JLOptions().can_inline == 1 && !isCI
    @info "running benchmarks"
    include("benchmarks.jl")  # these fail if inlining is off
end

end

ext: jl


text:
using Distributed
addprocs(4)

# using Revise

##

# for drawing Bayes tree with images (see debug tricks below)
# using Cairo, Fontconfig
# using Gadfly

using Dates
using RoME
@everywhere using RoME

#  Do some plotting
using RoMEPlotting
Gadfly.set_default_plot_size(35cm,25cm)

function driveHex(fgl::G,
                  posecount::Int;
                  graphinit::Bool=false ) where G <: AbstractDFG
    #
    # Drive around in a hexagon
    for i in (posecount-1):(posecount-1+5)
        psym = Symbol("x$i")
        posecount += 1
        nsym = Symbol("x$(i+1)")
        addVariable!(fgl, nsym, Pose2)
        pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
        addFactor!(fgl, [psym;nsym], pp, autoinit=graphinit )
    end

    return posecount
end

function offsetHexLeg(fgl::G,
                      posecount::Int; direction=:right,
                      graphinit::Bool=false ) where G <: AbstractDFG
    #
    psym = Symbol("x$(posecount-1)")
    nsym = Symbol("x$(posecount)")
    posecount += 1
    addVariable!(fgl, nsym, Pose2)
    pp = nothing
    if direction == :right
        pp = Pose2Pose2(MvNormal([10.0;0;-pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    elseif direction == :left
        pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    end
    addFactor!(fgl, [psym; nsym], pp, autoinit=graphinit )
    return posecount
end

## start with an empty factor graph object
fg = initfg()
getSolverParams(fg).isfixedlag = true
getSolverParams(fg).qfl = 20
posecount = 0

# Add the first pose :x0
addVariable!(fg, :x0, Pose2)
posecount += 1

# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor!(fg, [:x0], PriorPose2( MvNormal([0.0; 0.0; 0.0],
                                           Matrix(Diagonal([0.1;0.1;0.05].^2))) ), autoinit=false )

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, labels=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x0; :l1], p2br, autoinit=false )

## hex 1

posecount = driveHex(fg, posecount)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [Symbol("x$(posecount-1)"); :l1], p2br2, autoinit=false )

printGraphSummary(fg)

# draw figures for debugging
getSolverParams(fg).drawtree = true
# getSolverParams(fg).showtree = false

# burn out USB drive instead of solid storage
# getSolverParams(fg).logpath="/media/dehann/temp/caesar/$(now())"

# debugging options
getSolverParams(fg).multiproc = true
# getSolverParams(fg).async = true
# getSolverParams(fg).dbg = true
# getSolverParams(fg).limititers = 50
# getSolverParams(fg).downsolve = false

tree, smt, hist = solveTree!(fg) # , recordcliqs=ls(fg))

printGraphSummary(fg)

pl = plotBeehive_6(fg, meanmax=:mean)
# pl = plotBeehive_6(fg, meanmax=:max)
# drawPosesLandms(fg, meanmax=:max) |> SVG("/tmp/test.svg"); @async run(`eog /tmp/test.svg`)

## fun additions
# fetchCliqTaskHistoryAll!(smt, hist)
# assignTreeHistory!(tree, hist)
# printCliqHistorySummary(hist, tree, :l0)
# makeCsmMovie(fg, tree)

## hex 2

posecount = offsetHexLeg(fg, posecount, direction=:right)

# Add landmarks with Bearing range measurements
addVariable!(fg, :l2, Point2, labels=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [Symbol("x$(posecount-1)"); :l2], p2br, autoinit=false )

posecount = driveHex(fg, posecount)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [Symbol("x$(posecount-1)"); :l2], p2br2, autoinit=false )

# writeGraphPdf(fg,show=true)

# getSolverParams(fg).downsolve = false

tree, smt, hist = solveTree!(fg, tree) #, recordcliqs=ls(fg))

# drawPosesLandms(fg, meanmax=:max) |> SVG("/tmp/test.svg") # || @async run(`eog /tmp/test.svg`)
pl = plotBeehive_6(fg, meanmax=:mean)

# new sighting on :l0

addVariable!(fg, :l0, Point2, labels=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x5; :l0], p2br, autoinit=false )

p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x12; :l0], p2br2, autoinit=false )

# and here it goes totally wacko

# getSolverParams(fg).dbg = true

tree, smt, hist = solveTree!(fg, tree) #, recordcliqs=ls(fg))

pl = plotBeehive_6(fg, meanmax=:mean)

# getSolverParams(fg).dbg = false

## hex 3

posecount = offsetHexLeg(fg, posecount, direction=:right, graphinit=true)

posecount = driveHex(fg, posecount, graphinit=true)

# getSolverParams(fg).downsolve=false
tree, smt, hist = solveTree!(fg, tree) #, recordcliqs=ls(fg))

pl = plotBeehive_6(fg, meanmax=:mean)

p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x19; :l0], p2br2, autoinit=false )

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x18; :l1], p2br2, autoinit=false )

# Add landmarks with Bearing range measurements
addVariable!(fg, :l3, Point2, labels=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x14; :l3], p2br, autoinit=false )

p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x20; :l3], p2br, autoinit=false )

# getSolverParams(fg).dbg = true

tree, smt, hist = solveTree!(fg, tree)

pl = plotBeehive_6(fg, meanmax=:mean)
# drawPosesLandms(fg, meanmax=:max) |> SVG("/tmp/test.svg") # || @async run(`eog /tmp/test.svg`)

## hex 4

posecount = offsetHexLeg(fg, posecount, direction=:right)

posecount = driveHex(fg, posecount)

# Add landmarks with Bearing range measurements
addVariable!(fg, :l4, Point2, labels=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x21; :l4], p2br, autoinit=false )

p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x27; :l4], p2br, autoinit=false )

# writeGraphPdf(fg)

tree, smt, hist = solveTree!(fg, tree)

pl = plotBeehive_6(fg, meanmax=:mean)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x26; :l0], p2br2, autoinit=false )

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x25; :l2], p2br2, autoinit=false )

tree, smt, hist = solveTree!(fg, tree, recordcliqs=ls(fg))

pl = plotBeehive_6(fg, meanmax=:mean)
# drawPosesLandms(fg, meanmax=:max) |> SVG("/tmp/test.svg") # || @async run(`eog /tmp/test.svg`)

## hex 5

posecount = offsetHexLeg(fg, posecount, direction=:right)

# Add landmarks with Bearing range measurements
addVariable!(fg, :l5, Point2, labels=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x28; :l5], p2br, autoinit=false )

p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x4; :l5], p2br, autoinit=false )

tree, smt, hist = solveTree!(fg, tree)

pl = plotBeehive_6(fg, meanmax=:mean)

## Weird initialization issues happen here

posecount = driveHex(fg, posecount, graphinit=true)

# # Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x34; :l5], p2br2, autoinit=false )

p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x32; :l3], p2br2, autoinit=false )

p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x33; :l0], p2br2, autoinit=false )

# writeGraphPdf(fg)

tree, smt, hist = solveTree!(fg, tree)

pl = plotBeehive_6(fg, meanmax=:mean)

## hex 6

posecount = offsetHexLeg(fg, posecount, direction=:right)

posecount = driveHex(fg, posecount, graphinit=true)

# l6 from x35, x41, x11
# Add landmarks with Bearing range measurements
addVariable!(fg, :l6, Point2, labels=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x35; :l6], p2br, autoinit=false )

p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x41; :l6], p2br, autoinit=false )

p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x11; :l6], p2br, autoinit=false )

# # Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x40; :l0], p2br2, autoinit=false )

p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x39; :l4], p2br2, autoinit=false )

# writeGraphPdf(fg)

tree, smt, hist = solveTree!(fg, tree)

pl = plotBeehive_6(fg, meanmax=:mean)

## complete solve

# dontMarginalizeVariablesAll!(fg)
#
# # getSolverParams(fg).dbg = false
# # getSolverParams(fg).multiproc = false
#
# tree, smt, hist = solveTree!(fg) # , recordcliqs=ls(fg)
#
# pl = plotBeehive_6(fg, meanmax=:mean)
#
#
# # plotLocalProduct(fg, :l6)
# getSolverParams(fg).isfixedlag = true
# getSolverParams(fg).qfl = 20

## hex 7

posecount = offsetHexLeg(fg, posecount, direction=:left)
posecount = offsetHexLeg(fg, posecount, direction=:right)

posecount = driveHex(fg, posecount, graphinit=true)

#
# # # Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x47; :l5], p2br2, autoinit=false )
#
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x48; :l6], p2br2, autoinit=false )

tree, smt, hist = solveTree!(fg, tree)

pl = plotBeehive_6(fg, meanmax=:mean)
# resetVariableAllInitializations!(fg)

# ## manual
# pts = predictbelief(fg, :x49, [:x48x49f1;])
# setValKDE!(fg, :x49, manikde!(pts[1], Pose2().manifolds), false, pts[2])

## hex 8

posecount = offsetHexLeg(fg, posecount, direction=:right)

posecount = driveHex(fg, posecount)

# # Add landmarks with Bearing range measurements
# p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
# addFactor!(fg, [Symbol("x$(posecount-1)"); :l1], p2br2, autoinit=false )

# writeGraphPdf(fg)

tree, smt, hist = solveTree!(fg, tree)

## hex 9

posecount = offsetHexLeg(fg, posecount, direction=:right)

posecount = driveHex(fg, posecount)

# # Add landmarks with Bearing range measurements
# p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
# addFactor!(fg, [Symbol("x$(posecount-1)"); :l1], p2br2, autoinit=false )

# writeGraphPdf(fg)

tree, smt, hist = solveTree!(fg, tree)

## hex 10

posecount = offsetHexLeg(fg, posecount, direction=:left)
posecount = offsetHexLeg(fg, posecount, direction=:right)

posecount = driveHex(fg, posecount)

# # Add landmarks with Bearing range measurements
# p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
# addFactor!(fg, [Symbol("x$(posecount-1)"); :l1], p2br2, autoinit=false )

writeGraphPdf(fg)

tree, smt, hist = solveTree!(fg, tree)

drawPosesLandms(fg, meanmax=:max) |> SVG("/tmp/test.svg") || @async run(`eog /tmp/test.svg`)

## additional dev code

# from hex1

getSolverParams(fg).upsolve = false
getSolverParams(fg).downsolve = true

stuff = solveCliq!( fg, tree, :x3 )
stuff = solveCliq!( fg, tree, :x1 )
stuff = solveCliq!( fg, tree, :x6 )
stuff = solveCliq!( fg, tree, :x4 )
stuff = solveCliq!( fg, tree, :x2 )
stuff = solveCliq!( fg, tree, :x0 )

fg2bd = loadDFG("/tmp/caesar/cliqSubFgs/cliq2/fg_beforedownsolve", Main)

writeGraphPdf(fg2bd, show=true)

## dev auto init after hex 3

printCliqHistorySummary(tree, :x16)

cliq = getCliq(tree, :x16)
cid = cliq.index
ste = 15
sfg = hist[cid][ste][4].cliqSubFg
cliqh = hist[cid][ste][4].cliq

printCliqSummary(sfg, cliq)

pl = plotBeehive_6(sfg)

stuff = sandboxCliqResolveStep(tree, :x16, 15)

cond = getSolveCondition(cliqh)
notify(cond)

cond.waitq = Any[]

isready(cond)

wait(getSolveCondition(cliqh))

notify(getSolveCondition(cliq))
IIF.notifyCliqUpInitStatus!(cliq, :downsolved)

getCliqStatus(cliq)

writeGraphPdf(sfg, show=true)

pl = plotBeehive_6(fg, to=15)

# leaf cliq

printCliqHistorySummary(tree, :l0)

cliq = getCliq(tree, :l0)
cid = cliq.index
ste = 9
sfg = hist[cid][ste][4].cliqSubFg
cliqh = hist[cid][ste][4].cliq

printCliqSummary(sfg, cliqh)

writeGraphPdf(sfg, show=true)

plotBeehive_6(sfg)

ls(sfg, :x19)

pts = approxConv(sfg, :x19l0f1, :x19)

X19 = manikde!(pts, Pose2().manifolds)

plotPose(Pose2(), X19)

pts[3,:] = 2*pi*rand(100).-pi
setValKDE!(getVariable(sfg, :x19), X19)

plotTreeProductDown(fg,tree,:l2, :l2)

## debug down on hex1, x4

drawTree(tree)

printCliqHistorySummary(tree,:x4)

csmcStep19 = getCliqSolveHistory(tree,:x4)[19][4]

drawGraph(csmcStep19.cliqSubFg, show=true)

getSolverParams(csmcStep19.dfg).dbg = true

stuff = sandboxCliqResolveStep(tree,:x4,19)

plotKDE(map(x->getKDE(x, :x4),[fg; csmcStep19.cliqSubFg; stuff[4].dfg]), levels=2 )

drawPosesLandms(fg,meanmax=:mean)

plotKDE([getKDE(csmcStep19.cliqSubFg, :x3);getKDE(csmcStep19.cliqSubFg, :x4); getKDE(csmcStep19.cliqSubFg, :x5)], levels=1)

plotLocalProduct(csmcStep19.cliqSubFg)

# what are down messages into (:x4|)

plotCliqDownMsgs(tree,:x3, existing=pl, dims=[1;2]);

plotTreeProductDown(csmcStep19.cliqSubFg, csmcStep19.tree,:x4,:x5)

## what is going on with x3 and x5 in root

printCliqHistorySummary(tree, :x3)
csmc = getCliqSolveHistory(tree, :x3)[9][4]

getSolverParams(csmc.dfg).dbg = true

stuff = sandboxCliqResolveStep(tree,:x3,9)

plotTreeUpMsgs(fg, tree, :x5) #plotCliqUpMsgs

plotCliqDownMsgs(tree,:x3, existing=pl, dims=[1;2]);

plotTreeProductUp(fg, tree, :x5)

hist = getCliqSolveHistory(tree, :x6)

dfg = hist[11][4].dfg
subfg = hist[11][4].cliqSubFg
cliq = hist[11][4].cliq

printCliqHistorySummary(hist)

## DEBUG x19

plotKDE(fg, :x19, levels=3)

getCliq(tree, :x19)

printCliqHistorySummary(tree, :x19)

sfg_ad = hist[4][12][4].cliqSubFg

drawGraph(sfg_ad, show=true)

plotLocalProduct(sfg_ad, :x19)

getVariableInferredDim(fg, :x19)
getVariableInferredDim(fg, :l0)

## END DEBUG

"""
    $SIGNATURES

Related

plotLocalProduct
"""
function plotLocalInverseProduct()

end

plotKDE(fg, :x41, levels=2)

plotLocalProduct(fg, :x41)

## DEBUG weird hex5

dontMarginalizeVariablesAll!(fg)
# getSolverParams(fg).isfixedlag = false
# getSolverParams(fg).qfl = 99999999

getSolverParams(fg).dbg = true
# getSolverParams(fg).multiproc = false

tree, smt, hist = solveTree!(fg, recordcliqs=ls(fg))

pl = plotBeehive_6(fg, meanmax=:mean)

plotLocalProduct(fg, :l5)

stuff = localProduct(fg,:l5)

setValKDE!(fg, :l5, stuff[1], false, stuff[5])

ls(fg, :x33)

saveDFG(fg, "/home/dehann/Documents/beehive_hex5_fail_9_13")

printCliqHistorySummary(tree,:l0)

l0bdfg = getCliqSolveHistory(tree,:l0)[11][4].dfg
l0adfg = getCliqSolveHistory(tree,:l0)[12][4].dfg

plotKDE(l0bdfg, :l0, levels=3)
plotKDE(l0adfg, :l0, levels=3)

getPoints(getKDE(l0bdfg, :l0))

getPoints(getKDE(l0adfg, :l0))

## DEBUG END

## DEBUG post 6

pl = plotKDE(fg, ls(fg, r"l"), levels=1)
pl = plotKDE(fg, ls(fg, r"x"), levels=1)

dontMarginalizeVariablesAll!(fg)
getSolverParams(fg).dbg=true
tree, smt, hist = solveTree!(fg, recordcliqs=ls(fg)) #, tree)

pl = plotBeehive_6(fg, meanmax=:mean)

cliq = getCliq(tree, :l0)

printCliqHistorySummary(tree, :l1)

sfg31_ad = getCliqSolveHistory(tree, :l0)[12][4].cliqSubFg
fg31_11 = getCliqSolveHistory(tree, :l0)[12][4].dfg
fg31_12 = getCliqSolveHistory(tree, :l0)[12][4].dfg

plotKDE(sfg31_ad, :l0, levels=1)

plotLocalProduct(sfg31_ad, :l0)
plotLocalProduct(fg, :l0)

drawGraph(sfg31_ad)

plotKDE(sfg31_ad, sortVarNested(ls(sfg31_ad)), levels=1)

plotKDE(fg31_11, sortVarNested(ls(sfg31_ad)), levels=1)
plotKDE(fg31_12, sortVarNested(ls(sfg31_ad)), levels=1)

plotCliqUpMsgs(fg, tree, :x5)

plotCliqDownMsgs(tree, :x5)

getDwnMsgs(tree, :x9)

ls(fg, :l6)

plotKDE(fg, :l6, levels=2)

plotLocalProduct(fg, :l6)

getSolverParams(fg).isfixedlag = true
getSolverParams(fg).qfl = 15

## DEBUG END

#


text:
function compute_S!(I::IntegralHelper{T, <:AbstractERI, O}, ints::IntegralHelper{T, <:AbstractERI, AtomicOrbitals}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}
    # Compute AO overlap
    Sao = ints["S"]

    # Convert AO overlap to MO overlap
    C = I.orbitals.C
    if eltype(I.orbitals.C) !== T
        C = T.(C)
    end
    @tensoropt S[i,j] := Sao[μ, ν]*C[μ,i]*C[ν,j]
    I.cache["S"] = S
end

function compute_T!(I::IntegralHelper{T, <:AbstractERI, O}, ints::IntegralHelper{T, <:AbstractERI, AtomicOrbitals}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}
    # Compute AO overlap
    Tao = ints["T"]

    # Convert AO overlap to MO overlap
    C = I.orbitals.C
    if eltype(I.orbitals.C) !== T
        C = T.(C)
    end
    @tensoropt Tmo[i,j] := Tao[μ, ν]*C[μ,i]*C[ν,j]
    I.cache["T"] = Tmo
end

function compute_V!(I::IntegralHelper{T, <:AbstractERI, O}, ints::IntegralHelper{T, <:AbstractERI, AtomicOrbitals}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}
    # Compute AO overlap
    Vao = ints["V"]

    # Convert AO overlap to MO overlap
    C = I.orbitals.C
    if eltype(I.orbitals.C) !== T
        C = T.(C)
    end
    @tensoropt V[i,j] := Vao[μ, ν]*C[μ,i]*C[ν,j]
    I.cache["V"] = V
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
#=--- TerminalGraphics / terminal.jl --------------------------------------=#

"""
Module `Terminal`, a submodule to `TerminalGraphics` provides functions
to query and manipulate features of the terminal emulator displaying the
current Julia session, in particular the Julia REPL.
"""
module Terminal

using .._default_sixel_io

"Helper function to send similar requests to a terminal emulator."
function _commandeer_terminal(io::Base.Terminals.TTYTerminal,
                             command::AbstractString)
    Base.Terminals.raw!(io, true)
    write(io, string(command))
    Base.Terminals.raw!(io, false)
end

"""
Helper function to send similar requests to a terminal emulator,
and expect an answer which is analyzed using a regular expression.
"""
function _query_terminal(io::Base.Terminals.TTYTerminal,
                        question::AbstractString, answer::Regex)
    Base.Terminals.raw!(io, true)
    write(io, string(question))
    response = transcode(String, readavailable(io))
    Base.Terminals.raw!(io, false)
    m = match(answer, response)
    (map(parse, m.captures)...)
end

"""
    windowtitle([io])

Get the title of the current terminal window.
"""
function windowtitle end

windowtitle() = windowtitle(_default_sixel_io())
windowtitle(io::IO) = _query_terminal(io, "\033[21t", r"\033\]l(.*)\033\\")

"""
    window_icon_title([io])

Get the current terminal window icon title.
"""
function window_icon_title end

window_icon_title() = window_icon_title(_default_sixel_io())
window_icon_title(io::IO) = _query_terminal(io, "\033[20t", r"\033\]L(.*)\033\\")

"""
    windowsize([io])

Determine the size of the terminal window in pixels, (width, height).
"""
function windowsize end

windowsize() = windowsize(_default_sixel_io())
windowsize(io::IO) = _query_terminal(io, "\033[14t", r"\033\[4;([0-9]*);([0-9]*)t")

"""
    window_position([io])

Determine the current terminal window position in pixels.
"""
function window_position end

window_position() = window_position(_default_sixel_io())
window_position(io::IO) = _query_terminal(io, "\033[13t", r"\033\[3;([0-9]*);([0-9]*)t")

"""
    cursorposition([io])

Determine the current position of the curser within the textarea in
characters, (row, col).
"""
function cursorposition end

cursorposition() = cursorposition(_default_sixel_io())
cursorposition(io::IO) = _query_terminal(io, "\033[6n", r"\033\[([0-9]*);([0-9]*)R")
        # an alternative would be "\033[?6n"

"""
    textarea([io])

Determine the current terminal textarea size in characters, (cols, rows).
"""
function textarea end

textarea() = textarea(_default_sixel_io())
textarea(io::IO) = _query_terminal(io, "\033[18t", r"\033\[8;([0-9]*);([0-9]*)t")

"""
    screensize([io])

Determine the screen size the terminal window is displayed on, in characters,
(maxcols, maxrows).
"""
function screensize end

screensize() = screensize(_default_sixel_io())
screensize(io::IO) = _query_terminal(io, "\033[19t", r"\033\[9;([0-9]*);([0-9]*)t")

"""
    colors_reverse([io])

Set the terminal to reversed color mode, typically black-on-white.
"""
function colors_reverse end

colors_reverse() = colors_reverse(_default_sixel_io())
colors_reverse(io::IO) = _commandeer_terminal(io, "\033[?5h")

"""
    colors_normal([io])

Set the terminal to normal color mode, typically white-on-black.
"""
function colors_normal end

colors_normal() = colors_normal(_default_sixel_io())
colors_normal(io::IO) = _commandeer_terminal(io, "\033[?5l")

"""
    hassixel([io])

Determine if the terminal emulator supports Sixel graphics.
"""
function hassixel end

hassixel() = hassixel(_default_sixel_io())

hassixel(repl::Base.REPL.AbstractREPL) = hassixel(outstream(repl))

# CSI Ps c  Send Device Attributes (Primary DA).
#     Ps = 0  or omitted -> request attributes from terminal.
#   The response depends on the decTerminalID resource setting.
#     -> CSI ? 1 ; 2 c  ("VT100 with Advanced Video Option")
#     -> CSI ? 1 ; 0 c  ("VT101 with No Options")
#     -> CSI ? 6 c  ("VT102")
#     -> CSI ? 6 2 ; Psc  ("VT220")
#     -> CSI ? 6 3 ; Psc  ("VT320")
#     -> CSI ? 6 4 ; Psc  ("VT420")
#   The VT100-style response parameters do not mean anything by
#   themselves.  VT220 (and higher) parameters do, telling the
#   host what features the terminal supports:
#     Ps = 1  -> 132-columns.
#     Ps = 2  -> Printer.
#     Ps = 3  -> ReGIS graphics.
#     Ps = 4  -> Sixel graphics.
#     Ps = 6  -> Selective erase.
#     Ps = 8  -> User-defined keys.
#     Ps = 9  -> National Replacement Character sets.
#     Ps = 1 5  -> Technical characters.
#     Ps = 1 8  -> User windows.
#     Ps = 2 1  -> Horizontal scrolling.
#     Ps = 2 2  -> ANSI color, e.g., VT525.
#     Ps = 2 9  -> ANSI text locator (i.e., DEC Locator mode).

# Note: it appears e.g. mlterm does not strictly comply.

function hassixel(io::IO)
    Base.Terminals.raw!(io, true)
    write(io, "\033[0c")
    r = transcode(String, Base.readavailable(io))
    Base.Terminals.raw!(io, false)

    # first entry is terminal generation; must be larger than 6 (VT102)
    # Note: it seems MLterm doesn't comply, but a '4' is nowhere else to be
    # found anyway, so just skip this.
    #first(m) <= 6 && return false

    m = map(parse, matchall(r"[0-9]+", r))
    return 4 in m
end

end  # submodule Terminal

#=--- end of file ---------------------------------------------------------=#


text:
using PyPlot
using JLD
nSteps = 100001
X = load("../data/sens/solenoid/metrics_$nSteps.jld")
dJds_ust = X["dJds_ust"]
dJds_st = X["dJds_st"]
x_trj = X["x_trj"]
t_trj = (atan.(x_trj[2,:],x_trj[1,:]) .+
		 2π) .% 2π
g = X["g"]
a = X["a"]
b = X["b"]
nv = X["normvs"]

runup = 5000
t_trj = t_trj[runup:end]
g = g[runup:end]
a = a[runup:end]
b = b[runup:end]
nv = nv[runup:end]

fig, ax = subplots()
ax.plot(t_trj, g, "r.", ms=0.5, label=L"g")
ax.axis("scaled")
ax.xaxis.set_tick_params(labelsize=28)
ax.yaxis.set_tick_params(labelsize=28)
ax.grid(true)
leg = ax.legend(fontsize=28)
leg.legendHandles[1]._legmarker.set_markersize(10)

fig, ax = subplots()
ax.plot(t_trj, a, "r.", ms=0.5, label=L"a")
ax.axis("scaled")
ax.xaxis.set_tick_params(labelsize=28)
ax.yaxis.set_tick_params(labelsize=28)
ax.grid(true)
leg = ax.legend(fontsize=28)
leg.legendHandles[1]._legmarker.set_markersize(10)

fig, ax = subplots()
ax.plot(t_trj, b, "r.", ms=0.5, label=L"b")
ax.axis("scaled")
ax.xaxis.set_tick_params(labelsize=28)
ax.yaxis.set_tick_params(labelsize=28)
ax.grid(true)
ax.set_xlabel(L"\theta", fontsize=28)
leg = ax.legend(fontsize=28)
leg.legendHandles[1]._legmarker.set_markersize(10)

fig, ax = subplots()
ax.plot(t_trj, nv, "b.", ms=0.5, label=L"\||v\||")
ax.set_xlabel(L"\theta", fontsize=28)
ax.xaxis.set_tick_params(labelsize=28)
ax.yaxis.set_tick_params(labelsize=28)
ax.grid(true)
leg = ax.legend(fontsize=28)
leg.legendHandles[1]._legmarker.set_markersize(10)

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
# using Plots
#
# function nchanges(c1, c2)
#   out = Int64(0)
#   n = length(c1)
#   for i in 1:(n - 1)
#     j = i + 1
#     while j <= n
#       if (c1[i] == c1[j]) != (c2[i] == c2[j])
#         out += Int64(1)
#         j = n + 1
#       else
#         j += 1
#       end
#     end
#   end
#   return out
# end
#
# function clustacf(output::Array, maxlag = 20)
#   n = size(output, 2)
#   acfpoints = Vector{Float64}(maxlag + 1)
#   acfpoints[1] = Float64(1)
#   for lag in 1:maxlag
#     tmp = Float64(1)
#     for i in (lag + 1):size(output, 1)
#       tmp -= nchanges(output[i, :], output[i - lag, :]) / ((n - 1) * (size(output, 1) - lag))
#     end
#     acfpoints[lag + 1] = tmp
#   end
#   plot(x=0:maxlag, y = acfpoints, yend = acfpoints * 0, xend = 0:maxlag,
#         Geom.segment,
#         Geom.point,
#         Guide.xticks(ticks = collect(0:maxlag)),
#         Guide.xlabel("Lag"),
#         Guide.ylabel("ACF"),
#         Coord.cartesian(ymin = 0, ymax = 1))
# end
#
#
# function clustrand(output::Array, maxlag = 20)
#   n = size(output, 2)
#   acfpoints = Vector{Float64}(maxlag + 1)
#   acfpoints[1] = Float64(1)
#   for lag in 1:maxlag
#     tmp = Float64(0)
#     for i in (lag + 1):size(output, 1)
#       tmp += randindex(output[i, :], output[i - lag, :])[1] / (size(output, 1) - lag)
#     end
#     acfpoints[lag + 1] = tmp
#   end
#   plot(x=0:maxlag, y = acfpoints, yend = acfpoints * 0, xend = 0:maxlag,
#         Geom.segment,
#         Geom.point,
#         Guide.xticks(ticks = collect(0:maxlag)),
#         Guide.xlabel("Lag"),
#         Guide.ylabel("ACF"),
#         Coord.cartesian(ymin = 0, ymax = 1))
# end


text:
"""
Partitioned Gauss-Legendre Runge-Kutta tableau with s stages

```julia
PartitionedTableauGauss(::Type{T}, s)
PartitionedTableauGauss(s) = PartitionedTableauGauss(Float64, s)
```

The constructor takes the number of stages `s` and optionally the element type `T` of the tableau.

This partitioned tableau uses [`TableauGauss`](@ref) for both coefficients `a` and `ā`.
"""
function PartitionedTableauGauss(s::Int)
    PartitionedTableau(Symbol("PGauss", s), TableauGauss(s))
end

"""
Partitioned Gauss-Lobatto IIIA-IIIB tableau with s stages

```julia
TableauLobattoIIIAIIIB(::Type{T}, s)
TableauLobattoIIIAIIIB(s) = TableauLobattoIIIAIIIB(Float64, s)
```

The constructor takes the number of stages `s` and optionally the element type `T` of the tableau.

This partitioned tableau uses [`TableauLobattoIIIA`](@ref) for `a` and [`TableauLobattoIIIB`](@ref) for `ā`.
"""
function TableauLobattoIIIAIIIB(s)
    PartitionedTableau(Symbol("LobattoIIIAIIIB", s), TableauLobattoIIIA(s), TableauLobattoIIIB(s))
end

"""
Tableau for Gauss-Lobatto IIIB-IIIA method with s stages

```julia
TableauLobattoIIIBIIIA(::Type{T}, s)
TableauLobattoIIIBIIIA(s) = TableauLobattoIIIBIIIA(Float64, s)
```

The constructor takes the number of stages `s` and optionally the element type `T` of the tableau.

This partitioned tableau uses [`TableauLobattoIIIB`](@ref) for `a` and [`TableauLobattoIIIA`](@ref) for `ā`.
"""
function TableauLobattoIIIBIIIA(s)
    PartitionedTableau(Symbol("LobattoIIIBIIIA", s), TableauLobattoIIIB(s), TableauLobattoIIIA(s))
end

"""
Tableau for Gauss-Lobatto IIIC-IIIC̄ method with s stages

```julia
TableauLobattoIIICIIIC̄(::Type{T}, s)
TableauLobattoIIICIIIC̄(s) = TableauLobattoIIICIIIC̄(Float64, s)
```

The constructor takes the number of stages `s` and optionally the element type `T` of the tableau.

This partitioned tableau uses [`TableauLobattoIIIC`](@ref) for `a` and [`TableauLobattoIIIC̄`](@ref) for `ā`.
"""
function TableauLobattoIIICIIIC̄(s)
    PartitionedTableau(Symbol("LobattoIIICIIIC̄", s), TableauLobattoIIIC(s), TableauLobattoIIIC̄(s))
end

"""
Tableau for Gauss-Lobatto IIIC̄-IIIC method with s stages

```julia
TableauLobattoIIIC̄IIIC(::Type{T}, s)
TableauLobattoIIIC̄IIIC(s) = TableauLobattoIIIC̄IIIC(Float64, s)
```

The constructor takes the number of stages `s` and optionally the element type `T` of the tableau.

This partitioned tableau uses [`TableauLobattoIIIC̄`](@ref) for `a` and [`TableauLobattoIIIC`](@ref) for `ā`.
"""
function TableauLobattoIIIC̄IIIC(s)
    PartitionedTableau(Symbol("LobattoIIIC̄IIIC", s), TableauLobattoIIIC̄(s), TableauLobattoIIIC(s))
end


text:
module CLIMA

include(joinpath("Utilities", "TicToc", "TicToc.jl"))
include(joinpath("Utilities", "ArtifactWrappers", "ArtifactWrappers.jl"))
include(joinpath("Utilities", "ParametersType", "ParametersType.jl"))
include(joinpath("InputOutput", "Writers", "Writers.jl"))
include(joinpath("Common", "ConfigTypes", "ConfigTypes.jl"))
include(joinpath("Common", "PlanetParameters", "UniversalConstants.jl"))
include(joinpath("Common", "PlanetParameters", "PlanetaryConstants.jl"))
include(joinpath("Common", "PlanetParameters", "PlanetParameters.jl"))
include(joinpath("Utilities", "VariableTemplates", "VariableTemplates.jl"))
include(joinpath("Common", "MoistThermodynamics", "MoistThermodynamics.jl"))
include(joinpath(
    "Atmos",
    "Parameterizations",
    "CloudPhysics",
    "MicrophysicsParameters.jl",
))
include(joinpath(
    "Atmos",
    "Parameterizations",
    "CloudPhysics",
    "Microphysics.jl",
))
include(joinpath(
    "Atmos",
    "Parameterizations",
    "SurfaceFluxes",
    "SurfaceFluxes.jl",
))
include(joinpath(
    "Atmos",
    "Parameterizations",
    "SubgridScaleTurbulence",
    "SubgridScaleParameters.jl",
))
include(joinpath("Arrays", "MPIStateArrays.jl"))
include(joinpath("Mesh", "Mesh.jl"))
include(joinpath("DGmethods", "Courant.jl"))
include(joinpath("DGmethods", "SpaceMethods.jl"))
include(joinpath("DGmethods", "DGmethods.jl"))
include(joinpath("Ocean", "ShallowWater", "ShallowWaterModel.jl"))
include(joinpath(
    "Ocean",
    "HydrostaticBoussinesq",
    "HydrostaticBoussinesqModel.jl",
))
include(joinpath("DGmethods_old", "DGBalanceLawDiscretizations.jl"))
include(joinpath("LinearSolvers", "LinearSolvers.jl"))
include(joinpath("LinearSolvers", "GeneralizedConjugateResidualSolver.jl"))
include(joinpath("LinearSolvers", "GeneralizedMinimalResidualSolver.jl"))
include(joinpath("LinearSolvers", "ColumnwiseLUSolver.jl"))
include(joinpath("ODESolvers", "ODESolvers.jl"))
include(joinpath("ODESolvers", "GenericCallbacks.jl"))
include(joinpath("Utilities", "Callbacks", "Callbacks.jl"))
include(joinpath("Atmos", "Model", "AtmosModel.jl"))
include(joinpath("InputOutput", "VTK", "VTK.jl"))
include(joinpath("Diagnostics", "Diagnostics.jl"))
include(joinpath("Driver", "Driver.jl"))

end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
# run examples/numbering3d.jl to check node numbers
# using Revise
using Test
using Logging
using LarSurf
using SparseArrays
# using LinearAlgebraicRepresentation
# Lar = LinearAlgebraicRepresentation
# Logging.configure(level==Logging.Debug)

# include("../src/LarSurf.jl")
# include("../src/block.jl")

@testset "Block basic function Tests" begin
    slides0 = LarSurf.data_sub_from_block_sub([5, 5, 5], 0, [1,2,3])
    @test all(slides0 .== [1, 5, 6, 10, 11, 15])

    # slides0 = LarSurf.data_sub_from_block_sub([5, 5, 5], 1, [1,2,3])
    # @test all(slides0 .== [0, 6, 5, 11, 10, 16])

end

@testset "Block complex function Tests" begin
    block_size = [5, 5, 5]
    margin_size = 0

    ## Artifical data
    data3d = ones(Int16, 19, 15, 31)
    data3d[2:14,8:13,9:11] .= 10
    # data3d

    voxelsize_mm = [0.5, 0.9, 0.8]
    threshold=0

    blocks_number, blocks_number_axis = LarSurf.number_of_blocks_per_axis(
        size(data3d), block_size)
    @test blocks_number == 84
    @test blocks_number_axis == [4, 3, 7]

    block1, offset1, size1 = LarSurf.get_block(1,
        data3d, block_size, margin_size, blocks_number_axis
    )

    block4, offset4, size4 = LarSurf.get_block(4,
        data3d, block_size, margin_size, blocks_number_axis
    )
    # println("block1:", block1, offset1, size1)
    # println("block4:", block4, offset4, size4)
    @debug "block1" block1
    @debug "block1" block4
    @test size(block1)[1] == 5
    @test size(block1)[3] == 5
    @test size(block1)[3] == 5
    @test size(block4)[1] == 4
    @test size(block4)[2] == 5
    @test size(block4)[3] == 5
    # block12 = LarSurf.get_block(
    #     data3d, block_size, 1, blocks_number_axis, 2
    # )
    LarSurf.debug_array_stats(block1)
    LarSurf.debug_array_stats(block4)
    LarSurf.debug_array_stats(data3d)

end

@testset "Get face ID test" begin
    faces= LarSurf.get_face_ids_from_cube_in_grid([1,2,3], [1,1,1], false)
    @test collect(faces) == [1, 13, 22]
    # print(faces, "\n")
    faces = LarSurf.get_face_ids_from_cube_in_grid([1,2,3], [1,2,1], false)
    @test collect(faces) == [4, 16, 26]
    # print(faces, "\n")
    faces = LarSurf.get_face_ids_from_cube_in_grid([1,2,3], [1,1,2], false)
    @test collect(faces) == [2, 14, 23]
    # print(faces, "\n")

    faces = LarSurf.get_face_ids_from_cube_in_grid([2,3,4], [2,2,3], false)
    # print(faces, "\n")
    @test collect(faces) == [19, 59, 91]

end

@testset "Decode face ID in sub block test" begin
    data_size = [2,3,4]
    block_size = [2,2,2]
    offset = [0,2,2]
    # block_size [2,2,2]
    #  1 - first face on axis1
    # 12 - last face on axis1
    # 13 - first face axis 2
    # 24 - last face axis 2
    # 25 - first face axis 3
    # 36 - last face axis 3

    fid_big, cart = LarSurf.sub_grid_face_id_to_orig_grid_face_id(data_size, block_size, [0,0,0], 1)
    @test fid_big == 1

    fid_big, cart = LarSurf.sub_grid_face_id_to_orig_grid_face_id(data_size, block_size, [0,0,0], 2)
    @test fid_big == 2

    fid_big, cart = LarSurf.sub_grid_face_id_to_orig_grid_face_id(data_size, block_size, [0,0,0], 3)
    @test fid_big == 5

    fid_big, cart = LarSurf.sub_grid_face_id_to_orig_grid_face_id(data_size, block_size, [0,0,0], 12)
    @test fid_big == 30

    fid_big, cart = LarSurf.sub_grid_face_id_to_orig_grid_face_id(data_size, block_size, [0,0,0], 13)
    @test fid_big == 37

    fid_big, cart = LarSurf.sub_grid_face_id_to_orig_grid_face_id(data_size, block_size, [0,0,0], 24)
    @test fid_big == 62

    # moved offset [0,0,2]
    fid_big, cart = LarSurf.sub_grid_face_id_to_orig_grid_face_id(data_size, block_size, [0,0,2], 24)
    @test fid_big == 64

    # offset [0,2,2]
    fid_big, cart = LarSurf.sub_grid_face_id_to_orig_grid_face_id(data_size, block_size, [0,2,2], 22)
    @test fid_big == 68
    fid_big, cart = LarSurf.sub_grid_face_id_to_orig_grid_face_id(data_size, block_size, [0,2,2], 27)
    @test fid_big == 83
    # faces= LarSurf.get_face_ids_from_cube_in_grid([1,2,3], [1,1,1], false)
    # @test collect(faces) == [1, 13, 22]
    # # print(faces, "\n")
    # faces = LarSurf.get_face_ids_from_cube_in_grid([1,2,3], [1,2,1], false)
    # @test collect(faces) == [4, 16, 26]
    # # print(faces, "\n")
    # faces = LarSurf.get_face_ids_from_cube_in_grid([1,2,3], [1,1,2], false)
    # @test collect(faces) == [2, 14, 23]
    # # print(faces, "\n")
    #
    #
    # faces = LarSurf.get_face_ids_from_cube_in_grid([2,3,4], [2,2,3], false)
    # # print(faces, "\n")
    # @test collect(faces) == [19, 59, 91]

end

@testset "Get inner small block" begin
    faces = LarSurf.cube_in_block_surface([2,3,4], [1,1,1], [1,1,1])
    @test sort(collect(faces)) == sort([1, 13, 37, 41, 69, 70])
    # print(faces, "\n")
end

@testset "Get inner big block" begin
    faces = LarSurf.cube_in_block_surface([2,3,4], [1,1,1], [1,2,2])
    @test sort(collect(faces)) == sort([1, 13, 2, 14, 5, 17, 6, 18, 37, 45, 38, 46, 69, 71, 74, 76])
    # print(faces, "\n")
end

# TODO test subgrid face id based on small_block_face_id...

@testset "Get inner block slope" begin

    segmentation = LarSurf.generate_slope([11,12,13])
    data_size = LarSurf.size_as_array(size(segmentation))
    data3d = segmentation
    block_size = [2,3,4]
    block_i = 91
    margin_size = 0

block_number, blocks_number_axis = LarSurf.number_of_blocks_per_axis(data_size, block_size)
    a = Array{Int}(
        undef,
        blocks_number_axis[1],
        blocks_number_axis[2],
        blocks_number_axis[3]
    )

    bsub = CartesianIndices(a)[block_i]
#     faces = LarSurf.cube_in_block_surface([1,2,3], [1,1,1], [1,2,2])
    bsub_arr = [bsub[1], bsub[2], bsub[3]]

    first = (bsub_arr .== [1, 1, 1])
    last = (bsub_arr .== blocks_number_axis)
#     print(" bsub :", bsub, "block number axis", blocks_number_axis, " first last ", first, last, "\n")
    # if ----------
    # print(" end of col, row or slice ", bsub, "\n")

    xst, xsp, yst, ysp, zst, zsp = LarSurf.data_sub_from_block_sub(
        block_size, margin_size, bsub
    )
    # print("input block size ")
    LarSurf.print_slice3(xst, xsp, yst, ysp, zst, zsp)
    xst, oxst, xsh = LarSurf.get_start_and_outstart_ind(xst, margin_size)
    yst, oyst, ysh = LarSurf.get_start_and_outstart_ind(yst, margin_size)
    zst, ozst, zsh = LarSurf.get_start_and_outstart_ind(zst, margin_size)
#         print("[", xsh, ", ", ysh, ",", zsh, "]")
#     szx, szy, szz = size(data3d)
    szx, szy, szz = data_size
    bszx, bszy, bszz = block_size
    xsp, oxsp = LarSurf.get_end_and_outend_ind(xst, xsp, szx, xsh)
    ysp, oysp = LarSurf.get_end_and_outend_ind(yst, ysp, szy, ysh)
    zsp, ozsp = LarSurf.get_end_and_outend_ind(zst, zsp, szz, zsh)
    # print("postprocessing input")
    LarSurf.debug_slice3(xst, xsp, yst, ysp, zst, zsp)
    # print("postprocessing output")
    LarSurf.debug_slice3(oxst, oxsp, oyst, oysp, ozst, ozsp)

    outdata = zeros(
        eltype(data3d),
        xsp,
        ysp,
        zsp
#         block_size[1] + margin_size,
#         block_size[2] + margin_size,
#         block_size[3] + margin_size
    )
    outdata[oxst:oxsp, oyst:oysp, ozst:ozsp] = data3d[
        xst:xsp, yst:ysp, zst:zsp
    ]

#     @test collect(faces) == [1, 13, 22]
    # print(faces, "\n")
end

@testset "node id in grid" begin
    grid_size = [2, 3, 4]
    @test 27 == LarSurf.grid_voxel_cart_to_node_id(grid_size, [2, 2, 2])
    @test 41 == LarSurf.grid_voxel_cart_to_node_id(grid_size, [3, 1, 1])
    @test 60 == LarSurf.grid_voxel_cart_to_node_id(grid_size, [3, 4, 5])

end

@testset "number faces in grid" begin
    # lmodel::Lar.LARmodel = Lar.cuboidGrid(block_size, true)
    data_size = [3,4,5]
    # bigV, (bigVV, bigEV, bigFV, bigCV) = Lar.cuboidGrid(data_size, true)
    bigV, (bigVV, bigEV, bigFV, bigCV) = LarSurf.Lar.cuboidGrid(data_size, true)
    @test size(bigFV)[1] == LarSurf.grid_number_of_faces(data_size)

end

@testset "node ids in grid from face id - no correspondence check" begin
    data_size = [2, 3, 4]
    nodes_ids, nodes_carts = LarSurf.grid_face_id_to_node_ids(data_size, 20)
    @test sort(nodes_ids) == [29, 30, 34, 35]
    @test nodes_carts[1] == [2,2,4]
    @test nodes_carts[2] == [2,3,4]
    @test nodes_carts[3] == [2,3,5]
    @test nodes_carts[4] == [2,2,5]
    # test one of last faces, there was a problem
    nodes_ids, nodes_carts = LarSurf.grid_face_id_to_node_ids(data_size, 98)
    @test sort(nodes_ids) == [35, 40, 55, 60]
    nodes_carts_sorted = sort(nodes_carts)
    @test nodes_carts_sorted[1] == [2,3,5]
    @test nodes_carts_sorted[2] == [2,4,5]
    @test nodes_carts_sorted[3] == [3,3,5]
    @test nodes_carts_sorted[4] == [3,4,5]

    nodes_ids, nodes_carts = LarSurf.grid_face_id_to_node_ids(data_size, 67)
    @test sort(nodes_ids) == [38, 39, 58, 59]
    nodes_carts_sorted = sort(nodes_carts)
    @test nodes_carts_sorted[1] == [2,4,3]
    @test nodes_carts_sorted[2] == [2,4,4]
    @test nodes_carts_sorted[3] == [3,4,3]
    @test nodes_carts_sorted[4] == [3,4,4]

    nodes_ids, nodes_carts = LarSurf.grid_face_id_to_node_ids(data_size, 52)
    @test sort(nodes_ids) == [19,20,39,40]
    nodes_carts_sorted = sort(nodes_carts)
    @test nodes_carts_sorted[1] == [1,4,4]
    @test nodes_carts_sorted[2] == [1,4,5]
    @test nodes_carts_sorted[3] == [2,4,4]
    @test nodes_carts_sorted[4] == [2,4,5]
end

@testset "node ids in grid from face id with checking correspondence of id" begin
    # here testing also the positions
    data_size = [2,3,4]
    # data_size = [3,4,5]
    nodes_ids, nodes_carts = LarSurf.grid_face_id_to_node_ids(data_size, 83)
    # display(nodes_ids)
    @test sort(nodes_ids) == [15,20,35,40]
    @test LarSurf.check_faces_equal(nodes_ids, [20,15,35,40])
    @test LarSurf.array_equal_roll_invariant(nodes_ids, [20,15,35,40])
    nodes_carts_sorted = sort(nodes_carts)
    @test nodes_carts[findall(x->x==15, nodes_ids)[1]] == [1,3,5]
    @test nodes_carts[findall(x->x==20, nodes_ids)[1]] == [1,4,5]
    @test nodes_carts[findall(x->x==35, nodes_ids)[1]] == [2,3,5]
    @test nodes_carts[findall(x->x==40, nodes_ids)[1]] == [2,4,5]
end

@testset "get block test with fixed block size" begin
    segmentation = zeros(Int8, 5, 6, 7)
    segmentation[2:5,2:5,2:6] .= 1
    block_size = [2,2,2]
    data_size = LarSurf.size_as_array(size(segmentation))
    block_number, blocks_number_axis = LarSurf.number_of_blocks_per_axis(
        data_size, block_size)
    fixed_block_size=true
    block1, offset1, block_size1 = LarSurf.get_block(3,
        segmentation, block_size, 0,
        blocks_number_axis, fixed_block_size
        )
    @test size(block1,1) == 2
    @test size(block1,2) == 2
    @test size(block1,3) == 2
end

@testset "Block getter fixed size" begin

    segmentation = LarSurf.data234()
    block_size = [2,2,2]
    n, bgetter = LarSurf.block_getter(segmentation, block_size; fixed_block_size=true)
    @test n==4
    block1, offset1, block_size1 = LarSurf.get_block(2, bgetter...)
    @test size(block1,1) == 2
    @test size(block1,2) == 2
    @test size(block1,3) == 2
end

@testset "Block getter not fixed size" begin
    segmentation = LarSurf.data234()
    block_size = [2,2,2]
    n, bgetter = LarSurf.block_getter(segmentation, block_size; fixed_block_size=false)
    @test n==4
    block1, offset1, block_size1 = LarSurf.get_block(2, bgetter...)
    @test size(block1,1) == 2
    @test size(block1,2) == 1
    @test size(block1,3) == 2
end


text:
using SMC, Base.Test

include("_legacy_.jl")

srand(125)

## Linear Gaussian

dx, dy = 5, 3
A = randn(dx,dx); A /= 0.9norm(A)
B = randn(dy,dx); B /= 1.1norm(B)
Q = randn(dx,dx); Q *= Q'; Q += 0.5*eye(dx); Q += Q'; Q /= 5
R = randn(dy,dy); R *= R'; R += 0.5*eye(dy); R += R'; R /= 20

lg  = LinearGaussian(A,B,Q,R)
hmm = HMM(lg)
x0  = randn(dx)

K = 100
(states, observations) = generate(lg, x0, K)

x00 = x0+randn(dx)/5

srand(12)
kf = kalmanfilter(lg, observations, x00, eye(dx))
srand(12)
(kfm_leg, kfc_leg, kfm__leg, kfc__leg) = kf_legacy(observations, A, B, Q,
                                                    R, K, x00, eye(dx))

@test isapprox(kf.means, kfm_leg)
@test isapprox(kf.covs, kfc_leg)
@test isapprox(kf.means_, kfm__leg)
@test isapprox(kf.covs_, kfc__leg)

# fragile test
@test norm(kf.means-states)/norm(states) < 0.25

srand(32)
ks = kalmansmoother(lg, observations, kf)
srand(32)
(ksm_leg, ksc_leg) = ks_legacy(observations, A, B, Q, R, K,
                                kfm__leg, kfc__leg)

@test isapprox(ks.means, ksm_leg)
@test isapprox(ks.covs, ksc_leg)

# fragile test
@test norm(ks.means-states)/norm(states) < 0.2


text:
module SuboptimalLocalAlignments

# package code goes here

    include("AminoAcids.jl")
    using Gadfly
    #### Parameters ####
    indelPenalty = -3
    affineStart = -100
    affineExtend = -1
    threshold = 30
    width = 3
    ####################

    function plotAllSuboptimalAlignment(alignmentListFilt, xAminoAcid, yAminoAcid, path)
        xAxis = copy(xAminoAcid.content)
        unshift!(xAxis, "0")
        yAxis = copy(yAminoAcid.content)
        unshift!(yAxis, "0")

        xLen = length(xAminoAcid.content)+1
        yLen = length(yAminoAcid.content)+1

        allShowPathMatrix = getAllShowPathMatrix(alignmentListFilt, xLen, yLen, path)
        draw(SVG("allPath.svg", 20inch, 20inch),
        spy(allShowPathMatrix, Scale.y_discrete(labels = i->yAxis[i]), Scale.x_discrete(labels = i->xAxis[i])))
    end

    function plotSuboptimalAlignmentGroup(res, xAminoAcid, yAminoAcid, path)
        xAxis = copy(xAminoAcid.content)
        unshift!(xAxis, "0")
        yAxis = copy(yAminoAcid.content)
        unshift!(yAxis, "0")

        xLen = length(xAminoAcid.content)+1
        yLen = length(yAminoAcid.content)+1

        for group = res[1]
            pathMark = 0
            if group.memberNum > 1
                pathMark += 1
                showPathMatrix = getGroupShowPathMatrix(group.members, pathMark, xLen, yLen, path)
                fileName = string("GroupX", dec(pathMark))
                fileName = string(fileName, ".svg")
                draw(SVG(fileName, 20inch, 20inch),
                spy(showPathMatrix, Scale.y_discrete(labels = i->yAxis[i]), Scale.x_discrete(labels = i->xAxis[i])))
            end
        end

        for group = res[2]
            pathMark = 0
            if group.memberNum > 1
                pathMark += 1
                showPathMatrix = getGroupShowPathMatrix(group.members, pathMark, xLen, yLen, path)
                fileName = string("GroupY", dec(pathMark))
                fileName = string(fileName, ".svg")
                draw(SVG(fileName, 20inch, 20inch),
                spy(showPathMatrix, Scale.y_discrete(labels = i->yAxis[i]), Scale.x_discrete(labels = i->xAxis[i])))
            end
        end

    end

    function getGroupShowPathMatrix(alignmentList, pathMark, xLen, yLen, path)

        showPath = zeros(yLen, xLen)

        for a in alignmentList
            row = a.endPos[1]
            col = a.endPos[2]

            while path[row, col] != 0
                showPath[row, col] = pathMark

                if path[row, col] == 1
                    row -= 1
                    col -= 1
                elseif path[row, col] == 2 # indel from left
                    col -= 1
                elseif path[row, col] == 3 # indel from right
                    row -= 1
                end
            end
        end

        return showPath
    end

    function getAllShowPathMatrix(alignmentList, xLen, yLen, path)
        pathMark = 0
        showPath = zeros(yLen, xLen)

        for a in alignmentList
            pathMark += 1
            row = a.endPos[1]
            col = a.endPos[2]

            while path[row, col] != 0
                showPath[row, col] = pathMark

                if path[row, col] == 1
                    row -= 1
                    col -= 1
                elseif path[row, col] == 2 # indel from left
                    col -= 1
                elseif path[row, col] == 3 # indel from up
                    row -= 1
                end
            end
        end
        return showPath
    end

    function reportGroup(res)
        for group = res[1]
            pathMark = 0
            if group.memberNum > 1
                pathMark += 1
                println("X axis Alignment Group$(pathMark)")
                for member in group.members
                    println("Score = $(member.score), start = $(member.startPos), end = $(member.endPos):")
                    println(member.seqX)
                    println(member.seqY)
                    println("start position = $(member.startPos)")
                    println("end position = $(member.endPos)")
                    println()
                end
            end
        end

        for group = res[2]
            pathMark = 0
            if group.memberNum > 1
                pathMark += 1
                println("Y axis Alignment Group$(pathMark):")
                for member in group.members
                    println("Score = $(member.score), start = $(member.startPos), end = $(member.endPos):")
                    println(member.seqX)
                    println(member.seqY)
                    println("start position = $(member.startPos)")
                    println("end position = $(member.endPos)")
                    println()
                end
            end
        end
    end

    function getFilteredAllSuboptimalAlignment(dataPath; indelPenalty = -3, affineStart = -12, affineExtend = -1, threshold = 30, width = 3, useAffine = true, plotAllAlignments = false, plotGroupAlignments = false)

        list = loadData(dataPath)
        blosum = getBlosum()
        xAminoAcid = list[1]
        yAminoAcid = list[2]
        xLen = length(xAminoAcid.content)+1
        yLen = length(yAminoAcid.content)+1
        matrix = zeros(yLen, xLen)

        path = zeros(yLen, xLen)

        for row = 2:yLen
            for col = 2:xLen

                #=
                = In path matrix, 1 stands for match/mismatch,
                = 2 stands for indel from left
                = 3 stands for indel from up
                =#

                # if we chose the path with match/mismatch
                a = matrix[row-1, col-1] + getMappingScore(yAminoAcid.content[row-1], xAminoAcid.content[col-1], blosum)
                # calculate indel

                b = 0
                if useAffine
                    if path[row-1, col] == 2 # There already indels from left at previous left cell
                        b = matrix[row-1, col] + affineExtend # 2: indel from left
                    else
                        b = matrix[row-1, col] + affineStart
                    end
                else
                    b = matrix[row-1, col] + indelPenalty
                end

                c = 0
                if useAffine
                    if path[row, col-1] == 3 # There already indels from up at previous up cell
                        c = matrix[row, col-1] + affineExtend
                    else
                        c = matrix[row, col-1] + affineStart # 3: indel from up
                    end
                else
                    c = matrix[row, col-1] + indelPenalty
                end

                max = 0
                if (a > max)
                    max = a
                    path[row, col] = 1
                end
                if (b > max)
                    max = b
                    path[row, col] = 2
                end
                if (c > max)
                    max = c
                    path[row, col] = 3
                end
                matrix[row, col] = max
            end
        end

        ## Get all of the alignments above the threshold score in a list.
        alignmentList = []
        while getAlignment(matrix, path, xLen, yLen, alignmentList, threshold, width, xAminoAcid, yAminoAcid)
        end

        ## Filter out overlapping suboptimal alignments
        alignmentListFilt = filterAlignments(alignmentList)

        if plotAllAlignments
            plotAllSuboptimalAlignment(alignmentListFilt, xAminoAcid, yAminoAcid, path)
        end

        if plotGroupAlignments
            res = mergeAlignment(alignmentListFilt)
            reportGroup(res)
            plotSuboptimalAlignmentGroup(res, xAminoAcid, yAminoAcid, path)
        end

        return alignmentListFilt

    end

    #=@time alignmentListFilt = getFilteredAllSuboptimalAlignment("../data/test2.fa", plotAllAlignments = true, plotGroupAlignments = true)
    ## Print all of the alignments
    for a in alignmentListFilt
        println("Score = $(a.score), start = $(a.startPos), end = $(a.endPos):")
        println(a.seqX)
        println(a.seqY)
        println("start position = $(a.startPos)")
        println("end position = $(a.endPos)")
        println()
    end=#

end # module

max_stars_repo_name: yguoogo/Suboptimal-Local-Alignment

max_issues_repo_name: yguoogo/Suboptimal-Local-Alignment

max_forks_repo_name: yguoogo/Suboptimal-Local-Alignment


text:
mutable struct Box{T}
  val::T
end

@inline val(box::Box) = box.val
@inline val(ref::Ref) = ref.x


text:
"""
    AbstractLayeredArray{layer, T, N}

Supertype for `N`-dimensional layered arrays with elements of type `T`.
"""
abstract type AbstractLayeredArray{layer, T, N} <: AbstractArray{T, N} end

const AbstractLayeredVector{layer, T} = AbstractLayeredArray{layer, T, 1}
const AbstractLayeredMatrix{layer, T} = AbstractLayeredArray{layer, T, 2}

layerof(::AbstractLayeredArray{layer}) where {layer} = layer
layerof(::Type{<: AbstractLayeredArray{layer}}) where {layer} = layer
layerof(::Any) = 0 # bottom layer

# getindex for slice
# cannot overload `getindex` since `getindex(A::AbstractLayeredArray, i::Int...)` also calls it.
# so following unexported `_getindex` is overloaded
@inline function Base._getindex(::IndexCartesian, x::AbstractLayeredArray{layer}, i::Union{Real, AbstractArray}...) where {layer}
    @boundscheck checkbounds(x, i...)
    @inbounds LayeredArray{layer}(view(x, i...))
end
# following function fixes ambiguity
@inline function Base._getindex(::IndexCartesian, A::AbstractLayeredArray, I::Int...)
    @boundscheck checkbounds(A, I...) # generally _to_subscript_indices requires bounds checking
    @inbounds r = getindex(A, Base._to_subscript_indices(A, I...)...)
    r
end

function set!(dest::AbstractLayeredArray, src::AbstractLayeredArray)
    layerof(dest) == layerof(src) || throw(ArgumentError("layers must match in setting values, tried layer $(layerof(src)) -> $(layerof(dest))"))
    @simd for i in eachindex(dest, src)
        @inbounds dest[i] = src[i]
    end
    dest
end

function Base.similar(dest::AbstractLayeredArray{layer}, ::Type{ElType}, dims::Tuple{Vararg{Int}}) where {layer, ElType}
    LayeredArray{layer}(Array{ElType}(undef, dims))
end


text:
mutable struct CachedArray{T,N,NC,AA<:AbstractArray} <: AbstractCachedArray{T,N,T,NC,AA}
    parent::AA
    cached::Array{T,NC}
    current_I::Tuple
end

#not type stable, but as a constructor it's not performance-critical
function CachedArray(A::AbstractArray{T,N}, nd_cache::Int) where {T,N}
    @assert nd_cache < N
    cache_rngs = axes(A)[nd_cache+1:N]
    ci = map(first, cache_rngs)
    cached = similar(view(A, axes(A)[1:nd_cache]..., ci...))
    s = CachedArray{T,N,nd_cache,typeof(A)}(A, cached, ci)
    update_cache!(s, ci)
    return s
end


text:
COUNT=0

function great_painter(r::Robot)
    global COUNT=0
    overal_temp=0
    side=Ost
    while isborder(r,Nord)==false
        overal_temp+=row_temp(r,side)
        move!(r,Nord)
        side=inverse(side)
    end
    overal_temp+=row_temp(r,side)
    average_temp=overal_temp/COUNT
    return average_temp
end

function row_temp(r::Robot, side::HorizonSide)
    row_temp=0
    while isborder(r,side)==false
        row_temp+=check_temp(r)
        move!(r,side)
    end
    row_temp+=check_temp(r)
    return row_temp
end

function check_temp(r::Robot)
    global COUNT
    if ismarker(r)==true
        COUNT+=1
        return temperature(r)
    else
        return 0
    end
end

inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4))


text:
global active_repl
global active_repl_backend

type JuliaView <: TextView
  handle :: Ptr{GObject}
  function JuliaView()
    view = @TextView()
    Gtk.gobject_move_ref(new(view.handle), view)
  end
end
immutable REPLManager <: ClusterManager
  np :: Integer
end

function launch(manager :: REPLManager, params :: Dict, launched :: Array, c :: Condition)

end

function testREPELInit()
  global active_repl
  global active_repl_backend
  sin = Base.TTY(fd(STDIN))
  sout = Base.TTY(fd(STDOUT))
  serr = Base.TTY(fd(STDERR))
  term = Base.Terminals.TTYTerminal("", sin, sout, serr)
  Base.REPL.banner(term, term)
  active_repl = Base.REPL.LineEditREPL(term, true)
  active_repl.hascolor = false
  Base._atreplinit(active_repl)
  Base.REPL.run_repl(active_repl, backend -> (global active_repl_backend = backend))
end

max_stars_repo_name: barcharcraz/eve-julia-gui

max_issues_repo_name: barcharcraz/eve-julia-gui

max_forks_repo_name: barcharcraz/eve-julia-gui


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

#=
The intersection type can be one of three types:

1. intersect at one point
2. overlap at more than one point
3. do not overlap nor intersect
=#
function intersection(f, l1::Line{3,T}, l2::Line{3,T}) where {T}
  a, b = l1(0), l1(1)
  c, d = l2(0), l2(1)

  if measure(Tetrahedron(a, b, c, d)) > 0 # not in same plane
    return @IT NoIntersection nothing f
  elseif isapprox(norm((b - a) × (c - d)), zero(T), atol=atol(T)^2)
    if a in l2
      return @IT OverlappingLines l1 f
    else # parallel lines
      return @IT NoIntersection nothing f
    end
  else
    return @IT CrossingLines intersectpoint(l1, l2) f
  end
end

function intersection(f, l1::Line{2,T}, l2::Line{2,T}) where {T}
  a, b = l1(0), l1(1)
  c, d = l2(0), l2(1)

  if !isapprox(abs((b - a) × (c - d)), zero(T), atol=atol(T)^2)
    return @IT CrossingLines intersectpoint(l1, l2) f
  elseif isapprox(measure(Triangle(a, b, c)), zero(T), atol=atol(T)^2)
    return @IT OverlappingLines l1 f
  else
    return @IT NoIntersection nothing f
  end
end

# compute the intersection of two lines assuming that it is a point
function intersectpoint(l1::Line, l2::Line)
  a, b = l1(0), l1(1)
  c, d = l2(0), l2(1)

  v1  = a - b
  v2  = c - d
  v12 = a - c

  # the intersection point lies in between a and b at a fraction s
  # (https://en.wikipedia.org/wiki/Line-line_intersection#Formulas)
  s = (v12[1] * v2[2] - v12[2] * v2[1]) / (v1[1] * v2[2] - v1[2] * v2[1])

  a - s*v1
end


text:
# NOTE: these tests were taken from the scipy.optimize test suite

#=
- f1 is a simple quadratic with roots 0 and 1
- f2 is a symmetric parabola, x^2 - 1
- f3 is a quartic polynomial with large hump in interval
- f4 is step function with a discontinuity at 1
- f5 is a hyperbola with vertical asymptote at 1
- f6 has random values positive to left of 1, negative to right

Of course these are not real problems. They just test how the 'good' solvers
behave in bad circumstances where bisection is really the best. A good solver
should not be much worse than bisection in such circumstance, while being faster
for smooth monotone sorts of functions.
=#

f1(x) = x*(x-1)
f2(x) = x^2 - 1
f3(x) = x * (x-1) * (x-2) * (x-3.)

f4(x) = x > 1 ? 1.0 + .1*x :
        x < 1 ? -1.0 + .1*x :
        0.0

f5(x) = x != 1 ? 1.0/(1. - x) : 0.0

f6(x) = x > 1 ? rand() :
        x < 1 ? -rand() :
        0.0

funcs = [f1, f2, f3, f4, f5, f6]

@testset "Testing univariate root finders" begin

    a = 0.5
    b = sqrt(3)

    for root_finder in [bisect, brent, brenth, ridder]
        nm = string(root_finder)

        @testset "testing $root_finder" begin
            for f in funcs
                root = root_finder(f, a, b; xtol=0.1e-12)
                @test abs(root - 1.0) < 1e-11
            end
        end
    end

end


text:
export filter_events_by_duration

# ============================================================================ #
function filter_events_by_duration(evt::Vector{<:Real}, dur::Real; blank::Real=NaN)
    df = diff(evt)
    evt_ts = Vector{Float64}(undef, 0)
    evt_off = Vector{Float64}(undef, 0)
    nblank = 0
    warning = false
    for k in eachindex(df)
        if dur * 0.95 <= df[k] <= dur * 1.05
            push!(evt_ts, evt[k])
        elseif blank * 0.95 <= df[k] <= blank * 1.05
            nblank += 1
            push!(evt_off, evt[k])
        elseif !isempty(evt_ts) || nblank > 0
            # we have seen correct events and now we got one of odd duration
            if k < length(df)
                @warn("Some events that *should* be valid have invalid duration")
                warning = true
            end
        end
    end
    return evt_ts, warning, evt_off
end
# ============================================================================ #
function filter_events_by_duration(evt::Vector{<:Real},
            dur::AbstractFloat, nevt::Integer, limits::Tuple{<:Real,<:Real}=(-0.05, 0.05))

    warning = false

    if length(evt) < nevt
        @warn("Number of event timestamps is < the given number of events")
        warning = true
    end

    pad = 1.0 .+ limits
    df = diff(evt)
    # bgood = (df .>= dur*pad[1]) .& (df .<= dur*pad[2])
    bgood = (dur*pad[1]) .<= df .<= (dur*pad[2])
    ngood = sum(bgood)

    evt_off = Vector{Float64}(undef, 0)

    if ngood == nevt-1
        evt_ts = evt[[bgood; true]]

    elseif ngood == nevt
        evt_ts = evt[[bgood; false]]

    else
        # proportion of events of duration dur
        pdur = ngood ./ length(df)

        if pdur <= .1
            # few / no events are of length dur: events *MUST* only mark stim
            # onset
            total_dur = nanmedian(df)
            # buse = (df .>= total_dur*pad[1]) .& (df .<= total_dur*pad[2])
            buse = (total_dur*pad[1]) .<= df .<= (total_dur*pad[2])
            nuse = sum(buse)

            if nuse == nevt
                buse = [buse; false]
            elseif nuse == nevt-1
                buse = [buse; true]
            elseif nuse < nevt-1
                buse = [buse; true]
                @warn("Failed to locate the correct number of events")
                @show(length(evt), nevt, sum(buse))
                warning = true
            end

            evt_ts = evt[buse]

        elseif pdur >= .9
            if ngood + 1 >= nevt*2
                # b/c of the number of events, how we mark the final event
                # makes no difference
                push!(bgood, true)

                # most / all events are of length dur: events mark on and offset
                # of stim and blank_dur == stim_dur
                evt_ts = evt[bgood]

                #index of first stim on event (by definition precedes stim off)
                # buse = falses(length(evt_ts))
                # buse[1:2:(nevt*2)] .= true
                # evt_ts = evt_ts[buse]
                kon = 1:2:(nevt*2)
                evt_ts = evt_ts[kon]
                evt_off = evt_ts[kon .+ 1]

            elseif ngood > nevt
                # b/c of the number of events, how we mark the final event
                # makes no difference
                push!(bgood, true)

                # events mark stim on but we have extra events of approx. the
                # correct duration so just take the first nevt...
                @warn("Correct events could not be unambiguously identified...")
                @warn("falling back to heuristic method 1. Check the results!")
                evt_ts = evt[bgood]
                evt_ts = evt_ts[1:nevt]
                warning = true
            else
                error("Failed to locate the correct number of events")
            end

        elseif 0.4 <= pdur <= 0.6
            # approx half of events are of the correct length, so just go with that...
            @warn("Correct events could not be unambiguously identified...")
            @warn("falling back to heuristic method 2. Check the results!")
            warning = true

            # event duration is our only diagnostic info, so the final event
            # cannot mark a valid event onset
            push!(bgood, false)

            @show(length(evt), length(bgood), sum(bgood), nevt)
            evt_ts = evt[bgood]
            if length(evt_ts) > nevt
                evt_ts = evt_ts[1:nevt]
            end
        else
            error("Event durations are too inconsistent to estimate blank duration")
        end
    end

    # if warning
    #     @info("Requested dur: $(dur), median(diff(evt_ts)) = $(median(diff(evt_ts)))")
    # end

    return evt_ts, warning, evt_off
end
# ============================================================================ #
nanmedian(x::Vector{<:AbstractFloat}) = median(filter(!isnan, x))
# ============================================================================= #


text:
# Load Julia packages (libraries)

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

df = CSV.read(sr_datadir("chimpanzees.csv"), DataFrame);

# Define the Stan language model

stan10_1s = "
data{
    int N;
    int pulled_left[N];
}
parameters{
    real a;
}
model{
    real p;
    a ~ normal( 0 , 10 );
    pulled_left ~ binomial( 1 , inv_logit(a) );
}
";

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

m10_1s = SampleModel("m10.1s", stan10_1s);

# Input data for cmdstan

m10_1_data = Dict("N" => size(df, 1),
  "pulled_left" => df[!, :pulled_left]);

# Sample using cmdstan

rc10_1s = stan_sample(m10_1s, data=m10_1_data);

# Result rethinking

rethinking = "
  mean   sd  5.5% 94.5% n_eff Rhat
a 0.32 0.09 0.18  0.46   166    1
";

# Describe the draws

if success(rc10_1s)
    part10_1s = read_samples(m10_1s, :particles)
    part10_1s |> display
end


text:
module RingLists
export RingList, insertafter!, previous

import Base: ==, length, getindex, keys, haskey, insert!, eltype
import Base: Vector, show, hash, reverse, first, delete!, Set

struct RingList{T}
    data::Dict{T,T}
    revdata::Dict{T,T}
end

"""
`RingList{T}()` creates a new, empty `RingList` holding elements of type `T`.
Example: `RingList{Int}()`.

`RingList(list)` creates a new `RingList` from the elements in the
one-dimensional array `list`. Example: `RingList([1,2,3])`.

`RingList(x...)` creates a new `RingList` from the arguments.
Example: `RingList(1,2,3)`.
"""
function RingList(T::Type=Any)
    return RingList(Dict{T,T}(),Dict{T,T}())
end

function _rev(d::Dict{T,T}) where T
    rd = Dict{T,T}()
    for k in keys(d)
        v = d[k]
        rd[v] = k
    end
    return rd
end

function RingList(vals::Vector{T}) where T
    d = Dict{T,T}()
    n = length(vals)
    if n != length(unique(vals))  # check there are no repeats
        error("List of values may not have a repeat")
    end
    if n == 0
        return RingList(d,d)
    end
    if n==1
        a = vals[1]
        d[a] = a
        return RingList(d,d)
    end
    for i=1:n-1
        a = vals[i]
        b = vals[i+1]
        d[a] = b
    end
    a = vals[end]
    b = vals[1]
    d[a] = b
    return RingList(d,_rev(d))
end

function RingList(x...)
    return RingList(collect(x))
end

function RingList{T}() where T
    d = Dict{T,T}()
    return RingList(d,d)
end

==(a::RingList,b::RingList) = a.data == b.data
length(a::RingList) = length(a.data)
keys(a::RingList) = keys(a.data)
getindex(a::RingList, x) = a.data[x]

"""
`previous(a,x)` returns the element `y` so that `a[y]==x`.
"""
previous(a::RingList, x) = a.revdata[x]

haskey(a::RingList,x) = haskey(a.data,x)
eltype(a::RingList{T}) where T = T

"""
`first(a::RingList, true_first::Bool = false)` returns
an element of `a`. If `true_first` is `true`, then try to
return the smallest element held in `a`.
"""
function first(a::RingList, true_first::Bool = false)
    if true_first
        try
            x = minimum(keys(a))
            return x
        catch
        end
    end
    return first(a.data)[1]
end

"""
`insert!(a,x)` adds the element `x` into `a`.
No guarantee where it will go. See also
`insertafter!`
"""
function insert!(a::RingList{T},x::T) where T
    if length(a) == 0
        a.data[x] = x
        a.revdata[x] = x
        return nothing
    end
    if haskey(a,x)
        error("$x already in this RingList")
    end
    y = first(a)  # get the other elements
    if length(a) == 1
        a.data[x] = y
        a.data[y] = x
        a.revdata[x] = y
        a.revdata[y] = x
        return nothing
    end
    # y-->z   becomes y --> x --> z
    z = a[y]
    a.data[y] = x
    a.data[x] = z

    a.revdata[z] = x
    a.revdata[x] = y

    nothing
end

"""
`insertafter!(a::RingList,x,y)` inserts `x` into `a` after `y`.
"""
function insertafter!(a::RingList, x, y)
    if haskey(a,x)
        error("Element $x alread in this RingList")
    end
    if !haskey(a,y)
        error("Element $y not in this RingList, cannot insert after")
    end
    # who is currently after y?
    z = a[y]
    # we have y --> z
    # and change to y --> x --> z

    a.data[y] = x
    a.data[x] = z

    a.revdata[x] = y
    a.revdata[z] = x

    nothing
end

"""
`delete!(a,x)` removes `x` from the `RingList`.
"""
function delete!(a,x)
    if !haskey(a,x)
        error("$x not in the RingList")
    end
    if length(a)==1
        delete!(a.data,x)
        return
    end

    # prev --> x --> next
    next = a[x]
    prev = previous(a,x)
    delete!(a.data,x)
    delete!(a.revdata,x)
    a.data[prev] = next
    a.revdata[next] = prev
    return
end

function Vector(a::RingList{T}) where T
    n = length(a)
    if n == 0
        return T[]
    end
    K = keys(a)
    result = Vector{T}(undef,n)

    result[1] = first(a,true)

    i = 1
    while i<n
        i += 1
        result[i] = a[result[i-1]]
    end

    return result
end

"""
`reverse(a::RingList)` returns a new `RingList` containing the
same elements as `a` but in reverse order.
"""
reverse(a::RingList{T}) where T = RingList{T}(a.revdata,a.data)

Set(a::RingList{T}) where T = Set{T}(keys(a))

function show(io::IO, a::RingList{T}) where T
    v = Vector(a)
    result = "RingList{$T}("
    n = length(a)
    for i=1:n
        result *= "$(v[i])"
        if i<n
            result *= ","
        end
    end
    result *= ")"
    print(io, result)
end

hash(a::RingList) = hash(a.data)
hash(a::RingList, h::UInt64) = hash(a.data, h)

end # module


text:
@debug("test-proof-of-concept-linearmodel-readme.jl")
include("test-proof-of-concept-linearmodel-readme.jl")

@debug("test-proof-of-concept-linearmodel-extra.jl")
include("test-proof-of-concept-linearmodel-extra.jl")


text:
using SIMD

function chipsort_large(input::AbstractVector{T}, ::Val{V}, ::Val{K}) where {T,V,K}
    @inbounds begin
    J = div(length(input), V*K)
    chunks = reshape(valloc(T, div(32, sizeof(T)), V*J*K), V*J, K)
    output = valloc(T, div(32, sizeof(T)), V*J*K)
    chunks[:] .= input
    for k in 1:K
        chipsort_medium!((@view chunks[:,k]), Val(V), Val(V), Val(div(J,V)))
    end

    ds = VecMergeTree(chunks.parent, Val(V), Val(K))

    pout = pointer(output, 1)
    out = fill_tree!(ds, 1)
    vstorea(out, pout)
    pout += V*sizeof(T)
    for it in 1:(J*K-1)
        out = pop!(ds)
        vstorea(out, pout)
        pout += V*sizeof(T)
    end
    end
    output
end

struct VecMergeTree{A<:AbstractVector, B<:AbstractMatrix, C<:AbstractMatrix, T,V,K}   #{T,V,K,I}
    data::A
    cache::B
    ranges::C
    over::Vector{Bool}

    function VecMergeTree(data::AbstractVector{T}, ::Val{V},::Val{K}) where {T,V,K,I}
        cache = reshape(valloc(T, div(32, sizeof(T)), V*(2*K-1)), V, 2*K-1)
        J = div(length(data), V*K)
        ranges = copy(reshape([1:J:K*J; J:J:K*J+J-1],:,2)')
        over = zeros(Bool, 2*K-1)
        new{typeof(data), typeof(cache), typeof(ranges),T,V,K}(data, cache, ranges, over)
    end
end

function fill_tree!(ds::VecMergeTree{A,B,C,T,V,K}, k::Int64) where {A,B,C,T,V,K}
    if k < K
        ka = k<<1
        kb = ka+1
        va=fill_tree!(ds, ka)
        vb=fill_tree!(ds, kb)
        out,state = bitonic_merge(va,vb)
        vstore(state, pointer(ds.cache, 1+(k-1)*size(ds.cache,1)))
        out
    else
        datak = k-K+1
        vout = vloada(Vec{V,T}, pointer(ds.data, 1+(ds.ranges[1,datak]-1)*V))
        vstate = vloada(Vec{V,T}, pointer(ds.data, 1+(ds.ranges[1,datak])*V))
        vstorea(vstate, pointer(ds.cache, 1+(k-1)*size(ds.cache,1)))
        ds.ranges[1,datak] += 1
        vout
    end
end

@inline function pop!(ds::VecMergeTree{A,B,C,T,V,K}) where {A,B,C,T,V,K}
    @inbounds begin
    ## Find next leaf to merge and update
    k = 1
    ka = k<<1
    kb = ka+1
    while k < K && !(ds.over[ka] && ds.over[kb])
        k = if ds.over[kb] || !ds.over[ka] && ds.cache[1+V*(ka-1)] <= ds.cache[1+V*(kb-1)] ka else kb end
        ka = k<<1
        kb = ka+1
    end

    ## First vector to be merged descending the tree, from the leaves to root.
    out = vloada(Vec{V,T}, pointer(ds.cache, 1+(k-1)*V))

    @inbounds if k >= K && ds.ranges[1,k-K+1] < ds.ranges[2,k-K+1]
        datak = k-K+1
        vnew_state = vloada(Vec{V,T}, pointer(ds.data, 1+ds.ranges[1,datak]*V))
        ds.ranges[1,datak] += 1
        vstorea(vnew_state, pointer(ds.cache, 1+(k-1)*V))
    else  # implies `ds.over[ka] && ds.over[kb]` or ds.ranges[1,k-K+1] >= ds.ranges[2,k-K+1]
        ds.over[k] = true
    end

    k = k>>1
    @inbounds while k>0
        state = vloada(Vec{V,T}, pointer(ds.cache, 1+(k-1)*V))
        out, new_state = bitonic_merge(state, out)
        vstorea(new_state, pointer(ds.cache, 1+(k-1)*V))
        k = k>>1
    end
    end
    out
end


text:
################################################################################
# Functions required to generate data based on different models

"""
	genData(m::Model,nCons::Int,nProdTot::Int,nProdCons::Int; fname::String = "" )

Returns a simulated dataset for specified model.

### Arguments:
`m`:		Model specification. See `Model`. \\
`nCons`: 	# of consumers to simulate \\
`nProdTot`: 	# of alternatives in total (from which randomly chosen to display to consumer) \\
`nProdCons`: 	# of alternatives displayed to each consumer (incl. outside option). \\

### Optional keyword arguments:
`fname`: 	Filename to save in BSON format (can use directory path).  \\
"""
function genData(m::FullInfo,nCons::Int,nProdTot::Int,nProdCons::Int; fname::String = "" )
	Random.seed!(m.seed)

	# Consumer index into full data
	seli = [ 1+ii*nProdCons:(ii+1)*nProdCons for ii = 0:nCons-1]

	# Generate products
	pid,chars = genProducts(m::Model,seli,nProdTot,nProdCons)

	# Gather data in NamedTuple for convenience
	d = (seli = seli,pid = pid,chars = chars)

	# Get purchase indicators
	xβv = d.chars * m.β .+ rand(m.dV,size(d.chars,1))
	u   = xβv .+ rand(m.dE,size(d.chars,1))
	purch = genSearchPaths(m,d,u)

	# Gather again data
	d = (seli = seli,pid = pid,chars = chars, purch = purch)

	if !isempty(fname)
		@info "Saved data in bson format to $fname."
		@save fname d
	end

	return d
end

###################################################
function genSearchPaths(m::FullInfo,d,u)
	n = length(d.seli)

	# Same for all consumers
	purch = fill(0,n)

  @inbounds Threads.@threads for k = 1:length(d.seli)
		for ii in getrange(n) # Loop over consumers allocated to specific thread
			purch[ii] = genSearchPathi!(m::FullInfo,d.seli[ii],u)
		end
	end
	return purch
end

function genSearchPathi!(m::FullInfo,seli,u)
	_,purch = findmax_range(u,seli)
	return purch -seli[1] + 1  # adjustment as findmax_range gives index in whole array
end

###################################################
function negLogLik(m::FullInfo,β::Vector{T},d; e = [], v = [] ) where T<:Real

	# Draw random numbers if not provided
	if isempty(e)
		Random.seed!(m.seed+1)
		e = rand(m.dE,size(d.chars,1),m.options.nDraws)
	end
	if isempty(v)
		Random.seed!(m.seed+298)
		v = rand(m.dV,size(d.chars,1),m.options.nDraws)
	end

	# Form utility
	u = zeros(T,size(d.chars,1),m.options.nDraws)
	@views fillU!(u,d,β,e,v)

	# Loop through consumers (multi-threaded)
	logLik = zeros(eltype(β), Threads.nthreads())

	@inbounds Threads.@threads for k = 1:Threads.nthreads()
		local locSumLL = 0.0
		for i in getrange(length(d.seli))
			r = d.seli[i] # selects consumers
			locV = @views calcVFI(u,r,r[d.purch[i]],m.smo)
			locSumLL += log(locV)
		end
		logLik[k] = locSumLL
	end
	sum_logLik = abs(sum(logLik)) == Inf ? -1e10 : sum(logLik)
	return -sum_logLik
end

function fillU!(u,d,β,e,v)
	nChar = size(d.chars,2)
	@inbounds Threads.@threads for i in 1:size(u,1)
		for k in 1: size(u,2)
		u[i,k] = e[i,k] + v[i,k]
			for j in 1:nChar
				u[i,k] += d.chars[i,j] * β[j]
			end
		end
	end
end

function calcVFI(u,r,ind_p,smo)
	v = zero(eltype(u))
	n = size(u,2)
	@inbounds for ee in 1:n
		v2 = zero(eltype(u))
		for j in r
				v2 += exp(-smo*(u[ind_p,ee]-u[j,ee]))
		end
		v += 1 / (1 + v2) / n
	end
	return v
end

###################################################

"""
	estimateModel(m::Model,d; startvals = [], pars_true = [])

Estimates model `m`. Returns tuple containing parameter estimates and full output from optimizer.

### Arguments:
`m`:		Model specification. See `Model`.  \\
`d`: 		Data (see structure in a dataset generated through `genData`.) \\

### Optional keyword arguments:
`startvals`: 	Vector of starting values for estimation. \\
`pars_true`: 	Only used for printing (useful when known).  \\
"""
function estimateModel(m::FullInfo,d; startvals = [] , pars_true = []  )

	# Preliminaries
	nBeta = size(d.chars,2)

	if isempty(pars_true)
		pars_true = fill(NaN,nBeta)
	end

	if isempty(startvals)
		startvals = rand(Uniform(),nBeta)
	end

	# Random draws, pre-allocate
	Random.seed!(m.seed+651168) # ensures draws not same as in dataGen
	e = rand(m.dE,size(d.chars,1),m.options.nDraws)
	v = rand(m.dV,size(d.chars,1),m.options.nDraws)

		############################################################################
	# Optimization

		# Define objective function
	obj_fun(θ,p) = negLogLik(m,θ,d;v=v,e=e)
	obj = OptimizationFunction(obj_fun,m.options.autodiff)
	prob = OptimizationProblem(obj,startvals)

		# Run optimization
	res = solve(prob,m.options.algo;
				show_trace = true , show_every=m.options.showEvery,
				maxiters = m.options.maxIter,
				x_tol = m.options.x_tol,
				f_tol = m.options.f_tol,
				g_tol = m.options.g_tol)

	pars_opt = res.minimizer
	fval_opt = res.minimum

	# Calculate obj. fun at start and true val
	startVal = obj_fun(startvals,[])
	trueVal = isnan(pars_true[1]) ? NaN : obj_fun(pars_true,[])

	# Extract results
	β = pars_opt
	finalVal = fval_opt

	infotxt = "
	####################################
	# Results of Estimation:
	Estimation model:	$(typeof(m))
	Final Objective Value: 		$finalVal
	Start Objective Value:		$startVal
	True Objective Value:		$trueVal
	Parameters:
		β  = $(round.(β,digits=2))
		β0 = $(round.(pars_true[1:nBeta],digits=2))
	#################################### "

	@info infotxt

	return β,finalVal

end

"""
	calcWelfareDemand(m::Model,d; e = [] , v = [] )

Simulates welfare and demand for a model given data.

### Arguments:
`m`:		Model specification. See `Model`.  \\
`d`: 		Data (see structure in a dataset generated through `genData`.)
### Optional keyword arguments:
`e`: 	Matrix of ϵ error terms (draws along columns). \\
`v`: 	Matrix of ν error terms (draws along columns).   \\
"""
function calcWelfareDemand(m::FullInfo,d; e = [] , v = [] )

	# Draw random numbers if not provided
	if isempty(e)
		Random.seed!(m.seed+43)
		e = rand(m.dE,size(d.chars,1),m.options.nDraws)
	end
	if isempty(v)
		Random.seed!(m.seed+323)
		v = rand(m.dV,size(d.chars,1),m.options.nDraws)
	end

	# Derived values
	nCons = length(d.seli)

	# Form utility
	xβv = d.chars * m.β .+ v
	u   = xβv .+ e

	# Pre allocate
	sumU = zeros(Threads.nthreads())
	demand = zeros(Int64,length(d.pid[1]),Threads.nthreads())
	@inbounds Threads.@threads for k = 1:Threads.nthreads()
		threadRange = getrange(length(d.seli))
		locU = 0.0
		locD = zeros(Int64,length(d.pid[1]))
		for ii in threadRange, dd in 1:size(e,2)
			@views r 	= d.seli[ii]
			@views ud 	= u[:,dd]
			purchi 	= genSearchPathi!(m,r,ud)
			locU 	+= ud[r[purchi]]
			locD[purchi] += 1
		end
		sumU[k] = locU / size(e,2)
		demand[:,k] .= locD
	end
	return sum(sumU), sum(sumU), 0.0, sum(demand,dims=2) ./ size(e,2) ./ nCons
end


text:
"""
    PELT(segment_cost, n[, β = log(n)])

Runs the PELT algorithm using cost function `segment_cost` for a time series of length `n` and a penalty `β` and returns the position of changepoints and optimal segmentation cost.

See also: [`@PELT`](@ref), [`@segment_cost`](@ref), [`CROPS`](@ref)

# Returns
* `(CP::Vector{Int}, cost::Float64)`:
  * `CP::Vector{Int}`: Vector of indices of detected changepoints
  * `cost::Float64`: Cost of optimal segmentation

# Example
```julia-repl
n, λ = 1000, 100
μ, σ = Normal(0.0, 10.0), 1.0
sample, cps = @changepoint_sampler n λ Normal(μ, σ)  # Sample changepoints
seg_cost = NormalMeanChange(sample, σ)               # Create segment cost function
pelt_cps, pelt_cost = PELT(seg_cost, n)              # Run PELT
```

# References
Killick, R., Fearnhead, P. and Eckley, I.A. (2012) Optimal detection of changepoints with a linear computational cost, JASA 107(500), 1590-1598
"""
function PELT( segment_cost::Function , n::Int; pen::Float64 = log(n) )

    # F[t] is optimal cost of segmentation upto time t
    F = Array{Float64}(undef, n+1)
    F[1] = -pen
    F[2] = 0

    # last chpt prior to time t
    chpts =  Array{Int64}(undef, n)
    chpts[1] = 0

    # vector of candidate chpts at t
    R = Int64[0]

    for t in 2:n
        cpt_cands = R
        seg_costs = Array{Float64}(undef, length(cpt_cands))
        for i in 1:length(cpt_cands)
            seg_costs[i] = segment_cost(cpt_cands[i], t)
        end

        F[t+1] , tau = findmin( F[cpt_cands .+ 1] .+ seg_costs .+ pen )
        chpts[t] = cpt_cands[tau]

        # pruning step
        ineq_prune = (F[cpt_cands .+ 1] .+ seg_costs) .< F[t+1]
        R = push!(cpt_cands[ineq_prune], t - 1)
    end

    # get changepoints
    CP = Array{Int64}(undef,0)
    last = chpts[n]
    push!(CP,last)
    while last > 0
        last = chpts[last]
        push!(CP,last)
    end
    sort!(CP)

    return CP , F[n+1]

end


text:
using LinearAlgebra
"""
Quasi-Newton Methods for Solving Systems of Nonlinear Equations
@author Oscar Veliz
"""

"""
    F(X)

System of nonlinear equations on vector X ([x;y]).
Overwrite for your own testing, currently set to:
x^2 - y - 1
x - y^2 + 1

# Example
```julia-repl
julia> F([1.0; 2.0])
[-2.0, -2.0]
```
"""
F(X) = [X[1]^2 - X[2] - 1; X[1] - X[2]^2 + 1]

"""
    Secant(X, eps=10^-6)

Wolfe-Bittner Secant Method for the system of nonlinear equations in [`F`](@ref).
Requires `d+1` starting vectors in X.
# Example
```julia-repl
julia> Secant(transpose([1.0 1.0;1.0 2.0;1.5 2.0]))
2-element Array{Float64,1}:
 1.6180339736070495
 1.6180339507485015
```
"""
function Secant(X, eps::Float64 = 10^-6)
    n = size(X,1)
    FN = mapslices(F, X; dims=1)
    FX = FN[:,end]
    #println(X) #uncomment to see starting points
    while norm(FX) > eps
        A = [ones(1,n+1); FN]
        b = zeros(n+1,1)
        b[1,1] = 1
        p = A \ b
        Xbar = X * p
        #print(Xbar) #uncomment to see each itr
        X = hcat(X, Xbar)
        X = X[1:end,2:end]
        FX = F(Xbar)
        FN = hcat(FN,FX)
        FN = FN[1:end,2:end]
    end
    X[:,end]
end

"""
    steffJ(X, FX = F(X), n = size(X,1))

Traub-Steffensen approximation of the Jacobian.
Requires starting vector X.
FX is optional in case F(X) was already computed for efficiency.
# Example
```julia-repl
julia> steffJ([1.5;2.0])
2×2 Array{Float64,2}:
 2.25  -1.0
 1.0   -2.5
```
"""
function steffJ(X,FX = F(X),n = size(X,1))
    H = Diagonal(FX)
    J = X #placeholder to set dims
    for i = 1:n
        J = hcat(J, F(X + H[:,i]) - FX)
    end
    J = J[1:end, 2:end] * inv(H) #removes placeholder before mult
end

"""
    Steffensen(X, eps=10^-6)

Traub-Steffensen Method for a system of nonlinear equations in [`F`](@ref).
Requires starting vector X.
# Example
```julia-repl
julia> Steffensen([1.5;2.0])
2-element Array{Float64,1}:
 1.6180339736216265
 1.6180339387805316
```
"""
function Steffensen(X, eps::Float64 = 10^-6)
    n = size(X,1)
    #println(X) #uncomment to see itr
    FX = F(X)
    while(norm(FX) > eps)
        J = steffJ(X,FX,n)
        X = X - J \ FX
        #println(X) #uncomment to see itr
        FX = F(X)
    end
    X
end

"""
    fdJ(X, FX = F(X), h = 10^-4, n = size(X,1))

(Forward) Finite Difference approximation of the Jacobian.
Requires starting vector X. Optional h and n.
FX is optional in case F(X) was already computed for efficiency.
# Example
```julia-repl
julia> fdJ([1.5;2.0])
2×2 Array{Float64,2}:
 3.0001  -1.0
 1.0     -4.0001
```
"""
function fdJ(X,FX = F(X),h::Float64 = 10^-4,n = size(X,1))
    H = h * I(n)
    J = X #placeholder to set dims
    for i = 1:n
        J = hcat(J, F(X + H[:,i]) - FX)
    end
    J = h^-1 * J[1:end, 2:end] #removes placeholder and divide by h
end

"""
    FiniteDiff(X, h = 10^-4, eps = 10^-6)

Generalized Finite-Difference Method for a system of nonlinear equations in [`F`](@ref).
Requires starting vector X. Optional h and eps.
# Example
```julia-repl
julia> FiniteDiff([1.5;2.0])
2-element Array{Float64,1}:
 1.6180340699246842
 1.6180341416981485
```
"""
function FiniteDiff(X,h::Float64 = 10^-4, eps::Float64 = 10^-6)
    n = size(X,1)
    #println(X) #uncomment to see itr
    FX = F(X)
    while(norm(FX) > eps)
        J = fdJ(X,FX,h,n)
        X = X - J \ FX
        #println(X) #uncomment to see itr
        FX = F(X)
    end
    X
end

"""
    Broyden(X, eps = 10^-6)

"Good" Broyden Method for a system of nonlinear equations in [`F`](@ref).
Requires starting vector X. Optional eps. Initial J approximated by differencing.
# Example
```julia-repl
julia> Broyden([1.5;2.0])
2-element Array{Float64,1}:
 1.6180340721138766
 1.6180340507162783
```
"""
function Broyden(X, eps::Float64 = 10^-6)
    #println(X) #uncomment to see itr
    FX = F(X)
    invJ = inv(fdJ(X,FX))
    while(norm(FX) > eps)
        Xold = X
        X = X - invJ*FX
        #println(X) #uncomment to see itr
        deltaX = X - Xold
        FXold = FX
        FX = F(X)
        deltaF = FX - FXold
        trans = transpose(deltaX)
        invJ = invJ + (deltaX - invJ*deltaF)/(trans*invJ*deltaF)*trans*invJ
    end
    X
end

"""
Main
"""

println(Secant(transpose([1.0 1.0;1.0 2.0;1.5 2.0])))
println(Steffensen([1.5;2.0]))
println(FiniteDiff([1.5;2.0]))
println(Broyden([1.5;2.0]))

max_stars_repo_name: deboradeben/numerical-veliz

max_issues_repo_name: deboradeben/numerical-veliz

max_forks_repo_name: deboradeben/numerical-veliz


text:
using ExamplePackage
using Base.Test

# The testing setup
# Here we split different categories into tests
# Each category is a separate file
# By placing `@testset` on the file, it will run all of the tests there,
# and report back when the entire file is complete

# Optional: If you @testset begin ... end around the entire testing setup
# then all of your tests will run, and the failures will be counted
# and printed at the end. There is a balance between the completeness of
# the test results and the time it takes to run tests!

# Optional: tic() ... toc() around the tests as a quick way to keep track of
# large performance regressions.

# @testset begin
@testset "Output String Tests" begin include("output_string_tests.jl") end
# end


text:
"""
    secant(x0, x1, f, e, n)

Compute the root of the given function using Secant method, an Extension of Newton's method.
Insted of using a tangent line, which requres knoledge of the function's derivative, Secant method uses two values of
the initial approximation `x0` and `x1` to find secant line and uses the intersection of the secant line and x-axis
to find the next approximation.

    secant(x0::Float64, x1::Float64, f::Function, e::Float64, n::Integer)

### OUTPUT:
  * `p` - approximation of the root of the function `f` (Float64) or,
  * `-1` - if `n` the maximum number of iterations was exceeded.
"""
function secant(x0::Float64, x1::Float64, f::Function, e::Float64, n::Integer)

    f0 = f(x0)
    f1 = f(x1)

    for i in 1:n

        x = x1 - f1*(x1-x0)/(f1 - f0)
        fx = f(x)
        x0 = x1
        x1 = x
        f0 = f1
        f1 = fx

        if(abs(x1 - x0) < e)
            return x1
        end
    end

    return -1
end


text:
using Statistics
using LinearAlgebra

function ES(layers::Array, num_generations::Int, f::Function, npop::Int=50, sigma=0.5, alpha=0.01)
    µ = Individual(layers::Array) # Center
    n = length(µ.genes)
    evaluate!(µ, f)

    expert = µ # Best individual so far

    fits = zeros(num_generations)

    for i in eachindex(fits)
        population = Array{Individual}(undef, 𝜆)
        best = 1
        for j in eachindex(population)
            new_genes = µ.genes + randn(n) .* sigma
            population[j] = Individual(layers, new_genes)
            evaluate!(population[j], f)
            if population[j].fitness > population[best].fitness
                best = j
            end
        end

        if population[best].fitness >= expert.fitness
            expert = population[best]
        end

        fits[i] = expert.fitness

        u = max.(0, log(n/2+1) .- log.(1:n))
        u .= u ./ sum(u) .- 1/n

        # Sort from highest to lowest fitness
        idx = sortperm(population, by=x -> x.fitness, rev=true)

        step = zeros(n)
        for i in eachindex(population)
            step = step.+ u[i] .* population[idx[i]].genes
        end

        µ.genes = µ.genes .+ step .* alpha

        print(i, "/",  num_generations, " > ", expert.fitness, "\t\t\r")
    end
    fits, expert
end