text:
using PSCOPF

using Test
using Dates
using DataStructures

@testset "test_verify_firmness" begin

    hollow_uncertain_value = PSCOPF.UncertainValue{Int64}(missing,
                                                SortedDict("S1"=>10,"S2"=>missing))
    non_firm_uncertain_value = PSCOPF.UncertainValue{Int64}(missing,
                                                SortedDict("S1"=>10,"S2"=>15))
    same_uncertain_value = PSCOPF.UncertainValue{Int64}(missing,
                                                SortedDict("S1"=>10,"S2"=>10))
    definitive_uncertain_value = PSCOPF.UncertainValue{Int64}(10,
                                                SortedDict("S1"=>10,"S2"=>10))

    @testset "verifies_on_generator_schedule" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE,
            Dates.DateTime("2015-01-01T11:30:00") => PSCOPF.FREE,
            Dates.DateTime("2015-01-01T12:00:00") => PSCOPF.TO_DECIDE,
            Dates.DateTime("2015-01-01T13:00:00") => PSCOPF.DECIDED )

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => non_firm_uncertain_value,
            Dates.DateTime("2015-01-01T11:30:00") => same_uncertain_value,
            Dates.DateTime("2015-01-01T12:00:00") => definitive_uncertain_value,
            Dates.DateTime("2015-01-01T13:00:00") => definitive_uncertain_value )

        @test PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "violates_on_generator_schedule" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE,
            Dates.DateTime("2015-01-01T11:30:00") => PSCOPF.FREE,
            Dates.DateTime("2015-01-01T12:00:00") => PSCOPF.TO_DECIDE,
            Dates.DateTime("2015-01-01T13:00:00") => PSCOPF.DECIDED )

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => non_firm_uncertain_value,
            Dates.DateTime("2015-01-01T11:30:00") => same_uncertain_value,
            Dates.DateTime("2015-01-01T12:00:00") => non_firm_uncertain_value, # the problematic value
            Dates.DateTime("2015-01-01T13:00:00") => definitive_uncertain_value )

        @test !PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "missing_schedule_value" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T10:00:00") => PSCOPF.FREE, #No corresponding value in schedule
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => hollow_uncertain_value)

        @test !PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "ignored_schedule_values" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T10:00:00") => hollow_uncertain_value, # No firmness constraint
            Dates.DateTime("2015-01-01T11:00:00") => non_firm_uncertain_value)

        @test PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "free_with_hollow" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => hollow_uncertain_value)

        @test !PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "free_with_non_firm" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => non_firm_uncertain_value)

        @test PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "free_with_same" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => same_uncertain_value)

        @test PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "free_with_definitive" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => definitive_uncertain_value)

        @test PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "decided_with_hollow" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => hollow_uncertain_value)

        @test !PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "decided_with_non_firm" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => non_firm_uncertain_value)

        @test !PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "decided_with_same" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => same_uncertain_value)

        @test !PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "decided_with_definitive" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => definitive_uncertain_value)

        @test PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "to_decided_with_hollow" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.TO_DECIDE)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => hollow_uncertain_value)

        @test !PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "to_decided_with_non_firm" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.TO_DECIDE)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => non_firm_uncertain_value)

        @test !PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "to_decided_with_same" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.TO_DECIDE)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => same_uncertain_value)

        @test !PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

    @testset "to_decided_with_definitive" begin
        generator_firmness = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.TO_DECIDE)

        scheduled_values = SortedDict(
            Dates.DateTime("2015-01-01T11:00:00") => definitive_uncertain_value)

        @test PSCOPF.verify_firmness(generator_firmness, scheduled_values)
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
"""
    struct sDAMMmat
Fields required to plot surface of DAMM

# Fields
$(TYPEDFIELDS)
"""
struct sDAMMmat
  porosity::Float64
  params::NTuple{6, Float64}
  x::Array{Float64, 1}
  y::Array{Float64, 1}
  DAMM_Matrix::Array{Float64, 2}
end

"""
    struct sDAMMmatq
Fields required to plot surface of DAMM

# Fields
$(TYPEDFIELDS)
"""
struct sDAMMmatq
  porosity::Float64
  Tmed::Array{Float64, 1}
  θmed::Array{Float64, 1}
  Rmed::Array{Float64, 1}
  params::NTuple{6, Float64}
  x::Array{Float64, 1}
  y::Array{Float64, 1}
  DAMM_Matrix::Array{Float64, 2}
end

#= IF LsqFit.jl allows params to be a struct, then use:
struct DAMMparams
  αₛₓ::Float64 # Pre-exponential factor
  Eaₛₓ::Float64 # Activation energy
  kMₛₓ::Float64 # Michaelis constant for substrate
  kMₒ₂::Float64 # Michaelis constant for O₂
  porosity::Float64 # Soil bulk density / soil particle density
  Sxₜₒₜ::Float64 # Total soil carbon
end

#IF LsqFit.jl allows params to be a Dict, then use:
function DAMMparams(p::NTuple{6, Float64})
  p = Dict([:αₛₓ, :Eaₛₓ, :kMₛₓ, :kMₒ₂, :porosity, :Sxₜₒₜ]
	   .=> [p[1], p[2], p[3], p[4], p[5], p[6]])
end

function DAMMparams_u(p::Dict{Symbol, Float64}) # give unit to params
  αₛₓ = p[:αₛₓ]mg*cm^-3*hr^-1
  Eaₛₓ = p[:Eaₛₓ]kJ*mol^-1
  kMₛₓ = p[:kMₛₓ]g*cm^-3
  kMₒ₂ = p[:kMₒ₂]L*L^-1
  porosity = p[:porosity]
  Sxₜₒₜ = p[:Sxₜₒₜ]g*cm^-3
  P = Dict([:αₛₓ, :Eaₛₓ, :kMₛₓ, :kMₒ₂, :porosity, :Sxₜₒₜ]
	   .=> [αₛₓ, Eaₛₓ, kMₛₓ, kMₒ₂, porosity, Sxₜₒₜ])
  return P
end
=#

ext: jl


text:
# Double Integrator
T = Float64
model = TrajectoryOptimization.Dynamics.doubleintegrator
model_d = rk3(model)
n = model.n; m = model.m

# costs
Q = 1.0*Diagonal(I,n)
Qf = 1.0*Diagonal(I,n)
R = 1.0e-1*Diagonal(I,m)
x0 = [0; 0.]
xf = [1.; 0] # (ie, swing up)

N = 21
dt = 0.1
U0 = [0.001*rand(m) for k = 1:N-1]
obj = TrajectoryOptimization.LQRObjective(Q,R,Qf,xf,N)

u_max = 1.5
u_min = -1.5

bnd = BoundConstraint(n,m,u_max=u_max, u_min=u_min,trim=true)
goal = goal_constraint(xf)
constraints = Constraints(N)
for k = 1:N-1
    constraints[k] += bnd
end
constraints[N] += goal

doubleintegrator = TrajectoryOptimization.Problem(model_d, obj, constraints=constraints, x0=x0, xf=xf, N=N, dt=dt)
initial_controls!(doubleintegrator, U0)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
#!/usr/bin/env julia

# print out a page full of the same glyph in different fonts
# resizing each one to be the same width/height

# hardcoded list of fonts, handpicked for ampersandy goodness
fontslist = [
"Aachen-Bold", "AcademyEngravedLetPlain", "ACaslonPro-Semibold", "AGaramondPro-Bold", "Agenda-Medium", "AkzidenzGroteskBQ-BdCnd", "AlbertusMTStd", "AlgerianCondensedLET", "Americana", "AmericanTypewriterStd-Light", "AmplitudeCond-Light", "AnnaStd", "AnonymousPro", "AntiqueOliveStd-Roman", "ArcherPro-Bold", "ArialUnicodeMS", "ArnoldBoecklin", "AvantGarde-Book", "Avenir-Roman", "BankGothic-Light", "BaseNineB", "BaseTwelveSansB", "Baskerville", "BebasNeue", "BellCentennialStd-BdListing", "BellGothicStd-Black", "Bello-Script", "BelweStd-Medium", "BemboStd", "BenguiatStd-Book", "BenguiatGothic-BookOblique", "BernhardModernStd-Bold", "BickhamScriptPro-Regular", "BigCaslon-Medium", "BitstreamVeraSansMono-Bold", "Bodoni-Poster", "Bodoni-PosterCompressed", "BodoniStd", "BrandonGrotesque-Thin", "BrushScript", "CaeciliaLTStd-Bold", "Campanile", "CasablancaLHF", "Caslon3LTStd-Italic", "Caslon224Std-Black", "CaslonTwoTwentyFour-Black", "CenturyGothic", "CenturyOldStyle-Bold", "Chalkboard", "Champion-HTF-Bantamweight", "Charcoal", "CheltenhamStd-Book", "ChocICG", "ChopinScript", "ChronicleTextG1-Roman", "CircularStd-Medium", "Cochin-BoldItalic", "Code2000", "CooperBlackStd", "CopperplateGothic-Bold", "CorpidOffice-Bold", "Courier", "CutiveMono-Regular", "DejaVuSans", "DINNextLTPro-Regular", "DINPro-Regular", "DomCasualStd", "DorchesterScriptMTStd", "DTLFleischmann-DRegular", "DTLNobelT-Bold", "DTLProkyonSTBold", "Duvall", "EncorpadaClassic-Black", "EnvyCodeR", "EurostileLTStd", "EversonMono", "FagoOfficeSerif-Bold", "FantasqueSansMono-RegItalic", "FilosofiaGrandOT-Bold", "FiraMonoOT", "FiraSansOT", "Formata-Bold", "Forza-Medium", "FranklinGothicStd-Roman", "FrizQuadrataStd", "FrutigerLTStd-Roman", "FuturaStd-Book", "FuturaStd-Condensed", "GadgetRegular", "Geneva", "Georgia", "GillSansMT", "GillSansStd-UltraBold", "GothamBlack", "GothamBook", "GothamCondensed-Book", "GoudyModernMTStd", "GoudySansStd-Black", "GoudyStd", "GoudyStout", "Helvetica", "Helvetica-Neue", "Hermit-bold", "HoboStd", "HouseGothicHG23Ext-BOLD1", "Impact", "Inconsolata-Bold", "IndustriaLTStd-Solid", "InputMono-Black", "InputSans-Black", "InputSerif-Black", "Interstate", "ITCAvantGardeStd-Bk", "ITCFranklinGothicStd-Med", "ITCKabelStd-Book", "JoannaMTStd-Bold", "Knockout-HTF51-Middleweight", "Knox", "Korinna-Bold", "Kosmik-BoldOne", "LaserICG", "Lato-Regular", "Lavoisier", "LeagueGothic", "LegacySansStd-Bold", "LetterGothicMTStd", "LexiconNo1ItalicD-Txt", "LHFAmbrosia", "LHFLarcherRoman", "LiberationMono", "LucidaCalligraphy-Italic", "LucidaSansUnicode", "Luculent", "Lust-Regular", "MatrixInlineExtraBold", "MeridienLTStd-Bold", "MetaPro-Book", "Miller-Text", "MiltonMN", "MinionPro-Regular", "MissionScript", "MistralStd", "Monaco", "Monofonto", "Monofur", "Motter-Ombra", "MotterCorpusStd-Condensed", "MrsEavesOT-Roman", "MuseoSlab-500", "NeutraText-Book", "NeuzeitSLTStd-Book", "NewsGothicStd", "NovaMono", "OCRA", "OCRB", "OpenSans", "Optima-Regular", "Osaka", "Osaka-Mono", "OxygenMono-Regular", "P22BayerUniversal", "P22CezannePro", "P22JohnstonUnderground", "PalaceScriptMT-Std", "Palatino-Roman", "PalatinoLinotype-Bold", "PeignotLTStd-Demi", "Perla", "PerpetuaStd", "Pistilli", "PosterCut-Regular", "PrestigeEliteStd", "ProggySmallTT", "ProggyTinyTT", "ProximaNova-Regular", "PTMono-Bold", "Reina12Standard", "ReinaEngravedPro", "Roboto-Medium", "RockwellStd-Bold", "RomicStd-Medium", "RunicMT-Condensed", "Scala-Italic", "SegoeUI", "Sentinel-Black", "SerifGothicStd", "Share-TechMono", "SinhalaMN-Bold", "SnellRoundhand-BlackScript", "Sofachrome", "SourceCodePro-Black", "SouvenirStd-Light", "SteepTypewriter-Medium", "Stilla", "StoneSansStd-Medium", "StoneSerifStd-Medium", "Superclarendon-Black", "SweetSansPro-Heavy", "TechnoRegular", "Telex-Regular", "Textile", "TiffanyStd-Demi", "TiffanyStd-Heavy", "TimesNewRomanMTStd", "TradeGothicLTStd", "TrajanPro-Regular", "TriniteNo1-Bold", "Triplex-Bold", "TriplexSerif-Bold", "Trixie-Text", "Tungsten-Medium", "Ubuntu", "UbuntuMono-Bold", "Unispace", "UniversityRomanStd", "UniversLTStd", "UtopiaStd-Regular", "VAGRoundedStd-Bold", "VendomeICG", "Verdana", "Verlag-Book", "Verve", "Vitesse-Book", "Vitesse-Medium", "Whitney-Light", "WilhelmKlingsporGotLTStd", "WindsorBT-Roman", "Zapfino"
]

