text:
import TimeZones: Transition
import Base.Dates: Hour, Second, UTM

# Constructor FixedTimeZone from a string.
@test FixedTimeZone("0123") == FixedTimeZone("UTC+01:23", 4980)
@test FixedTimeZone("+0123") == FixedTimeZone("UTC+01:23", 4980)
@test FixedTimeZone("-0123") == FixedTimeZone("UTC-01:23", -4980)
@test FixedTimeZone("01:23") == FixedTimeZone("UTC+01:23", 4980)
@test FixedTimeZone("+01:23") == FixedTimeZone("UTC+01:23", 4980)
@test FixedTimeZone("-01:23") == FixedTimeZone("UTC-01:23", -4980)
@test FixedTimeZone("01:23:45") == FixedTimeZone("UTC+01:23:45", 5025)
@test FixedTimeZone("+01:23:45") == FixedTimeZone("UTC+01:23:45", 5025)
@test FixedTimeZone("-01:23:45") == FixedTimeZone("UTC-01:23:45", -5025)
@test FixedTimeZone("99:99:99") == FixedTimeZone("UTC+99:99:99", 362439)
@test FixedTimeZone("UTC") == FixedTimeZone("UTC", 0)
@test FixedTimeZone("UTC+00") == FixedTimeZone("UTC", 0)
@test FixedTimeZone("UTC+1") == FixedTimeZone("UTC+01:00", 3600)
@test FixedTimeZone("UTC-1") == FixedTimeZone("UTC-01:00", -3600)
@test FixedTimeZone("UTC+01") == FixedTimeZone("UTC+01:00", 3600)
@test FixedTimeZone("UTC-01") == FixedTimeZone("UTC-01:00", -3600)
@test FixedTimeZone("UTC+0123") == FixedTimeZone("UTC+01:23", 4980)
@test FixedTimeZone("UTC-0123") == FixedTimeZone("UTC-01:23", -4980)

@test_throws ArgumentError FixedTimeZone("1")
@test_throws ArgumentError FixedTimeZone("01")
@test_throws ArgumentError FixedTimeZone("123")
@test_throws ArgumentError FixedTimeZone("012345")
@test_throws ArgumentError FixedTimeZone("0123:45")
@test_throws ArgumentError FixedTimeZone("01:2345")
@test_throws ArgumentError FixedTimeZone("01:-23:45")
@test_throws ArgumentError FixedTimeZone("01:23:-45")
@test_throws ArgumentError FixedTimeZone("01:23:45:67")
@test_throws ArgumentError FixedTimeZone("UTC1")

# Test exception messages
tz = FixedTimeZone("Imaginary/Zone", 0, 0)

buffer = IOBuffer()
showerror(buffer, AmbiguousTimeError(DateTime(2015,1,1), tz))
@test takebuf_string(buffer) == "Local DateTime 2015-01-01T00:00:00 is ambiguious"

buffer = IOBuffer()
showerror(buffer, NonExistentTimeError(DateTime(2015,1,1), tz))
@test takebuf_string(buffer) == "DateTime 2015-01-01T00:00:00 does not exist within Imaginary/Zone"

warsaw = resolve("Europe/Warsaw", tzdata["europe"]...)

# Standard time behaviour
local_dt = DateTime(1916, 2, 1, 0)
utc_dt = DateTime(1916, 1, 31, 23)

# Disambiguating parameters ignored when there is no ambiguity.
@test ZonedDateTime(local_dt, warsaw).zone.name == :CET
@test ZonedDateTime(local_dt, warsaw, 1).zone.name == :CET
@test ZonedDateTime(local_dt, warsaw, 2).zone.name == :CET
@test ZonedDateTime(local_dt, warsaw, true).zone.name == :CET
@test ZonedDateTime(local_dt, warsaw, false).zone.name == :CET
@test ZonedDateTime(utc_dt, warsaw, from_utc=true).zone.name == :CET

@test ZonedDateTime(local_dt, warsaw).utc_datetime == utc_dt
@test ZonedDateTime(local_dt, warsaw, 1).utc_datetime == utc_dt
@test ZonedDateTime(local_dt, warsaw, 2).utc_datetime == utc_dt
@test ZonedDateTime(local_dt, warsaw, true).utc_datetime == utc_dt
@test ZonedDateTime(local_dt, warsaw, false).utc_datetime == utc_dt
@test ZonedDateTime(utc_dt, warsaw, from_utc=true).utc_datetime == utc_dt

# Daylight saving time behaviour
local_dt = DateTime(1916, 6, 1, 0)
utc_dt = DateTime(1916, 5, 31, 22)

# Disambiguating parameters ignored when there is no ambiguity.
@test ZonedDateTime(local_dt, warsaw).zone.name == :CEST
@test ZonedDateTime(local_dt, warsaw, 1).zone.name == :CEST
@test ZonedDateTime(local_dt, warsaw, 2).zone.name == :CEST
@test ZonedDateTime(local_dt, warsaw, true).zone.name == :CEST
@test ZonedDateTime(local_dt, warsaw, false).zone.name == :CEST
@test ZonedDateTime(utc_dt, warsaw, from_utc=true).zone.name == :CEST

@test ZonedDateTime(local_dt, warsaw).utc_datetime == utc_dt
@test ZonedDateTime(local_dt, warsaw, 1).utc_datetime == utc_dt
@test ZonedDateTime(local_dt, warsaw, 2).utc_datetime == utc_dt
@test ZonedDateTime(local_dt, warsaw, true).utc_datetime == utc_dt
@test ZonedDateTime(local_dt, warsaw, false).utc_datetime == utc_dt
@test ZonedDateTime(utc_dt, warsaw, from_utc=true).utc_datetime == utc_dt

# Typical "spring-forward" behaviour
local_dts = (
    DateTime(1916,4,30,22),
    DateTime(1916,4,30,23),
    DateTime(1916,5,1,0),
)
utc_dts = (
    DateTime(1916,4,30,21),
    DateTime(1916,4,30,22),
)
@test_throws NonExistentTimeError ZonedDateTime(local_dts[2], warsaw)
@test_throws NonExistentTimeError ZonedDateTime(local_dts[2], warsaw, 1)
@test_throws NonExistentTimeError ZonedDateTime(local_dts[2], warsaw, 2)
@test_throws NonExistentTimeError ZonedDateTime(local_dts[2], warsaw, true)
@test_throws NonExistentTimeError ZonedDateTime(local_dts[2], warsaw, false)

@test ZonedDateTime(local_dts[1], warsaw).zone.name == :CET
@test ZonedDateTime(local_dts[3], warsaw).zone.name == :CEST
@test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).zone.name == :CET
@test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).zone.name == :CEST

@test ZonedDateTime(local_dts[1], warsaw).utc_datetime == utc_dts[1]
@test ZonedDateTime(local_dts[3], warsaw).utc_datetime == utc_dts[2]
@test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).utc_datetime == utc_dts[1]
@test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).utc_datetime == utc_dts[2]

# Typical "fall-back" behaviour
local_dt = DateTime(1916, 10, 1, 0)
utc_dts = (DateTime(1916, 9, 30, 22), DateTime(1916, 9, 30, 23))
@test_throws AmbiguousTimeError ZonedDateTime(local_dt, warsaw)

@test ZonedDateTime(local_dt, warsaw, 1).zone.name == :CEST
@test ZonedDateTime(local_dt, warsaw, 2).zone.name == :CET
@test ZonedDateTime(local_dt, warsaw, true).zone.name == :CEST
@test ZonedDateTime(local_dt, warsaw, false).zone.name == :CET
@test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).zone.name == :CEST
@test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).zone.name == :CET

@test ZonedDateTime(local_dt, warsaw, 1).utc_datetime == utc_dts[1]
@test ZonedDateTime(local_dt, warsaw, 2).utc_datetime == utc_dts[2]
@test ZonedDateTime(local_dt, warsaw, true).utc_datetime == utc_dts[1]
@test ZonedDateTime(local_dt, warsaw, false).utc_datetime == utc_dts[2]
@test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).utc_datetime == utc_dts[1]
@test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).utc_datetime == utc_dts[2]

# Zone offset reduced creating an ambigious hour
local_dt = DateTime(1922,5,31,23)
utc_dts = (DateTime(1922, 5, 31, 21), DateTime(1922, 5, 31, 22))
@test_throws AmbiguousTimeError ZonedDateTime(local_dt, warsaw)

@test ZonedDateTime(local_dt, warsaw, 1).zone.name == :EET
@test ZonedDateTime(local_dt, warsaw, 2).zone.name == :CET
@test_throws AmbiguousTimeError ZonedDateTime(local_dt, warsaw, true)
@test_throws AmbiguousTimeError ZonedDateTime(local_dt, warsaw, false)
@test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).zone.name == :EET
@test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).zone.name == :CET

@test ZonedDateTime(local_dt, warsaw, 1).utc_datetime == utc_dts[1]
@test ZonedDateTime(local_dt, warsaw, 2).utc_datetime == utc_dts[2]
@test ZonedDateTime(utc_dts[1], warsaw, from_utc=true).utc_datetime == utc_dts[1]
@test ZonedDateTime(utc_dts[2], warsaw, from_utc=true).utc_datetime == utc_dts[2]

# Check behaviour when save is larger than an hour.
paris = resolve("Europe/Paris", tzdata["europe"]...)

@test ZonedDateTime(DateTime(1945,4,2,1), paris).zone == FixedTimeZone("WEST", 0, 3600)
@test_throws NonExistentTimeError ZonedDateTime(DateTime(1945,4,2,2), paris)
@test ZonedDateTime(DateTime(1945,4,2,3), paris).zone == FixedTimeZone("WEMT", 0, 7200)

@test_throws AmbiguousTimeError ZonedDateTime(DateTime(1945,9,16,2), paris)
@test ZonedDateTime(DateTime(1945,9,16,2), paris, 1).zone == FixedTimeZone("WEMT", 0, 7200)
@test ZonedDateTime(DateTime(1945,9,16,2), paris, 2).zone == FixedTimeZone("CET", 3600, 0)

# Ensure that dates are continuous when both a UTC offset and the DST offset change.
@test ZonedDateTime(DateTime(1945,9,16,1), paris).utc_datetime == DateTime(1945,9,15,23)
@test ZonedDateTime(DateTime(1945,9,16,2), paris, 1).utc_datetime == DateTime(1945,9,16,0)
@test ZonedDateTime(DateTime(1945,9,16,2), paris, 2).utc_datetime == DateTime(1945,9,16,1)
@test ZonedDateTime(DateTime(1945,9,16,3), paris).utc_datetime == DateTime(1945,9,16,2)

# Transitions changes that exceed an hour.
t = VariableTimeZone("Testing", [
    Transition(DateTime(1800,1,1), FixedTimeZone("TST",0,0)),
    Transition(DateTime(1950,4,1), FixedTimeZone("TDT",0,7200)),
    Transition(DateTime(1950,9,1), FixedTimeZone("TST",0,0)),
])

# A "spring forward" where 2 hours are skipped.
@test ZonedDateTime(DateTime(1950,3,31,23), t).zone == FixedTimeZone("TST",0,0)
@test_throws NonExistentTimeError ZonedDateTime(DateTime(1950,4,1,0), t)
@test_throws NonExistentTimeError ZonedDateTime(DateTime(1950,4,1,1), t)
@test ZonedDateTime(DateTime(1950,4,1,2), t).zone == FixedTimeZone("TDT",0,7200)

# A "fall back" where 2 hours are duplicated. Never appears to occur in reality.
@test ZonedDateTime(DateTime(1950,8,31,23), t).utc_datetime == DateTime(1950,8,31,21)  # TDT

# First occurrences of duplicated hours.
@test ZonedDateTime(DateTime(1950,9,1,0), t, 1).utc_datetime == DateTime(1950,8,31,22) # TST
@test ZonedDateTime(DateTime(1950,9,1,1), t, 1).utc_datetime == DateTime(1950,8,31,23) # TST

# Second occurrences of duplicated hours.
@test ZonedDateTime(DateTime(1950,9,1,0), t, 2).utc_datetime == DateTime(1950,9,1,0)   # TDT
@test ZonedDateTime(DateTime(1950,9,1,1), t, 2).utc_datetime == DateTime(1950,9,1,1)   # TDT

@test ZonedDateTime(DateTime(1950,9,1,2), t).utc_datetime == DateTime(1950,9,1,2)      # TDT

# Ambiguous local DateTime that has more than 2 solutions. Never occurs in reality.
t = VariableTimeZone("Testing", [
    Transition(DateTime(1800,1,1), FixedTimeZone("TST",0,0)),
    Transition(DateTime(1960,4,1), FixedTimeZone("TDT",0,7200)),
    Transition(DateTime(1960,8,31,23), FixedTimeZone("TXT",0,3600)),
    Transition(DateTime(1960,9,1), FixedTimeZone("TST",0,0)),
])

@test ZonedDateTime(DateTime(1960,8,31,23), t).utc_datetime == DateTime(1960,8,31,21)  # TDT
@test ZonedDateTime(DateTime(1960,9,1,0), t, 1).utc_datetime == DateTime(1960,8,31,22) # TDT
@test ZonedDateTime(DateTime(1960,9,1,0), t, 2).utc_datetime == DateTime(1960,8,31,23) # TXT
@test ZonedDateTime(DateTime(1960,9,1,0), t, 3).utc_datetime == DateTime(1960,9,1,0)   # TST
@test ZonedDateTime(DateTime(1960,9,1,1), t).utc_datetime == DateTime(1960,9,1,1)      # TST

@test_throws AmbiguousTimeError ZonedDateTime(DateTime(1960,9,1,0), t, true)
@test_throws AmbiguousTimeError ZonedDateTime(DateTime(1960,9,1,0), t, false)

# Significant offset change: -11:00 -> 13:00.
apia = resolve("Pacific/Apia", tzdata["australasia"]...)

# Skips an entire day.
@test ZonedDateTime(DateTime(2011,12,29,23),apia).utc_datetime == DateTime(2011,12,30,9)
@test_throws NonExistentTimeError ZonedDateTime(DateTime(2011,12,30,0),apia)
@test_throws NonExistentTimeError ZonedDateTime(DateTime(2011,12,30,23),apia)
@test ZonedDateTime(DateTime(2011,12,31,0),apia).utc_datetime == DateTime(2011,12,30,10)

# Redundant transitions should be ignored.
# Note: that this can occur in reality if the TZ database parse has a Zone that ends at
# the same time a Rule starts. When this occurs the duplicates always in standard time
# with the same abbreviation.
zone = Dict{AbstractString,FixedTimeZone}()
zone["DTST"] = FixedTimeZone("DTST", 0, 0)
zone["DTDT-1"] = FixedTimeZone("DTDT-1", 0, 3600)
zone["DTDT-2"] = FixedTimeZone("DTDT-2", 0, 3600)