using Luxor, Colors

function editfontname(ff1)
    ff2 = replace(ff1, "Std", "")
    ff2 = replace(ff2, "-", " ")
    ff2 = replace(ff2, r"MT|LT|TT", "")
    ff2 = replace(ff2, "Bd", "Bold")
    ff2 = replace(ff2, "Bk", "Book")
    ff2 = replace(ff2, r"([A-z])(Pro)", s"\1 \2")
    ff2 = replace(ff2, "Regular", "")
    return ff2
end

function heart()
    move(127,1) # bother, it's offset from 0/0
    curve(134.2, -6.5, 134.2, -6.5, 156.1, -29.6)
    curve(185.8, -60.5, 198.1, -74.3, 213.7, -95.3)
    curve(240.4, -131, 253.3, -163.7, 253.3, -194.9)
    curve(253.3, -218, 246.4, -237.8, 232.6, -253.7)
    curve(219.1, -268.7, 204.1, -275.3, 181.9, -275.3)
    curve(154, -275.3, 136.3, -265.1, 127, -243.8)
    curve(124, -252.5, 120.4, -257.6, 112.9, -263.6)
    curve(103.6, -270.8, 88.3, -275.3, 73.3, -275.3)
    curve(59.2, -275.3, 46, -271.4, 35.2, -264.5)
    curve(14.5, -250.7, 1, -223.4, 1, -194.6)
    curve(1, -167.3, 13, -136.4, 37.3, -101)
    curve(53.8, -77, 86.5, -39.8, 127, 1)
    closepath()
end

function placefonts(fontlist, str, background="white")
    global x, y, xstep, ystep, hmargin, vmargin, rowspacing, oversizefonts
    for font in fontlist
        fontface(font)
        # work out size at which every letter is the same height
        fsize = 50
        te = []
        while true
            fontsize(fsize)
            te = textextents(str) # [xb, yb, width, height, xadvance, yadvance]
            if te[3] == 0 break end # some fonts don't  have certain characters...
            if te[4] > ystep
                break
            else
                fsize += 1
            end
        end
        # draw character at the size we found
        fontsize(fsize)
        textcentred(str, x, y)
        # underneath, draw its name
        fontface("BellCentennialStd-Address")
        fontsize(6)
        fontnamemodified = editfontname(font)
        # in case glyph overlaps with text, put a semitransparent box behind the text (yuk)
        gsave()
        sethue(background)
        setopacity(0.5)
        te = textextents(fontnamemodified)
        rect(x - te[3]/2 - te[1], y + te[4], te[3] + te[1], -te[4] - te[2] + 12, :fill)
        grestore()
        # show font name
        textcentred(fontnamemodified, x, y + 12)
        # room for another one?
        nextxpos = x + te[3]/2 + xstep
        if nextxpos > (pagewidth - hmargin - xstep - xstep)
            x = hmargin
            y += ystep + rowspacing
        else
            x += xstep + rowspacing
        end
    end
end

const pagewidth = 1190 # A2 size paper, points
const pageheight = 1684 # A2 size paper, points

function main(str)
    global x, y, xstep, ystep, hmargin, vmargin, rowspacing, oversizefonts
    filename = replace("i-love-" * str * "-fonts.pdf", r"[^A-z0-9.]", "")
    Drawing(pagewidth, pageheight, "/tmp/$(filename)")
    # background and title
    background = "mintcream"
    sethue(background)
    rect(0,0, pagewidth, pageheight, :fill)
    sethue(0,0,0)
    fontface("Perla")
    fontsize(60)
    text("I ", (pagewidth/2) - 295, 115)
    text(" ampersands!", (pagewidth/2) - 210, 115)
    gsave()
    translate((pagewidth/2) - 255, 115)
    scale(.175, .175)
    heart()
    sethue("red")
    fillstroke()
    grestore()

    sethue(0,0,0)
    gsave()
    translate(70,200) # start below top, point size is up

    # row is y, increasing down from the top
    # col is x, increasing right from the left

    xstep     = 60
    ystep     = 46
    hmargin   = 20
    vmargin   = 20
    x         = hmargin
    y         = vmargin
    rowspacing = 30 # horizontal column gap
    placefonts(fontslist, str, background)
    grestore()
    # decoration at the bottom
    gsave()
    translate(pagewidth/2, pageheight - vmargin - 12)
    setline(0.2)
    x_vals = -14 * 2 * pi:pi/200: 14 * 2 * pi
    pl = [Point(d * 2 * pi , -sin(d) * cos(12 * d) * 8 * sin(d/10)) for d in x_vals]
    poly(simplify(pl, 0.05), close=false, :stroke)
    grestore()
    # border decoration around the edges
    setline(1.6)
    translate((pagewidth/2),(pageheight/2))
    rect(-(pagewidth/2) + 6, -(pageheight/2) + 6, pagewidth - 12, pageheight - 12, :stroke)
    setline(.5)
    rect(-(pagewidth/2) + 10, -(pageheight/2) + 10, pagewidth - 20, pageheight - 20, :stroke)
    finish()
    preview()
end

main("&")


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

using Markdown
using InteractiveUtils

# ╔═╡ 0d696b90-14b9-11eb-19c6-7110e9eded9e
using PotentialFlow

# ╔═╡ 40e53300-14b9-11eb-28ca-c39ecf8c05ef
using Plots

# ╔═╡ 59e4bc40-14b9-11eb-3670-cdc8fa8f0ded
begin
	X = range(-1, 1, length=21)
	Y = range(-1, 1, length=21)
	Z = [x + im*y for y in Y, x in X]
end

# ╔═╡ 7eee64f0-14b9-11eb-2c53-650b51c7db37
begin
	freestream = Freestreams.Freestream(rand(ComplexF64))
	u = induce_velocity(Z, freestream, 0.0);
	quiver(real.(Z)[:], imag.(Z)[:], quiver = (real.(u[:]), imag.(u[:])), ratio = 1)
end

# ╔═╡ ed7e6140-14b9-11eb-0c7b-594e8b4fd1ed
begin
#	z = [x + im*y for y in linspace(-1, 1, 21), x in linspace(-1, 1, 21)][:],
    filter(X->abs(X)>0.25, Z)

	source = Source.Point(0.0, 1);
    usource = induce_velocity(Z, source, 0.0);

    quiver(real.(Z)[:], imag.(Z)[:], quiver = 0.2.*(real.(usource[:]), imag.(usource[:])), ratio = 1)
end

# ╔═╡ 911afe2e-14ba-11eb-3c75-a97deaefba0c
begin
	filter(X->abs(X)>0.25, Z)

	vortex = Vortex.Point(0.0, 2);
    uvortex = induce_velocity(Z, vortex, 0.0);
    quiver(real.(Z)[:], imag.(Z)[:], quiver = 0.2.*(real.(uvortex[:]), imag.(uvortex[:])), ratio = 1)
end

# ╔═╡ 06a79d20-14bb-11eb-3030-cdd4d53dfac5
begin
	filter(X->abs(X)>0.25, Z)
	f1 = Freestreams.Freestream(1.0)
    s1 = Source.Point(0.0, 2)

	u1 = induce_velocity(Z, (f1,s1), 0.0);

	quiver(real.(Z)[:], imag.(Z)[:], quiver = 0.2.*(real.(u1[:]), imag.(u1[:])), ratio = 1)
end

# ╔═╡ afbebea2-14bd-11eb-257a-63c212245316
begin

	vs = Vortex.Point.([-0.5, 0.5], [1, 1])
    ss = Source.Point.([-0.5, 0.5], [-1, 1])

	us = induce_velocity(Z, (vs,ss), 0.0);

	quiver(real.(Z)[:], imag.(Z)[:], quiver = 0.2.*(real.(us[:]), imag.(us[:])), ratio = 1)
end

# ╔═╡ 27ebaaf0-14be-11eb-3b42-d3c031377d61
begin

	v2 = Vortex.Point(0.0, 2)

	ψ = streamfunction(Z,v2)

	contour(X, Y, ψ, color = cgrad([:black, :black]), grid = false)
end

# ╔═╡ Cell order:
# ╠═0d696b90-14b9-11eb-19c6-7110e9eded9e
# ╠═40e53300-14b9-11eb-28ca-c39ecf8c05ef
# ╠═59e4bc40-14b9-11eb-3670-cdc8fa8f0ded
# ╠═7eee64f0-14b9-11eb-2c53-650b51c7db37
# ╠═ed7e6140-14b9-11eb-0c7b-594e8b4fd1ed
# ╠═911afe2e-14ba-11eb-3c75-a97deaefba0c
# ╠═06a79d20-14bb-11eb-3030-cdd4d53dfac5
# ╠═afbebea2-14bd-11eb-257a-63c212245316
# ╠═27ebaaf0-14be-11eb-3b42-d3c031377d61


text:
#=

    discrepancies.jl

    Contains code to assess discrepancy between data and model simulation using Cramér–von Mises distance.

    Author:     Alexander P. Browning
                ======================
                School of Mathematical Sciences
                Queensland University of Technology
                ======================
                ap.browning@icloud.com
                alexbrowning.me

=#

"""
    makecdf(x::Vector;n=1000)

Creates a linear interpolation approximation of the distribution function of `x`.

"""
function makecdf(x::Vector;n=1000)
    n_x = length(x)
    xp = range(minimum(x),maximum(x) + 1e-3,length=n)
    p = [count(x .< xpᵢ) for xpᵢ in xp] / n_x
    itp = LinearInterpolation(xp, p, extrapolation_bc = Interpolations.Flat())
    return itp
end

"""
    cvmises(cdf::AbstractInterpolation,y::Vector)

Computes Cramér–von Mises distance between distribution defined by the distribution function `cdf`
(given as an interpolation) and the univariate sample `y`.

"""
function cvmises(cdf::AbstractInterpolation,y::Vector)
    n = length(y)
    cdfs = cdf(sort(y))
    return sqrt(sum( ((2(1:n) .- 1.0) / 2n - cdfs).^2 ) / n + 1 / 12n^2)
end

"""
    addist(cdf::AbstractInterpolation,y::Vector)

Computes Anderson-Darling distance between distribution defined by the distribution function `cdf`
(given as an interpolation) and the univariate sample `y`.

"""
function addist(cdf::AbstractInterpolation,y::Vector)
    n = length(y)
    cdfs = min.(max.(cdf(sort(y)),1e-9),1-1e-9)
    sqrt(-n - sum( (2(1:n) .- 1.0) / n .* (log.(cdfs) + log.(1.0 .- cdfs[end:-1:1]) )  ))
end

"""
    ksdist(cdf::AbstractInterpolation,y::Vector)

Computes Kolmogorov–Smirnov distance between distribution defined by the distribution function `cdf`
(given as an interpolation) and the univariate sample `y`.

"""
function ksdist(cdf::AbstractInterpolation,y::Vector)
    n = length(y)
    cdfs = cdf(sort(y))
    δp = maximum((1:n) / n - cdfs)
    δn = -minimum((0:n-1) / n - cdfs)
    return max(δn, δp)
end

"""
kuiperdist(cdf::AbstractInterpolation,y::Vector)

Computes Kuiper's distance between distribution defined by the distribution function `cdf`
(given as an interpolation) and the univariate sample `y`.

"""
function kuiperdist(cdf::AbstractInterpolation,y::Vector)
    n = length(y)
    cdfs = cdf(sort(y))
    δp = maximum((1:n) / n - cdfs)
    δn = -minimum((0:n-1) / n - cdfs)
    return δn + δp
end

"""
    makediscrepancy(X;w=[0.5,0.5])

Create discrepancy function `dist(Y)` using data `X` where:
    X[i] = [[Q̄,Ū],[Q,U]]
and
    dist(Y) = ΣᵢΣⱼ( 2w[2] * cor(X[i][j]...) + w[1] * Σₖ addist(X[i][j][k]) )

By default, `w = [0.5,0.5]` so Anderson-Darling distance is weighted equally with correlation.

"""
function makediscrepancy(X::Vector{Vector{Vector{Vector{Float64}}}},f=addist;w=[0.5,0.5])
    Xcdf = [[makecdf.(X[i][j]) for j = 1:length(X[1])] for i = 1:length(X)]
    Xcor = [[cor(X[i][j]...) for j = 1:length(X[1])] for i = 1:length(X)]
    function dist(Y::Vector{Vector{Vector{Vector{Float64}}}})
        d = 0.0
        for i = 1:length(X), j = 1:length(X[1])
            d += 2 * w[1] * f(Xcdf[i][j][1],Y[i][j][1]) +    # Quenchable
                 1 * w[1] * f(Xcdf[i][j][2],Y[i][j][2]) +    # Unquenchable
                 1 * w[2] * norm(Xcor[i][j] - cor(Y[i][j]...))  # Correlations
        end
        return d
    end
    return dist
end

"""
    makediscrepancy_univariate(X;w=[0.5,0.5])

Create discrepancy function `dist(Y)` using data `X` where:
    X[i] = [[Q̄,Ū],[Q,U]]
and
    dist(Y) = ΣᵢΣⱼ addist(X[i][j][1])

Compared to `makediscrepancy`, only information relating to the quenchable probe is considered.

"""
function makediscrepancy_univariate(X::Vector{Vector{Vector{Vector{Float64}}}},f=addist)
    Xcdf = [[makecdf.(X[i][j]) for j = 1:length(X[1])] for i = 1:length(X)]
    function dist(Y::Vector{Vector{Vector{Vector{Float64}}}})
        d = 0.0
        for i = 1:length(X), j = 1:length(X[1])
            d += w[1] * f(Xcdf[i][j][1],Y[i][j][1])   # Quenchable
        end
        return d
    end
    return dist
end

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
"""
    QUERYNODE

Internal map from manipulation verb names (as `Symbol`s) to tuples
`(T<:QueryNode, H<:QueryHelper)`.
"""
const VERB = Dict{Symbol, DataType}(
    :select => Select,
    :filter => Filter,
    :groupby => GroupBy,
    :orderby => OrderBy,
    :summarize => Summarize,
    :summarise => Summarize,

    # <: JoinNode
    :leftjoin => LeftJoin,
    :outerjoin => OuterJoin,
    :innerjoin => InnerJoin,
    :join => InnerJoin,
    :crossjoin => CrossJoin,
)

_leaf(x) = DataNode(x)
_leaf(c::Cursor) = c.graph

function _graph_ex(args, index::Dict{Symbol, Symbol}, primary)
    stash = Dict{Set{Symbol}, Expr}(
        Set([src]) => Expr(:call, :(StructuredQueries._leaf), esc(src)) for src in values(index)
    )
    # for recording which manipulation tracks have merged via joins
    join_hist = Dict{Symbol, Set{Symbol}}(
        src => Set([src]) for src in values(index)
    )
    graph_ex = traverse!(args, stash, join_hist, index, primary)
    return graph_ex
end

"""

Stores the extant graph in `stash`.
"""
function traverse!(exs, stash, join_hist, index, primary)
    ex = exs[1]
    if ex.head == :line # skip line number Exprs
        return traverse!(exs[2:end], stash, join_hist, index, primary)
    elseif ex.head == :call
        verb, args = ex.args[1], ex.args[2:end]
        V = VERB[verb]
        node_ex, srcs_used = build_node_ex!(stash, join_hist, V, args, index, primary)
        # Have any such sources been joined *together* in previous verbs?
        # If so, they will have matching join_hist values
        srcs_hist = unique([join_hist[src] for src in srcs_used])
        # If we have two or more sources with different join histories,
        # record the join
        joined_srcs = union(srcs_hist...)
        foreach(src->(join_hist[src] = joined_srcs), joined_srcs)
        stash[joined_srcs] = node_ex
        if length(exs) > 1 # more verbs to go
            return traverse!(exs[2:end], stash, join_hist, index, primary)
        else
            return node_ex # no more verbs, return graph expression
        end
    else # Each following case is parsed as a select
        args = Vector{Expr}()
        if ex.head == :return
            _args = ex.args[1]
            if isa(_args, Expr)
                if _args.head == :tuple # multiple column selections
                    append!(args, _args.args)
                else
                    # assume a single column selection
                    # TODO: handle things
                    push!(args, _args)
                end
            else # ???
            end
        elseif ex.head == :tuple # multiple column selections
            append!(args, ex.args) # TODO: handle assignments...
        # TODO: handle assignment, but don't expect a row token
        # elseif ex.head == :(=)
        else
            # TODO: check for validity ...
            push!(args, ex)
        end
        node_ex, srcs_used = build_node_ex!(stash, join_hist, Select, args, index, primary)
        return node_ex
    end
end

function build_node_ex!{V}(stash, join_hist, ::Type{V}, args, index, primary)
    # process args into Verb instances
    # which sources are pertinent to the present verb invocation?
    srcs_used, helpers_ex = process_args(V, args, index, primary)
    # If no declared sources are used in any verbs, throw an error
    isempty(srcs_used) && throw(ArgumentError("No declared sources used in verbs."))
    # build expression to instantiate node for present verb
    src_nodes_ex = Expr(:tuple)
    unique_src_hists = unique([ join_hist[src] for src in srcs_used ])
    append!(src_nodes_ex.args, [ stash[src_hist] for src_hist in unique_src_hists ])
    node_ex = Expr(:call, Node{V}, src_nodes_ex, args, helpers_ex)
    # node_ex = Expr(:call, :(StructuredQueries._with), Expr(:call, Q, src_nodes_ex, args, helpers_ex))
    return node_ex, srcs_used
end


text:
mutable struct Poly <: Shape
    position::Vector{Luxor.Point}
    color::Union{Blending,Coloring}
    style::Union{Symbol, Vector{Float64}}
end

function Poly( ; position, color, style=:fill, kwargs...)
    return Poly(position, color, style)
end


text:
### Functions to Generate Pose Graphs ###
### Ted Steiner ###
### October 2014 ###

# Generate pose graph measurements from a route, set of poses, and list of loop-closure sites.
function makePoseGraph( nodes::Dict,
                        anchors::Array{Int,1},
                        route::Array{Int,1},
                        poses::Vector{Pose2D},
                        filename::String="posegraph.txt" )

    anchor_sets = Array{Int,1}[]
    push!(anchor_sets,anchors)

    filenames = String[]
    push!(filenames,filename)

    return makePoseGraph(nodes,anchor_sets,route,poses,filenames)
end

# Generate pose graphs from a single route and set of poses and loop-closures to multiple anchor sets
# Array of filenames should align with array of anchor_sets
function makePoseGraph( nodes::Dict,
                        anchor_sets::Array{Array{Int,1},1},
                        route::Array{Int,1},
                        poses::Vector{Pose2D},
                        filenames::Array{String,1} )

    odom_inf = diagm([5,5,50])
    loop_inf = diagm([50,50,100])
    odom_cov = inv(odom_inf)
    loop_cov = inv(loop_inf)

    odometry = getPoseOdometry(route,poses,odom_cov)
    meas = getMeasurements(route,poses,anchor_sets,odometry,loop_cov)

    for k = 1:length(anchor_sets)
        writeMeasurements(filenames[k],meas[k])
    end

    return nothing
end

function makePoseGraph( nodes::Dict,
                        anchor_sets::Array{Array{Int,1},1},
                        route::Array{Int,1},
                        poses::Vector{Pose2D},
                        odometry::Vector{BetweenFactor},
                        filenames::Array{String,1} )

    loop_inf = diagm([50,50,100])
    loop_cov = inv(loop_inf)

    meas = getMeasurements(route,poses,anchor_sets,odometry,loop_cov)

    for k = 1:length(anchor_sets)
        writeMeasurements(filenames[k],meas[k])
    end

    return nothing
end

function makePoseGraph( anchor_sets::Array{Array{Int,1},1},
                        route::Vector{Int},
                        poses::Vector{Pose2D},
                        odometry_factors::Vector{BetweenFactor},
                        loop_factors::Dict{(Int,Int),BetweenFactor},
                        prior_meas::Dict{Int,Meas2D},
                        prior_inf::Inf2D,
                        filenames::Vector{String} )

    meas = getMeasurements(route,poses,anchor_sets,odometry_factors,loop_factors,prior_meas,prior_inf)

    for k = 1:length(anchor_sets)
        writeMeasurements(filenames[k],meas[k])
    end

    return nothing
end


text:
# Load input
masses = open("day01_input.txt") do f
    [parse(Int32, line) for line in eachline(f)]
end

# Part 1
fuel_requirement = sum([Int32(floor(mass / 3)) - 2 for mass in masses])
println("Part 1: Fuel requirement: $fuel_requirement")

# Part 2
function compute_mass_recursive(fuel_mass::Int32)::Int32
    extra_fuel_mass = Int32(floor(fuel_mass / 3)) - 2
    if extra_fuel_mass <= 0
        return fuel_mass
    else
        return fuel_mass + compute_mass_recursive(extra_fuel_mass)
    end
end

fuel_requirement = sum([compute_mass_recursive(mass) - mass for mass in masses])
println("Part 2: Fuel requirement: $fuel_requirement")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
module TestExamplesThreaded

using Test
using Trixi

include("test_trixi.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
Trixi.mpi_isroot() && isdir(outdir) && rm(outdir, recursive=true)

@testset "Threaded tests" begin
  @testset "TreeMesh" begin
    @trixi_testset "elixir_advection_restart.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_restart.jl"),
        # Expected errors are exactly the same as in the serial test!
        l2   = [7.81674284320524e-6],
        linf = [6.314906965243505e-5])
    end

    @trixi_testset "elixir_advection_amr_refine_twice.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_amr_refine_twice.jl"),
        l2   = [0.00020547512522578292],
        linf = [0.007831753383083506])
    end

    @trixi_testset "elixir_advection_amr_coarsen_twice.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_amr_coarsen_twice.jl"),
        l2   = [0.0014321062757891826],
        linf = [0.0253454486893413])
    end

    @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_source_terms_nonperiodic.jl"),
        l2   = [2.259440511766445e-6, 2.318888155713922e-6, 2.3188881557894307e-6, 6.3327863238858925e-6],
        linf = [1.498738264560373e-5, 1.9182011928187137e-5, 1.918201192685487e-5, 6.0526717141407005e-5],
        rtol = 0.001)
    end

    @trixi_testset "elixir_euler_ec.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_ec.jl"),
        l2   = [0.061751715597716854, 0.05018223615408711, 0.05018989446443463, 0.225871559730513],
        linf = [0.29347582879608825, 0.31081249232844693, 0.3107380389947736, 1.0540358049885143])
    end
  end

  @testset "StructuredMesh" begin
    @trixi_testset "elixir_advection_restart.jl with waving flag mesh" begin
      @test_trixi_include(joinpath(examples_dir(), "structured_2d_dgsem", "elixir_advection_restart.jl"),
        l2   = [0.00016265538283059892],
        linf = [0.0015194298895079283],
        rtol = 7e-6, # Higher tolerance to make tests pass in CI with macOS
        elixir_file="elixir_advection_waving_flag.jl",
        restart_file="restart_000021.h5")
    end

    @trixi_testset "elixir_mhd_ec.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "structured_2d_dgsem", "elixir_mhd_ec.jl"),
        l2   = [0.04937480811868297, 0.06117033019988596, 0.060998028674664716, 0.03155145889799417,
                0.2319175391388658, 0.02476283192966346, 0.024483244374818587, 0.035439957899127385,
                0.0016022148194667542],
        linf = [0.24749024430983746, 0.2990608279625713, 0.3966937932860247, 0.22265033744519683,
                0.9757376320946505, 0.12123736788315098, 0.12837436699267113, 0.17793825293524734,
                0.03460761690059514],
        tspan = (0.0, 0.3))
    end
  end

  @testset "UnstructuredMesh" begin
    @trixi_testset "elixir_acoustics_gauss_wall.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "unstructured_2d_dgsem", "elixir_acoustics_gauss_wall.jl"),
      l2   = [0.029330394861252995, 0.029345079728907965, 0.03803795043486467, 0.0,
              7.175152371650832e-16, 1.4350304743301665e-15, 1.4350304743301665e-15],
      linf = [0.36236334472179443, 0.3690785638275256, 0.8475748723784078, 0.0,
              8.881784197001252e-16, 1.7763568394002505e-15, 1.7763568394002505e-15],
        tspan = (0.0, 5.0))
    end
  end

  @testset "P4estMesh" begin
    @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"),
        l2   = [0.0034516244508588046, 0.0023420334036925493, 0.0024261923964557187, 0.004731710454271893],
        linf = [0.04155789011775046, 0.024772109862748914, 0.03759938693042297, 0.08039824959535657])
    end

    @trixi_testset "elixir_eulergravity_convergence.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_eulergravity_convergence.jl"),
        l2   = [0.00024871265138964204, 0.0003370077102132591, 0.0003370077102131964, 0.0007231525513793697],
        linf = [0.0015813032944647087, 0.0020494288423820173, 0.0020494288423824614, 0.004793821195083758],
        tspan = (0.0, 0.1))
    end
  end

  @testset "DGMulti" begin
    @trixi_testset "elixir_euler_weakform.jl (SBP, EC)" begin
      @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_euler_weakform.jl"),
        cells_per_dimension = (4, 4),
        volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
        surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
        approximation_type = SBP(),
        l2 = [0.01280067571168776, 0.010607599608273158, 0.010607599608239976, 0.026408338014056502],
        linf = [0.03798302318566282, 0.05321027922532284, 0.05321027922605448, 0.13392025411839015],
      )
    end

    @trixi_testset "elixir_euler_triangulate_pkg_mesh.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_euler_triangulate_pkg_mesh.jl"),
        l2 = [4.664661209491976e-6, 3.7033509525940745e-6, 4.794877426562555e-6, 1.2682723101532175e-5],
        linf = [2.5099852761334418e-5, 2.2683684021362893e-5, 2.6180448559287584e-5, 5.5752932611508044e-5]
      )
    end
  end
end

# Clean up afterwards: delete Trixi output directory
Trixi.mpi_isroot() && isdir(outdir) && @test_nowarn rm(outdir, recursive=true)