dup = VariableTimeZone("DuplicateTest", [
    Transition(DateTime(1800,1,1), zone["DTST"])
    Transition(DateTime(1935,4,1), zone["DTDT-1"])  # Ignored
    Transition(DateTime(1935,4,1), zone["DTDT-2"])
    Transition(DateTime(1935,9,1), zone["DTST"])
])

# Make sure that the duplicated hour only doesn't contain an additional entry.
@test_throws AmbiguousTimeError ZonedDateTime(DateTime(1935,9,1), dup)
@test ZonedDateTime(DateTime(1935,9,1), dup, 1).zone.name == Symbol("DTDT-2")
@test ZonedDateTime(DateTime(1935,9,1), dup, 2).zone.name == :DTST
@test_throws BoundsError ZonedDateTime(DateTime(1935,9,1), dup, 3)

# Ensure that DTDT-1 is completely ignored.
@test_throws NonExistentTimeError ZonedDateTime(DateTime(1935,4,1), dup)
@test ZonedDateTime(DateTime(1935,4,1,1), dup).zone.name == Symbol("DTDT-2")
@test ZonedDateTime(DateTime(1935,8,31,23), dup).zone.name == Symbol("DTDT-2")

# Check equality between ZonedDateTimes
utc = FixedTimeZone("UTC", 0, 0)

spring_utc = ZonedDateTime(DateTime(2010, 5, 1, 12), utc)
spring_apia = ZonedDateTime(DateTime(2010, 5, 1, 1), apia)

# The absolutely min DateTime you can create. Even smaller than typemin(DateTime)
early_utc = ZonedDateTime(DateTime(UTM(typemin(Int64))), utc)

@test spring_utc.zone == FixedTimeZone("UTC", 0, 0)
@test spring_apia.zone == FixedTimeZone("SST", -39600, 0)
@test spring_utc == spring_apia
@test spring_utc !== spring_apia
@test !isequal(spring_utc, spring_apia)
@test hash(spring_utc) != hash(spring_apia)
@test astimezone(spring_utc, apia) === spring_apia  # Since ZonedDateTime is immutable
@test astimezone(spring_apia, utc) === spring_utc
@test isequal(astimezone(spring_utc, apia), spring_apia)
@test hash(astimezone(spring_utc, apia)) == hash(spring_apia)

fall_utc = ZonedDateTime(DateTime(2010, 10, 1, 12), utc)
fall_apia = ZonedDateTime(DateTime(2010, 10, 1, 2), apia)

@test fall_utc.zone == FixedTimeZone("UTC", 0, 0)
@test fall_apia.zone == FixedTimeZone("SDT", -39600, 3600)
@test fall_utc == fall_apia
@test fall_utc !== fall_apia
@test !(fall_utc < fall_apia)
@test !(fall_utc > fall_apia)
@test !isequal(fall_utc, fall_apia)
@test hash(fall_utc) != hash(fall_apia)
@test astimezone(fall_utc, apia) === fall_apia  # Since ZonedDateTime is immutable
@test astimezone(fall_apia, utc) === fall_utc
@test isequal(astimezone(fall_utc, apia), fall_apia)
@test hash(astimezone(fall_utc, apia)) == hash(fall_apia)

# A FixedTimeZone is effective for all of time where as a VariableTimeZone has as start.
@test TimeZones.utc(early_utc) < apia.transitions[1].utc_datetime
@test_throws NonExistentTimeError astimezone(early_utc, apia)

# ZonedDateTime constructor that takes any number of Period or TimeZone types
@test_throws ArgumentError ZonedDateTime(FixedTimeZone("UTC", 0, 0), FixedTimeZone("TMW", 86400, 0))

# Equality for VariableTimeZones
another_warsaw = resolve("Europe/Warsaw", tzdata["europe"]...)

@test warsaw == warsaw
@test warsaw === warsaw
@test warsaw == another_warsaw
@test warsaw !== another_warsaw

# VariableTimeZone with a cutoff set
cutoff_tz = VariableTimeZone(
    "cutoff", [Transition(DateTime(1970, 1, 1), utc)], DateTime(1988, 5, 6),
)

ZonedDateTime(DateTime(1970, 1, 1), cutoff_tz)  # pre cutoff
@test_throws UnhandledTimeError ZonedDateTime(DateTime(1988, 5, 6), cutoff_tz)  # on cutoff
@test_throws UnhandledTimeError ZonedDateTime(DateTime(1989, 5, 7), cutoff_tz)
@test_throws UnhandledTimeError ZonedDateTime(DateTime(1988, 5, 5), cutoff_tz) + Hour(24)

zdt = ZonedDateTime(DateTime(2038, 3, 28), warsaw, from_utc=true)
@test_throws UnhandledTimeError zdt + Hour(1)

# TimeZones that no longer have any transitions after the max_year shouldn't have a cutoff
# eg. Asia/Hong_Kong, Pacific/Honolulu, Australia/Perth
perth = resolve("Australia/Perth", tzdata["australasia"]...)
zdt = ZonedDateTime(DateTime(2200, 1, 1), perth, from_utc=true)

# Convenience constructors for making a DateTime on-the-fly
digits = [2010, 1, 2, 3, 4, 5, 6]
for i in eachindex(digits)
    @test ZonedDateTime(digits[1:i]..., warsaw) == ZonedDateTime(DateTime(digits[1:i]...), warsaw)
    @test ZonedDateTime(digits[1:i]..., utc) == ZonedDateTime(DateTime(digits[1:i]...), utc)
end

# Convenience constructor dealing with ambiguous time
digits = [1916, 10, 1, 0, 2, 3, 4]  # Fall DST transition in Europe/Warsaw
for i in eachindex(digits)
    expected = [
        ZonedDateTime(DateTime(digits[1:i]...), warsaw, 1)
        ZonedDateTime(DateTime(digits[1:i]...), warsaw, 2)
    ]

    if i > 1
        @test_throws AmbiguousTimeError ZonedDateTime(digits[1:i]..., warsaw)
    end

    @test ZonedDateTime(digits[1:i]..., warsaw, 1) == expected[1]
    @test ZonedDateTime(digits[1:i]..., warsaw, 2) == expected[2]
    @test ZonedDateTime(digits[1:i]..., warsaw, true) == expected[1]
    @test ZonedDateTime(digits[1:i]..., warsaw, false) == expected[2]
end


text:
######################################################################################################
# Merge from                                                                                #
# * Wouter Koolen: https://bitbucket.org/wmkoolen/tidnabbil/src/master/purex_games_paper/, and       #
# * Xuedong Shang: https://github.com/xuedong/LinBAI.jl                                              #
# We only modify the `dump_stats` with our added `oracle` for the linear and Lipschitz structure.    #
######################################################################################################
using Statistics;

function solve(pep, μ, δ, β)
    Tstar, wstar = oracle(pep, μ);
    ⋆ = istar(pep, μ);
    # lower bound
    kl = (1-2δ)*log((1-δ)/δ);
    lbd = Tstar*kl;
    # more practical lower bound with the employed threshold β
    practical = binary_search(t -> t-Tstar*β(t), max(1, lbd), 1e10);
    Tstar, wstar, ⋆, lbd, practical;
end