end # module


text:
function ComputeShift(s,l,m,n,a,Op::OperatorShift)
    ##make domain
    r₊ = ψ.R.r₊ ; r₋ = ψ.R.r₋ ; s = ψ.s ; Δr = 0.1*(r₊-r₋); ϵ = eps(0.1);

    point1 = r₊ + Δr - Δr*im
    point2 = r₊ - Δr - Δr*im

    radial1 = SemiInfiniteLine(point1 , point1 + Δr*im , false)
    angular = LineSegment(-1.0+100*ϵ , 1.0-100*ϵ , true) #to avoid the NaNs at the edges
    C1 = radial1 ⊗ angular

    radial2 = LineSegment(point1,point2,true)
    C2 = radial2 ⊗ angular

    radial3 = SemiInfiniteLine(point2 , point2 + Δr*im , true)
    C3 = radial3 ⊗ angular

    TheContour = C1⊕C2⊕C3

end

max_stars_repo_name: Potatoasad/KerrQNMShifts

max_issues_repo_name: Potatoasad/KerrQNMShifts

max_forks_repo_name: Potatoasad/KerrQNMShifts


text:
module LMDB_Common
    using LMDB
    using Base.Test
    import LMDB.MDBValue

    @test LMDB.version()[1] == v"0.9.15"

    # LMDBError
    ex = LMDBError(Cint(0))
    @test_throws LMDBError throw(ex)
    @test ex.code == 0

    # MDBValue
    val = "abcd"
    mdb_val_ref = Ref(MDBValue(val))
    mdb_val = mdb_val_ref[]
    @test val == bytestring(convert(Ptr{UInt8}, mdb_val.data), mdb_val.size)

    val = [1233]
    T = eltype(val)
    val_size = sizeof(val)
    mdb_val = MDBValue(val[1])
    @test val_size == mdb_val.size
    nvals = floor(Int, mdb_val.size/sizeof(T))
    value = pointer_to_array(convert(Ptr{T}, mdb_val.data), nvals)
    @test val == value

    val = [0x0003, 0xff45]
    val_size = sizeof(val)
    T = eltype(val)
    mdb_val = MDBValue(val)
    @test val_size == mdb_val.size
    nvals = floor(Int, mdb_val.size/sizeof(T))
    value = pointer_to_array(convert(Ptr{T}, mdb_val.data), nvals)
    @test val == value
end


text:
mutable struct MvTuner{F <: Real}
  current_mean::Vector{F}
  current_cov::Matrix{F}
  beta::F
  iter::Int
  d::Int
end

MvTuner(d::Integer, T::Type=Float64) = MvTuner(zeros(T, d), eye(T, d)*0.01/d,
                                               T(.05), 1, d)

function MvTuner(v::AbstractVector{<:Real})
  d = length(v)
  return MvTuner(v, eye(d)*0.01/d, .05, 1, d)
end

function update!(x::AbstractVector{<:Real}, tuner::MvTuner{<:Real})
  tuner.iter += 1
  update_mean!(tuner.current_mean, x, tuner.iter)
  update_cov!(tuner.current_cov, tuner.current_mean, x, tuner.iter)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
###############################################
# Finding Nash Equlibria through optimization #
###############################################

"""
vFunction(game::Dict{String,<:AbstractArray},
                s::Vector{Vector{T}}) where T<:Real

Computes numerically for a given game and a given profile V - function crucial
to implement finding NE through optimization, if it equal to 0 - given
profile is a Nash Equilibrium, if greater - not

**Input Parameters**
* `game::Dict{String,<:AbstractArray}` - dictionary containing a game
* `s::Vector{Vector{T}}` - strategy profile for which we want to compute V - function
"""
function vFunction(game::Dict{String,<:AbstractArray},
                s::Vector{Vector{T}}) where T<:Real

    if !all(@. round(sum(s), digits = 2) == 1)
        error("Provided vectors are not a strategy profile")
    end

v = 0.0

for i in 1:length(s)
    PlayerActions = Matrix(1I, length(s[i]), length(s[i]))
    for j in 1:length(s[i])
        temp = deepcopy(s)
        temp[i] = PlayerActions[j,:]
        for k in 1:length(s)
        v = v + max(0, get_payoff(game, temp)[string("player",k)] - get_payoff(game, s)[string("player",k)])^2
end
    end
end

return v

end

"""
cart_prod_simplices(s::Vector{Vector{T}}) where T<:Real

Creates cartesian product of probability simplices for a given vector of vectors (especially - for a strategy profile)

from Manifold.jl
function(ProbabilitySimplex(n)) creates n-probability simplex
× - creates product

**Input Parameters**
* `s::Vector{Vector{T}}` - vector of vectors for which we create cartesian product of probability simplices
"""

function cart_prod_simplices(s::Vector{Vector{T}}) where T<:Real

vec = ProbabilitySimplex(length(s[1]))

for i in 2:length(s)
vec = vec × ProbabilitySimplex(length(s[i]))
end

return vec

end

"""
It is not possible to fully implement finding NE through optimization
without heavily modifing Julia optimization packages

The most promising at the time being is Manopt.jl (Manifold Optimization)
https://github.com/JuliaManifolds/Manopt.jl using certain data types from
Manifolds.jl (https://github.com/JuliaManifolds/Manifolds.jl).

However, at the time being (06.09.2020) not every manifold type from Manifolds.jl
is supported/imported to that library
https://github.com/JuliaManifolds/Manopt.jl/blob/master/src/Manopt.jl
- especially not Probability Simplex
https://github.com/JuliaManifolds/Manifolds.jl/blob/master/src/manifolds/ProbabilitySimplex.jl
crucial to finding NE through vFunction optimization

Other packages are far less advanced - JuMP
https://jump.dev/JuMP.jl/v0.21.1/installation/index.html
is a Julia frontend for various solvers and solver syntax is the best for
that library.

Optim.jl at the time being supports only very simple manifold constraints and
only for first-order algorithms
https://julianlsolvers.github.io/Optim.jl/stable/#algo/manifolds/
"""


text:
# Autogenerated wrapper script for rocm_smi_lib_jll for x86_64-linux-musl-cxx03
export liboam, librocm_smi, rocm_smipy, rsmiBindingspy, rsmiBindingspylink

JLLWrappers.@generate_wrapper_header("rocm_smi_lib")
JLLWrappers.@declare_library_product(liboam, "liboam.so.1")
JLLWrappers.@declare_library_product(librocm_smi, "librocm_smi64.so.1")
JLLWrappers.@declare_file_product(rocm_smipy)
JLLWrappers.@declare_file_product(rsmiBindingspy)
JLLWrappers.@declare_file_product(rsmiBindingspylink)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        liboam,
        "oam/lib/liboam.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        librocm_smi,
        "rocm_smi/lib/librocm_smi64.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        rocm_smipy,
        "bin/rocm_smi.py",
    )

    JLLWrappers.@init_file_product(
        rsmiBindingspy,
        "rocm_smi/bindings/rsmiBindings.py",
    )

    JLLWrappers.@init_file_product(
        rsmiBindingspylink,
        "bin/rsmiBindings.py",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "Articulation" begin
    gint = Graph(13)
    add_edge!(gint, 1, 7)
    add_edge!(gint, 1, 2)
    add_edge!(gint, 1, 3)
    add_edge!(gint, 12, 13)
    add_edge!(gint, 10, 13)
    add_edge!(gint, 10, 12)
    add_edge!(gint, 12, 11)
    add_edge!(gint, 5, 4)
    add_edge!(gint, 6, 4)
    add_edge!(gint, 8, 9)
    add_edge!(gint, 6, 5)
    add_edge!(gint, 1, 6)
    add_edge!(gint, 7, 5)
    add_edge!(gint, 7, 3)
    add_edge!(gint, 7, 8)
    add_edge!(gint, 7, 10)
    add_edge!(gint, 7, 12)

    for g in testgraphs(gint)
      art = @inferred(articulation(g))
      ans = [1, 7, 8, 12]
      @test art == ans
    end
    for level in 1:6
        btree = LightGraphs.BinaryTree(level)
        for tree in [btree, Graph{UInt8}(btree), Graph{Int16}(btree)]
          artpts = @inferred(articulation(tree))
          @test artpts == collect(1:(2^(level - 1) - 1))
        end
    end

    hint = LightGraphs.blkdiag(WheelGraph(5), WheelGraph(5))
    add_edge!(hint, 5, 6)
    for h in (hint, Graph{UInt8}(hint), Graph{Int16}(hint))
      @test @inferred(articulation(h)) == [5, 6]
    end
end


text:
######################
## type declaration ##
######################

@doc doc"""
RVMatrix constructor. RVMatrix must be a lower triangular matrix with
certain properties.
"""->
type RVMatrix
    matrix::Array{Int, 2}
    function RVMatrix(x::Array{Int, 2})
        if !istril(x)
            error("matrix must be triangular")
        end
        chkColumnSubsets(x)
        return new(x)
    end
end

## contructor helper function
##---------------------------

function chkColumnSubsets(x::Array{Int, 2})
    nVars = size(x, 1)
    for ii=nVars:-1:2
        if !(Set(x[ii:nVars, ii]) == Set(x[ii:nVars, ii-1]))
            error("column entries must be contained in column to the left")
        end
    end
    for ii=1:nVars-1
        if issubset(x[ii, ii], Set(x[ii+1:end, ii]))
            error("diagonal entries may not occur elsewhere in column")
        end
    end
    return true
end

################
## conversion ##
################

import Base.convert

function convert(::Type{RVMatrix}, vn::Vine)
    ## get number of variables
    nVars = size(vn.trees, 1)

    ## init RVMatrix
    rmatr = zeros(Int, nVars, nVars)

    ## convert vn to array of CTreePaths
    paths = convert(Copulas.CTreePaths, vn)

    ## set current column
    for ii=1:(nVars-1)
        ## find longest conditioning tree
        indMaxDepth = findMaxDepth(paths)

        ## current root
        currTree = paths[indMaxDepth]
        rt = currTree.root

        ## write down root
        rmatr[ii, ii] = rt

        ## get longest path of current tree
        longestPath = getLongestPath(currTree)
        rmatr[ii+1:end, ii] = flipud(longestPath)

        ## eliminate current root
        paths = allExcept(paths, indMaxDepth)

        ## eliminate root from each tree
        for jj=1:length(paths)
            for ll=1:length(dg(paths[jj]))
                if !isempty(dg(paths[jj])[ll])
                    if dg(paths[jj])[ll][end] == rt
                        pop!(dg(paths[jj])[ll])
                    end
                end
            end
        end
    end
    rmatr[end, end] = paths[1].root
    return RVMatrix(rmatr)
end

function convert(::Type{Copulas.Vine}, rmat::Copulas.RVMatrix)

    rmatr = dg(rmat)
    nVars = size(rmatr, 1)
    ## paths = Array(Array{Array{Int, 1}, 1}, nVars)
    paths = Array(Copulas.CTreePaths, nVars);
    for ii=1:(nVars-1)
        rt = rmatr[ii, ii]
        kk = Array(Array{Int, 1}, 0)
        push!(kk, flipud(rmatr[(ii+1):end, ii]))
        paths[rt] = Copulas.CTreePaths(rt, kk)
    end

    ## already include unconditional link for last node
    kk = Array(Array{Int, 1}, 0)
    push!(kk, [rmatr[nVars-1, nVars-1]])
    paths[rmatr[nVars, nVars]] = Copulas.CTreePaths(rmatr[nVars, nVars], kk)

    ## horizontal values

    ## unconditional links
    for ll=1:(nVars-2)
        rt = rmatr[nVars, ll]
        otherNode = rmatr[ll, ll]
        ## push!(paths[rt], [otherNode])
        Copulas.attach!(paths[rt], otherNode, [])
    end

    for jj=(nVars-1):-1:2
        for ll=1:(jj-1)
            rt = rmatr[jj, ll]
            otherNode = rmatr[ll, ll]
            condSet = rmatr[(jj+1):end, ll]
            Copulas.attach!(paths[rt], otherNode, condSet)
        end
    end

    display(paths)

    vnParRef = Array(Int, nVars, nVars)
    for ii=1:nVars
        kk = convert(Copulas.CTreeParRef, paths[ii])
        vnParRef[:, ii] = dg(kk)
    end
    vn = Copulas.Vine(vnParRef)
    return vn
end

@doc doc"""
Eliminate a single entry from an array of CTreePaths instances.
"""->
function allExcept(paths::Array{Copulas.CTreePaths,1}, ind)
    nPaths = length(paths)
    takeElems = trues(nPaths)
    takeElems[ind] = false
    return paths[takeElems]
end

@doc doc"""
For an array of conditioning trees find the first tree which has
maximum depth.
"""->
function findMaxDepth(ctrs::Array{Copulas.CTreePaths,1})
    depths = Int[Copulas.maxDepth(ctrs[ii]) for ii=1:size(ctrs, 1)]
    return indmax(depths)
end

@doc doc"""
For a given conditioning tree get the first path which has the maximum
path length.
"""->
function getLongestPath(ctr::Copulas.CTreePaths)
    maxDep = Copulas.maxDepth(ctr)
    isMaxDep = Bool[length(path) == maxDep for path in dg(ctr)]
    ind = find(isMaxDep)[1]
    return dg(ctr)[ind]
end

@doc doc"""
Delegate functions to single field of composite type.
"""->
function dg(rvm::RVMatrix)
    return rvm.matrix
end


text:
############################################################
# make_unique
############################################################
function make_unique!(s)
	count=1
	sp=fsortperm(s)
	current=s[sp[1]]
	for index in 2:length(s)
		if s[sp[index]]==current
			s[sp[index]]=s[sp[index]]*"."*string(count)
			count+=1
		else
			current=s[sp[index]]
			count=1
		end
	end
end

function make_unique(s)
	r=similar(s)
	count=1
	sp=fsortperm(s)
	r[sp[1]]=s[sp[1]]
	for index in 2:length(s)
		if s[sp[index]]==s[sp[index-1]]
			r[sp[index]]=s[sp[index]]*"."*string(count)
			count+=1
		else
			r[sp[index]]=s[sp[index]]
			count=1
		end
	end
	r
end

############################################################
# read_csv
############################################################
function read_csv(path::String)
    df = GZip.open(path, "r") do io
        CSV.read(io, DataFrame, header = false)
    end
    return(df)
end

############################################################
# read_mtx
############################################################
function read_mtx(filename::String)
    _parseint(x) = parse(Int, x)
    GZip.open(filename, "r") do mmfile
        # Read first line
        firstline = chomp(readline(mmfile))
        tokens = split(firstline)
        if length(tokens) != 5
            throw(string("Not enough words on first line: ", firstline))
        end
        if tokens[1] != "%%MatrixMarket"
            throw(string("Not a valid MatrixMarket header:", firstline))
        end
        (head1, rep, field, symm) = map(lowercase, tokens[2:5])
        if head1 != "matrix"
            throw("Unknown MatrixMarket data type: $head1 (only \"matrix\" is supported)")
        end

        eltypes = field == "real" ? Float32 :
                 field == "integer" ? Int32 :
                 throw(ParseError("Unsupported field $field (only real and complex are supported)"))

        symlabel = symm == "general" ? identity :
                   throw("Unknown matrix symmetry: $symm (only general, symmetric, skew-symmetric and hermitian are supported)")

        # Skip all comments and empty lines
        ll   = readline(mmfile)
        while length(chomp(ll))==0 || (length(ll) > 0 && ll[1] == '%')
            ll = readline(mmfile)
        end
        # Read matrix dimensions (and number of entries) from first non-comment line
        dd = map(_parseint, split(ll))
        if length(dd) < (rep == "coordinate" ? 3 : 2)
            throw(ParseError(string("Could not read in matrix dimensions from line: ", ll)))
        end
        rows = dd[1]
        cols = dd[2]
        dataf = CSV.read(mmfile, DataFrame, header = false, type = eltypes, delim = " ");
        return symlabel(sparse(dataf.Column1, dataf.Column2, dataf.Column3, rows, cols))
    end
end

############################################################
# sum_same_genes
############################################################
function sum_same_genes(raw_matrix::Array{Int32,2}, features::Vector{String}, barcodes::Vector{String})::NamedArray
    if (features |> unique |> length) == (features |> length)
        FinalMat = NamedArray(raw_matrix)
        setnames!(FinalMat, features, 1)
        setnames!(FinalMat, barcodes, 2)
    else
            # Count occurence of each genes
            countmap_features =  countmap(features)
            #Separate non unique genes with unique ones
            non_unique_features = (collect(keys(countmap_features))[values(countmap_features) .> 1])
            unique_features = collect(keys(countmap_features))[values(countmap_features) .== 1]
            #split matrix with those genes
            sub_X_unique_features = raw_matrix[in.(df_features.Column2, (unique_features,)),:]
            features1 = df_features.Column2[in.(df_features.Column2, (unique_features,))]
            sub_X_non_unique_features::SparseMatrixCSC{Int32,Int32} = reduce(vcat, map(x -> sum(raw_matrix[df_features.Column2 .== x, :], dims = 1), non_unique_features))
            FinalMat = vcat(sub_X_unique_features, sub_X_non_unique_features)
            ordered_features = append!(features1, non_unique_features)
            FinalMat = NamedArray(FinalMat)
            setnames!(FinalMat, ordered_features, 1)
            setnames!(FinalMat, barcodes[1], 2)
    end
    return(FinalMat)
end

############################################################
# read_10X
############################################################
function read_10X(path::String; features_col::Int = 2, barcodes_col::Int = 1, unique_features::Bool = true)
    # Check directory files
    dir_file = readdir(path)

    # Get files
    barc_path = dir_file[occursin.(r"^barcode", dir_file)]
    feat_path = dir_file[occursin.(r"^gene|^feature", dir_file)]
    rmat_path = dir_file[occursin.(r"^matrix", dir_file)]

    # Check if file exists and get full path
    isempty(barc_path) ? throw("$dir_file doesn't contain barcode file") : barc_path = join([path, barc_path[1]])
    isempty(feat_path) ? throw("$dir_file doesn't contain feature file") : feat_path = join([path, feat_path[1]])
    isempty(rmat_path) ? throw("$dir_file doesn't contain matrix file") : rmat_path = join([path, rmat_path[1]])

    # Read different components of teh 10X folder
    println("\nReading Matrix... \n")
    rmat = read_mtx(rmat_path);
    println("\nReading Features... \n")
    feat_df = read_csv(feat_path);
    println("\nReading Cell Barcodes... \n")
    barc_df = read_csv(barc_path);

    # Construct Named Arrays
    println("\nConstructing Named Arrays... \n")
    feat = feat_df[!, features_col]
    barc = barc_df[!, barcodes_col]
    feat_df_end = ncol(feat_df)
    if  ncol(feat_df) > 2
        list_assay = Dict{String,NamedArray{Int32,2}}()
        type_vec =  feat_df[!, feat_df_end]
        type =  type_vec |> unique
        for i in type
            feat_ind = feat_df[!,feat_df_end] .== i
            feat_split = feat_df[feat_ind, features_col];
            rmat_split = rmat[feat_ind,:];
                if unique_features
                    list_assay[i] = NamedArray(rmat_split, (make_unique(feat_split), barc))
                else
                    list_assay[i] = sum_same_genes(rmat_split, make_unique(feat_split), barc)
                end
        end
    else
        if unique_features
            list_assay = NamedArray(rmat, (make_unique(feat), barc))
        else
            list_assay = sum_same_genes(rmat, make_unique(feat), barc)
        end
    end
    println("\nDONE! \n")
    return(list_assay)
end

############################################################
# init_pissaro
############################################################
function init_pissaro(matrix::AbstractArray; min_cells::Int = 5, min_features::Int = 0, assay_name::String = "RNA")::Pissaro
    # Assay
    nzm =  matrix.array .!= 0
    n_cell = Vector{Int32}((sum(nzm .!= 0, dims = 2))[:,1])
    n_feat = Vector{Int32}((sum(nzm .!= 0, dims = 1))[1,:])
    min_cell_vec = n_cell .>= min_cells
    min_feat_vec = n_feat .>= min_features
    matrix = matrix[min_cell_vec, min_feat_vec]
    features_data = DataFrame(features = names(matrix,1), nCells = n_cell[min_cell_vec])
    cell_data = DataFrame(cells = names(matrix, 2), nFeatures = n_feat[min_feat_vec], nCount = Vector{Int32}(sum(matrix.array, dims = 1)[1,:]))
    sp_matrix::SparseMatrixCSC{Float32,Int32} = spzeros(size(matrix)[1], size(matrix)[2])
    empty_mat = NamedArray(sp_matrix, (names(matrix,1), names(matrix,2)))
    dict_assay = Dict(assay_name => Assays(matrix, deepcopy(empty_mat), NamedArray{Float32}(1,1), features_data, Array{String,1}(), assay_name))
    dim_red = Dict{String,Reductions}()
    pissaro_object = Pissaro(dict_assay, dim_red, cell_data, assay_name, "")
    return(pissaro_object)
end

############################################################
# FastExpMean
############################################################
function FastExpMean(matrix::AbstractArray)
    a = mean(expm1.(matrix), dims = 2)
    res = log1p.(a)[:,1]
    return res
end
############################################################
# FastLogVMR
############################################################
function FastLogVMR(matrix::AbstractArray)
    rmat = expm1.(matrix)
    a = mean(rmat, dims = 2)
    b = var(rmat , dims = 2)
    res = log1p.(a ./ b)[:,1]
    return res
end

max_stars_repo_name: Cortalak/Pissaro

max_issues_repo_name: Cortalak/Pissaro

max_forks_repo_name: Cortalak/Pissaro


text:
@doc """
  lift(a::padic) -> fmpz

  Returns the positive canonical representative in Z. a needs
  to be integral.
""" ->
function lift(a::padic)
  b = fmpz()
  R = parent(a)
  ccall((:padic_get_fmpz, :libflint), Void, (Ptr{fmpz}, Ptr{padic}, Ptr{FlintPadicField}), &b, &a, &R)
  return b
end


text:
using GeoInterface
using Documenter
using Test

include("test_primitives.jl")
doctest(GeoInterface)


text:
module Rsvg

using Cairo
using Glib_jll
using Librsvg_jll

include("gerror.jl")
include("gio.jl")
include("types.jl")
include("calls.jl")

export RsvgDimensionData, RsvgHandle

function __init__()
	Rsvg.set_default_dpi(72.0)
end

end # module


text:
struct Scan
  reactors::Vector{Reactor}
end

default_B_0_list = Vector{Number}(range(1, stop=20, length=41))

function Scan(reactor::Symbol, B_0_list::Vector{Number}=default_B_0_list; R_0_guess=R_0_max, kwargs...)
  input_dict = Dict()

  default = getfield(FusionSystems, Symbol(reactor))()

  for (cur_key, cur_value) in kwargs
    input_dict[cur_key] = getfield(default, cur_key) * cur_value
  end

  _tmp_solve(tmp_R_0, tmp_B_0) = solve!(
    getfield(FusionSystems, reactor)(; B_0=tmp_B_0, input_dict...), tmp_R_0
  )

  @assert length(B_0_list) > 1

  beg_reactor = _tmp_solve(R_0_guess, B_0_list[1])
  end_reactor = _tmp_solve(R_0_guess, B_0_list[end])

  if beg_reactor.valid && end_reactor.valid
    scan_reactors = [beg_reactor]
    R_0_guess = beg_reactor.R_0

    for B_0 in B_0_list[2:end-1]
      tmp_reactor = _tmp_solve(R_0_guess, B_0)
      @assert tmp_reactor.valid

      R_0_guess = tmp_reactor.R_0
      push!(scan_reactors, tmp_reactor)
    end

    push!(scan_reactors, end_reactor)
    return Scan(scan_reactors)
  end

  fast_scan(
    reactor, B_0_list, R_0_guess, input_dict,
    beg_reactor, end_reactor, _tmp_solve,
  )
end

function plot(scan::Scan, x::Symbol, y::Symbol)
  _filter_function(cur_value) = (
    isnothing(cur_value) || isnan(cur_value) ||
    isa(cur_value, SymEngine.Basic)
  )

  xx = [ getfield(reactor, x) for reactor in scan.reactors ]
  yy = [ getfield(reactor, y) for reactor in scan.reactors ]

  delete_indices = []
  append!(delete_indices, findall(_filter_function, xx))
  append!(delete_indices, findall(_filter_function, yy))

  delete_indices = sort(unique(delete_indices))
  deleteat!(xx, delete_indices)
  deleteat!(yy, delete_indices)

  @assert length(xx) == length(yy)
  SimplePlots.plot(xx, yy, label="")

  xlabel!(String(x))
  ylabel!(String(y))

  iszero(length(xx)) && return

  xlims!(0, 1.1*maximum(xx))
  ylims!(0, 1.1*maximum(yy))
end

function fast_scan(reactor, B_0_list, R_0_guess, input_dict, beg_reactor, end_reactor, _tmp_solve)
  scan_reactors = Vector{Reactor}(undef, length(B_0_list))

  scan_reactors[1] = beg_reactor
  scan_reactors[end] = end_reactor

  if beg_reactor.valid
    R_0_guess = beg_reactor.R_0
  elseif end_reactor.valid
    R_0_guess = end_reactor.R_0
  end

  used_index = NaN
  for cur_index in bisection_indices(length(B_0_list))
    B_0 = B_0_list[cur_index]
    tmp_reactor = _tmp_solve(R_0_guess, B_0)

    scan_reactors[cur_index] = tmp_reactor
    tmp_reactor.valid || continue

    used_index = cur_index
    R_0_guess = tmp_reactor.R_0

    break
  end

  @assert !isnan(used_index)

  init_R_0_guess = R_0_guess
  for cur_index in (used_index-1):-1:2
    B_0 = B_0_list[cur_index]
    tmp_reactor = _tmp_solve(R_0_guess, B_0)

    scan_reactors[cur_index] = tmp_reactor

    if beg_reactor.valid
      @assert tmp_reactor.valid
    else
      tmp_reactor.valid || break
    end

    R_0_guess = tmp_reactor.R_0
  end

  R_0_guess = init_R_0_guess
  for cur_index in (used_index+1):length(B_0_list)
    B_0 = B_0_list[cur_index]
    tmp_reactor = _tmp_solve(R_0_guess, B_0)

    scan_reactors[cur_index] = tmp_reactor

    if end_reactor.valid
      @assert tmp_reactor.valid
    else
      tmp_reactor.valid || break
    end

    R_0_guess = tmp_reactor.R_0
  end

  for cur_index in 1:length(B_0_list)
    isassigned(scan_reactors, cur_index) && continue
    B_0 = B_0_list[cur_index]

    tmp_reactor = getfield(FusionSystems, reactor)(; B_0=B_0, input_dict...)
    scan_reactors[cur_index] = tmp_reactor

    tmp_reactor.valid = false
    tmp_reactor.solved = true
  end

  Scan(scan_reactors)
end

export Scan


text:
using LightGraphs

sizesnbm1 = Int64[@allocated non_backtracking_matrix(CycleGraph(2^i))for i in 4:10]
sizesnbm2 = Int64[@allocated Nonbacktracking(CycleGraph(2^i)) for i in 4:10]

println("Relative Sizes:\n $(float(sizesnbm1./sizesnbm2))")

macro storetime(name, expression)
    ex = quote
        val = $expression;
        timeinfo[$name] = @elapsed $expression;
        val
    end
    return ex
end

function bench(g)
    nbt = @storetime :Construction_Nonbacktracking  Nonbacktracking(g)
    x  = ones(Float64, size(nbt)[1])
    info("Cycle with $n vertices has nbt in R^$(size(nbt))")

    B, nmap = @storetime :Construction_Dense        non_backtracking_matrix(g)
    y  = @storetime :Multiplication_Nonbacktracking nbt*x
    z  = @storetime :Multiplication_Dense           B*x;
    S  = @storetime :Construction_DS                sparse(B)
    z  = @storetime :Multiplication_Sparse          z = S*x;
    Sp = @storetime :Construction_Sparse            sparse(nbt)
    z  = @storetime :Multiplication_Sparse          Sp*x
end

function report(timeinfo)
    info("Times")
    println("Function\t Constructors\t Multiplication")
    println("Dense   \t $(timeinfo[:Construction_Dense])\t $(timeinfo[:Multiplication_Dense])")
    println("Sparse  \t $(timeinfo[:Construction_Sparse])\t $(timeinfo[:Multiplication_Sparse])")
    println("Implicit\t $(timeinfo[:Construction_Nonbacktracking])\t $(timeinfo[:Multiplication_Nonbacktracking])")
    info("Implicit Multiplication is $(timeinfo[:Multiplication_Dense]/timeinfo[:Multiplication_Nonbacktracking]) faster than dense.")
    info("Sparse Multiplication is $(timeinfo[:Multiplication_Nonbacktracking]/timeinfo[:Multiplication_Sparse]) faster than implicit.")
    info("Direct Sparse Construction took $(timeinfo[:Construction_Sparse])
    Dense to Sparse took: $(timeinfo[:Construction_DS])")
end

n = 2^13
C = CycleGraph(n)
timeinfo = Dict{Symbol, Float64}()
bench(C)
report(timeinfo)


text:
import LinearAlgebra._chol!
using LinearAlgebra: checksquare

const chol_mt_threshold = Ref(64.0)
mt_thresholds[:chol] = chol_mt_threshold;
const chol_dot_threshold = Ref(48)

function _chol!(A::AbstractMatrix{DT}, ::Type{UpperTriangular}
      ) where {DT <: Union{DoubleFloat{T},Complex{DoubleFloat{T}}}} where T
    @assert !has_offset_axes(A)
    n = checksquare(A)
    use_threads = (nthreads() > 1) &&
        (Float64(n)^2 > chol_mt_threshold[])
    use_views = n > chol_dot_threshold[]
    @inbounds begin
        for k = 1:n
            if use_views
                vk = uview(A,1:k-1,k)
                A[k,k] -= dot(vk,vk)
            else
                for i = 1:k - 1
                    A[k,k] -= A[i,k]'A[i,k]
                end
            end
            Akk, info = _chol!(A[k,k], UpperTriangular)
            if info != 0
                return UpperTriangular(A), info
            end
            A[k,k] = Akk
            AkkInv = inv(copy(Akk'))
            if use_threads
                _mt_uchol_loop(n,k,A,AkkInv,Vec{Npref,T})
            else
                if n < chol_dot_threshold[]
                    for j = k + 1:n
                        for i = 1:k - 1
                            A[k,j] -= A[i,k]'A[i,j]
                        end
                        A[k,j] = AkkInv*A[k,j]
                    end
                else
                    liA = LinearIndices(A)
                    for j = k + 1:n
                    #for i = 1:k - 1
                    #    A[k,j] -= A[i,k]'A[i,j]
                    #end
                    # uview should work here since A isn't yet UT
                    # but it is currently not inferred
                    # @inbounds A[k,j] -= dot(uview(A,1:k-1,k),uview(A,1:k-1,j))
                        @inbounds A[k,j] -= _dot(k-1,A,liA[1,k],A,liA[1,j],Vec{Npref,T})
                        A[k,j] = AkkInv*A[k,j]
                    end
                end
            end
        end
    end
    return UpperTriangular(A), convert(BlasInt, 0)
end

# function barrier helps inference (even for non-threaded path)
function _mt_uchol_loop(n,k,A,AkkInv,::Type{Vec{N,T}}) where {N,T}
    liA = LinearIndices(A)
    @threads for j = k + 1:n
#        @inbounds A[k,j] -= dot(uview(A,1:k-1,k),uview(A,1:k-1,j))
        @inbounds A[k,j] -= _dot(k-1,A,liA[1,k],A,liA[1,j],Vec{N,T})
        @inbounds A[k,j] = AkkInv*A[k,j]
    end
end

function _chol!(A::AbstractMatrix{DoubleFloat{T}}, ::Type{LowerTriangular}) where T
    @assert !has_offset_axes(A)
    n = checksquare(A)
    use_threads = (nthreads() > 1) &&
        (Float64(n)^2 > chol_mt_threshold[])
    @inbounds begin
        for k = 1:n
            #for i = 1:k - 1
            #    A[k,k] -= A[k,i]*A[k,i]'
            #end
            vk = view(A,k,1:k-1)
            A[k,k] -= dot(vk,vk)
            Akk, info = _chol!(A[k,k], LowerTriangular)
            if info != 0
                return LowerTriangular(A), info
            end
            A[k,k] = Akk
            AkkInv = inv(Akk)
            if k<n
                if use_threads
                    _mt_lchol_loop(n,k,A,AkkInv)
                else
                    for j = 1:k - 1
                        #@simd for i = k + 1:n
                        #    A[i,k] -= A[i,j]*A[k,j]'
                        #end
                        axpy!(-A[k,j]',uview(A,k+1:n,j),uview(A,k+1:n,k))
                    end
                    for i = k + 1:n
                        A[i,k] *= AkkInv'
                    end
                end
            end
        end
     end
    return LowerTriangular(A), convert(BlasInt, 0)
end

# can't just thread the j-loop because of dependence
#    @threads for j = 1:k - 1
#        @inbounds _axpy!(n-k,-A[k,j]',A,liA[k+1,j],A,liA[k+1,k],VT)
#    end

function _mt_lchol_loop(n,k,A,AkkInv)
    # liA = LinearIndices(A)
    @threads for i = k + 1:n
        # for j = 1:k - 1
            # A[i,k] -= A[i,j]*A[k,j]'
        # end
        # can't use LI until we have strided version of dot
        @inbounds A[i,k] -= dot(uview(A,k,1:k-1),uview(A,i,1:k-1))
    end
    @threads for i = k + 1:n
        @inbounds A[i,k] *= AkkInv'
    end
end


text:
const CURLOPT_URL = 10002
const CURLOPT_FOLLOWLOCATION = 52;
const CURLE_OK = 0
jlo = "http://julialang.org";
curl = ccall( (:curl_easy_init, "libcurl"), Ptr{Uint8}, ())
ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8},
   (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_URL, jlo.data)
ccall((:curl_easy_perform,"libcurl"),
 				Ptr{Uint8}, (Ptr{Uint8},), curl)
ccall((:curl_easy_cleanup,"libcurl"),
Ptr{Uint8},(Ptr{Uint8},), curl);

max_stars_repo_path: Module 3/Chapter04/call_curl.jl

max_issues_repo_path: Module 3/Chapter04/call_curl.jl

max_forks_repo_path: Module 3/Chapter04/call_curl.jl


text:
#=
This file contains structs and functions used the branch and bound
implementation.
=#

###
### Structs to hold upper and lower bounds on the treewidth of a graph.
###

"""A struct to hold the current upper bound of the branch and bound search."""
mutable struct UpperBound
    tw::UInt16            # The best treewidth found so far.
    order::Vector{UInt16} # The best elimination order found so far.
    lock::SpinLock        # A lock to prevent multiple threads updating the tw at once.
end

"""Use a heuristic to get an initial upper bound on the treewidth."""
function initialise_upper_bound(G::lg.AbstractGraph, seed::Integer)
    initial_ub_order, initial_ub_tw = min_fill(G; seed=seed)
    initial_ub = UpperBound(initial_ub_tw, initial_ub_order, SpinLock())
    @info "The initial treewidth is $initial_ub_tw"
    initial_ub
end

"""A struct to store best intermediate treewidths. (see FPBB Yaun 2011)"""
struct LowerBounds
    tabel::Dict{BitVector, UInt16}
    lock::SpinLock

    LowerBounds() = new(Dict{BitVector, UInt16}(), SpinLock())
end

###
### A struct to maintain the search state of a depth first search through elimination orders.
###

"""Represents the search state of a single instance of the branch and bound search."""
mutable struct DFSState
    graph::Graph                       # The graph corresponding to the current state.
    N::UInt16                          # The size of the original graph.

    intermediate_graph_key::BitVector  # A bit vector indicating which vertices of the original graph
                                       # have been eliminated.

    curr_order::Vector{UInt16}         # The current elimination order being constructed.
    ub::UpperBound                     # A struct holding the best elimination order found so far.
    lbs::LowerBounds                   # A struct holding the lower bounds found so far (See Yaun 2011)
    lb_ub::UInt16

    branches::Vector{Vector{UInt16}}   # A vector of vetrices to search for each node depth
    pq::MMDQueue                       # A pre-allocated priority queue to be used for mmd calls.

    # Some varibales to record where pruning happends (for introspection)
    lbs_prune_at_depth::Vector{UInt64}
    mmd_prune_at_depth::Vector{UInt64}
    ub_prune_at_depth::Vector{UInt64}

    # Variables to store various quantities related to the search.
    finish_time::Float64
    timed_out::Bool
    space_covered::Float64
    nodes_visited::UInt64
    depth::UInt16
    rng::MersenneTwister
end

"""DFSState constructor."""
function DFSState(g::Graph, ub::UpperBound, lbs::LowerBounds, finish_time::Float64, seed::Int=42)
    N = g.num_vertices
    rng = MersenneTwister(seed)
    initial_curr_order = collect(0x0001:N)

    branches = [Vector{UInt16}(undef, N-i) for i = 0:N-2]

    pq = MMDQueue(g.vertices, g.degree)

    DFSState(g, N, falses(N),
             initial_curr_order, ub, lbs, 0x0000, branches, pq,
             zeros(UInt64, N), zeros(UInt64, N), zeros(UInt64, N),
             finish_time, false, 0.0, 0, 1, rng)
end

function initialise_dfsstates(G::lg.AbstractGraph,
                                ub::UpperBound,
                                finish_time::Float64,
                                seed::Integer,
                                n::Integer)
    # Create a search state for each thread.
    lbs = LowerBounds()
    state = DFSState(Graph(G), ub, lbs, finish_time, seed)
    [copy(state, seed + t) for t = 1:n]
end

function prep_state!(state::DFSState, run_time::AbstractFloat)::Nothing
    state.finish_time = time() + run_time
    state.space_covered = 0.0
    state.timed_out = false
    nothing
end

"""Return a copy of the given DFSState."""
function Base.copy(s::DFSState, seed::Int=42)
    DFSState(deepcopy(s.graph), s.N, copy(s.intermediate_graph_key),
             copy(s.curr_order), s.ub, s.lbs, 0x0000, deepcopy(s.branches), deepcopy(s.pq),
             zeros(UInt64, s.N), zeros(UInt64, s.N), zeros(UInt64, s.N),
             s.finish_time, false, 0.0, 0, s.depth, MersenneTwister(seed))
end

"""Prints some info about the given DFSState."""
function Base.show(io::IO, s::DFSState)
    compact = get(io, :compact, false)
    if !compact
        println(io, "The graph size is {nodes, edges} = ", s.graph)
        println(io, "The best treewidth found was ", s.ub.tw)
        println(io, "The number of nodes visted was ", s.nodes_visited)
    else
        show(io, (s.ub.tw, s.ub.order))
    end
end

###
### A struct to hold the results of a depth first search.
###

"""Holds the results of a depth first search"""
struct DFSReport
    states::Vector{DFSState}    # A vector containing the states returned from each thread used.
    tw::UInt16                  # The best treewidth found
    order::Vector{UInt16}       # The best elimination order found

    # Variables to indicate how much of the
    # search space was covered.
    nodes_visited::UInt64
    space_covered::AbstractFloat

    # Time measurements.
    duration::AbstractFloat
    actual_time::AbstractFloat
    heuristic_time::AbstractFloat
    dfs_time::AbstractFloat

    # Vectors indicating where pruning happened.
    ub_pruned::Vector{UInt64}
    lbs_pruned::Vector{UInt64}
    mmd_pruned::Vector{UInt64}
end

"""Constructor for a DFSReport"""
function DFSReport(
                states::Vector{DFSState},
                allocated_time,
                start_time,
                heurisitc_finish_time,
                actual_finish_time
                )
    actual_time = actual_finish_time - start_time
    heuristic_time = heurisitc_finish_time - start_time
    dfs_time = actual_finish_time - heurisitc_finish_time

    tw = states[1].ub.tw
    order = copy(states[1].ub.order)

    total_nodes_visited = sum([state.nodes_visited for state in states])
    total_space_covered = sum([s.timed_out ? s.space_covered : 1/length(states) for s in states])

    ub_pruned = sum([state.ub_prune_at_depth for state in states])
    lbs_pruned = sum([state.lbs_prune_at_depth for state in states])
    mmd_pruned = sum([state.mmd_prune_at_depth for state in states])

    DFSReport(
            states,
            tw,
            order,
            total_nodes_visited,
            total_space_covered,
            allocated_time,
            actual_time,
            heuristic_time,
            dfs_time,
            ub_pruned,
            lbs_pruned,
            mmd_pruned
            )
end

"""Prints the results of a branch-and-bound depth first search."""
function Base.show(io::IO, r::DFSReport)
    compact = get(io, :compact, false)
    if !compact
        println(io, "")
        println(io, r.states[1].graph)
        println(io, "The best treewidth found was ", r.tw)
        println(io, "The number of nodes visted was ", r.nodes_visited)
        println(io, "The fraction of the search space covered was (to machine precision) ", r.space_covered)
        println(io, "Time (seconds):")
        println(io, " allocated   = ", r.duration)
        println(io, " actual      = ", r.actual_time)
        println(io, " heuristic   = ", r.heuristic_time)
        println(io, " b&b dfs     = ", r.dfs_time)
    else
        show(io, (r.best_tw, r.best_order))
    end
end

###
### functions to alter the graph while searching through elimination orders
###

"""
    eliminate!(state::DFSState, v::Int)

Eliminates vertex 'v' from the given state and update all its relevant variables.

Returns a vector of neighbours of v and the edges added when v was eliminated. These can
be used to restore the eliminated vertex.
"""
function Graphs.eliminate!(state::DFSState, v::UInt16)
    state.intermediate_graph_key[v] = true
    state.lb_ub += eliminate!(state.graph, v)
    nothing
end

"""
    restore_last_eliminated!(state::DFSState,
                            edges_to_remove::Vector{Tuple{Int, Int}})

Restores a vertex 'v' which was eliminated to form the given 'state'.
"""
function Graphs.restore_last_eliminated!(state::DFSState)
    restore_last_eliminated!(state.graph)
    v = state.graph.vertices[state.graph.num_vertices]
    state.intermediate_graph_key[v] = false
    nothing
end

###
### Forwarding heuristic functions to a state's graph.
###

"""Use the mmd heuristic on the given state graph"""
function mmd(state::DFSState)
    g = state.graph
    verts = vertices(g)

    initialise_mmdqueue!(state.pq, verts, g.degree)
    mmd(state.graph, state.pq)
end


text:
import Term:
    remove_markup,
    has_markup,
    remove_ansi,
    has_ansi,
    get_last_ANSI_code,
    cleantext,
    textlen,
    escape_brackets,
    unescape_brackets,
    replace_text,
    nospaces,
    remove_brackets,
    unspace_commas,
    chars,
    join_lines,
    split_lines,
    textwidth,
    fillin,
    truncate,
    reshape_text
import Term.Style: apply_style
import Term.Measures: width as get_width

@testset "TU_markup" begin
    strings = [
        (
            "this is {red} some {blue} text {/blue} that I like{/red}",
            "this is  some  text  that I like",
        ),
        (
            "{bold underline} this is {red on_green} text {/red on_green} I like {/bold underline}",
            " this is  text  I like ",
        ),
    ]

    for (s1, s2) in strings
        @test has_markup(s1)
        @test remove_markup(s1) == s2
        @test cleantext(s1) == s2
        @test textlen(s1) == textwidth(s2)
    end

    @test remove_markup("text with {{double}} squares") == "text with {{double}} squares"
    @test !has_markup("text with {{double}} squares")
end

@testset "TU_ansi" begin
    apply_style("test{(.2, .5, .6)}coloooor{/(.2, .5, .6)}")
    strings = [
        (
            "this is \e[31m some \e[34m text \e[39m\e[31m that I like\e[39m",
            "this is  some  text  that I like",
            "\e[39m",
        ),
        (
            "\e[1m\e[4m this is \e[31m\e[42m text \e[39m\e[49m\e[4m I like \e[22m\e[24m",
            " this is  text  I like ",
            "\e[24m",
        ),
        (
            "test\e[38;2;51;128;153mcoloooor\e[39m and white",
            "testcoloooor and white",
            "\e[39m",
        ),
    ]

    for (s1, s2, ltag) in strings
        @test has_ansi(s1)
        @test remove_ansi(s1) == s2
        @test get_last_ANSI_code(s1) == ltag
    end
end

@testset "TU_replace_text" begin
    text = "abcdefghilmnopqrstuvz"

    @test replace_text(text, 0, 5, "aaa") == "aaafghilmnopqrstuvz"
    @test replace_text(text, 0, 5, ',') == ",,,,,fghilmnopqrstuvz"

    @test replace_text(text, 18, 21, "aaa") == "abcdefghilmnopqrstaaa"
    @test replace_text(text, 18, 21, ',') == "abcdefghilmnopqrst,,,"

    @test replace_text(text, 10, 15, "aaa") == "abcdefghilaaarstuvz"
    @test replace_text(text, 10, 15, ',') == "abcdefghil,,,,,rstuvz"

    @test nospaces("a (1, 2, 3) 4") == "a(1,2,3)4"
    @test remove_brackets("aaa (asdsd) BB") == "aaa asdsd BB"

    @test unspace_commas("a, 2, 3") == "a,2,3"
end

@testset "TU_misc" begin
    @test chars("abcd") == ['a', 'b', 'c', 'd']

    strings = [
        "aaa\nadasda\nasdasda",
        """
        asdasd
adsada
asddsa""",
    ]
    for str in strings
        @test join_lines(split_lines(str)) == str
    end
end

@testset "TU_reshape" begin
    str = """
Lorem ipsum {bold}dolor sit{/bold} amet, consectetur adipiscing elit,
ed do e{red}iusmod tempor incididunt{/red} ut {bold}labore et {underline}dolore{/underline} magna aliqua.{/bold} Ut enim ad minim
veniam, quis{green} nostrud exercitation {on_black}ullamco laboris nisi ut aliquip ex {/on_black}
ea commodo consequat.{blue} Duis aute irure dolor in{/blue} reprehenderit
in voluptate velit{/green} esse {italic}cillum dolore{/italic}{red} eu{/red}{italic green} fugiat {/italic green}nulla
pariatur. Excepteur{red} sint{/red}{blue} occaecat cupidatat {/blue}non proident,
sunt in culpa qui {italic}officia{/italic} deserunt mollit anim
id est laborum."""

    str_reshaped = "Lorem ipsum \e[1mdolor sit\e[22m amet, \nconsectetur adipiscing elit,\ned do e\e[31miusmod tempor incididunt\e[39m \nut \e[1mlabore et \e[4mdolore\e[24m\e[1m magna \naliqua.\e[22m\e[39m Ut enim ad minim\nveniam, quis\e[32m nostrud \n\e[32mexercitation \e[40mullamco laboris \n\e[32mnisi ut aliquip ex \e[49m\e[32m\e[39m\e[39m\e[32m\e[39m\e[40m\nea commodo consequat.\e[34m Duis aute \nirure dolor in\e[39m reprehenderit \nin voluptate velit\e[39m\e[3m esse \e[3mcillum \ndolore\e[23m\e[31m eu\e[39m\e[23m\e[3m\e[32m fugiat \e[23m{/green}nulla \npariatur. Excepteur\e[31m sint\e[39m\e[34m \noccaecat cupidatat \e[39mnon proident, \nsunt in culpa qui \e[3mofficia\e[23m \ndeserunt mollit anim \nid est laborum."

    logo_str = """Term.jl is a {#9558B2}Julia{/#9558B2} package for creating styled terminal outputs.

    Term provides a simple {italic green4 bold}markup language{/italic green4 bold} to add {bold bright_blue}color{/bold bright_blue} and {bold underline}styles{/bold underline} to your text.
    More complicated text layout can be created using {red}"Renderable"{/red} objects such
    as {red}"Panel"{/red} and {red}"TextBox"{/red}.
    These can also be nested and stacked to create {italic pink3}fancy{/italic pink3} and {underline}informative{/underline} terminal ouputs for your Julia code"""

    logo_str_reshaped = "Term.jl is a \e[38;2;149;88;178mJulia\e[39m package for \ncreating styled terminal outputs.\n\nTerm provides a simple \e[3m\e[38;5;28m\e[1mmarkup \nlanguage\e[23m\e[39m\e[3m\e[22m\e[38;5;28m to add \e[1m\e[38;5;12mcolor\e[22m\e[38;5;28m\e[39m\e[1m and \e[1m\e[4mstyles\e[22m\e[1m\e[24m\e[1m \nto your text.\nMore complicated text layout can \nbe created using \e[31m\"Renderable\"\e[39m \nobjects such \nas {red}\"Panel\"{/red} and {red}\"TextBox\"{/red}.\nThese can also be nested and \nstacked to create \e[3m\e[38;5;175mfancy\e[23m\e[39m\e[3m and \n\e[4minformative\e[24m\e[3m terminal ouputs for \nyour Julia code"

    strings = [
        (
            "Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua.",
            "Lorem ipsum dolor sit amet, \nconsectetur adipiscing elit, sed \ndo eiusmod tempor incididunt ut \nlabore et dolore magna aliqua.",
        ),
        (
            "Lorem {red}ipsum dolor sit {underline}amet, consectetur{/underline} adipiscing elit, {/red}{blue}sed do eiusmod tempor incididunt{/blue} ut labore et dolore magna aliqua.",
            "Lorem \e[31mipsum dolor sit \e[4mamet, \nconsectetur\e[24m\e[31m adipiscing elit, \e[39m\e[34msed \ndo eiusmod tempor incididunt\e[39m ut \nlabore et dolore magna aliqua.",
        ),
        (
            "Lorem{red}ipsumdolorsit{underline}amet, consectetur{/underline} adipiscing elit, {/red}seddoeiusmo{blue}dtemporincididunt{/blue}ut labore et dolore magna aliqua.",
            "Lorem\e[31mipsumdolorsit\e[4mamet, \nconsectetur\e[24m\e[31m adipiscing elit, \n\e[39mseddoeiusmo\e[34mdtemporincididunt\e[39mut \nlabore et dolore magna aliqua.",
        ),
        (
            "ต้าอ่วยวาทกรรมอาว์เซี้ยว กระดี๊กระด๊า ช็อปซาดิสต์โมจิดีพาร์ตเมนต์ อินดอร์วิว สี่แยกมาร์กจ๊อกกี้ โซนี่บัตเตอร์ฮันนีมูน ยาวีแพลนหงวนสคริปต์ แจ็กพ็อตต่อรองโทรโข่งยากูซ่ารุมบ้า บอมบ์เบอร์รีวีเจดีพาร์ทเมนท์ บอยคอตต์เฟอร์รี่บึมมาราธอน ",
            "ต้าอ่วยวาทกรรมอาว์เซี้ยว กระดี๊กระด๊า \nช็อปซาดิสต์โมจิดีพาร์ตเมนต์ อินดอร์วิว \nสี่แยกมาร์กจ๊อกกี้ โซนี่บัตเตอร์ฮันนีมูน \nยาวีแพลนหงวนสคริปต์ \nแจ็กพ็อตต่อรองโทรโข่งยากูซ่ารุมบ้า \nบอมบ์เบอร์รีวีเจดีพาร์ทเมนท์ \nบอยคอตต์เฟอร์รี่บึมมาราธอน ",
        ),
        (
            "ต้าอ่วยวาท{red}กรรมอาว์เซี้ยว กระดี๊กระด๊า {/red}ช็อปซาดิสต์โมจิดีพาร์ตเม{blue underline}นต์ อินดอร์วิว สี่แยกมาร์กจ๊อกกี้ โซนี่บัตเต{/blue underline}อร์ฮันนีมูน ยาวีแพลนหงวนสคริปต์ แจ็กพ็อตต่อรองโทรโข่งยากูซ่ารุมบ้า บอมบ์เบอร์รีวีเจดีพาร์ทเมนท์ บอยคอตต์เฟอร์รี่บึมมาราธอน ",
            "ต้าอ่วยวาท\e[31mกรรมอาว์เซี้ยว กระดี๊กระด๊า \n\e[39mช็อปซาดิสต์โมจิดีพาร์ตเม\e[34m\e[4mต์ อินดอร์วิว \nสี่แยกมาร์กจ๊อกกี้ โซนี่บัตเต\e[39m\e[24m\e[34m}อร์ฮันนีมูน \nยาวีแพลนหงวนสคริปต์ \nแจ็กพ็อตต่อรองโทรโข่งยากูซ่ารุมบ้า \nบอมบ์เบอร์รีวีเจดีพาร์ทเมนท์ \nบอยคอตต์เฟอร์รี่บึมมาราธอน ",
        ),
        (
            "국가유공자·상이군경 및 전몰군경의 유가족은 법률이 정하는 바에 의하여",
            "국가유공자·상이군경 및 \n전몰군경의 유가족은 법률이 \n정하는 바에 의하여",
        ),
        (
            "국{red}가유공자·상이군{bold}경 및 전{/bold}몰군경의 유{/red}가족은 법률이 정하는 바에 의하여",
            "국\e[31m가유공자·상이군\e[1m경 및 \n전\e[22m\e[31m몰군경의 유\e[39m가족은 법률이 \n정하는 바에 의하여",
        ),
        (
            "朗眠裕安無際集正聞進士健音社野件草売規作独特認権価官家複入豚末告設悟自職遠氷育教載最週場仕踪持白炎組特曲強真雅立覧自価宰身訴側善論住理案者券真犯著避銀楽験館稿告",
            "朗眠裕安無際集正聞進士健音社野件\n草売規作独特認権価官家複入豚末告\n設悟自職遠氷育教載最週場仕踪持白\n炎組特曲強真雅立覧自価宰身訴側善\n論住理案者券真犯著避銀楽験館稿告",
        ),
        (
            "┌────────────────┬────────────────┬────────────────┬────────────────┬──────────────",
            "┌────────────────┬───────────────\n─┬────────────────┬──────────────\n──┬──────────────",
        ),
        (
            "┌────────────abcde────┬──────────── ────┬────────abcde fghi────────┬────────────────┬──────────────",
            "┌────────────abcde────┬──────────\n── ────┬────────abcde \nfghi────────┬────────────────┬───\n───────────",
        ),
        (
            "┌─────────{red}───ab{/red}cde────┬──────{green}────── ────┬────────abcde fghi{/green}────────┬────────────────┬──────────────",
            "┌─────────\e[31m───ab\e[39mcde────\n┬──────\e[32m────── \n────┬────────abcde \nfghi\e[39m────────┬────────────\n────┬──────────────",
        ),
        (
            "┌──────────{red}────{/red}──┬{blue bold}────────────────┬──{/blue bold}──────────────┬────────────────┬──────────────end",
            "┌──────────\e[31m────\e[39m──┬\e[1m───────────────┬──{/blue}\n\e[22m\e[39m}──────────────┬──────────\n──────┬──────────────end\e[39m",
        ),
        (
            "."^100,
            ".................................\n.................................\n.................................\n.",
        ),
        (
            ".{red}|||{/red}...."^10,
            ".\e[31m|||\e[39m.....\e[31m|||{/red\n}.....\e[31m|||\e[39m.....\e[31m|||\n\e[39m\e[31m.....\e[31m|||\e[39m.....|||\e[39m\e[31m.....\e[31m|||\e[39m\e[31m...\n..\e[31m|||\e[39m\e[31m.....\e[31m|||\e[39m.....\e[31m|||\e[39m\e[31m....\e[39m",
        ),
        (
            ".|||...."^10,
            ".|||.....|||.....|||.....|||.....\n|||.....|||.....|||.....|||.....|\n||.....|||....",
        ),
        (str, str_reshaped),
        (logo_str, logo_str_reshaped),
    ]

    width = 33
    debug = false
    for (i, (input, expected)) in enumerate(strings)
        reshaped = reshape_text(input, width)
        reshaped_no_ansi = remove_ansi(reshaped)
        lens = length.(split(reshaped_no_ansi, '\n'))
        if debug && reshaped != expected
            println("== reshaped == ")
            println(reshaped)
            println(repr(reshaped))
            println("\n== reshaped no ansi == ")
            println(reshaped_no_ansi)
            println("\n== expected == ")
            println(expected)
        end
        # FIXME: should work when `length(input) != ncodeunits(input)` using non unit byte characters: see docs.julialang.org/en/v1/manual/strings/#Unicode-and-UTF-8
        if length(input) == ncodeunits(input) && !occursin('\n', input)
            (debug && any(lens .> width)) && println(lens)
            @test all(lens .≤ width)
        end
        @test reshaped == expected
    end

    for width in (40, 60, 99)
        rh = reshape_text(str, width)
        @test all(textlen.(split(rh, '\n')) .≤ width)
    end

    for i in 5:10
        width = 2^i
        for offset in (-(width ÷ 2)):(width ÷ 2)
            txt = '.'^(2width + offset)
            rt = reshape_text(txt, width)
            len = length.(split(rt, '\n'))
            # @show length(txt) width rt len
            @test all(len[1:(end - 1)] .== width)
            @test len[end] == (offset > 0 ? offset : width + offset)
        end
    end
end


text:
using Base.StackTraces

"Executes a feature or scenario, and presents the results in real time."
struct Executor
    stepdefmatcher::StepDefinitionMatcher
    presenter::RealTimePresenter
    executionenv::ExecutionEnvironment

    Executor(matcher::StepDefinitionMatcher,
             presenter::RealTimePresenter = QuietRealTimePresenter();
             executionenv::ExecutionEnvironment = NoExecutionEnvironment()) = new(matcher, presenter, executionenv)
end

"Abstract type for a result when executing a scenario step."
abstract type StepExecutionResult end

"No matching step definition was found."
struct NoStepDefinitionFound <: StepExecutionResult
    step::Gherkin.ScenarioStep
end

"More than one matching step definition was found."
struct NonUniqueMatch <: StepExecutionResult
    locations::Vector{StepDefinitionLocation}
end

"Successfully executed a step."
struct SuccessfulStepExecution <: StepExecutionResult end

"An assert failed in the step."
struct StepFailed <: StepExecutionResult
    assertion::String
end

"An unexpected error or exception occurred during the step execution."
struct UnexpectedStepError <: StepExecutionResult
    ex::Exception
    stack::StackTrace
end

"The step was not executed."
struct SkippedStep <: StepExecutionResult end

"""
    issuccess(::StepExecutionResult)

True for successful step executions, and false otherwise.
"""
issuccess(::SuccessfulStepExecution) = true
issuccess(::StepExecutionResult) = false

"A result for a complete scenario."
struct ScenarioResult
    steps::Vector{StepExecutionResult}
    scenario::Scenario
end

"""
    executescenario(::Executor, ::Gherkin.Scenario)

Execute each step in a `Scenario`. Present the results in real time.
Returns a `ScenarioResult`.
"""
function executescenario(executor::Executor, scenario::Gherkin.Scenario)
    # Tell the presenter of the start of a scenario execution.
    present(executor.presenter, scenario)

    # The `context` object for a scenario execution is provided to each step,
    # so they may store intermediate values.
    context = StepDefinitionContext()

    beforescenario(executor.executionenv, context, scenario)

    # The `steps` vector contains the results for all steps. At initialization,
    # they are all `Skipped`, because if one step fails then we stop the execution of the following
    # steps.
    steps = Vector{StepExecutionResult}(undef, length(scenario.steps))
    fill!(steps, SkippedStep())

    # Keep track of which the last executed step is. This could be removed if the code was shuffled
    # a bit below.
    lastexecutedstep = 0

    # Find a unique step definition for each step and execute it.
    # Present each step, first as it is about to be executed, and then once more with the result.
    for i = 1:length(scenario.steps)
        step = scenario.steps[i]
        # Tell the presenter that this step is about to be executed.
        present(executor.presenter, step)

        steps[i] = try
            # Find a step definition matching the step text.
            stepdefinition = findstepdefinition(executor.stepdefmatcher, step)

            # The block text is provided to the step definition via the context.
            context[:block_text] = step.block_text

            try
                # Execute the step definition. Note that it's important to use `Base.invokelatest` here,
                # because otherwise it might not find that function defined yet.
                Base.invokelatest(stepdefinition.definition, context)
            catch ex
                UnexpectedStepError(ex, stacktrace(catch_backtrace()))
            end
        catch ex
            if ex isa NoMatchingStepDefinition
                NoStepDefinitionFound(step)
            elseif ex isa NonUniqueStepDefinition
                NonUniqueMatch(ex.locations)
            else
                rethrow(ex)
            end
        end
        # Present the result after execution.
        present(executor.presenter, step, steps[i])
        lastexecutedstep = i

        # If this step failed, then we skip any remaining steps.
        if !issuccess(steps[i])
            break
        end
    end

    # Present any remaining steps as skipped.
    for k = lastexecutedstep + 1:length(scenario.steps)
        present(executor.presenter, scenario.steps[k])
        present(executor.presenter, scenario.steps[k], steps[k])
    end

    afterscenario(executor.executionenv, context, scenario)

    ScenarioResult(steps, scenario)
end

"""
    executescenario(::Executor, ::Gherkin.ScenarioOutline)

Execute a `Scenario Outline`, which contains one or more examples. Each example is transformed into
a regular `Scenario`, and it's executed.
Reeturns a list of `ScenarioResult`s.
"""
function executescenario(executor::Executor, outline::Gherkin.ScenarioOutline)
    scenarios = transformoutline(outline)
    [executescenario(executor, scenario)
     for scenario in scenarios]
end

"The execution result for a feature, containing one or more scenarios."
struct FeatureResult
    feature::Feature
    scenarioresults::Vector{ScenarioResult}
end

"""
    executefeature(::Executor, ::Gherkin.Feature)

Execute all scenarios and scenario outlines in a feature.
"""
function executefeature(executor::Executor, feature::Gherkin.Feature)
    # Present that a new feature is about to be executed.
    present(executor.presenter, feature)

    # Execute each scenario and scenario outline in the feature.
    scenarioresults = [executescenario(executor, s) for s in feature.scenarios]

    # Return a list of `ScenarioResults`.
    # Since regular scenarios return a `ScenarioResult` directly, and scenario outlines return lists
    # of `ScenarioResult`s, we have to flatten that list.
    FeatureResult(feature, reduce(vcat, scenarioresults, init=[]))
end

#
# Interpolation of Scenario Outlines
#

function transformoutline(outline::ScenarioOutline)
    # The examples are in a multidimensional array.
    # The size of dimension 1 is the number of placeholders.
    # The size of dimension 2 is the number of examples.
    [interpolatexample(outline, outline.examples[:,exampleindex])
     for exampleindex in 1:size(outline.examples, 2)]
end

function interpolatexample(outline::ScenarioOutline, example::Vector{T}) where {T <: AbstractString}
    # Create placeholders on the form `<foo>` for each placeholder `foo`, and map it against the
    # value for this particular example.
    #
    # An example Scenario Outline could have two placeholders `foo` and `bar`:
    #    Scenario Outline: An example
    #       When some action <foo>
    #       Then some postcondition <bar>
    #
    #    Examples:
    #    | foo | bar |
    #    | 1   | 2   |
    #    | 17  | 42  |
    #
    # The parsed Scenario Outline object already has a list of the placeholders `foo` and `bar`.
    # This method is called twice, once with example vector `[1, 2]` and once with `[17, 42]`.
    #
    placeholders_kv = ["<$(outline.placeholders[i])>" => example[i] for i in 1:length(example)]
    placeholders = Dict{String, T}(placeholders_kv...)

    # `interpolatestep` just creates a new scenario step of the same type, but with all occurrences
    # of placeholders replace with the example value.
    fromplaceholders = x -> placeholders[x]
    steps = [interpolatestep(step, fromplaceholders) for step in outline.steps]

    Scenario(outline.description, outline.tags, steps)
end

interpolatestep(step::Given, fromplaceholders::Function) = Given(interpolatesteptext(step.text, fromplaceholders);
                                                                 block_text=interpolatesteptext(step.block_text, fromplaceholders))
interpolatestep(step::When, fromplaceholders::Function) = When(interpolatesteptext(step.text, fromplaceholders);
                                                               block_text=interpolatesteptext(step.block_text, fromplaceholders))
interpolatestep(step::Then, fromplaceholders::Function) = Then(interpolatesteptext(step.text, fromplaceholders);
                                                               block_text=interpolatesteptext(step.block_text, fromplaceholders))

interpolatesteptext(text::String, fromplaceholders::Function) = replace(text, r"<[^>]*>" => fromplaceholders)


text:
using CSV, DataFrames, Statistics

results_dir = ARGS[1]
results = CSV.read("$results_dir/results.csv")

by(results, [:sampler, :Nd], :B0_ess=>mean, :sigma_ess=>mean, Symbol("B[1]_ess")=>mean, Symbol("B[2]_ess")=>mean, Symbol("B[3]_ess")=>mean, :tree_depth=>mean, :epsilon=>mean) |> println