function dump_stats(pep, μ, δs, βs, srs, datas, repeats)
    best_arm = typeof(pep) == BestArm;
    linear_bai = typeof(pep) == LinearBestArm;
    linear_threshold = typeof(pep) == LinearThreshold;
    linear = linear_bai || linear_threshold;
    lipschitz = typeof(pep) == LipschitzBestArm;
    linear ? K=length(pep.arms) : K=length(μ);
    if lipschitz
        value_lipschitz, wstar_lipschitz = oracle_lipschitz_bai(pep);
        value_standard, wstar_standard = oracle(pep, μ);
        println("[optimal value] standard=$(value_standard), Lipschitz=$(value_lipschitz)");
    elseif linear
        if linear_bai
            value_linear, wstar_linear = oracle_linear_bai(pep, μ);
            value_standard, wstar_standard = oracle(pep, [μ'pep.arms[k] for k=1:K]);
            println("[optimal value] standard=$(value_standard), linear=$(value_linear)");
        elseif linear_threshold
            value_linear, wstar_linear = oracle_linear_threshold(pep, μ);
            println("[optimal value] linear=$(value_linear)");
        end
    end

    for i in 1:length(δs)
        δ = δs[i];
        β = βs[i];
        data = getindex.(datas, i);
        ⋆ = istar(pep, µ);
        if best_arm
            Tstar, wstar, ⋆, lbd, practical = solve(pep, μ, δ, β);
        end
        rule = repeat("-", 60);
        println("");
        println(rule);
        println("$pep at δ = $δ");
        println(@sprintf("%27s", "samples"), " ",
                @sprintf("%6s", "err"), " ",
                @sprintf("%5s", "time"), " ",
                join(map(k -> @sprintf("%15s", k), 1:length(μ))),
        );
        if best_arm
            println(@sprintf("%-42s", "μ"), join(map(x -> @sprintf("%0.4f      ", x), μ)));
            println(@sprintf("%-42s", "w⋆"), join(map(x -> @sprintf("%0.4f      ", x), wstar)));
            println(rule);
            println(@sprintf("%-20s", "oracle"), @sprintf("%7.2f", lbd), @sprintf("%-13s", " "), join(map(w -> @sprintf("%6.0f(%.2f)", lbd*w, w), wstar)), " ");
            println(@sprintf("%-20s", "practical"), @sprintf("%7.2f", practical), @sprintf("%-13s", " "), join(map(w -> @sprintf("%6.0f(%.2f)", practical*w, w), wstar)), " ");
        elseif linear
            println(@sprintf("%-42s", "a\'μ"), join(map(x -> @sprintf("%9.4f", x'µ), pep.arms)));
            println(@sprintf("%-20s", "w⋆_linear"), " ",
                    @sprintf("%7.0f", (1-2δ)*log((1-δ)/δ)*value_linear), " ",
                    @sprintf("%-11s", "--"), " ",
                    join(map(x -> @sprintf("%0.4f      ", x), wstar_linear))
            );
            if linear_bai
                println(@sprintf("%-20s", "w⋆_standard"), " ",
                        @sprintf("%7.0f", (1-2δ)*log((1-δ)/δ)*value_standard), " ",
                        @sprintf("%-11s", "--"), " ",
                        join(map(x -> @sprintf("%0.4f      ", x), wstar_standard))
                );
            end
        elseif lipschitz
            println(@sprintf("%-42s", "μ"), join(map(x -> @sprintf("%0.8f      ", x), μ)));
            println(@sprintf("%-20s", "w⋆_lipschitz"), " ",
                    @sprintf("%7.0f", (1-2δ)*log((1-δ)/δ)*value_lipschitz), " ",
                    @sprintf("%-11s", "--"), " ",
                    join(map(x -> @sprintf("%0.4f      ", x), wstar_lipschitz))
            );
            println(@sprintf("%-20s", "w⋆_standard"), " ",
                    @sprintf("%7.0f", (1-2δ)*log((1-δ)/δ)*value_standard), " ",
                    @sprintf("%-11s", "--"), " ",
                    join(map(x -> @sprintf("%0.4f      ", x), wstar_standard))
            );
        end
        println(rule);

        for r in eachindex(srs)
            Eτ = sum(x->sum(x[2]), data[r,:])/repeats;
            Estd = std(map(x->sum(x[2]), data[r,:]));
            if linear_threshold
                answer = [((pep.arms[k]'µ > pep.threshold) ? 1 : 0) for k=1:length(pep.arms)];
                err = sum(x->(sum(x[1].!=answer)>0 ? 1 : 0), data[r,:])/repeats;
            else # bai
                err = sum(x->x[1].!=⋆, data[r,:])/repeats;
            end
            tim = sum(x->x[3],data[r,:])/repeats;
            println(@sprintf("%-20s", long(srs[r])),
                    @sprintf("%7.0f", Eτ), " ",
                    @sprintf("%0.4f", err), " ",
                    @sprintf("%3.1f", tim/1e6),
                    join(map(k -> @sprintf("%6.0f(%.4f)", sum(x->x[2][k], data[r,:])/repeats, sum(x->x[2][k], data[r,:])/(repeats*Eτ)),
                        1:K)), " "
                    );
            if err > δ
                @warn "too many errors for $(srs[r])";
            end
        end
        println(rule);
    end
end

function τhist(pep, μ, δ, β, srs, data)
    Tstar, wstar, ⋆, lbd, practical = solve(pep, μ, δ, β)
    stephist(map(x -> sum(x[2]), data)', label=permutedims(collect(abbrev.(srs))));
    vline!([lbd], label="lower bd");
    vline!([practical], label="practical");
end

function _boxes(pep, μ, δ, β, srs, data, repeats)
    xs = permutedims(collect(abbrev.(srs)));
    means = sum(sum.(getindex.(data,2)),dims=2)/repeats;
    boxplot(xs, map(x -> sum(x[2]), data)', label="", notch=true, outliers=false, xtickfontsize=15, ytickfontsize=15); # yaxis=:log) # xguidefontsize=30, yguidefontsize=30, legendfontsize=30
    plot!(xs, means', marker=(:star4,10,:black), label="");
end

ext: jl


text:
import Bibliography
using Test

for file in ["test.bib"]
    test_import = Bibliography.import_bibtex("../examples/$file")
    println(typeof(test_import))
    for e in test_import
        println(e)
    end
    println("publications: $(Bibliography.export_web(test_import))")
    println(Bibliography.export_bibtex("result.bib", test_import))
    rm("result.bib")
end


text:
@testset "149.max-points-on-a-line.jl" begin
    @test max_points([[1, 1], [3, 2], [5, 3], [4, 1], [2, 3], [1, 4]]) == 4
    @test max_points([[1, 1], [2, 2], [3, 3]]) == 3
    @test max_points([[0, 0], [94911150, 94911151], [94911151, 94911152]]) == 2
    @test max_points(Vector{Int}[]) == 0
end


text:
# wrappers of the low-level sparse CUSOLVER functionality

using SparseArrays

using LinearAlgebra

using ..CUSPARSE: CuSparseMatrixCSR, CuSparseMatrixCSC, cusparseindex,
                  CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER,
                  CUSPARSE_FILL_MODE_UPPER, CUSPARSE_DIAG_TYPE_NON_UNIT,
                  CUSPARSE_DIAG_TYPE_UNIT, cusparseMatDescr

#csrlsvlu
for (fname, elty, relty) in ((:cusolverSpScsrlsvluHost, :Float32, :Float32),
                             (:cusolverSpDcsrlsvluHost, :Float64, :Float64),
                             (:cusolverSpCcsrlsvluHost, :ComplexF32, :Float32),
                             (:cusolverSpZcsrlsvluHost, :ComplexF64, Float64))
    @eval begin
        function csrlsvlu!(A::SparseMatrixCSC{$elty},
                           b::Vector{$elty},
                           x::Vector{$elty},
                           tol::$relty,
                           reorder::Cint,
                           inda::Char)
            cuinda = cusparseindex(inda)
            n = size(A,1)
            if size(A,2) != n
                throw(DimensionMismatch("LU factorization is only possible for square matrices!"))
            end
            if size(A,2) != length(b)
                throw(DimensionMismatch("second dimension of A, $(size(A,2)), must match the length of b, $(length(b))"))
            end
            if length(x) != length(b)
                throw(DimensionMismatch("length of x, $(length(x)), must match the length of b, $(length(b))"))
            end

            Mat     = similar(A)
            transpose!(Mat, A)
            cudesca = cusparseMatDescr(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            rcudesca = Ref(cudesca)
            singularity = Ref{Cint}(1)
            $fname(sparse_handle(), n, length(A.nzval), rcudesca, Mat.nzval,
                   convert(Vector{Cint},Mat.colptr), convert(Vector{Cint},Mat.rowval), b,
                   tol, reorder, x, singularity)

            if singularity[] != -1
                throw(SingularException(singularity[]))
            end

            x
        end
    end
end

#csrlsvqr
for (fname, elty, relty) in ((:cusolverSpScsrlsvqr, :Float32, :Float32),
                             (:cusolverSpDcsrlsvqr, :Float64, :Float64),
                             (:cusolverSpCcsrlsvqr, :ComplexF32, :Float32),
                             (:cusolverSpZcsrlsvqr, :ComplexF64, Float64))
    @eval begin
        function csrlsvqr!(A::CuSparseMatrixCSR{$elty},
                           b::CuVector{$elty},
                           x::CuVector{$elty},
                           tol::$relty,
                           reorder::Cint,
                           inda::Char)
            cuinda = cusparseindex(inda)
            n = size(A,1)
            if size(A,2) != n
                throw(DimensionMismatch("QR factorization is only possible for square matrices!"))
            end
            if size(A,2) != length(b)
                throw(DimensionMismatch("second dimension of A, $(size(A,2)), must match the length of b, $(length(b))"))
            end
            if length(x) != length(b)
                throw(DimensionMismatch("length of x, $(length(x)), must match the length of b, $(length(b))"))
            end

            cudesca = cusparseMatDescr(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            rcudesca = Ref(cudesca)
            singularity = Ref{Cint}(1)
            $fname(sparse_handle(), n, A.nnz, rcudesca, A.nzVal, A.rowPtr, A.colVal, b, tol, reorder, x, singularity)

            if singularity[] != -1
                throw(SingularException(singularity[]))
            end

            x
        end
    end
end

#csrlsvchol
for (fname, elty, relty) in ((:cusolverSpScsrlsvchol, :Float32, :Float32),
                             (:cusolverSpDcsrlsvchol, :Float64, :Float64),
                             (:cusolverSpCcsrlsvchol, :ComplexF32, :Float32),
                             (:cusolverSpZcsrlsvchol, :ComplexF64, Float64))
    @eval begin
        function csrlsvchol!(A::CuSparseMatrixCSR{$elty},
                             b::CuVector{$elty},
                             x::CuVector{$elty},
                             tol::$relty,
                             reorder::Cint,
                             inda::Char)
            cuinda = cusparseindex(inda)
            n      = size(A,1)
            if size(A,2) != n
                throw(DimensionMismatch("Cholesky factorization is only possible for square matrices!"))
            end
            if size(A,2) != length(b)
                throw(DimensionMismatch("second dimension of A, $(size(A,2)), must match the length of b, $(length(b))"))
            end
            if length(x) != length(b)
                throw(DimensionMismatch("length of x, $(length(x)), must match the length of b, $(length(b))"))
            end

            cudesca = cusparseMatDescr(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            rcudesca = Ref(cudesca)
            singularity = zeros(Cint,1)
            $fname(sparse_handle(), n, A.nnz, rcudesca, A.nzVal, A.rowPtr, A.colVal, b, tol, reorder, x, singularity)

            if singularity[1] != -1
                throw(SingularException(singularity[1]))
            end

            x
        end
    end
end

#csrlsqvqr
for (fname, elty, relty) in ((:cusolverSpScsrlsqvqrHost, :Float32, :Float32),
                             (:cusolverSpDcsrlsqvqrHost, :Float64, :Float64),
                             (:cusolverSpCcsrlsqvqrHost, :ComplexF32, :Float32),
                             (:cusolverSpZcsrlsqvqrHost, :ComplexF64, Float64))
    @eval begin
        function csrlsqvqr!(A::SparseMatrixCSC{$elty},
                            b::Vector{$elty},
                            x::Vector{$elty},
                            tol::$relty,
                            inda::Char)
            cuinda = cusparseindex(inda)
            m,n    = size(A)
            if m < n
                throw(DimensionMismatch("csrlsqvqr only works when the first dimension of A, $m, is greater than or equal to the second dimension of A, $n"))
            end
            if size(A,2) != length(b)
                throw(DimensionMismatch("second dimension of A, $(size(A,2)), must match the length of b, $(length(b))"))
            end
            if length(x) != length(b)
                throw(DimensionMismatch("length of x, $(length(x)), must match the length of b, $(length(b))"))
            end

            cudesca  = cusparseMatDescr(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            rcudesca = Ref(cudesca)
            p        = zeros(Cint,n)
            min_norm = zeros($relty,1)
            rankA    = zeros(Cint,1)
            Mat      = similar(A)
            transpose!(Mat, A)
            $fname(sparse_handle(), m, n, length(A.nzval), rcudesca, Mat.nzval,
                   convert(Vector{Cint},Mat.colptr), convert(Vector{Cint},Mat.rowval),
                   b, tol, rankA, x, p, min_norm)

            x, rankA[1], p, min_norm[1]
        end
    end
end

#csreigvsi
for (fname, elty, relty) in ((:cusolverSpScsreigvsi, :Float32, :Float32),
                             (:cusolverSpDcsreigvsi, :Float64, :Float64),
                             (:cusolverSpCcsreigvsi, :ComplexF32, :Float32),
                             (:cusolverSpZcsreigvsi, :ComplexF64, Float64))
    @eval begin
        function csreigvsi(A::CuSparseMatrixCSR{$elty},
                           μ_0::$elty,
                           x_0::CuVector{$elty},
                           tol::$relty,
                           maxite::Cint,
                           inda::Char)
            cuinda = cusparseindex(inda)
            m,n    = size(A)
            if m != n
                throw(DimensionMismatch("A must be square!"))
            end
            if n != length(x_0)
                throw(DimensionMismatch("second dimension of A, $(size(A,2)), must match the length of x_0, $(length(x_0))"))
            end

            cudesca = cusparseMatDescr(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            rcudesca = Ref(cudesca)
            x       = copy(x_0)
            μ       = CuArrays.zeros($elty,1)
            $fname(sparse_handle(), n, A.nnz, rcudesca, A.nzVal, A.rowPtr, A.colVal,
                   μ_0, x_0, maxite, tol, μ, x)

            collect(μ)[1], x
        end
    end
end

#csreigs
for (fname, elty, relty) in ((:cusolverSpScsreigsHost, :ComplexF32, :Float32),
                             (:cusolverSpDcsreigsHost, :ComplexF64, :Float64),
                             (:cusolverSpCcsreigsHost, :ComplexF32, :ComplexF32),
                             (:cusolverSpZcsreigsHost, :ComplexF64, :ComplexF64))
    @eval begin
        function csreigs(A::SparseMatrixCSC{$relty},
                         lbc::$elty,
                         ruc::$elty,
                         inda::Char)
            cuinda = cusparseindex(inda)
            m,n    = size(A)
            if m != n
                throw(DimensionMismatch("A must be square!"))
            end

            cudesca = cusparseMatDescr(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
            rcudesca = Ref(cudesca)
            numeigs = zeros(Cint,1)
            Mat     = similar(A)
            transpose!(Mat, A)
            $fname(sparse_handle(), n, length(A.nzval), rcudesca, Mat.nzval,
                   convert(Vector{Cint},Mat.colptr), convert(Vector{Cint},Mat.rowval),
                   lbc, ruc, numeigs)

            numeigs[1]
        end
    end
end


text:
# basic integral struct with upper bound and lower bound.
struct Integral{X, T <: Symbolics.VarDomainPairing} <: Function
    x
    domain::T
    Integral(x,domain) = new{typeof(x),typeof(domain)}(Symbolics.value.(x), domain)
end

(I::Integral)(x) = Term{SymbolicUtils.symtype(x)}(I, [x])
(I::Integral)(x::Num) = Num(I(Symbolics.value(x)))
SymbolicUtils.promote_symtype(::Integral, x) = x

function Base.show(io::IO, I::Integral)
    print(io, "Integral(", I.x, ", ", I.domain, ")")
end

Base.:(==)(I1::Integral, I2::Integral) = convert(Bool, simplify(isequal(I1.x, I2.x) && isequal(I1.domain, I2.domain)))


text:
function vmap_quote(N, ::Type{T}) where {T}
    W, Wshift = VectorizationBase.pick_vector_width_shift(T)
    val = Expr(:curly, :Val, W)
    q = Expr(:block, Expr(:(=), :M, Expr(:call, :length, :dest)), Expr(:(=), :vdest, Expr(:call, :vectorizable, :dest)), Expr(:(=), :m, 0))
    fcall = Expr(:call, :f)
    loopbody = Expr(:block, Expr(:call, :vstore!, :vdest, fcall, :m), Expr(:(+=), :m, W))
    fcallmask = Expr(:call, :f)
    bodymask = Expr(:block, Expr(:(=), :__mask__, Expr(:call, :mask, val, Expr(:call, :&, :M, W-1))), Expr(:call, :vstore!, :vdest, fcallmask, :m, :__mask__))
    for n ∈ 1:N
        arg_n = Symbol(:varg_,n)
        push!(q.args, Expr(:(=), arg_n, Expr(:call, :vectorizable, Expr(:ref, :args, n))))
        push!(fcall.args, Expr(:call, :vload, Expr(:call, val), arg_n, :m))
        push!(fcallmask.args, Expr(:call, :vload, Expr(:call, val), arg_n, :m, :__mask__))
    end
    loop = Expr(:for, Expr(:(=), :_, Expr(:call, :(:), 0, Expr(:call, :-, Expr(:call, :(>>>), :M, Wshift), 1))), loopbody)
    push!(q.args, loop)
    ifmask = Expr(:if, Expr(:call, :(!=), :m, :M), bodymask)
    push!(q.args, ifmask)
    push!(q.args, :dest)
    q
end
@generated function vmap!(f::F, dest::AbstractArray{T}, args::Vararg{<:AbstractArray,N}) where {F,T,N}
    vmap_quote(N, T)
end
function vmap(f::F, args...) where {F}
    T = Base._return_type(f, Base.Broadcast.eltypes(args))
    dest = similar(first(args), T)
    vmap!(f, dest, args...)
end


text:
using Test
using Random

using MinimalRLCore

# function test_hashlesstilecoder()

#     @testset "Test HashlessTileCoder" begin
#         tc = MinimalRLCore.HashlessTileCoder([2,2], [0 0 ; 1 1], 4)
#         Random.seed!(1)
#         s = rand(2)
#         @test MinimalRLCore.feature_size(tc) == 36
#         @test all(MinimalRLCore.create_features(tc, s) .== [1, 13, 22, 29])
#     end

#     @testset "Test HashlessTileCoderFull" begin
#         tc = MinimalRLCore.HashlessTileCoderFull([2,2], [0 0 ; 1 1], 4)
#         Random.seed!(1)
#         s = rand(2)
#         @test MinimalRLCore.feature_size(tc) == 36
#         @test all(MinimalRLCore.create_features(tc, s) .== [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0])
#     end

# end


text:
#=
ellipsoid_im.jl
2019-07 Helena H
=#

export ellipsoid_im

using Plots
using Test: @test, @test_throws
#using MIRT: jim, MIRT_image_geom, image_geom

"""
    phantom = ellipsoid_im(ig, params;
	oversample=1, checkfov=false, how=:slow, showmem=false, hu_scale=1, return_params=false)`

generate ellipsoid phantom image from parameters:
	`[x_center, y_center, z_center, x_radius, y_radius, z_radius,
		xy_angle_degrees, z_angle_degrees, density]`

in
- `ig`			from `image_geom()`
- `params`		`[N 9]` ellipsoid parameters.

option
- `oversample::Int`		oversampling factor (default:1)
- `checkfov::Bool`		warn if any ellipsoid is out of fov
- `how::Symbol`				`:fast` does it fast -- to do, only works slow
							`:lowmem` uses less memory than :fast but slower
							`:slow` default
- `showmem::Bool`
- `hu_scale::Real`			use 1000 to scale shepp-logan to HU; default 1
- `return_params::Bool`		if true, return both phantom and params

out
- `phantom`		`[nx ny nz]` image
- `params`		`[N 9]` ellipsoid parameters (only if `return_params=true`)
"""
function ellipsoid_im(ig::MIRT_image_geom, params::AbstractMatrix{<:Real} ;
		oversample::Int = 1,
		checkfov::Bool = false,
		how::Symbol = :slow,
		showmem::Bool = false,
		hu_scale::Real = 1,
		return_params::Bool = false)

	size(params,2) != 9 && throw("bad cuboid parameter vector size")
	params[:,9] .*= hu_scale

	if oversample > 1
		ig = ig.over(oversample)
	end

	args = (ig.nx, ig.ny, ig.nz, params, ig.dx, ig.dy, ig.dz,
			ig.offset_x, ig.offset_y, ig.offset_z)
	phantom = zeros(Float32, ig.nx, ig.ny, ig.nz)

	checkfov && !ellipsoid_im_check_fov(args...) && throw("ellipsoid exceeds FOV")

	if how === :slow
		phantom += ellipsoid_im_slow(args..., showmem)
#=
	elseif how === :lowmem
		phantom += ellipsoid_im_lowmem(args..., showmem)
	elseif how === :fast
		phantom += ellipsoid_im_fast(args..., showmem)
=#
	else
		throw("bad how $how")
	end

	if oversample > 1
		phantom = downsample3(phantom, oversample)
	end

	if return_params
		return (phantom, params)
	end
	return phantom
end

"""
ellipsoid_im_slow()

brute force fine grid - can use lots of memory
"""
function ellipsoid_im_slow(nx, ny, nz, params, dx, dy, dz,
		offset_x, offset_y, offset_z, showmem)

	phantom = zeros(Float32, nx, ny, nz)

	wx = (nx - 1)/2 + offset_x
	wy = (ny - 1)/2 + offset_y
	wz = (nz - 1)/2 + offset_z
	x = ((0:(nx-1)) .- wx) * dx
	y = ((0:(ny-1)) .- wy) * dy
	z = ((0:(nz-1)) .- wz) * dz
	xmax = maximum(x)
	xmin = minimum(x)
	ymax = maximum(y)
	ymin = minimum(y)
	zmax = maximum(z)
	zmin = minimum(z)
	(xx, yy, zz) = ndgrid(x, y, z)

	#ticker reset
	np = size(params)[1]
	for ip in 1:np
		#ticker(ip, np)

		par = params[ip, :]
		cx = par[1]
		cy = par[2]
		cz = par[3]
		rx = par[4]
		ry = par[5]
		rz = par[6]

		azim = deg2rad(par[7])
		polar = deg2rad(par[8])
		(xr, yr, zr) = rot3(xx .- cx, yy .- cy, zz .- cz, azim, polar)
		tmp = ((xr / rx).^2 + (yr / ry).^2 + (zr / rz).^2) .<= 1
		value = Float32(par[9])
		phantom += value * tmp
	end
	return phantom
end

#=
"""
ellipsoid_im_fast()

currently not working
only slow option works
"""
function ellipsoid_im_fast(nx, ny, nz, params, dx, dy, dz,
	offset_x, offset_y, offset_z, over, showmem)

	phantom = zeros(Float32, nx, ny, nz)

	wx = (nx-1)/2 + offset_x
	wy = (ny-1)/2 + offset_y
	wz = (nz-1)/2 + offset_z
	x = ((0:nx-1) .- wx) * dx
	y = ((0:ny-1) .- wy) * dy
	z = ((0:nz-1) .- wz) * dz
	xmax = maximum(x)
	xmin = minimum(x)
	ymax = maximum(y)
	ymin = minimum(y)
	zmax = maximum(z)
	zmin = minimum(z)
	(xx, yy, zz) = ndgrid(x, y, z)

	if over > 1
		tmp = ((1:over) - (over+1)/2) / over
		(xf, yf, zf) = ndgrid(tmp*dx, tmp*dy, tmp*dz)
		xf = vec(xf)'
		yf = vec(yf)'
		zf = vec(zf)'
	end

	hx = abs(dx) / 2
	hy = abs(dy) / 2
	hz = abs(dz) / 2

	#ticker reset
	np = size(params)[1]
	for ip in 1:np
		#ticker(mfilename, ip, np)
		par = params[ip,:]
		cx = par[1]
		cy = par[2]
		cz = par[3]
		rx = par[4]
		ry = par[5]
		rz = par[6]

		azim = par[7] * (pi/180)
		polar = par[8] * (pi/180)

		xs = xx .- cx
		ys = yy .- cy
		zs = zz .- cz

		(xr, yr, zr) = rot3(xs, ys, zs, azim, polar)
		if over == 1
			vi = ((xr / rx).^2 + (yr / ry).^2 + (zr / rz).^2) .<= 1
			value = Float32(par[9])
			phantom += value * vi
		end

		#= check all 8 corners of the cube that bounds
		all possible 3D rotations of the voxel =#
		hh = sqrt(hx^2 + hy^2 + hz^2)
		vi = true # voxels entirely inside the ellipsoid
		#vo = true # voxels entirely outside the ellipsoid

		for ix in -1:1, iy in -1:1, iz in -1:1
			xo = xr .+ ix*hh
			yo = yr .+ iy*hh
			zo = zr .+ iz*hh
			vi = vi .& (((xo/rx).^2 + (yo/ry).^2 + (zo/rz).^2) .< 1)
		end

		vo = false #to do. for now, "outside" test is failing
		if ip == 1
			throw("todo:must debug this")
		end
		if any(vec(vi) & vec(vo))
			throw("bug")
		end
		#=
		if 0
			% coordinates of "outer" corner of each voxel, relative to ellipsoid center
			sign_p = @(x) (x >= 0) * sqrt(3); % modified sign()
			xo = xr + sign_p(xr) * hx;
			yo = yr + sign_p(yr) * hy;
			zo = zr + sign_p(zr) * hz;

			% voxels that are entirely inside the ellipsoid
			vi = (xo / rx).^2 + (yo / ry).^2 + (zo / rz).^2 <= 1;
		end

		if 0
			% coordinates of "inner" corner of each pixel, relative to ellipse center
			sign_n = @(x) (x > 0) * sqrt(3); % modified sign()
			xi = xr - sign_n(xs) * hx;
			yi = yr - sign_n(ys) * hy;
			zi = zr - sign_n(zs) * hz;

			% voxels that are entirely outside the ellipsoid
			vo = (max(abs(xi),0) / rx).^2 + (max(abs(yi),0) / ry).^2 ...
				+ (max(abs(zi),0) / rz).^2 >= 1;
		end
		=#

		edge = !vi & !vo
		if showmem
			println("edge fraction ",
				sum(vec(edge)) / prod(size(edge)), "=",
				sum(vec(edge)), "/", prod(size(edge)))
		end

		x = xx[edge] .- cx
		y = yy[edge] .- cy
		z = zz[edge] .- cz
		x = x .+ xf'
		y = y .+ yf'
		z = z .+ zf'
		(xr, yr, zr) = rot3(x, y, z, azim, polar)
		in = ((xr/rx).^2 + (yr/ry).^2 + (zr/rz).^2) .<= 1
		tmp = mean(in, 2)
		gray = Float32(vi)
		gray[edge] = tmp
		#=
		if 0 % todo: help debug
			ee = (gray > 0) & (gray < 1);
			im(ee)
			im(vi)
		prompt
			%minmax(tmp)
			%clf, im pl 1 2
			%im(tmp)
		end=#

		value = Float32(par[9])
		phantom += value*gray
	end
	return phantom
end

"""
ellipsoid_im_lowmem()

Do 'one slice at a time' to reduce memory
"""
function ellipsoid_im_lowmem(nx, ny, nz, params, dx, dy, dz,
		offset_x, offset_y, offset_z)

	phantom = zeros(Float32, nx, ny, nz)
	for iz in 1:nz
		offset_z_new = (nz-1)/2 + offset_z - (iz-1)
		phantom[:,:,iz] = ellipsoid_im_fast(nx, ny, 1, params, dx, dy, dz,
		offset_x, offset_y, offset_z_new, showmem && iz == 1)
	end
	return phantom
end
=#

"""
ellipsoid_im_check_fov()
"""
function ellipsoid_im_check_fov(nx, ny, nz, params,
		dx, dy, dz, offset_x, offset_y, offset_z)
	wx = (nx - 1)/2 + offset_x
	wy = (ny - 1)/2 + offset_y
	wz = (nz - 1)/2 + offset_z
	xx = ((0:nx - 1) .- wx) * dx
	yy = ((0:ny - 1) .- wy) * dy
	zz = ((0:nz - 1) .- wz) * dz

	xmax = maximum(xx)
	xmin = minimum(xx)
	ymax = maximum(yy)
	ymin = minimum(yy)
	zmax = maximum(zz)
	zmin = minimum(zz)

	for ip in 1:size(params)[1]
		par = params[ip, :]
		cx = par[1]
		cy = par[2]
		cz = par[3]
		rx = par[4]
		ry = par[5]
		rz = par[6]

		if (cx + rx > xmax) || (cx - rx < xmin)
			isinteractive() && @warn("fov: x range $xmin $xmax, cx=$cx, rx=$rx")
			return false
		end
		if (cy + ry > ymax) || (cy - ry < ymin)
			isinteractive() && @warn("fov: y range $ymin $ymax, cy=$cy, ry=$ry")
			return false
		end

		if (cz + rz > zmax) || (cz - rz < zmin)
			isinteractive() && @warn("fov: z range $zmin $zmax, cz=$cz, rz=$rz")
			return false
		end
	end
	return true
end

"""
    phantom = ellipsoid_im(ig, ptype ; args...)

`ptype = :zhu | :kak | :e3d | :spheroid`
"""
function ellipsoid_im(ig::MIRT_image_geom, ptype::Symbol ; args...)
	xfov = ig.fovs[1]
	yfov = ig.fovs[2]
	zfov = ig.fovs[3]

	if ptype === :zhu
		params = shepp_logan_3d_parameters(xfov, yfov, zfov, :zhu)
	elseif ptype === :kak
		params = shepp_logan_3d_parameters(xfov, yfov, zfov, :kak)
#=
	elseif ptype === :e3d
		params = shepp_logan_3d_parameters(xfov, yfov, zfov, :e3d)
=#
	elseif ptype === :spheroid
		params = spheroid_params(xfov, yfov, zfov, ig.dx, ig.dy, ig.dz)
	else
		throw("bad phantom symbol $ptype")
	end
	return ellipsoid_im(ig, params ; args...)
end

"""
    phantom = ellipsoid_im(ig ; args...)

`:zhu` (default) for given image geometry `ig`
"""
function ellipsoid_im(ig::MIRT_image_geom ; args...)
	return ellipsoid_im(ig, :zhu ; args...)
end

"""
    rot3
"""
function rot3(x, y, z, azim, polar)
	polar != 0 && throw("z (polar) rotation not done")
	xr = cos(azim) * x + sin(azim) * y
	yr = -sin(azim) * x + cos(azim) * y
	zr = z
	return (xr, yr, zr)
end

"""
shepp_logan_3d_parameters()

most of these values are unitless 'fractions of field of view'
"""
function shepp_logan_3d_parameters(xfov, yfov, zfov, ptype)
	# parameters from Kak and Slaney, typos ?
	ekak = Float32[
		0		0		0		0.69	0.92	0.9		0	0	2.0
		0		0		0		0.6624	0.874	0.88	0	0	-0.98
		-0.22	0		-0.25	0.41	0.16	0.21	0	0	-0.98
		0.22	0		-0.25	0.31	0.11	0.22	72	0	-0.02
		0		0.1		-0.25	0.046	0.046	0.046	0	0	0.02
		0		0.1		-0.25	0.046	0.046	0.046	0	0	0.02
		-0.8	-0.65	-0.25	0.046	0.023	0.02	0	0	0.01
		0.06	-0.065	 -0.25	0.046	0.023	0.02	90	0	0.01
		0.06	-0.105	 0.625	0.56	0.04	0.1		90	0	0.02
		0		0.1		-0.625	0.056	0.056	0.1		0	0	-0.02
	]

	#parameters from leizhu@standford.edu, says kak&slaney are incorrect
	ezhu = Float32[
	0		0		0		0.69	0.92	0.9		0	0	2.0
	0		-0.0184	0		0.6624	0.874	0.88	0	0	-0.98
	-0.22	0		-0.25	0.41	0.16	0.21	-72	0	-0.02
	0.22	0		-0.25	0.31	0.11	0.22	72	0	-0.02
	0		0.35	-0.25	0.21	0.25	0.35	0	0	0.01
	0		0.1		-0.25	0.046	0.046	0.046	0	0	0.01
	-0.08	-0.605	-0.25	0.046	0.023	0.02	0	0	0.01
	0		-0.1	-0.25	0.046	0.046	0.046	0	0	0.01
	0		-0.605	-0.25	0.023	0.023	0.023	0	0	0.01
	0.06	-0.605	-0.25	0.046	0.023	0.02	-90	0	0.01
	0.06	-0.105	0.0625	0.056	0.04	0.1		-90	0	0.02
	0		0.1		0.625	0.056	0.056	0.1		0	0	-0.02
	]

	if ptype === :zhu
		params = ezhu
	elseif ptype === :kak
		params = ekak
#=
	elseif ptype === :shepp_logan_e3d || ptype === :e3d
		params = e3d[:, [5:7 2:4 8 1]]
=#
	else
		throw("unknown parameter type $ptype")
	end

	params[:,[1,4]] .*= xfov/2
	params[:,[2,5]] .*= yfov/2
	params[:,[3,6]] .*= zfov/2
	return params
end

"""
spheroid_params()
"""
function spheroid_params(xfov, yfov, zfov, dx, dy, dz)
#	xfov = nx * dx, number * size
	xradius = (xfov/2) - dx
	yradius = (yfov/2) - dy
	zradius = (zfov/2) - dz
	return [ 0 0 0 xradius yradius zradius 0 0 1 ]
end

"""
    ellipsoid_im()

show docstring
"""
function ellipsoid_im()
	@doc ellipsoid_im
end

"""
ellipsoid_im_show()
"""
function ellipsoid_im_show()
	ig = image_geom(nx=512, nz=64, dz=0.625, fov=500)
	ig = ig.down(8)

	spheroid = ellipsoid_im(ig, :spheroid)
	p1 = jim(spheroid, title="spheroid")

	x = ellipsoid_im(ig, :zhu; hu_scale=1000)
	p2 = jim(x, title="zhu", clim=(900,1100))

	plot(p1, p2)
end

"""
ellipsoid_im_test()
"""
function ellipsoid_im_test()
	ig = image_geom(nx=512, nz=64*2, dz=0.625, fov=500)
	ig = ig.down(8)

	#test different ellipses
	ellipsoid_im(ig)
	ellipsoid_im(ig, :zhu)
	ellipsoid_im(ig, :kak, oversample=2)
	ellipsoid_im(ig, :spheroid, return_params=true)
	ellipsoid_im(ig; checkfov=true)
	ellipsoid_im(ig; showmem=true)
	@test_throws String ellipsoid_im(ig, :bad)
	@test_throws String ellipsoid_im(ig, :spheroid, how=:bad)
	@test_throws String shepp_logan_3d_parameters(0, 0, 0, :bad)
	tmp = [1, 1, 1, 0, 0, 0]
	@test !ellipsoid_im_check_fov(2, 9, 9, [3 0 0 1 1 1 0 0 1], tmp...)
	@test !ellipsoid_im_check_fov(9, 2, 9, [0 3 0 1 1 1 0 0 1], tmp...)
	@test !ellipsoid_im_check_fov(9, 9, 2, [0 0 3 1 1 1 0 0 1], tmp...)

	#ell1 = ellipsoid_im(ig, :zhu; how=:fast) # fast doesn't work
	#ell2 = ellipsoid_im(ig, :zhu; how=:lowmem) # lowmem calls fast - doesn't work
	# todo: test slow vs fast
	true
end

"""
    ellipsoid_im(:test)

`ellipsoid_im(:show)`

run tests
"""
function ellipsoid_im(test::Symbol)
	if test === :show
		return ellipsoid_im_show()
	end
	test != :test && throw(ArgumentError("test $test"))
	ellipsoid_im()
	ellipsoid_im(:show)
	@test ellipsoid_im_test()
	true
end

#ellipsoid_im(:show)
#ellipsoid_im(:test)


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    int{L} x = 7;
    int z = 0;
    void m() {
	final Condition c = new Condition();
	int y = ({L /c T}) 42;
	c.set();
	z = y;
    }

    public static void main(String[] args) {
	new C().m();
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


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

Makie.plottype(::AbstractVector{<:PointOrGeometry}) =
  Viz{<:Tuple{AbstractVector{<:PointOrGeometry}}}

function Makie.plot!(plot::Viz{<:Tuple{AbstractVector{<:PointOrGeometry}}})
  # retrieve items object
  items = plot[:object][]

  # fallback to collection recipe
  viz!(plot, Collection(items),
    size         = plot[:size],
    color        = plot[:color],
    colormap     = plot[:colormap],
    facetcolor   = plot[:boundarycolor],
    showfacets   = plot[:showboundary],
    decimation   = plot[:decimation],
  )
end

Makie.plottype(::PointOrGeometry) = Viz{<:Tuple{PointOrGeometry}}

function Makie.plot!(plot::Viz{<:Tuple{PointOrGeometry}})
  # retrieve item object
  item = plot[:object][]

  # fallback to vector recipe
  viz!(plot, [item],
    size          = plot[:size],
    color         = plot[:color],
    colormap      = plot[:colormap],
    boundarycolor = plot[:boundarycolor],
    showboundary  = plot[:showboundary],
    decimation    = plot[:decimation],
  )
end

Makie.plottype(::Domain) = Viz{<:Tuple{Domain}}

function Makie.plot!(plot::Viz{<:Tuple{Domain}})
  # retrieve domain object
  domain = plot[:object][]

  # fallback to vector recipe
  viz!(plot, collect(domain),
    size          = plot[:size],
    color         = plot[:color],
    colormap      = plot[:colormap],
    boundarycolor = plot[:facetcolor],
    showboundary  = plot[:showfacets],
    decimation    = plot[:decimation],
  )
end


text:
"""
Solving the assigment problem.
Introduce a costs square matrix.
The function give the minimun costs posible and the elements of the costs matrix that minimize the costs of the task assigment
"""
function solve_assigm_prob(costs_matrix)
    N=size(costs_matrix)[1];
    if N==size(costs_matrix)[2];
        modelo = JuMP.Model(JuMP.with_optimizer(GLPK.Optimizer))
        @variable(modelo, x[1:N,1:N] >=0 )

        @objective(modelo, Min, sum(costs_matrix.*x))
        @constraint(modelo, rows[i=1:N], sum(x[i,1:N])==1)
        @constraint(modelo, cols[j=1:N], sum(x[1:N,j])==1)
        optimize!(modelo)
        println("El coste optimo es: ")
        println(objective_value(modelo))
        for i=1:N
            for j=1:N
                if value(x[i,j])!=0
                    println("C[$i,$j]=",costs_matrix[i,j])
                end
            end
        end

    else println("La matriz debe ser cuadrada")
    end
end


text:
"A collection of functions and types for manipulating image data and doing computer vision."
module Bildefin

include("utils.jl")
include("opencv-wrapper.jl")

end


text:
using Plots, Random

function isInDisc(P, r)
    return norm(P) <= r
end

function disc_plot(x; color=:red)
    scatter(x[:,1], x[:,2]; c=color, aspect_ratio=1, legend=false, grid=false, mswidth=0, xlims=[-1.1, 1.1], ylims=[-1.1, 1.1])
end

function polar2cart(r, θ)
    return hcat(r .* cos.(θ), r .* sin.(θ))
end

function cart2polar(x)
    N = size(x,1)
    r = [norm(x[i,:]) for i in 1:N]
    θ = [atan(x[i,2], x[i,1]) for i in 1:N]
    return r, θ
end

## Annealed Networks on a Bounded Disc
Random.seed!(2020)

# rotation
N = 256; iter = 5000
θx = 2π*rand(N)
rx = sqrt.(rand(N))
x = polar2cart(rx, θx)
disc_plot(x)

# 1. Learning steps
α = 0.1
anim = @animate for n in 1:iter
    # update y
    global θx, rx, x
    θ = 2π*rand()
    θy = (θx .+ θ) .% 2π
    ry = rx
    y = polar2cart(ry, θy)
    # find Λⱼ(n)
    nnx = Dict()
    for j in 1:N
        nnx[j] = []
    end
    for i in 1:N
        arr = [norm(x[j,:] - y[i,:]) for j in 1:N]
        min_val = minimum(arr)
        min_ind = findall(arr .== min_val)
        for j in min_ind
            nnx[j] = push!(nnx[j], i)
        end
    end
    for j in 1:N
        if nnx[j] != []
            x[j,:] += α * sum([y[i,:]-x[j,:] for i in nnx[j]])
        end
    end
    rx, θx = cart2polar(x)
    if n % 50 == 0
        disc_plot(x)
        title!("iter=$(n)")
    end
end when n % 50 == 0

gif(anim, "selfSimNN_2D_rotation.gif", fps = 5)

## dilation
N = 256; iter = 100000
θx = 2π*rand(N)
rx = sqrt.(rand(N))
x = polar2cart(rx, θx)
disc_plot(x)

# 1. Learning steps
α = 0.01
anim = @animate for n in 1:iter
    global α
    if n > .99*iter
        α -= 0.01/(.01*iter)
    end
    # update y
    global θx, rx, x
    d = (rand()-0.5)*2*log(16)
    θy = θx .+ (d < 0)*π
    ry = abs.(d) .* rx
    y = polar2cart(ry, θy)
    # find Λⱼ(n)
    nnx = Dict()
    for j in 1:N
        nnx[j] = []
    end
    for i in 1:N
        if ry[i] > 1
            continue
        end
        arr = [norm(x[j,:] - y[i,:]) for j in 1:N]
        min_val = minimum(arr)
        min_ind = findall(arr .== min_val)
        for j in min_ind
            nnx[j] = push!(nnx[j], i)
        end
    end
    for j in 1:N
        if nnx[j] != []
            x[j,:] += α * sum([y[i,:]-x[j,:] for i in nnx[j]])
        end
    end
    rx, θx = cart2polar(x)
    rx = min.(rx, 1)
    x = polar2cart(rx, θx)
    if n % 1000 == 0
        disc_plot(x)
        title!("iter=$(n)")
    end
end when n % 1000 == 0

gif(anim, "selfSimNN_2D_dilation3.gif", fps = 5)

disc_plot(x)

## rotation + dilation
N = 256; iter = 100000
θx = 2π*rand(N)
rx = sqrt.(rand(N))
x = polar2cart(rx, θx)
disc_plot(x)

# 1. Learning steps
α = 0.01
anim = @animate for n in 1:iter
    global α
    if n > .99*iter
        α -= 0.01/(.01*iter)
    end
    # update y
    global θx, rx, x
    θ = 2π*rand()
    d = (rand()-0.5)*2*log(16)
    θy = (θx .+ θ .+ (d < 0)*π) .% 2π
    ry = abs.(d) .* rx
    y = polar2cart(ry, θy)
    # find Λⱼ(n)
    nnx = Dict()
    for j in 1:N
        nnx[j] = []
    end
    for i in 1:N
        if !isInDisc(y[i,:], 1)
            continue
        end
        arr = [norm(x[j,:] - y[i,:]) for j in 1:N]
        min_val = minimum(arr)
        min_ind = findall(arr .== min_val)
        for j in min_ind
            nnx[j] = push!(nnx[j], i)
        end
    end
    for j in 1:N
        if nnx[j] != []
            x[j,:] += α * sum([y[i,:]-x[j,:] for i in nnx[j]])
        end
    end
    rx, θx = cart2polar(x)
    rx = min.(rx, 1)
    x = polar2cart(rx, θx)
    if n % 1000 == 0
        disc_plot(x)
        title!("iter=$(n)")
    end
end when n % 1000 == 0

gif(anim, "selfSimNN_2D_rotdil2.gif", fps = 5)

plt = disc_plot(x)
savefig(plt, "figs/selfSimNN_2D_rotdil_final2.png")

using NGWP
_, _, X = SunFlowerGraph(; N=256); sf_plt = disc_plot(X)
savefig(sf_plt, "figs/sunflower_arangement.png")


text:
using SimpleTest
using YAML
println("Enter two numbers")
num1 = parse(Float64, readline())
num2 = parse(Float64, readline())
result = simple_operation(num1, num2)
println("The sum is ", result)
sep = "/"
working_path = pwd()
settings_path = joinpath(working_path, "settings.yml")
testpath = joinpath(pwd(), "src")
push!(LOAD_PATH, testpath)
settings = YAML.load(open(settings_path))
output_path = joinpath(working_path, "output")
start_of_opt(settings, sep, output_path)


text:
abstract type LenseFlowOp{I<:ODESolver,t₀,t₁,T} <: FlowOpWithAdjoint{I,t₀,t₁,T} end

# `L = LenseFlow(ϕ)` just creates a wrapper holding ϕ. Then when you do `L*f` or
# `cache(L,f)` we create a CachedLenseFlow object which holds all the
# precomputed quantities and preallocated memory needed to do the lense.

"""
    LenseFlow(ϕ, [n=7])

`LenseFlow` is the ODE-based lensing algorithm from [Millea, Anderes,
& Wandelt, 2019](https://arxiv.org/abs/1708.06753). The number of
steps in the ODE solver is controlled by `n`. The action of the
operator, as well as its adjoint, inverse, inverse-adjoint, and
gradient of any of these w.r.t. `ϕ` can all be computed. The
log-determinant of the operation is zero independent of `ϕ`, in the
limit of `n` high enough.
"""
struct LenseFlow{I<:ODESolver,t₀,t₁,T} <: LenseFlowOp{I,t₀,t₁,T}
    ϕ :: Field
end
LenseFlow(ϕ,n=7) = LenseFlow{RK4Solver{n}}(ϕ)
LenseFlow{I}(ϕ) where {I<:ODESolver} = LenseFlow{I,0,1}(ϕ)
LenseFlow{I,t₀,t₁}(ϕ) where {I,t₀,t₁} = LenseFlow{I,float(t₀),float(t₁),real(eltype(ϕ))}(ϕ)

struct CachedLenseFlow{N,t₀,t₁,ŁΦ<:Field,ÐΦ<:Field,ŁF<:Field,ÐF<:Field,T} <: LenseFlowOp{RK4Solver{N},t₀,t₁,T}

    # save ϕ to know when to trigger recaching
    ϕ :: Ref{Any}

    # p and M⁻¹ quantities precomputed at every time step
    p   :: Dict{Float16,SVector{2,Diagonal{T,ŁΦ}}}
    M⁻¹ :: Dict{Float16,SMatrix{2,2,Diagonal{T,ŁΦ},4}}

    # f type memory
    memŁf  :: ŁF
    memÐf  :: ÐF
    memŁvf :: FieldVector{ŁF}
    memÐvf :: FieldVector{ÐF}

    # ϕ type memory
    memŁϕ  :: ŁΦ
    memÐϕ  :: ÐΦ
    memŁvϕ :: FieldVector{ŁΦ}
    memÐvϕ :: FieldVector{ÐΦ}

end

### printing
typealias_def(::Type{<:RK4Solver{N}}) where {N} = "$N-step RK4"
typealias_def(::Type{<:CachedLenseFlow{N,t₀,t₁,ŁΦ,<:Any,ŁF}}) where {N,t₀,t₁,ŁΦ,ŁF} =
    "CachedLenseFlow{$t₀→$t₁, $N-step RK4}(ϕ::$(typealias(ŁΦ)), Łf::$(typealias(ŁF)))"
typealias_def(::Type{<:LenseFlow{I,t₀,t₁}}) where {I,t₀,t₁} =
    "LenseFlow{$t₀→$t₁, $(typealias(I))}(ϕ)"
size(L::CachedLenseFlow) = length(L.memŁf) .* (1,1)

# convenience for getting the actual ϕ map
getϕ(L::LenseFlow) = L.ϕ
getϕ(L::CachedLenseFlow) = L.ϕ[]

# if the type and ϕ are the same, its the same op
hash(L::LenseFlowOp, h::UInt64) = foldr(hash, (typeof(L), getϕ(L)), init=h)

### caching
τ(t) = Float16(t)
cache(cL::CachedLenseFlow, f) = cL
(cL::CachedLenseFlow)(ϕ::Field) = cache!(cL,ϕ)
function cache(L::LenseFlow, f)
    f′ = Ł(L.ϕ) .* Ł(f) # in case ϕ is batched but f is not, promote f to batched
    cache!(alloc_cache(L,f′), L, f′)
end
function cache!(cL::CachedLenseFlow{N,t₀,t₁}, ϕ) where {N,t₀,t₁}
    if cL.ϕ[] === ϕ
        cL
    else
        cache!(cL,LenseFlow{RK4Solver{N},t₀,t₁}(ϕ),cL.memŁf)
    end
end
function cache!(cL::CachedLenseFlow{N,t₀,t₁}, L::LenseFlow{RK4Solver{N},t₀,t₁}, f) where {N,t₀,t₁}
    ts = range(t₀,t₁,length=2N+1)
    ∇ϕ,∇∇ϕ = map(Ł, gradhess(L.ϕ))
    T = eltype(L.ϕ)
    for (t,τ) in zip(ts,τ.(ts))
        @! cL.M⁻¹[τ] = pinv(Diagonal.(I + T(t)*∇∇ϕ))
        @! cL.p[τ] = cL.M⁻¹[τ]' * Diagonal.(∇ϕ)
    end
    cL.ϕ[] = L.ϕ
    cL
end
function alloc_cache(L::LenseFlow{RK4Solver{N},t₀,t₁}, f) where {N,t₀,t₁}
    ts = range(t₀,t₁,length=2N+1)
    p, M⁻¹ = Dict(), Dict()
    Łf,Ðf = Ł(f),  Ð(f)
    Łϕ,Ðϕ = Ł(L.ϕ),Ð(L.ϕ)
    for (t,τ) in zip(ts,τ.(ts))
        M⁻¹[τ] = Diagonal.(similar.(@SMatrix[Łϕ Łϕ; Łϕ Łϕ]))
        p[τ]   = Diagonal.(similar.(@SVector[Łϕ,Łϕ]))
    end
    CachedLenseFlow{N,t₀,t₁,typeof(Łϕ),typeof(Ðϕ),typeof(Łf),typeof(Ðf),eltype(Łϕ)}(
        Ref{Any}(L.ϕ), p, M⁻¹,
        similar(Łf), similar(Ðf), similar.(@SVector[Łf,Łf]), similar.(@SVector[Ðf,Ðf]),
        similar(Łϕ), similar(Ðϕ), similar.(@SVector[Łϕ,Łϕ]), similar.(@SVector[Ðϕ,Ðϕ]),
    )
end

# the way these velocities work is that they unpack the preallocated fields
# stored in L.mem* into variables with more meaningful names, which are then
# used in a bunch of in-place (eg mul!, Ł!, etc...) functions. note the use of
# the @! macro, which just rewrites @! x = f(y) to x = f!(x,y) for easier
# reading.

function velocity(L::LenseFlowOp{<:RK4Solver}, f₀::Field)
    function v!(v::Field, t::Real, f::Field)
        Ðf, Ð∇f, Ł∇f = L.memÐf, L.memÐvf,  L.memŁvf
        p = L.p[τ(t)]

        @! Ðf  = Ð(f)
        @! Ð∇f = ∇ᵢ * Ðf
        @! Ł∇f = Ł(Ð∇f)
        @! v   = p' * Ł∇f
    end
    return (v!, L.memŁf .= Ł(f₀))
end

function velocityᴴ(L::LenseFlowOp{<:RK4Solver}, f₀::Field)
    function v!(v::Field, t::Real, f::Field)
        Łf, Łf_p, Ð_Łf_p = L.memŁf, L.memŁvf, L.memÐvf
        p = L.p[τ(t)]

        @! Łf = Ł(f)
        @! Łf_p = p * Łf
        @! Ð_Łf_p = Ð(Łf_p)
        @! v = -∇ᵢ' * Ð_Łf_p
    end
    return (v!, L.memÐf .= Ð(f₀))
end

function negδvelocityᴴ(L::LenseFlowOp{<:RK4Solver}, (f₀, δf₀)::FieldTuple)

    function v!((df_dt, dδf_dt, dδϕ_dt)::FieldTuple, t::Real, (f, δf, δϕ)::FieldTuple)

        p   = L.p[τ(t)]
        M⁻¹ = L.M⁻¹[τ(t)]

        # dδf/dt
        Łδf, Łδf_p, Ð_Łδf_p = L.memŁf, L.memŁvf, L.memÐvf
        @! Łδf     = Ł(δf)
        @! Łδf_p   = p * Łδf
        @! Ð_Łδf_p = Ð(Łδf_p)
        @! dδf_dt  = -∇ᵢ' * Ð_Łδf_p

        # df/dt
        Ðf, Ð∇f, Ł∇f = L.memÐf, L.memÐvf,  L.memŁvf
        @! Ðf     = Ð(f)
        @! Ð∇f    = ∇ᵢ * Ðf
        @! Ł∇f    = Ł(Ð∇f)
        @! df_dt  = p' * Ł∇f

        # dδϕ/dt
        δfᵀ_∇f, M⁻¹_δfᵀ_∇f, Ð_M⁻¹_δfᵀ_∇f = L.memŁvϕ, L.memŁvϕ, L.memÐvϕ
        @! δfᵀ_∇f       = spin_adjoint(Łδf) * Ł∇f
        @! M⁻¹_δfᵀ_∇f   = M⁻¹ * δfᵀ_∇f
        @! Ð_M⁻¹_δfᵀ_∇f = Ð(M⁻¹_δfᵀ_∇f)
        @! dδϕ_dt       = -∇ⁱ' * Ð_M⁻¹_δfᵀ_∇f
        memÐϕ = L.memÐϕ
        for i=1:2, j=1:2
            dδϕ_dt .+= (@! memÐϕ = ∇ⁱ[i]' * (@! memÐϕ = ∇ᵢ[j]' * (@! memÐϕ = Ð(@. L.memŁϕ = t * p[j].diag * M⁻¹_δfᵀ_∇f[i]))))
        end

        FieldTuple(df_dt, dδf_dt, dδϕ_dt)

    end

    return (v!, FieldTuple(Ł(f₀), Ð(δf₀), L.memÐϕ .= Ð(zero(getϕ(L)))))

end

# adapting storage
adapt_structure(storage, Lϕ::LenseFlow{I,t₀,t₁}) where {I<:ODESolver,t₀,t₁} = LenseFlow{I,t₀,t₁}(adapt(storage,Lϕ.ϕ))
function adapt_structure(storage, Lϕ::CachedLenseFlow{N,t₀,t₁}) where {N,t₀,t₁}
    _adapt(x) = adapt(storage, x)
    memŁf, memÐf, memŁϕ, memÐϕ = _adapt(Lϕ.memŁf), _adapt(Lϕ.memÐf), _adapt(Lϕ.memŁϕ), _adapt(Lϕ.memÐϕ)
    CachedLenseFlow{N,t₀,t₁,typeof(memŁϕ),typeof(memÐϕ),typeof(memŁf),typeof(memÐf),eltype(memŁϕ)}(
        Ref(_adapt(Lϕ.ϕ[])),
        Dict(t => _adapt.(x) for (t,x) in Lϕ.p),
        Dict(t => _adapt.(x) for (t,x) in Lϕ.M⁻¹),
        memŁf, memÐf, _adapt.(Lϕ.memŁvf), _adapt.(Lϕ.memÐvf),
        memŁϕ, memÐϕ, _adapt.(Lϕ.memŁvϕ), _adapt.(Lϕ.memÐvϕ)
    )
end

"""
Returns αmax such that 𝕀 + ∇∇(ϕ + α * η) has non-zero discriminant
(pixel-by-pixel) for all α values in [0, αmax].

This mean ϕ + αmax * η is the maximum step in the η direction which
can be added to ϕ and still yield a lensing potential in the
weak-lensing regime. This is important because it guarantees the
potential can be paseed to LenseFlow, which cannot handle the
strong-lensing / "shell-crossing" regime.
"""
function get_max_lensing_step(ϕ, η)

    ϕ₁₁, ϕ₁₂, ϕ₂₁, ϕ₂₂ = Map.(gradhess(ϕ)[2])
    η₁₁, η₁₂, η₂₁, η₂₂ = Map.(gradhess(η)[2])

    a = @. η₁₁*η₂₂ - η₁₂^2
    b = @. η₁₁*(1+ϕ₂₂) + η₂₂*(1+ϕ₁₁) - 2*η₁₂*ϕ₁₂
    c = @. (1+ϕ₁₁) * (1+ϕ₂₂) - ϕ₁₂^2

    α₁ = cpu(@. (-b + sqrt(b^2 - 4*a*c))/(2*a))
    α₂ = cpu(@. (-b - sqrt(b^2 - 4*a*c))/(2*a))

    αmax = min(minimum(α₁[α₁.>0]), minimum(α₂[α₂.>0]))

end


text:
julia>
julia> q = convert(Ptr{Float64}, 0x0000000113f70d68)
Ptr{Float64} @0x0000000113f70d68

julia> B = pointer_to_array(q, (2,))
2-element Array{Float64,1}:
 2.3
 3.14149

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Plug # module Bukdu

using ..Deps
using ..Bukdu: Assoc, AbstractPlug, ApplicationController, AbstractRender, Render, bukdu_env

function plug
end

include("plugs/conn.jl")

using Logging: AbstractLogger
include("plugs/Loggers.jl")

include("plugs/ContentParsers.jl")
include("plugs/parsers.jl")

include("plugs/static.jl")

include("plugs/head.jl")

end # module Bukdu.Plug

export Plug, Conn, ApplicationController, Render, plug

using .Plug: Conn, ApplicationController, AbstractPlug, AbstractRender, Render
import .Plug: plug

# module Bukdu


text:
"""
    couple_models(aero::QuasiSteady, stru::LiftingLineSection, flap::SimpleFlap,
        ctrl::LiftingLineSectionControl)

Create an aerostructural model using a quasi-steady aerodynamics model, a
lifting line aerodynamic model, and a linear, steady-state control surface model.
The existence of this coupling allows the [`QuasiSteady`](@ref) and [`SimpleFlap`](@ref)
to be used with [`LiftingLine`](@ref) and [`LiftingLineFlaps`](@ref). This model
introduces the freestream air density ``\\rho`` as an additional parameter.
"""
function couple_models(aero::QuasiSteady, stru::LiftingLineSection, flap::SimpleFlap,
    ctrl::LiftingLineSectionControl)

    return (aero, stru, flap, ctrl)
end

# --- Traits --- #

# steady
function number_of_additional_parameters(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return 1
end

function coupling_inplaceness(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return OutOfPlace()
end

function coupling_rate_jacobian_type(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Zeros()
end

function coupling_state_jacobian_type(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Nonlinear()
end

function coupling_parameter_jacobian_type(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Nonlinear()
end

function coupling_time_gradient_type(::Type{QuasiSteady{0}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Zeros()
end

# quasisteady without acceleration terms
function number_of_additional_parameters(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return 1
end

function coupling_inplaceness(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return OutOfPlace()
end

function coupling_rate_jacobian_type(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Zeros()
end

function coupling_state_jacobian_type(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Nonlinear()
end

function coupling_parameter_jacobian_type(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Nonlinear()
end

function coupling_time_gradient_type(::Type{QuasiSteady{1}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Zeros()
end

# quasisteady with acceleration terms
function number_of_additional_parameters(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})
    return 1
end

function coupling_inplaceness(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return OutOfPlace()
end

function coupling_rate_jacobian_type(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Linear()
end

function coupling_state_jacobian_type(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Nonlinear()
end

function coupling_parameter_jacobian_type(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Nonlinear()
end

function coupling_time_gradient_type(::Type{QuasiSteady{2}}, ::Type{LiftingLineSection},
    ::Type{SimpleFlap}, ::Type{LiftingLineSectionControl})

    return Zeros()
end

# --- Methods --- #

# steady
function get_coupling_inputs(aero::QuasiSteady{0}, stru::LiftingLineSection,
    flap::SimpleFlap, ctrl::LiftingLineSectionControl, dx, x, p, t)
    # extract rate variables
    dvx, dvy, dvz, dωx, dωy, dωz, dδ = dx
    # extract state variables
    vx, vy, vz, ωx, ωy, ωz, δ = x
    # extract parameters
    a, b, a0, α0, cd0, cm0, cnδ, caδ, cmδ, ρ = p
    # local freestream velocity components
    u, v, ω = liftingline_velocities(vx, vz, ωy)
    # calculate loads
    Na, Aa, Ma = quasisteady0_loads(a, b, ρ, a0, α0, cd0, cm0, u, v)
    # loads due to flap deflections
    Nf, Af, Mf = simpleflap_loads(b, u, ρ, cnδ, caδ, cmδ, δ)
    # loads per unit span
    f = SVector(Aa+Af, 0, Na+Nf)
    m = SVector(0, Ma+Mf, 0)
    # return inputs
    return SVector(f..., m...)
end

# quasisteady without acceleration terms
function get_coupling_inputs(aero::QuasiSteady{1}, stru::LiftingLineSection,
    flap::SimpleFlap, ctrl::LiftingLineSectionControl, dx, x, p, t)
    # extract rate variables
    dvx, dvy, dvz, dωx, dωy, dωz, dδ = dx
    # extract state variables
    vx, vy, vz, ωx, ωy, ωz, δ = x
    # extract parameters
    a, b, a0, α0, cd0, cm0, cnδ, caδ, cmδ, ρ = p
    # local freestream velocity components
    u, v, ω = liftingline_velocities(vx, vz, ωy)
    # calculate aerodynamic loads
    Na, Aa, Ma = quasisteady1_loads(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)
    # loads due to flap deflections
    Nf, Af, Mf = simpleflap_loads(b, u, ρ, cnδ, caδ, cmδ, δ)
    # loads per unit span
    f = SVector(Aa+Af, 0, Na+Nf)
    m = SVector(0, Ma+Mf, 0)
    # return inputs
    return SVector(f..., m...)
end

# quasisteady with acceleration terms
function get_coupling_inputs(aero::QuasiSteady{2}, stru::LiftingLineSection,
    flap::SimpleFlap, ctrl::LiftingLineSectionControl, dx, x, p, t)
    # extract rate variables
    dvx, dvy, dvz, dωx, dωy, dωz, dδ = dx
    # extract state variables
    vx, vy, vz, ωx, ωy, ωz, δ = x
    # extract parameters
    a, b, a0, α0, cd0, cm0, cnδ, caδ, cmδ, ρ = p
    # local freestream velocity components
    u, v, ω = liftingline_velocities(vx, vz, ωy)
    udot, vdot, ωdot = liftingline_accelerations(dvx, dvz, dωy)
    # calculate aerodynamic loads
    Na, Aa, Ma = quasisteady2_loads(a, b, ρ, a0, α0, cd0, cm0, u, v, ω, vdot, ωdot)
    # loads due to flap deflections
    Nf, Af, Mf = simpleflap_loads(b, u, ρ, cnδ, caδ, cmδ, δ)
    # loads per unit span
    f = SVector(Aa+Af, 0, Na+Nf)
    m = SVector(0, Ma+Mf, 0)
    # return inputs
    return SVector(f..., m...)
end

# --- Performance Overloads --- #

# TODO: State rate, state, and parameter jacobians

# --- Convenience Methods --- #

function set_additional_parameters!(padd, aero::QuasiSteady, stru::LiftingLineSection,
    flap::SimpleFlap, ctrl::LiftingLineSectionControl; rho)

    padd[1] = rho

    return padd
end

function separate_additional_parameters(aero::QuasiSteady, stru::LiftingLineSection,
    flap::SimpleFlap, ctrl::LiftingLineSectionControl, padd)

    return (rho = padd[1],)
end


text:
mutable struct SplashScene <: Ranger.AbstractScene
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}
    transitioning::Nodes.TransitionProperties

    replacement::Ranger.AbstractScene

    mesh::Geometry.Mesh

    function SplashScene(world::Ranger.World, name::String, replacement::Ranger.AbstractScene)
        obj = new()

        # We use "obj" to represent a lack of parent.
        obj.base = Nodes.NodeData(name, Nodes.NodeNil(), world)
        obj.replacement = replacement   # default to self/obj = No replacement present
        obj.transitioning = Nodes.TransitionProperties()
        obj.transform = Nodes.TransformProperties{Float64}()
        obj.mesh = Geometry.Mesh()

        obj
    end
end

function build(node::SplashScene, world::Ranger.World)
    # The splash scene doesn't not have a transform (aka translate + scale)
    # applied to it so we build the mesh using view-space coordinates
    y = -Float64(world.view_height) / 2.0
    w = 320

    # Construct a checker board grid of rectangles.
    while y <= world.view_height
        x = -Float64(world.view_width) / 2.0
        while x <= world.view_width
            Geometry.add_vertex!(node.mesh, x, y)  # top-left
            Geometry.add_vertex!(node.mesh, x + w, y + w)   # bottom-right
            x += w
        end
        y += w
    end

    Geometry.build!(node.mesh)
end

# --------------------------------------------------------
# Timing
# --------------------------------------------------------
function Nodes.update(node::SplashScene, dt::Float64)
    update(node.transitioning, dt);
end

# --------------------------------------------------------
# Visits
# --------------------------------------------------------
function Nodes.draw(node::SplashScene, context::Rendering.RenderContext)
    if Nodes.is_dirty(node)
        # Transform this node's vertices using the context
        Rendering.transform!(context, node.mesh)
        Nodes.set_dirty!(node, false)
    end

    # Draw background
    Rendering.render_checkerboard(context, node.mesh, RangerGame.darkgray, RangerGame.lightgray)

    Rendering.set_draw_color(context, RangerGame.white)
    Rendering.draw_text(context, 10, 10, node.base.name, 2, 2, false);
end

# --------------------------------------------------------
# Life cycle events
# --------------------------------------------------------
function Nodes.enter_node(node::SplashScene, man::Nodes.NodeManager)
    println("enter ", node);
    # Register node as a timing target in order to receive updates
    Nodes.register_target(man, node);
end

function Nodes.exit_node(node::SplashScene, man::Nodes.NodeManager)
    println("exit ", node);
    Nodes.unregister_target(man, node);
end

function Nodes.transition(node::SplashScene)
    if Nodes.ready(node.transitioning)
        Scenes.REPLACE_TAKE
    else
        Scenes.NO_ACTION
    end
end

function Nodes.get_replacement(node::SplashScene)
    node.replacement
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
struct RECA_TwoDim <: AbstractReca
    res_size::Int
    in_size::Int
    out_size::Int
    generations::Int
    permutations::Int
    train_data::AbstractArray{Float64}
    nla_type::ReservoirComputing.NonLinearAlgorithm
    states::AbstractArray{Float64}
    maps::AbstractArray{Int}
end

"""
    RECA_TwoDim(train_data, res_size, generations, permutations [, nla_type])

Create the RECA_TwoDim struct for two-dimensional Cellular Automata Reservoir Computing, as described in [1].

[1] Yilmaz, Ozgur. “Reservoir computing using cellular automata.” arXiv preprint arXiv:1410.0162 (2014).
"""
function RECA_TwoDim(train_data, res_size, generations, permutations; nla_type = NLADefault())

    in_size = size(train_data, 1)
    out_size = size(train_data, 1)

    initial_state = zeros(Bool, res_size, res_size)
    maps = single_mapping(in_size, res_size, permutations)
    states = harvest_states_standard(train_data, initial_state, generations, permutations, maps)
    train_data = convert(AbstractArray{Float64}, train_data)

    return RECA_TwoDim(res_size, in_size, out_size, generations, permutations, train_data, nla_type, states, maps)
end

"""
    RECATD_predict_discrete(reca, predict_len::Int, W_out::AbstractArray{Float64})

Return the prediction for a given length of the constructed RECA_TwoDim struct.
"""
function RECATD_predict_discrete(reca,
    predict_len::Int,
    W_out::AbstractArray{Float64})

    output = Array{Int}(undef, reca.in_size, predict_len)
    x = reca.states[:, end]
    last_states = reshape(x, reca.res_size, reca.res_size, reca.generations)
    last_state = convert(AbstractArray{Bool}, last_states[:,:, end])

    for i=1:predict_len
        out = convert(AbstractArray{Int}, W_out*x .> 0.5)
        output[:,i] = out#convert(AbstractArray{Int}, out .> 0.5)
        last_state = single_encoding(out, last_state, reca.maps)
        gol = GameOfLife(last_state, reca.generations)
        x = reshape(gol.all_runs, reca.res_size*reca.res_size*reca.generations)
        last_state = gol.all_runs[:, :, end]
    end
    return output
end

"""
    RECATDdirect_predict_discrete(reca::AbstractReca, W_out::AbstractArray{Float64}, test_data::AbstractArray{Int})

Given the input data return the corresponding predicted output, as described in [1].

[1] Yilmaz, Ozgur. “Reservoir computing using cellular automata.” arXiv preprint arXiv:1410.0162 (2014).
"""
function RECATDdirect_predict_discrete(reca::AbstractReca,
    W_out::AbstractArray{Float64},
    test_data::AbstractArray{Int})

    predict_len = size(test_data, 2)
    output = Array{Int}(undef, size(W_out, 1), predict_len)
    last_state = zeros(Bool, reca.res_size, reca.res_size)

    for i=1:predict_len

        last_state = single_encoding(test_data[:,i], last_state, reca.maps)
        gol = GameOfLife(last_state, reca.generations)
        x = reshape(gol.all_runs, reca.res_size*reca.res_size*reca.generations)
        out = W_out*x
        output[:,i] = convert(AbstractArray{Int}, out .> 0.5)
         last_state = gol.all_runs[:, :, end]

    end
    return output
end

function harvest_states_standard(input_data::AbstractArray{Int},
    initial_state::AbstractArray{T},
    generations::Int,
    permutations::Int,
    single_map::AbstractArray{Int}) where T<: Bool

    input_size = size(input_data, 1)
    train_time = size(input_data, 2)
    res_size = size(initial_state, 1)

    states = zeros(T, res_size*res_size*generations , train_time)
    #single_map = single_mapping(input_size, res_size, permutations)#outside

    for i=1:train_time
        initial_state = single_encoding(input_data[:,i], initial_state, single_map)
        gol = GameOfLife(initial_state, generations)#+1
        gol_states = copy(gol.all_runs)#2:end
        states[:, i] = reshape(gol_states, res_size*res_size*generations)
        initial_state = gol.all_runs[:, :, end]
    end
    return states
end

function single_encoding(input_vector::AbstractArray{Int},
    reservoir::AbstractArray{T}, #must be size(reservoir, 1) = size(reservoir, 2)
    single_map::AbstractArray{Int}) where T<:Bool

    map_size = size(single_map, 2)
    input_size = size(input_vector, 1)
    res_size = size(reservoir, 1)
    permutations = Int(map_size/input_size)
    new_reservoir = copy(reservoir)

    counter = 0
    for i=1:permutations
        for j=1:input_size
            counter += 1
            first_coordinate = mod1(single_map[counter], res_size)
            second_coordinate = cld(single_map[counter], res_size)

            new_reservoir[first_coordinate, second_coordinate] = input_vector[j]
        end
    end

    return new_reservoir
end

#binary input onto reservoir first approach: random initial conditions, single mapping
function single_mapping(input_size::T, res_size::T, permutations::T) where T<: Int

    dims = input_size*permutations
    first_coordinate = sample(1:res_size*res_size, input_size*permutations, replace=false)

    single_map = collect(transpose(first_coordinate))

    return single_map
end


text:
# ### m13.6nct.jl

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

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

# m13.6nc1 is equivalent to m13.6nc in following Turing model

@model ppl13_6nc(actor, block_id, condition, prosoc_left, pulled_left)
    # fixed priors
    Rho_block ~ LKJ(3, 4.)
    Rho_actor ~ LKJ(3, 4.)
    sigma_block ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 3)
    sigma_actor ~ filldist(truncated(Cauchy(0, 2), 0, Inf), 3)
    a ~ Normal(0, 1)
    bp ~ Normal(0, 1)
    bpc ~ Normal(0, 1)

    # adaptive NON-CENTERED priors
    Rho_block = (Rho_block' + Rho_block) / 2
    Rho_actor = (Rho_actor' + Rho_actor) / 2

    L_Rho_block = cholesky(Rho_block).L
    L_Rho_actor = cholesky(Rho_actor).L

    z_N_block ~ filldist(Normal(0, 1), 3, 6)
    z_N_actor ~ filldist(Normal(0, 1), 3, 7)

    # @show size(L_Rho_block) size(sigma_block) size(z_N_block_id)

    a_block_bp_block_bpc_block = sigma_block .* L_Rho_block * z_N_block
    a_actor_bp_actor_bpc_actor = sigma_actor .* L_Rho_actor * z_N_actor

    a_block = a_block_bp_block_bpc_block[1, :]
    bp_block = a_block_bp_block_bpc_block[2, :]
    bpc_block = a_block_bp_block_bpc_block[3, :]
    a_actor = a_actor_bp_actor_bpc_actor[1, :]
    bp_actor = a_actor_bp_actor_bpc_actor[2, :]
    bpc_actor = a_actor_bp_actor_bpc_actor[3, :]

    # linear models
    BPC = bpc .+ bpc_actor[actor] + bpc_block[block_id]
    BP = bp .+ bp_actor[actor] + bp_block[block_id]
    A = a .+ a_actor[actor] + a_block[block_id]
    logit_p = A + (BP + BPC .* condition) .* prosoc_left

    # likelihood
    pulled_left .~ BinomialLogit.(1, logit_p)
end

m13_6nct = ppl13_6nc(d.actor, d.block_id, d.condition,
    d.prosoc_left, d.pulled_left)
nchains = 4; sampler = NUTS(0.65); nsamples=2000
nchains = 1 # Fails with nchains=4
chns13_6nct = mapreduce(c -> sample(m13_6nct, sampler, nsamples), chainscat, 1:nchains)

m_13_6_ncs_results = """
Inference for Stan model: 78757717ec3ffa528c3819708907eb7f.
    3 chains, each with iter=5000; warmup=1000; thin=1;
    post-warmup draws per chain=4000, total post-warmup draws=12000.

                         mean se_mean   sd    2.5%     25%     50%     75%   97.5%
    z_N_block_id[1,1]   -0.44    0.01 0.90   -2.17   -1.03   -0.46    0.13    1.39
    z_N_block_id[1,2]    0.21    0.01 0.86   -1.50   -0.35    0.22    0.78    1.89
    z_N_block_id[1,3]    0.34    0.01 0.87   -1.42   -0.22    0.35    0.92    2.02
    z_N_block_id[1,4]   -0.02    0.01 0.84   -1.71   -0.56   -0.03    0.52    1.63
    z_N_block_id[1,5]   -0.09    0.01 0.84   -1.73   -0.64   -0.10    0.47    1.56
    z_N_block_id[1,6]    0.01    0.01 0.92   -1.85   -0.59    0.03    0.63    1.80
    z_N_block_id[2,1]   -0.37    0.01 0.82   -1.99   -0.90   -0.37    0.15    1.27
    z_N_block_id[2,2]   -0.18    0.01 0.78   -1.76   -0.66   -0.17    0.31    1.38
    z_N_block_id[2,3]   -0.36    0.01 0.79   -1.93   -0.87   -0.35    0.16    1.18
    z_N_block_id[2,4]    0.23    0.01 0.77   -1.33   -0.26    0.23    0.72    1.74
    z_N_block_id[2,5]    0.04    0.01 0.76   -1.49   -0.46    0.04    0.52    1.58
    z_N_block_id[2,6]    1.01    0.01 0.84   -0.76    0.49    1.04    1.56    2.58
    z_N_block_id[3,1]   -0.59    0.01 0.91   -2.32   -1.20   -0.61    0.00    1.28
    z_N_block_id[3,2]    0.28    0.01 0.88   -1.50   -0.29    0.28    0.86    2.04
    z_N_block_id[3,3]    0.32    0.01 0.89   -1.46   -0.26    0.34    0.91    2.09
    z_N_block_id[3,4]   -0.11    0.01 0.89   -1.87   -0.69   -0.10    0.46    1.64
    z_N_block_id[3,5]   -0.12    0.01 0.86   -1.83   -0.68   -0.12    0.43    1.64
    z_N_block_id[3,6]    0.20    0.01 0.91   -1.63   -0.40    0.20    0.81    1.99
    L_Rho_block[1,1]     1.00     NaN 0.00    1.00    1.00    1.00    1.00    1.00
    L_Rho_block[1,2]     0.00     NaN 0.00    0.00    0.00    0.00    0.00    0.00
    L_Rho_block[1,3]     0.00     NaN 0.00    0.00    0.00    0.00    0.00    0.00
    L_Rho_block[2,1]    -0.06    0.00 0.32   -0.65   -0.29   -0.06    0.17    0.57
    L_Rho_block[2,2]     0.94    0.00 0.07    0.73    0.92    0.97    0.99    1.00
    L_Rho_block[2,3]     0.00     NaN 0.00    0.00    0.00    0.00    0.00    0.00
    L_Rho_block[3,1]     0.04    0.00 0.32   -0.57   -0.19    0.04    0.27    0.63
    L_Rho_block[3,2]    -0.04    0.00 0.31   -0.62   -0.26   -0.04    0.18    0.57
    L_Rho_block[3,3]     0.89    0.00 0.10    0.64    0.84    0.92    0.97    1.00
    z_N_actor[1,1]      -0.49    0.01 0.36   -1.25   -0.71   -0.46   -0.24    0.15
    z_N_actor[1,2]       1.99    0.01 0.62    0.88    1.55    1.96    2.40    3.29
    z_N_actor[1,3]      -0.63    0.01 0.38   -1.44   -0.87   -0.60   -0.35    0.04
    z_N_actor[1,4]      -0.65    0.01 0.39   -1.48   -0.89   -0.61   -0.37    0.03
    z_N_actor[1,5]      -0.48    0.01 0.36   -1.25   -0.71   -0.45   -0.23    0.15
    z_N_actor[1,6]       0.09    0.00 0.33   -0.56   -0.13    0.09    0.30    0.75
    z_N_actor[1,7]       0.73    0.01 0.41    0.04    0.45    0.69    0.98    1.62
    z_N_actor[2,1]       0.06    0.01 0.83   -1.60   -0.49    0.07    0.59    1.71
    z_N_actor[2,2]       0.14    0.01 0.99   -1.80   -0.52    0.14    0.81    2.07
    z_N_actor[2,3]       0.38    0.01 0.86   -1.34   -0.18    0.38    0.95    2.04
    z_N_actor[2,4]       0.25    0.01 0.85   -1.46   -0.29    0.26    0.81    1.91
    z_N_actor[2,5]       0.13    0.01 0.84   -1.57   -0.41    0.13    0.68    1.79
    z_N_actor[2,6]      -0.62    0.01 0.87   -2.26   -1.20   -0.65   -0.09    1.20
    z_N_actor[2,7]      -0.04    0.01 0.92   -1.84   -0.66   -0.04    0.58    1.78
    z_N_actor[3,1]       0.10    0.01 0.86   -1.63   -0.46    0.10    0.67    1.81
    z_N_actor[3,2]       0.07    0.01 0.98   -1.84   -0.60    0.08    0.74    1.98
    z_N_actor[3,3]      -0.35    0.01 0.86   -2.02   -0.94   -0.37    0.22    1.38
    z_N_actor[3,4]       0.02    0.01 0.86   -1.70   -0.55    0.01    0.58    1.73
    z_N_actor[3,5]      -0.05    0.01 0.85   -1.73   -0.62   -0.06    0.50    1.65
    z_N_actor[3,6]      -0.25    0.01 0.89   -1.97   -0.83   -0.26    0.31    1.56
    z_N_actor[3,7]       0.50    0.01 0.97   -1.43   -0.13    0.50    1.16    2.38
    L_Rho_actor[1,1]     1.00     NaN 0.00    1.00    1.00    1.00    1.00    1.00
    L_Rho_actor[1,2]     0.00     NaN 0.00    0.00    0.00    0.00    0.00    0.00
    L_Rho_actor[1,3]     0.00     NaN 0.00    0.00    0.00    0.00    0.00    0.00
    L_Rho_actor[2,1]    -0.07    0.00 0.32   -0.65   -0.30   -0.07    0.15    0.56
    L_Rho_actor[2,2]     0.94    0.00 0.07    0.73    0.92    0.97    0.99    1.00
    L_Rho_actor[2,3]     0.00     NaN 0.00    0.00    0.00    0.00    0.00    0.00
    L_Rho_actor[3,1]     0.08    0.00 0.32   -0.55   -0.15    0.09    0.31    0.66
    L_Rho_actor[3,2]    -0.04    0.00 0.31   -0.64   -0.27   -0.05    0.18    0.57
    L_Rho_actor[3,3]     0.88    0.00 0.10    0.61    0.83    0.91    0.96    1.00
    a                    0.24    0.01 0.66   -1.06   -0.19    0.24    0.67    1.56
    bp                   0.72    0.00 0.40   -0.14    0.47    0.73    0.98    1.48
    bpc                 -0.03    0.00 0.44   -0.88   -0.31   -0.03    0.25    0.88
    sigma_actor[1]       2.35    0.01 0.90    1.20    1.73    2.16    2.75    4.67
    sigma_actor[2]       0.46    0.00 0.37    0.02    0.19    0.39    0.64    1.37
    sigma_actor[3]       0.53    0.01 0.49    0.02    0.19    0.40    0.72    1.81
    sigma_block[1]       0.23    0.00 0.20    0.01    0.09    0.18    0.31    0.74
    sigma_block[2]       0.57    0.01 0.41    0.03    0.29    0.50    0.77    1.56
    sigma_block[3]       0.52    0.01 0.43    0.02    0.20    0.42    0.72    1.59
    a_block[1]          -0.12    0.00 0.23   -0.68   -0.22   -0.06    0.01    0.23
    a_block[2]           0.06    0.00 0.21   -0.32   -0.04    0.02    0.14    0.57
    a_block[3]           0.10    0.00 0.22   -0.27   -0.02    0.04    0.19    0.65
    a_block[4]          -0.01    0.00 0.20   -0.44   -0.09    0.00    0.08    0.42
    a_block[5]          -0.02    0.00 0.20   -0.49   -0.11   -0.01    0.06    0.39
    a_block[6]           0.00    0.00 0.22   -0.49   -0.09    0.00    0.10    0.46
    bp_block[1]         -0.18    0.00 0.42   -1.12   -0.42   -0.14    0.05    0.64
    bp_block[2]         -0.11    0.00 0.41   -0.99   -0.32   -0.07    0.11    0.70
    bp_block[3]         -0.23    0.00 0.43   -1.24   -0.46   -0.16    0.03    0.54
    bp_block[4]          0.15    0.00 0.41   -0.59   -0.08    0.09    0.37    1.11
    bp_block[5]          0.04    0.00 0.40   -0.75   -0.18    0.01    0.24    0.92
    bp_block[6]          0.63    0.01 0.57   -0.14    0.17    0.53    0.96    1.98
    bpc_block[1]        -0.37    0.01 0.54   -1.75   -0.63   -0.22   -0.01    0.38
    bpc_block[2]         0.18    0.00 0.44   -0.57   -0.05    0.09    0.39    1.28
    bpc_block[3]         0.22    0.00 0.47   -0.56   -0.04    0.11    0.43    1.40
    bpc_block[4]        -0.08    0.00 0.45   -1.11   -0.28   -0.03    0.12    0.82
    bpc_block[5]        -0.07    0.00 0.42   -1.02   -0.26   -0.03    0.12    0.78
    bpc_block[6]         0.09    0.00 0.48   -0.88   -0.13    0.04    0.30    1.19
    Rho_block[1,1]       1.00     NaN 0.00    1.00    1.00    1.00    1.00    1.00
    Rho_block[1,2]      -0.06    0.00 0.32   -0.65   -0.29   -0.06    0.17    0.57
    Rho_block[1,3]       0.04    0.00 0.32   -0.57   -0.19    0.04    0.27    0.63
    Rho_block[2,1]      -0.06    0.00 0.32   -0.65   -0.29   -0.06    0.17    0.57
    Rho_block[2,2]       1.00    0.00 0.00    1.00    1.00    1.00    1.00    1.00
    Rho_block[2,3]      -0.04    0.00 0.31   -0.62   -0.27   -0.05    0.18    0.57
    Rho_block[3,1]       0.04    0.00 0.32   -0.57   -0.19    0.04    0.27    0.63
    Rho_block[3,2]      -0.04    0.00 0.31   -0.62   -0.27   -0.05    0.18    0.57
    Rho_block[3,3]       1.00    0.00 0.00    1.00    1.00    1.00    1.00    1.00
    a_actor[1]          -1.02    0.01 0.71   -2.43   -1.48   -1.01   -0.55    0.38
    a_actor[2]           4.44    0.02 1.63    2.09    3.36    4.17    5.21    8.34
    a_actor[3]          -1.31    0.01 0.72   -2.75   -1.77   -1.30   -0.83    0.11
    a_actor[4]          -1.35    0.01 0.72   -2.80   -1.82   -1.34   -0.88    0.10
    a_actor[5]          -1.00    0.01 0.71   -2.42   -1.48   -0.99   -0.53    0.40
    a_actor[6]           0.20    0.01 0.72   -1.24   -0.28    0.20    0.66    1.60
    a_actor[7]           1.55    0.01 0.76    0.10    1.05    1.53    2.04    3.08
    bp_actor[1]          0.05    0.00 0.37   -0.68   -0.13    0.02    0.22    0.89
    bp_actor[2]          0.03    0.01 0.64   -1.21   -0.24    0.00    0.25    1.46
    bp_actor[3]          0.22    0.00 0.42   -0.44   -0.02    0.12    0.41    1.25
    bp_actor[4]          0.16    0.00 0.39   -0.52   -0.05    0.08    0.34    1.09
    bp_actor[5]          0.10    0.00 0.38   -0.61   -0.09    0.04    0.26    0.99
    bp_actor[6]         -0.31    0.00 0.42   -1.35   -0.54   -0.21   -0.01    0.28
    bp_actor[7]         -0.04    0.00 0.43   -1.00   -0.23   -0.02    0.15    0.87
    bpc_actor[1]         0.02    0.00 0.42   -0.92   -0.15    0.01    0.20    0.94
    bpc_actor[2]         0.20    0.01 0.83   -1.18   -0.13    0.05    0.41    2.31
    bpc_actor[3]        -0.28    0.01 0.50   -1.58   -0.49   -0.14    0.02    0.46
    bpc_actor[4]        -0.05    0.00 0.43   -1.05   -0.22   -0.01    0.15    0.79
    bpc_actor[5]        -0.08    0.00 0.43   -1.13   -0.25   -0.02    0.12    0.73
    bpc_actor[6]        -0.14    0.00 0.45   -1.25   -0.32   -0.05    0.08    0.69
    bpc_actor[7]         0.41    0.01 0.73   -0.41   -0.01    0.18    0.64    2.39
    Rho_actor[1,1]       1.00     NaN 0.00    1.00    1.00    1.00    1.00    1.00
    Rho_actor[1,2]      -0.07    0.00 0.32   -0.65   -0.30   -0.07    0.15    0.56
    Rho_actor[1,3]       0.08    0.00 0.32   -0.55   -0.15    0.09    0.31    0.66
    Rho_actor[2,1]      -0.07    0.00 0.32   -0.65   -0.30   -0.07    0.15    0.56
    Rho_actor[2,2]       1.00    0.00 0.00    1.00    1.00    1.00    1.00    1.00
    Rho_actor[2,3]      -0.05    0.00 0.32   -0.63   -0.28   -0.05    0.17    0.56
    Rho_actor[3,1]       0.08    0.00 0.32   -0.55   -0.15    0.09    0.31    0.66
    Rho_actor[3,2]      -0.05    0.00 0.32   -0.63   -0.28   -0.05    0.17    0.56
    Rho_actor[3,3]       1.00    0.00 0.00    1.00    1.00    1.00    1.00    1.00
    lp__              -285.47    0.11 6.46 -299.01 -289.56 -285.23 -281.02 -273.65
                      n_eff Rhat
    z_N_block_id[1,1] 13367    1
    z_N_block_id[1,2] 16316    1
    z_N_block_id[1,3] 13212    1
    z_N_block_id[1,4] 16245    1
    z_N_block_id[1,5] 14391    1
    z_N_block_id[1,6] 13216    1
    z_N_block_id[2,1] 11937    1
    z_N_block_id[2,2] 11867    1
    z_N_block_id[2,3] 12192    1
    z_N_block_id[2,4] 14644    1
    z_N_block_id[2,5] 13497    1
    z_N_block_id[2,6] 10467    1
    z_N_block_id[3,1] 15608    1
    z_N_block_id[3,2] 15577    1
    z_N_block_id[3,3] 17322    1
    z_N_block_id[3,4] 16464    1
    z_N_block_id[3,5] 17740    1
    z_N_block_id[3,6] 15742    1
    L_Rho_block[1,1]    NaN  NaN
    L_Rho_block[1,2]    NaN  NaN
    L_Rho_block[1,3]    NaN  NaN
    L_Rho_block[2,1]  13964    1
    L_Rho_block[2,2]   6881    1
    L_Rho_block[2,3]    NaN  NaN
    L_Rho_block[3,1]  18184    1
    L_Rho_block[3,2]  17024    1
    L_Rho_block[3,3]   5530    1
    z_N_actor[1,1]     3997    1
    z_N_actor[1,2]     8100    1
    z_N_actor[1,3]     4054    1
    z_N_actor[1,4]     4042    1
    z_N_actor[1,5]     4067    1
    z_N_actor[1,6]     4812    1
    z_N_actor[1,7]     5978    1
    z_N_actor[2,1]    17688    1
    z_N_actor[2,2]    19821    1
    z_N_actor[2,3]    15620    1
    z_N_actor[2,4]    15806    1
    z_N_actor[2,5]    16515    1
    z_N_actor[2,6]    13494    1
    z_N_actor[2,7]    19772    1
    z_N_actor[3,1]    16698    1
    z_N_actor[3,2]    25570    1
    z_N_actor[3,3]    15609    1
    z_N_actor[3,4]    17747    1
    z_N_actor[3,5]    16452    1
    z_N_actor[3,6]    16474    1
    z_N_actor[3,7]    17924    1
    L_Rho_actor[1,1]    NaN  NaN
    L_Rho_actor[1,2]    NaN  NaN
    L_Rho_actor[1,3]    NaN  NaN
    L_Rho_actor[2,1]  20181    1
    L_Rho_actor[2,2]   5099    1
    L_Rho_actor[2,3]    NaN  NaN
    L_Rho_actor[3,1]  16979    1
    L_Rho_actor[3,2]  16551    1
    L_Rho_actor[3,3]   5465    1
    a                  4604    1
    bp                 9915    1
    bpc               10895    1
    sigma_actor[1]     4844    1
    sigma_actor[2]     7378    1
    sigma_actor[3]     6330    1
    sigma_block[1]     7461    1
    sigma_block[2]     5026    1
    sigma_block[3]     6372    1
    a_block[1]        11376    1
    a_block[2]        14094    1
    a_block[3]        11484    1
    a_block[4]        14446    1
    a_block[5]        12892    1
    a_block[6]        12830    1
    bp_block[1]       10766    1
    bp_block[2]       11464    1
    bp_block[3]       10088    1
    bp_block[4]       10004    1
    bp_block[5]       10791    1
    bp_block[6]        6266    1
    bpc_block[1]       9908    1
    bpc_block[2]      11090    1
    bpc_block[3]      10070    1
    bpc_block[4]      12510    1
    bpc_block[5]      12340    1
    bpc_block[6]      11472    1
    Rho_block[1,1]      NaN  NaN
    Rho_block[1,2]    13964    1
    Rho_block[1,3]    18184    1
    Rho_block[2,1]    13964    1
    Rho_block[2,2]    11984    1
    Rho_block[2,3]    15006    1
    Rho_block[3,1]    18184    1
    Rho_block[3,2]    15006    1
    Rho_block[3,3]    11786    1
    a_actor[1]         4930    1
    a_actor[2]         7208    1
    a_actor[3]         4880    1
    a_actor[4]         5074    1
    a_actor[5]         4983    1
    a_actor[6]         4827    1
    a_actor[7]         5096    1
    bp_actor[1]       12464    1
    bp_actor[2]       14817    1
    bp_actor[3]       10567    1
    bp_actor[4]       12646    1
    bp_actor[5]       13023    1
    bp_actor[6]       10818    1
    bp_actor[7]       14796    1
    bpc_actor[1]      12189    1
    bpc_actor[2]      12373    1
    bpc_actor[3]       8575    1
    bpc_actor[4]      11896    1
    bpc_actor[5]      10873    1
    bpc_actor[6]      10218    1
    bpc_actor[7]       9055    1
    Rho_actor[1,1]      NaN  NaN
    Rho_actor[1,2]    20181    1
    Rho_actor[1,3]    16979    1
    Rho_actor[2,1]    20181    1
    Rho_actor[2,2]    10875    1
    Rho_actor[2,3]    14880    1
    Rho_actor[3,1]    16979    1
    Rho_actor[3,2]    14880    1
    Rho_actor[3,3]    10080    1
    lp__               3172    1
"""

# End of m13.6nct.jl


text:
module BenchGEMM

using BenchmarkTools

include("gemm.jl")

let
    A = randn(3, 3)
    B = randn(3, 3)
    desired = A * B
    C = zero(desired)

    @assert ndrmul!(C, A, B) ≈ desired
    @assert manmul!(C, A, B) ≈ desired
end

suite = BenchmarkGroup()

let sub = suite["mul"] = BenchmarkGroup()
    for impl in ["ndr", "man"]
        sub[impl] = BenchmarkGroup()
    end
end

# for n in [8, 32, 256]
for n in [32]
    A = randn(n, n)
    B = randn(n, n)
    desired = A * B
    C = zero(desired)

    sub = suite["mul"]
    for (impl, f) in [("ndr", ndrmul!), ("man", manmul!)]
        sub[impl][n] = @benchmarkable($f($C, $A, $B))
    end
end

end  # module
BenchGEMM.suite


text:
mutable struct Handle
    ptr::Ptr{Cvoid}

    function Handle()
        ref = Ref{Ptr{Cvoid}}()
        @cudnn :cudnnCreate (Ptr{Ptr{Cvoid}},) ref
        h = new(ref[])
        finalizer(h) do x
            @cudnn :cudnnDestroy (Ptr{Cvoid},) x
        end
        h
    end
end

Base.cconvert(::Type{Ptr{Cvoid}}, h::Handle) = h.ptr
Base.isequal(x::Handle, y::Handle) = isequal(x.ptr, y.ptr)

const HANDLES = Array{Handle}(undef, ndevices())

function gethandle()
    dev = getdevice()
    isassigned(HANDLES,dev+1) || (HANDLES[dev+1] = Handle())
    HANDLES[dev+1]
end