text:
# License is MIT: LICENSE.md

using ModuleInterfaceTools

@api extend StrRegex

using Random, Unicode, Test

const IndexError = isdefined(Base, :StringIndexError) ? StringIndexError : UnicodeError

cvt(v::Integer, m::Nothing) = v
cvt(v::Integer, m::Vector{Int}) = v < 1 ? v : m[v]
cvt(v::UnitRange, m) = cvt(v.start,m):cvt(v.stop,m)

# Should test GenericString also, once overthing else is working
const UnicodeStringTypes = (String, UTF8Str, UTF16Str, UTF32Str, UCS2Str)

const ASCIIStringTypes = (ASCIIStr, LatinStr, UTF8Str, UTF16Str, UTF32Str, UCS2Str, String)

function test2(str, list, umap::Any=nothing)
    mymap = (encoding(str) === encoding(UTF8Str) ? umap : nothing)
    for (pat, res) in list
        cvtres = cvt(res, mymap)
        (r = fnd(First, pat, str)) == cvtres ||
            println("fnd(First, $(typeof(pat)):\"$pat\", $(typeof(str)):\"$str\") => $r != $cvtres ($res)")
        @test fnd(First, pat, str) == cvtres
    end
end

function test3(str, list, umap::Any=nothing)
    mymap = encoding(str) == encoding(UTF8Str) ? umap : nothing
    for (pat, beg, res) in list
        cvtbeg = cvt(beg, mymap)
        cvtres = cvt(res, mymap)
        (r = fnd(Fwd, pat, str, cvtbeg)) == cvtres ||
            println("fnd(Fwd, $(typeof(pat)):\"$pat\", $(typeof(str)):\"$str\", $cvtbeg ($beg)) => $r != $cvtres ($res)")
        @test fnd(Fwd, pat, str, cvtbeg) == cvtres
    end
end

const fbbstr = "foo,bar,baz"
const astr = "Hello, world.\n"
const u8str = "∀ ε > 0, ∃ δ > 0: |x-y| < δ ⇒ |f(x)-f(y)| < ε"

# v0.6.2 returns a Vector{Any} for some reason, make sure it is converted to Vector{Int}
const utf8map = convert(Vector{Int}, collect(eachindex(u8str)))

@testset "ASCII Regex Tests" begin
    for T in ASCIIStringTypes
        fbb = T(fbbstr)
        str = T(astr)
        @testset "str: $T" begin
            # string forward search with a single-char regex
            let pats = (R"x", R"H", R"l", R"\n"),
                res  = (0:-1, 1:1, 3:3, 14:14)
                test2(str, zip(pats, res))
            end
            let pats = (R"H", R"l", R"l", R"l", R"\n"),
                pos  = (  2,    4,    5,   12,    14), # Was 15 for fndnext
                res  = (0:-1, 4:4,11:11, 0:-1, 14:14)
                test3(str, zip(pats, pos, res))
            end
            i = 1
            while i <= ncodeunits(str)
                @test fnd(Fwd, R"."s, str, i) == i:i
                # string forward search with a zero-char regex
                @test fnd(Fwd, R"", str, i) == i:i-1
                i = nextind(str, i)
            end
            let pats = (R"xx", R"fo", R"oo", R"o,", R",b", R"az"),
                res  = ( 0:-1,   1:2,   2:3,   3:4,   4:5, 10:11)
                test2(fbb, zip(pats, res))
            end
            let pats = (R"fo", R"oo", R"o,", R",b", R",b", R"az"),
                pos  = (    3,     4,     5,     6,    10,    11),
                res  = ( 0:-1,  0:-1,  0:-1,   8:9,  0:-1,  0:-1) # was 12 for fndnext
                test3(fbb, zip(pats, pos, res))
            end
        end
        rs = T("foo123 foo456 foo789")
        @testset "(starts|ends)_with: $T" begin
            @test startswith(rs, R"fo+\d+")
            @test endswith(rs, R"fo+\d+")
            @test !startswith(rs, R"f\d+")
            @test !endswith(rs, R"f\d+")
        end
    end
end

@testset "Unicode Regex Tests" begin
    for T in UnicodeStringTypes
        @testset "str: $T" begin
            str = T(u8str)
            @testset "Regex" begin
                let pats = (R"z", R"∄", R"∀", R"∃", R"x", R"ε"),
                    res  = (0:-1, 0:-1, 1:1, 10:10,20:20,  3:3)
                    test2(str, zip(pats, res), utf8map)
                end
                let pats = (R"∀", R"∃", R"x", R"x", R"ε", R"ε"),
                    pos  = (   2,   11,   21,   35,    4,   45), # was 56 for fndnext
                    res  = (0:-1, 0:-1,34:34, 0:-1,45:45,45:45)  # was 0:-1 for last
                    test3(str, zip(pats, pos, res), utf8map)
                end
                v = cvt(2, encoding(T) == encoding(UTF8Str) ? utf8map : nothing)
                @test fnd(First, R"∀", str)  == fnd(First, R"\u2200", str)
                @test fnd(Fwd, R"∀", str, v) == fnd(Fwd, R"\u2200", str, v)
                i = 1
                while i <= ncodeunits(str)
                    @test fnd(Fwd, R"."s, str, i) == i:i
                    # string forward search with a zero-char regex
                    @test fnd(Fwd, R"", str, i) == i:i-1
                    i = nextind(str, i)
                end
            end
        end
    end
end

const RegexStrings = (ASCIIStr, BinaryStr, Text1Str, LatinStr, _LatinStr, UTF8Str)

collect_eachmatch(re, str; overlap=false) =
    [m.match for m in collect(eachmatch(re, str, overlap = overlap))]

@testset "UTF8Str Regex" begin
    # Proper unicode handling
    @test match(R"∀∀", UTF8Str("∀x∀∀∀")).match == "∀∀"
    @test collect_eachmatch(R".\s", UTF8Str("x \u2200 x \u2203 y")) == ["x ", "∀ ", "x ", "∃ "]
end

@testset "Regex" begin
    for T in RegexStrings
        @test collect_eachmatch(R"a?b?", T("asbd")) == ["a","","b","",""] ==
            collect_eachmatch(R"""a?b?""", T("asbd"))
        @test collect_eachmatch(R"a?b?", T("asbd"), overlap=true) == ["a","","b","",""]
        @test collect_eachmatch(R"\w+", T("hello"), overlap=true) ==
            ["hello","ello","llo","lo","o"]
        @test collect_eachmatch(R"(\w+)(\s*)", T("The dark side of the moon")) ==
            ["The ", "dark ", "side ", "of ", "the ", "moon"]
        @test collect_eachmatch(R"", T("")) == [""]
        @test collect_eachmatch(R"", T(""), overlap=true) == [""]
        @test collect_eachmatch(R"aa", T("aaaa")) == ["aa", "aa"]
        @test collect_eachmatch(R"aa", T("aaaa"), overlap=true) == ["aa", "aa", "aa"]
        @test collect_eachmatch(R"", T("aaa")) == ["", "", "", ""]
        @test collect_eachmatch(R"", T("aaa"), overlap=true) == ["", "", "", ""]
        @test collect_eachmatch(R"GCG", T("GCGCG")) == ["GCG"]
        @test collect_eachmatch(R"GCG", T("GCGCG"),overlap=true) == ["GCG","GCG"]

# Issue 8278
target = """71.163.72.113 - - [30/Jul/2014:16:40:55 -0700] "GET emptymind.org/thevacantwall/wp-content/uploads/2013/02/DSC_006421.jpg HTTP/1.1" 200 492513 "http://images.search.yahoo.com/images/view;_ylt=AwrB8py9gdlTGEwADcSjzbkF;_ylu=X3oDMTI2cGZrZTA5BHNlYwNmcC1leHAEc2xrA2V4cARvaWQDNTA3NTRiMzYzY2E5OTEwNjBiMjc2YWJhMjkxMTEzY2MEZ3BvcwM0BGl0A2Jpbmc-?back=http%3A%2F%2Fus.yhs4.search.yahoo.com%2Fyhs%2Fsearch%3Fei%3DUTF-8%26p%3Dapartheid%2Bwall%2Bin%2Bpalestine%26type%3Dgrvydef%26param1%3D1%26param2%3Dsid%253Db01676f9c26355f014f8a9db87545d61%2526b%253DChrome%2526ip%253D71.163.72.113%2526p%253Dgroovorio%2526x%253DAC811262A746D3CD%2526dt%253DS940%2526f%253D7%2526a%253Dgrv_tuto1_14_30%26hsimp%3Dyhs-fullyhosted_003%26hspart%3Dironsource&w=588&h=387&imgurl=occupiedpalestine.files.wordpress.com%2F2012%2F08%2F5-peeking-through-the-wall.jpg%3Fw%3D588%26h%3D387&rurl=http%3A%2F%2Fwww.stopdebezetting.com%2Fwereldpers%2Fcompare-the-berlin-wall-vs-israel-s-apartheid-wall-in-palestine.html&size=49.0KB&name=...+%3Cb%3EApartheid+wall+in+Palestine%3C%2Fb%3E...+%7C+Or+you+go+peeking+through+the+%3Cb%3Ewall%3C%2Fb%3E&p=apartheid+wall+in+palestine&oid=50754b363ca991060b276aba291113cc&fr2=&fr=&tt=...+%3Cb%3EApartheid+wall+in+Palestine%3C%2Fb%3E...+%7C+Or+you+go+peeking+through+the+%3Cb%3Ewall%3C%2Fb%3E&b=0&ni=21&no=4&ts=&tab=organic&sigr=13evdtqdq&sigb=19k7nsjvb&sigi=12o2la1db&sigt=12lia2m0j&sign=12lia2m0j&.crumb=.yUtKgFI6DE&hsimp=yhs-fullyhosted_003&hspart=ironsource" "Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/36.0.1985.125 Safari/537.36"""

        pat = R"""([\d\.]+) ([\w.-]+) ([\w.-]+) (\[.+\]) "([^"\r\n]*|[^"\r\n\[]*\[.+\][^"]+|[^"\r\n]+.[^"]+)" (\d{3}) (\d+|-) ("(?:[^"]|\")+)"? ("(?:[^"]|\")+)"?"""

        match(pat, T(target))

        # issue #26829
        @test map(m -> m.match, eachmatch(R"^$|\S", "ö")) == ["ö"]

        # issue #26199
        @test map(m -> m.match, eachmatch(R"(\p{L}+)", "Tú")) == ["Tú"]
        @test map(m -> m.match, eachmatch(R"(\p{L}+)", "Tú lees.")) == ["Tú", "lees"]
        @test map(m -> m.match, eachmatch(R"(\p{L}+)", "¿Cuál es tu pregunta?")) == ["Cuál", "es", "tu", "pregunta"]

        # Issue 9545 (32 bit)
        buf = PipeBuffer()
        show(buf, R"")
        @test read(buf, String) == "R\"\""

        # see #10994, #11447: PCRE2 allows NUL chars in the pattern
        @test occurs_in(RegexStr(T("^a\0b\$")), T("a\0b"))

        # regex match / search string must be a String
        @test_throws ArgumentError match(R"test", GenericString("this is a test"))
        @test_throws ArgumentError fnd(First, R"test", GenericString("this is a test"))

        # Issue 27125
        @static if VERSION >= v"1.3"
            msg = "#Hello# from Julia"
            re = R"#(.+)# from (?<name>\w+)"
            subst = s"FROM: \g<name>\n MESSAGE: \1"
            @test replace(msg, re => subst) == "FROM: Julia\n MESSAGE: Hello"
        end

        # findall
        @test findall(R"\w+", "foo bar") == [1:3, 5:7]
        @test findall(R"\w+", "foo bar", overlap=true) == [1:3, 2:3, 3:3, 5:7, 6:7, 7:7]
        @test findall(R"\w*", "foo bar") == [1:3, 4:3, 5:7, 8:7]
        @test findall(R"\b", "foo bar") == [1:0, 4:3, 5:4, 8:7]

        # count
        @test count(R"\w+", "foo bar") == 2
        @test count(R"\w+", "foo bar", overlap=true) == 6
        @test count(R"\w*", "foo bar") == 4
        @test count(R"\b", "foo bar") == 4

        # Named subpatterns
        let m = match(R"(?<a>.)(.)(?<b>.)", T("xyz"))
            @test (m[:a], m[2], m["b"]) == ("x", "y", "z")
            typ = T === _LatinStr ? ASCIIStr : T
            @test sprint(show, m) == "RegexStrMatch{$typ}(\"xyz\", a=\"x\", 2=\"y\", b=\"z\")"
        end

        # Backcapture reference in substitution string
        @test replace(T("abcde"), R"(..)(?P<byname>d)" => s"\g<byname>xy\\\1") == "adxy\\bce"
        @test_throws ErrorException replace(T("a"), R"(?P<x>)" => s"\g<y>")
    end
end

for ST in UnicodeStringTypes
    @testset "$ST Regex Utility functions" begin
        foobarbaz = ST("foo,bar,baz")
        foo = ST("foo")
        bar = ST("bar")
        baz = ST("baz")
        abc = ST("abc")
        abcd = ST("abcd")
        @testset "rsplit/split" begin
            @test split(foobarbaz, R",") == [foo,bar,baz]

            let str = ST("a.:.ba..:..cba.:.:.dcba.:.")
                @test split(str, R"\.(:\.)+") == ["a","ba.",".cba","dcba",""]
                @test split(str, R"\.(:\.)+"; keepempty=false) == ["a","ba.",".cba","dcba"]
                @test split(str, R"\.+:\.+") == ["a","ba","cba",":.dcba",""]
                @test split(str, R"\.+:\.+"; keepempty=false) == ["a","ba","cba",":.dcba"]
            end

            # zero-width splits

            @test split(ST(""), R"") == [""]
            @test split(abc,  R"") == ["a","b","c"]
            @test split(abcd, R"b?") == ["a","c","d"]
            @test split(abcd, R"b*") == ["a","c","d"]
            @test split(abcd, R"b+") == ["a","cd"]
            @test split(abcd, R"b?c?") == ["a","d"]
            @test split(abcd, R"[bc]?") == ["a","","d"]
            @test split(abcd, R"a*") == ["","b","c","d"]
            @test split(abcd, R"a+") == ["","bcd"]
            @test split(abcd, R"d*") == ["a","b","c",""]
            @test split(abcd, R"d+") == [abc,""]
            @test split(abcd, R"[ad]?") == ["","b","c",""]
        end

        @testset "replace" begin
            @test replace(abcd, R"b?" => "^") == "^a^c^d^"
            @test replace(abcd, R"b+" => "^") == "a^cd"
            @test replace(abcd, R"b?c?" => "^") == "^a^d^"
            @test replace(abcd, R"[bc]?" => "^") == "^a^^d^"

            @test replace("foobarfoo", R"(fo|ba)" => "xx") == "xxoxxrxxo"
            @test replace("foobarfoo", R"(foo|ba)" => bar) == "barbarrbar"

            @test replace(ST("äƀçđ"), R"ƀ?" => "π") == "πäπçπđπ"
            @test replace(ST("äƀçđ"), R"ƀ+" => "π") == "äπçđ"
            @test replace(ST("äƀçđ"), R"ƀ?ç?" => "π") == "πäπđπ"
            @test replace(ST("äƀçđ"), R"[ƀç]?" => "π") == "πäππđπ"

            @test replace(ST("foobarfoo"), R"(fo|ba)" => "ẍẍ") == "ẍẍoẍẍrẍẍo"

            @test replace(ST("ḟøøbarḟøø"), R"(ḟø|ba)" => "xx") == "xxøxxrxxø"
            @test replace(ST("ḟøøbarḟøø"), R"(ḟøø|ba)" => bar) == "barbarrbar"

            @test replace(ST("fooƀäṙfoo"), R"(fo|ƀä)" => "xx") == "xxoxxṙxxo"
            @test replace(ST("fooƀäṙfoo"), R"(foo|ƀä)" => "ƀäṙ") == "ƀäṙƀäṙṙƀäṙ"

            @test replace(ST("ḟøøƀäṙḟøø"), R"(ḟø|ƀä)" => "xx") == "xxøxxṙxxø"
            @test replace(ST("ḟøøƀäṙḟøø"), R"(ḟøø|ƀä)" => "ƀäṙ") == "ƀäṙƀäṙṙƀäṙ"

            # for Char pattern call Char replacement function
            @test replace(ST("a"), R"a" => typeof) == "SubString{$ST}"
        end
    end
end


text:
using Nemo, GibbsTypePriors, JLD

grid_k = [collect(1:25:1000); 1000]

Cnk_numerical_accuracy = GibbsTypePriors.Cnk.(1000, grid_k, 0.6) |> x -> accuracy_bits.(x)
save("test/graphical_tests/saves_for_graphical_tests/accuracy_Cnk_1000.jld", "Cnk_numerical_accuracy", Cnk_numerical_accuracy)

ext: jl

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
using LinearAlgebra,Statistics

struct Regression
    α
    β
end

function Regression(x::AbstractVector,y::AbstractVector)
    x̄ = mean(x)
    ȳ = mean(y)
    sx = x .- x̄
    sy = y .- ȳ

    β̂ = dot(sx,sy) / dot(sx,sx)
    α̂ = ȳ - β̂*x̄

    return Regression(α̂,β̂)
end

model(p::Regression,ϵ,x) = p.α .+ (p.β .* x) .+ ϵ

max_forks_repo_name: langfzac/nbg-papers


text:
include(joinpath("gcn", "helpers.jl"))
include(joinpath("gcn", "math.jl"))
include(joinpath("gcn", "wavefront.jl"))
include(joinpath("gcn", "indexing.jl"))
include(joinpath("gcn", "assertion.jl"))
include(joinpath("gcn", "synchronization.jl"))
include(joinpath("gcn", "memory_static.jl"))
include(joinpath("gcn", "hostcall.jl"))
include(joinpath("gcn", "output.jl"))
include(joinpath("gcn", "memory_dynamic.jl"))
include(joinpath("gcn", "execution_control.jl"))
include(joinpath("gcn", "atomics.jl"))


text:
@testset "Maximum size rectangle binary sub-matrix with all 1s" begin

    matrix = [
        0  1  1  0  1
        1  1  0  1  0
        0  1  1  1  0
        1  1  1  1  0
        1  1  1  1  1
        0  0  0  0  0
    ]

    @test maximum_rectagular_submatrix(matrix) == (9, [(2, 4)], [(3, 5)])
    @test maximum_rectagular_submatrix(BitArray(matrix)) == (9, [(2, 4)], [(3, 5)])

    matrix = [
        0  1  1  0
        1  1  1  1
        1  1  1  1
        1  1  0  0
    ]

    @test maximum_rectagular_submatrix(matrix) == (8, [(1, 4)], [(2, 3)])
    @test maximum_rectagular_submatrix(BitArray(matrix)) == (8, [(1, 4)], [(2, 3)])
end


text:
src = [rand(10_000) for i in 1:100]

@testset "bubblesort" begin
    println("--------")
    println("tmap!")
    dst = similar(src)
    @time begin
        tmap!(dst, src) do x
            y = copy(x)
            n = length(x)
            for i in 1:n, j in 2:n
                if y[j-1] > y[j]
                    y[j-1], y[j] = y[j], y[j-1]
                end
            end
            y
        end
    end
    @test all(issorted.(dst))

    println("simple")
    function simple(src, dst)
        @Threads.threads for k in 1:length(src)
            x = src[k]
            y = copy(x)
            n = length(x)
            for i in 1:n, j in 2:n
                if y[j-1] > y[j]
                    y[j-1], y[j] = y[j], y[j-1]
                end
            end
            dst[k] = y
        end
    end

    dst = similar(src)
    @time simple(src, dst)
    @test all(issorted.(dst))
end


text:
module MyinterpolationA
export myinterpolationA

function myinterpolationA(grid, vals)
    function func(x)
        if x < grid[1] || x > grid[length(grid)]
            return "error"
        end

        num_2 = searchsortedfirst(grid, x)
        if num_2 == 1
            return vals[num_2]
        end

        num_1 = num_2 - 1
        x_1 = grid[num_1]
        x_2 = grid[num_2]
        y_1 = vals[num_1]
        y_2 = vals[num_2]

        y = ((y_2- y_1)/(x_2 - x_1))*(x - x_1) + y_1
        return y
    end
end

end


text:
using BinaryBuilder

cuda_version = v"10.1.243"  # NOTE: could be less specific

sources_linux = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/10.1_20191031/cudnn-10.1-linux-x64-v7.6.5.32.tgz",
                  "7eaec8039a2c30ab0bc758d303588767693def6bf49b22485a2c00bf2e136cb3")
]
sources_macos = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/10.1_20191031/cudnn-10.1-osx-x64-v7.6.5.32.tgz",
                  "8ecce28a5ed388a2b9b2d239e08d7c550f53b79288e6d9e5eb4c152bfc711aff")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/10.1_20191031/cudnn-10.1-windows10-x64-v7.6.5.32.zip",
                  "5e4275d738cc3a105cf6558b70b8a2ff514989ca1cd17bc8515086e20561a652")
]

include("../common.jl")


text:
using Plots	# ここの読み込みにとても時間がかかる
gr()		# バックエンドの設定

# ２要素スカラーのプロット
function plot_2scalar(time,value,filename)
	plot(time,value)
	savefig("./figs/" * filename * ".png")

end

function plot_2scalarp(time,value,filename)
	plot(time,value)
	savefig("./figs/" * filename * ".png", marker=:auto)

end

function plot_2scalar_range(time,value,xrange, yrange,filename)
	plot(time,value, xlims=(xrange[1], xrange[2]), ylims=(yrange[1], yrange[2]), marker = :auto)
	savefig("./figs/" * filename * ".png")

end

# ベクトルを3D上にプロット
function plot_vec(value,filename)
	plot(value[:,1], value[:, 2], value[:, 3])
	xlabel!("x")
	ylabel!("y")
	#zlabel!("z")

	savefig("./figs/" * filename * ".png")
end

function plot_vecs(value1, value2, valuename, filename)
	plot(value1[:,1], value1[:, 2], value1[:, 3], label = valuename[1])
	plot!(value2[:,1], value2[:, 2], value2[:, 3], label = valuename[2])
	xlabel!("x")
	ylabel!("y")
	#zlabel!("z")

	savefig("./figs/" * filename * ".png")
end

function plot_vec_range(value, xrange, yrange,zrange, filename)
	plot(value[:,1], value[:, 2], value[:, 3], xlims=(xrange[1], xrange[2]), ylims=(yrange[1], yrange[2]), zlims=(zrange[1], zrange[2]))
	xlabel!("x")
	ylabel!("y")
	#zlabel!("z")

	savefig("./figs/" * filename * ".png")
end

function plot_2vec_time(time,value,filename)
	plot(time, value[:, 1], value[:, 2])
	xlabel!("time")
	ylabel!("y")
	#zlabel!("z")

	savefig("./figs/" * filename * ".png")
end

function plot_ecef(x_ecef_log)
	plot(x_ecef_log[:,1], x_ecef_log[:, 2], x_ecef_log[:, 3])
	xlabel!("x")
	ylabel!("y")
	#zlabel!("z")

	savefig("./figs/" * "x_cecf.png")
end

function plot_3scalar(time, value1, value2, value3, valuename, filename)
	plot(time, value1, label = valuename[1])
	plot!(time, value2, label = valuename[2])
	plot!(time, value3, label = valuename[3])

	savefig("./figs/" * filename * ".png")
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
#include(joinpath(BridgeLandmarks.dir(),"scripts", "postprocessing.jl"))

"""
    landmarksmatching(
        xobs0::Array{PointF},xobsT::Array{PointF};
        pars = Pars_ms(),
        updatescheme = [:innov, :mala_mom],
        ITER = 100,
        outdir=@__DIR__,
        Σobs = nothing,
        ainit = nothing,
        printskip=20
    )

## Arguments
- `xobs0`: Array of PointF (coordinates of inital shape)
- `xobsT`: Array of PointF (coordinats of shape at time T)

## Optional arguments
- `pars`: either `pars_ms() or `pars_ahs()` (this selects the model and default parameter values)
- `updatescheme`: array of mcmc updates
- `ITER`: number of iterations
- `outdir`: path of directory to which output is written
- `Σobs`: If provided, a length two array, where the first (second) element gives an Array{Unc} elements, where each element specifies
    the covariance of the extrinsic noise for each landmark at time 0 (final time T). If not provided, Σobs is constructed using
    the value of pars.σobs
- `anit`: Hamiltonian kernel parameter. If not provided, defaults to setting
    `ainit = mean(norm.([xobs0[i]-xobs0[i-1] for i in 2:n]))`
- `printskip`: skip every printskip observations in writing output to console

## Example:
```
    # prepare some test data
    n = 7
    xobs0 = [PointF(2.0cos(t), sin(t))/4.0  for t in collect(0:(2pi/n):2pi)[2:end]]
    θ, ψ =  π/6, 0.1
    rot =  SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), cos(θ))
    stretch = SMatrix{2,2}(1.0 + ψ, 0.0, 0.0, 1.0 - ψ)
    shift = [0.5, -1.0]
    xobsT = [rot * stretch * xobs0[i] + shift  for i in 1:n]

    landmarksmatching(xobs0,xobsT; ITER=10)
    landmarksmatching(xobs0,xobsT; ITER=10, pars=BL.Pars_ahs())
```
"""
function landmarksmatching(
    xobs0::Array{PointF},xobsT::Array{PointF};
    pars = Pars_ms(),
    updatescheme = [:innov, :mala_mom],
    ITER = 100,
    outdir=@__DIR__,
    Σobs = nothing,
    ainit = nothing,
    printskip=20
    )

    model = pars.model
    n = length(xobs0)
    nshapes = 1
    @assert length(xobs0)==length(xobsT) "The two given landmark configurations do not have the same number of landmarks."

    if isnothing(Σobs)
        Σobs = fill([pars.σobs^2 * one(UncF) for i in 1:n],2) # noise on observations
    end
    T = 1.0; t = 0.0:pars.dt:T; tt = tc(t,T)
    obs_atzero = true
    fixinitmomentato0 = false
    subsamples = 0:pars.skip_saveITER:ITER
    obsinfo = set_obsinfo(xobs0,[xobsT],Σobs, obs_atzero,fixinitmomentato0)

    ################################# initialise P #################################
    if isnothing(ainit)
        ainit = 0.5*mean(norm.(diff(xobs0)))
    end
    cinit = pars.cinit
    γinit = pars.γinit
    if model == :ms
        P = MarslandShardlow(ainit, cinit, γinit, 0.0, n)
    elseif model == :ahs
        nfsinit = construct_nfs(pars.db, pars.stdev, γinit)
        P = Landmarks(ainit, cinit, n, pars.db , pars.stdev, nfsinit)
    end

    ################## prior specification with θ = (a, c, γ) ########################
    priorθ = product_distribution([pars.aprior, pars.γprior])
    priormom = MvNormalCanon(gramkernel(xobs0,P)/pars.κ)

    xinit = State(xobs0, zeros(PointF,P.n))
    mT = zeros(PointF, n)     # randn(PointF, n)

    start = time()
          Xsave, parsave, initendstates_save, accinfo, δ, ρ, δa =
                lm_mcmc(tt, obsinfo, mT, P, ITER, subsamples, xinit, pars, priorθ, priormom, updatescheme, outdir,printskip)
    elapsed = time() - start

    write_output(obsinfo.xobs0, obsinfo.xobsT, parsave, Xsave, initendstates_save, elapsed, accinfo, tt, n,nshapes,subsamples,ITER, updatescheme, Σobs, pars, ρ, δ, P, outdir)
    nothing
end

"""
    landmarksmatching(
        landmarks0::Matrix{Float64},landmarksT::Matrix{Float64};
        pars = Pars_ms(),
        updatescheme = [:innov, :mala_mom],
        ITER = 100,
        outdir=@__DIR__,
        Σobs = nothing,
        ainit = nothing
    )

landmarksmatching, where the input does not consist of arrays of points, but both the initial shape and final shape
are represented by a matrix. Each row of the matrix gives the coordinates of a landmark.
"""
function landmarksmatching(
    landmarks0::Matrix{Float64},landmarksT::Matrix{Float64};
    pars = Pars_ms(),
    updatescheme = [:innov, :mala_mom],
    ITER = 100,
    outdir=@__DIR__,
    Σobs = nothing,
    ainit = nothing,
    printskip=20
    )

    @assert size(landmarks0)==size(landmarksT)  "landmarks0 and landmarksT should have the same dimensions."
    # convert landmark coordinates to arrays of PointF
    xobs0 = [PointF(r...) for r in eachrow(landmarks0)]
    xobsT = [PointF(r...) for r in eachrow(landmarksT)]
    landmarksmatching(xobs0, xobsT; pars=pars,updatescheme=updatescheme,
        ITER=ITER, outdir=outdir, Σobs=Σobs, ainit=ainit)
    nothing
end


text:
using LinearAlgebra
using SparseArrays
using DifferentialEquations
using Statistics
using Sundials
using DataStructures
using Plots
using Plots.PlotMeasures
include("MC_functions.jl")

model_name = "US_Bank_Callcenter"
p = Dict("forwards" => true)
S = 470
t = 1440.0
u_0 = bd_stat(US_Bank_Q(0))
prob = ODEProblem(US_Bank_dudt, u_0, (0.0, 1440.0),p)
sol = solve(prob, CVODE_BDF(), reltol=1e-14, abstol=1e-14, saveat=5.0)

NS_perf_mes = zeros(288)
rhos = zeros(288)
lzs = zeros(288)
for i in range(1,stop=288)
    if ~ is_stochastic_vector(sol[i])
        sol[i][sol[i].<0].=0
        if ~ is_stochastic_vector(sol[i])
            print("Not stochastic after shift.")
        end
    end
    NS_perf_mes[i] = sol[i]' * US_Bank_Number_Waiting_r(sol.t[i])
    rhos[i] = US_Bank_Lambda(sol.t[i])/(1/3*US_Bank_Servers(sol.t[i]))
    lzs[i] = US_Bank_Lambda(sol.t[i])
end

perf_mes = zeros(288)
for i in range(1,stop=288)
    # perf_mes[i] = get_stat_dist_Q(US_Bank_Q(sol.t[i]))' * US_Bank_Number_Waiting_r(sol.t[i])
    print(i)
    perf_mes[i] = bd_stat_mode(US_Bank_Q(sol.t[i]))' * US_Bank_Number_Waiting_r(sol.t[i])
end

rel_error = (perf_mes.-NS_perf_mes)./NS_perf_mes

plot(sol.t[1:288], rel_error,
    title="PSA versus NS Relative Error", legend=:topleft,
    label="Abs (PSA-NS)/NS", xlabel="Time (Minutes)",
    ylabel="Abs (PSA-NS)/NS",right_margin=20mm)
plot!(NaN.*[1:288], label="rho(t)", grid=false,
    legend=:topleft,linecolor=:green)
plot!(twinx(), sol.t[1:288], rhos, label="rho(t)", legend=false,linecolor=:green, ylabel="rho(t)")
savefig(string(model_name,"PSA_vs_NS_relative_error_USBank.png"))

plot(sol.t[1:288], perf_mes, title="PSA and NSA Actual",
    label="PSA", xlabel="Time (Minutes)",yscale =:log,
        ylabel="Expected Number Waiting",
        right_margin=20mm, legend =:topleft)
plot!(sol.t[1:288], NS_perf_mes, label="NS", xlabel="Time (Minutes)",right_margin=20mm)
plot!(NaN.*[1:288], label="rho(t)", grid=false, legend=:left,linecolor=:green)
plot!(twinx(), sol.t[1:288], rhos, label="rho(t)", legend=false,linecolor=:green, ylabel="rho(t)")
savefig(string(model_name,"PSA_vs_NS_Actual_Values"))

# Also plot the different in measures in 2-norm.

# Same thing but for a simple sin M/M/s queue
# We look at the expected delay on such a thing.

model_name = "Queue_Saturated"

p = Dict("forwards" => true, "S"=> 600, "s"=>330, "mu"=> .333)
p["sin_base"] = .8*p["mu"]*p["s"]
p["sin_rel_amp"] = .1
p["phase_shift"] = 0
p["sin_rate"] = 2*pi/1440.0

u_0 = get_stat_dist_Q(Finite_Server_Q(0,p))
prob = ODEProblem(Finite_Server_dudt, u_0, (0.0, 1440.0),p)
sol = solve(prob, CVODE_BDF(), reltol=1e-12, abstol=1e-12, saveat=5.0)

NS_perf_mes = zeros(288)
rhos = zeros(288)
for i in range(1,stop=288)
    if ~ is_stochastic_vector(sol[i])
        sol[i][sol[i].<0].=0
        if ~ is_stochastic_vector(sol[i])
            print("Not stochastic after shift.")
        end
    end
    NS_perf_mes[i] = sol[i]' * Finite_Server_Number_Waiting_r(p)
    rhos[i] = lambda_func(sol.t[i], p)/(p["mu"]*p["s"])
end

perf_mes = zeros(288)
for i in range(1,stop=288)
    perf_mes[i] = (bd_stat(Finite_Server_Q(sol.t[i],p)))' * Finite_Server_Number_Waiting_r(p)
end

rel_error = (perf_mes.-NS_perf_mes)./NS_perf_mes

plot(sol.t[1:288], rel_error, title="PSA versus NS Relative Error", label="(PSA-NS)/NS",
        xlabel="Time (Minutes)", ylabel="(PSA-NS)/NS",right_margin=20mm)
plot!(NaN.*[1:288], label="rho(t)", grid=false, legend=:left,linecolor=:green)
plot!(twinx(), sol.t[1:288], rhos, label="rho(t)", legend=false,linecolor=:green, ylabel="rho(t)")
savefig(string(model_name,"PSA_vs_NS_relative_error_large_rho_USBankLevel.png"))

plot(sol.t[1:288], perf_mes, title="PSA and NS Actual",
 label="PSA", xlabel="Time (Minutes)", ylabel="Expected Number Waiting",
 right_margin=20mm)
plot!(sol.t[1:288], NS_perf_mes, label="NS")
plot!(NaN.*[1:288], label="rho(t)", grid=false,
    legend=:topleft, linecolor=:green)
plot!(twinx(), sol.t[1:288], rhos,
        label="rho(t)", legend=false,linecolor=:green, ylabel="rho(t)")
savefig(string(model_name,"PSA_vs_NS_Actual_Values_large_rho_USBankLevel.png"))

max_stars_repo_name: alexinfanger/blogfiles

max_issues_repo_name: alexinfanger/blogfiles

max_forks_repo_name: alexinfanger/blogfiles


text:
using GitLab, GitLab.name
using Base.Test

# The below tests are network-dependent, and actually make calls to GitLab's API

## testuser = Owner("julia-gitlab-test-bot")
testuser = Owner("mdpradeep")
julweb = Owner("TestProject1", true)
testcommit = Commit("93add52417601de893c3db3b30276c576933bf33")

hasghobj(obj, items) = any(x -> name(x) == name(obj), items)

auth = GitLab.authenticate(ENV["GITLAB_AUTH"]) # don't hardcode access tokens!
println("Authentication successful")
options = Dict("private_token" => auth.token)
ghjl = GitLab.repo_by_name("Calculus"; headers=options)

## @test rate_limit(; headers = options)["rate"]["limit"] == 5000

##########
# Owners #
##########

# test GitLab.owner
@test name(owner(testuser; headers = options)) == name(testuser)
## TODO CHECK @test name(owner(julweb; headers = options)) == name(julweb)

# test GitLab.orgs
## TODO CHECK @test hasghobj("TestProject1", orgs("TestProject1"; headers = options))

#= No equivalent APIs
# test GitLab.followers, GitLab.following
@test hasghobj("mdpradeep", first(followers(testuser; headers = options)))
@test hasghobj("mdpradeep", first(following(testuser; headers = options)))
=#

# test GitLab.repos
@test hasghobj(ghjl, first(repos(julweb; headers = options)))

################
# Repositories #
################

# test GitLab.repo
@test name(repo(get(ghjl.id); headers = options)) == name(ghjl)

#= No equivalent API
# test GitLab.forks
@test length(first(forks(ghjl; headers = options))) > 0
=#

# test GitLab.contributors
@test hasghobj("Pradeep Mudlapur", map(x->x["contributor"], first(contributors(ghjl; headers = options))))

# test GitLab.stats
## TODO - Enable @test stats(ghjl, "compound_metrics"; headers = options).status < 300

# test GitLab.branch, GitLab.branches
@test name(branch(ghjl, "master"; headers = options)) == "master"
@test hasghobj("master", first(branches(ghjl; headers = options)))

# test GitLab.commit, GitLab.commits
@test name(commit(ghjl, testcommit; headers = options)) == name(testcommit)
@test hasghobj(testcommit, first(commits(ghjl; headers = options)))

# test GitLab.file, GitLab.directory, GitLab.readme, GitLab.permalink
readme_file = file(ghjl, "README.md", "master"; headers = options)
#= No equivalent API - directory()
src_dir = first(directory(ghjl, "src"; headers = options))
owners_dir = src_dir[findfirst(c -> get(c.path) == "src/owners", src_dir)]
test_sha = "eab14e1ab7b4de848ef6390101b6d40b489d5d08"
readme_permalink = string(permalink(readme_file, test_sha))
owners_permalink = string(permalink(owners_dir, test_sha))
@test readme_permalink == "https://github.com/JuliaComputing/GitLab.jl/blob/$(test_sha)/README.md"
@test owners_permalink == "https://github.com/JuliaComputing/GitLab.jl/tree/$(test_sha)/src/owners"
@test hasghobj("src/GitLab.jl", src_dir)
=#

@test readme_file == readme(ghjl; headers = options)
# test GitLab.status, GitLab.statuses
## TODO No equivalent API @test get(status(ghjl, testcommit; headers = options).sha) == name(testcommit)
@test !(isempty(first(statuses(ghjl, testcommit; headers = options))))

test_repo = GitLab.repo_by_name("TestProject1"; headers=options)
# test GitLab.comment, GitLab.comments
## TODO - CHECK  @test name(comment(test_repo, 9, :commit; headers = options)) == 9
@test !(isempty(first(comments(test_repo, "5c35ae1de7f6d6bfadf0186e165f7af6537e7da8", :commit; headers = options))))

# These require `auth` to have push-access
@test hasghobj("mdpradeep", collaborators(ghjl; headers = options))
@test iscollaborator(ghjl, "mdpradeep"; headers = options)

##########
# Issues #
##########

## TODO - check
state_param = Dict("state" => "all")

# test GitLab.pull_request, GitLab.pull_requests
@test get(pull_request(test_repo, 3; headers = options).title) == "edit"
@test hasghobj(3, first(pull_requests(test_repo; headers = options, params = state_param)))

# test GitLab.issue, GitLab.issues
@test get(issue(test_repo, 6; headers = options).title) == "This is a test"
@test hasghobj(6, first(issues(test_repo; headers = options, params = state_param)))

############
# Activity #
############

# test GitLab.stargazers, GitLab.starred
## TODO - NO equivalent API @test length(first(stargazers(ghjl; headers = options))) > 10 # every package should fail tests if it's not popular enough :p
@test hasghobj(GitLab.repo_by_name("admin-project"; headers=options), first(starred(testuser; headers = options)))

#= TODO no equivalent API
# test GitLab.watched, GitLab.watched
@test hasghobj(testuser, first(watchers(test_repo; headers = options)))
@test hasghobj(ghjl, first(watched(testuser; headers = options)))
=#


text:
get_terms(e::E) where E <: Expression = (terms=Expression[];_get_terms!(e,terms);terms)
_get_terms!(e::Expression2{typeof(+),E1,E2},terms) where {E1,E2} = (_get_terms!(e.e1,terms);_get_terms!(e.e2,terms))
function _get_terms!(e::Expression2{typeof(*),E1,E2},terms,mul=1) where {E1,E2 <: ExpressionSum}
    _get_terms!(e.e2,terms,mul*e.e1)
end
function _get_terms!(e::Expression2{typeof(*),E1,E2},terms,mul=1) where {E1 <: ExpressionSum, E2}
    _get_terms!(e.e2,terms,mul*e.e2)
end
function _get_terms!(e::Expression2{typeof(*),E1,E2},terms,mul=1) where {E1 <: ExpressionSum, E2 <: ExpressionSum}
    error("not supported")
end
function _get_terms!(e::ExpressionSum{E,I},terms,mul=1) where {E,I}
    _get_terms!(e.inner,terms)
    for e in e.es
        _get_terms!(e,terms,mul)
    end
end
function _get_terms!(e::ExpressionSum{E,Nothing},terms,mul=1) where E
    for e in e.es
        _get_terms!(e,terms,mul)
    end
end
_get_terms!(e::E,terms,mul=1) where E = push!(terms,mul ==1 ? e : e*mul)

get_entries(f::Sink{Field}) = get_entries(inner(f))
get_entries(f::Field1{E,I}) where {E,I} = (entries=[];_get_entries!(f,entries);entries)
function _get_entries!(f::Field1{E,I},entries) where {E,I}
    _get_entries!(inner(f),entries)
    for e in f.es
        _get_entries!(e,entries)
    end
end
function _get_entries!(f::Field1{E,Nothing},entries) where E
    for e in f.es
        _get_entries!(e,entries)
    end
end
_get_entries!(e::E,entries) where E = push!(entries,e)

get_entries_expr(f::Sink{Field}) = get_entries_expr(inner(f))
function get_entries_expr(f::Field1{E,I}) where {E,I}
    entries = get_entries(f)
    fs = Vector{Expression}(undef,length(entries))
    for c in entries
        fs[index(c)] = c.e
    end
    return fs
end

sparsity(e::E) where E = (indices=Int[];_sparsity!(e::E,indices);indices)
_sparsity!(e::Variable,indices) = @inbounds union!(indices,index(e))
_sparsity!(e::Parameter,indices) = nothing
_sparsity!(e::Constant,indices) = nothing
_sparsity!(e::Expression1{F,F1},indices) where {F,F1} = _sparsity!(e.e1,indices)
_sparsity!(e::Expression2{F,F1,F2},indices) where {F,F1,F2} = (_sparsity!(e.e1,indices);_sparsity!(e.e2,indices))
_sparsity!(e::Expression2{F,F1,F2},indices) where {F,F1<:Real,F2} = _sparsity!(e.e2,indices)
_sparsity!(e::Expression2{F,F1,F2},indices) where {F,F1,F2<:Real} = _sparsity!(e.e1,indices)
function _sparsity!(e::ExpressionSum{E,I},indices) where {E,I}
    _sparsity!(e.inner,indices)
    @simd for i in eachindex(e.es)
        @inbounds  _sparsity!(e.es[i],indices)
    end
end
function _sparsity!(e::ExpressionSum{E,Nothing},indices) where E
    @simd for i in eachindex(e.es)
        @inbounds _sparsity!(e.es[i],indices)
    end
end


text:
using Flux
using Flux: onecold, onehotbatch
using Flux.Losses: logitcrossentropy
using GraphNeuralNetworks
using MLDatasets: Cora
using Statistics, Random
using CUDA
CUDA.allowscalar(false)

function eval_loss_accuracy(X, y, ids, model, g)
    ŷ = model(g, X)
    l = logitcrossentropy(ŷ[:,ids], y[:,ids])
    acc = mean(onecold(ŷ[:,ids]) .== onecold(y[:,ids]))
    return (loss = round(l, digits=4), acc = round(acc*100, digits=2))
end

# arguments for the `train` function
Base.@kwdef mutable struct Args
    η = 5f-3             # learning rate
    epochs = 10         # number of epochs
    seed = 17             # set seed > 0 for reproducibility
    usecuda = false      # if true use cuda (if available)
    nhidden = 64        # dimension of hidden features
end

function train(Layer; verbose=false, kws...)
    args = Args(; kws...)
    args.seed > 0 && Random.seed!(args.seed)

    if args.usecuda && CUDA.functional()
        device = Flux.gpu
        args.seed > 0 && CUDA.seed!(args.seed)
    else
        device = Flux.cpu
    end

    # LOAD DATA
    data = Cora.dataset()
    g = GNNGraph(data.adjacency_list, graph_type=GRAPH_T) |> device
    X = data.node_features |> device
    y = onehotbatch(data.node_labels, 1:data.num_classes) |> device
    train_ids = data.train_indices |> device
    val_ids = data.val_indices |> device
    test_ids = data.test_indices |> device
    ytrain = y[:,train_ids]

    nin, nhidden, nout = size(X,1), args.nhidden, data.num_classes

    ## DEFINE MODEL
    model = GNNChain(Layer(nin, nhidden),
                    #  Dropout(0.5),
                     Layer(nhidden, nhidden),
                     Dense(nhidden, nout))  |> device

    ps = Flux.params(model)
    opt = ADAM(args.η)

    ## TRAINING
    function report(epoch)
        train = eval_loss_accuracy(X, y, train_ids, model, g)
        test = eval_loss_accuracy(X, y, test_ids, model, g)
        println("Epoch: $epoch   Train: $(train)   Test: $(test)")
    end

    verbose && report(0)
    @time for epoch in 1:args.epochs
        gs = Flux.gradient(ps) do
            ŷ = model(g, X)
            logitcrossentropy(ŷ[:,train_ids], ytrain)
        end
        Flux.Optimise.update!(opt, ps, gs)
        verbose && report(epoch)
    end

    train_res = eval_loss_accuracy(X, y, train_ids, model, g)
    test_res = eval_loss_accuracy(X, y, test_ids, model, g)
    return train_res, test_res
end

function train_many(; usecuda=false)
    for (layer, Layer) in [
                ("GCNConv", (nin, nout) -> GCNConv(nin => nout, relu)),
                ("ResGatedGraphConv", (nin, nout) -> ResGatedGraphConv(nin => nout, relu)),
                ("GraphConv", (nin, nout) -> GraphConv(nin => nout, relu, aggr=mean)),
                ("SAGEConv", (nin, nout) -> SAGEConv(nin => nout, relu)),
                ("GATConv", (nin, nout) -> GATConv(nin => nout, relu)),
                ("GINConv", (nin, nout) -> GINConv(Dense(nin, nout, relu), 0.01, aggr=mean)),
                ## ("ChebConv", (nin, nout) -> ChebConv(nin => nout, 2)), # not working on gpu
                ## ("NNConv", (nin, nout) -> NNConv(nin => nout)),  # needs edge features
                ## ("GatedGraphConv", (nin, nout) -> GatedGraphConv(nout, 2)),  # needs nin = nout
                ## ("EdgeConv",(nin, nout) -> EdgeConv(Dense(2nin, nout, relu))), # Fits the traning set but does not generalize well
                ]

        @show layer
        @time train_res, test_res = train(Layer; usecuda, verbose=false)
        @test train_res.acc > 94
        @test test_res.acc > 70
    end
end

## if GRAPH_T != :dense # some erratic errors with :dense
train_many(usecuda=false)
if TEST_GPU
    train_many(usecuda=true)
end


text:
using DTCInitium
using Test

@testset "DTCInitium.jl" begin
    # Write your tests here.
end


text:
using .Samsara
using Conception
# using Conception: SAT, activate!

"""
mutable struct CircularSys <: AbstractSystem
This system is an implementation of circular double linked list,
where one can step both up and down:
Deterministic transmission :up => +1 and :down => -1

Constructor takes argument of type: Tuple{Any}, writing to _all_nodes::Tuple
Private functions _step_up! and _step_down! are used by step_system_mechanics!(CircularSys)
"""
mutable struct CircularSys <: Samsara.AbstractSystem
    _all_nodes
    _current_index::Int
    _all_actions      # key => value     value is [:up, :down, :nothing]
    _latent_variables
    function CircularSys(;nodes =(:A, :B, :C), actions = (:up, :down, nothing) )
        if length(actions) < 1
            throw(ArgumentError("Wrong number of actions: 3 actions are required for CircularSys, "*
                                "in an order representing the following:   :up, :down, :no-op"))
        end
        retval = new(nodes, 1, actions, nothing)
        #step_system_mechanics!(retval)
        return retval
    end
end

# Conception.activate!(everything_else) = return

"""
system_state(sys::CircularSys)
returnerer full system state for CircularSys -- a Tuple of observable_param
"""
function system_state(sys::CircularSys)
    (sys._all_nodes[sys._current_index], )
end

"""
current_state(sys::CircularSys)
Returns the current state, where we are currently located, only.
"""
current_state(sys::CircularSys) = sys._all_nodes[current_index(sys)]
current_index(sys::CircularSys) = sys._current_index
current_action(sys::CircularSys)= sys._latent_variables

function set_action_in!(sys::CircularSys, action)
    if action in sys._all_actions
        sys._latent_variables = action
    else
        throw(ArgumentError("set_action_in funker ikkje slik for det enkle circularsys."))
        #sys._latent_variables = nothing
    end
end

" Dispatch Base.length for this struct: "
function Base.length(sys::CircularSys)
    length(sys._all_nodes)
end

"""
    reset!(sys::CircularSys)
Reset current index to be the first: sys._current_index = 1
"""
function reset!(sys::CircularSys)
    return sys._current_index = 1
    #return sys._current_index = rand(1:length(sys._all_nodes))
end

"""
function step_up!(sys::CircularSys)
Step one up.
Return next, i.e. index of stepping one step "up".
When at max-position, go around to beginning
"""
function _step_up!(sys::CircularSys)
    (sys._current_index == length(sys)) && (return sys._current_index = 1)
    sys._current_index+=1
end

"""
function step_down!(sys::CircularSys)
Step one down.
Return next index, i.e. index of stepping one step "down"
When at min-position, go around to the end
"""
function _step_down!(sys::CircularSys)
    (sys._current_index == 1) && (return sys._current_index = length(sys))
    sys._current_index-=1
end

""" ======== Definerer system mechanics ===========
Gjennom function step_system_mechanics!(CircularSys)
kan vi definere korleis systemet funker.
"""
function step_system_mechanics!(sys::CircularSys)
    next_action = sys._latent_variables
    next_action == sys._all_actions[1] && (_step_up!(sys))
    next_action == sys._all_actions[2] && (_step_down!(sys))
    if (length(sys._all_actions) > 3) && (next_action == sys._all_actions[4])
        reset!(sys)
    end
    #:else:              && *stay where you are*

    # Fra Conception: TemporalType:
    # I tilfelle selv-transisjon er  det viktig å også deaktivere!
    activate!(next_action) # For å aktivere læringa!
    activate!(current_state(sys)) # For å aktivere læringa!

    sys._latent_variables = nothing

    current_state(sys)
end

export step_system_mechanics!, system_state, current_state, set_action_in!, current_action

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
#==============================================================================#
# AWSS3/test/runtests.jl
#
# Copyright Sam O'Connor 2014 - All rights reserved
#==============================================================================#

using AWSS3
using Base.Test
using Retry

using Compat.readstring

AWSCore.set_debug_level(1)

function test_without_catch(f)

    @protected try
        f()
    catch e
        @ignore if isa(e, Base.Test.Error)
            e = e.err
            rethrow(e)
        end
    end
end

#-------------------------------------------------------------------------------
# Load credentials...
#-------------------------------------------------------------------------------

aws = AWSCore.aws_config(region="ap-southeast-2")

#-------------------------------------------------------------------------------
# S3 tests
#-------------------------------------------------------------------------------

# Delete old test files...

for b in s3_list_buckets(aws)

    if ismatch(r"^ocaws.jl.test", b)
        @protected try
            println("Cleaning up old test bucket: " * b)
            for v in s3_list_versions(aws, b)
                s3_delete(aws, b, v["Key"]; version = v["VersionId"])
            end
            s3_delete_bucket(aws, b)
        catch e
            @ignore if e.code == "NoSuchBucket" end
        end
    end
end

# Temporary bucket name...

bucket_name = "ocaws.jl.test." * lowercase(Dates.format(now(Dates.UTC),
                                                        "yyyymmddTHHMMSSZ"))

# Test exception code for deleting non existand bucket...

@protected try

    s3_delete_bucket(aws, bucket_name)

catch e
     @ignore if e.code == "NoSuchBucket" end
end

# Create bucket...

s3_create_bucket(aws, bucket_name)

@repeat 4 try

    # Turn on object versioning for this bucket...

    s3_enable_versioning(aws, bucket_name)

    # Check that the new bucket is returned in the list of buckets...

    @test bucket_name in s3_list_buckets(aws)

    test_without_catch() do
        # Check that our test keys do not exist yet...
        @test !s3_exists(aws, bucket_name, "key1")
        @test !s3_exists(aws, bucket_name, "key2")
        @test !s3_exists(aws, bucket_name, "key3")
    end

catch e
    @delay_retry if e.code == "NoSuchBucket" end
end

# Create test objects...

s3_put(aws, bucket_name, "key1", "data1.v1")
s3_put(aws, bucket_name, "key2", "data2.v1")
s3_put(aws, bucket_name, "key3", "data3.v1")
s3_put(aws, bucket_name, "key3", "data3.v2")
s3_put(aws, bucket_name, "key3", "data3.v3")

# Check that test objects have expected content...

@test s3_get(aws, bucket_name, "key1") == b"data1.v1"
@test s3_get(aws, bucket_name, "key2") == b"data2.v1"
@test s3_get(aws, bucket_name, "key3") == b"data3.v3"

# Check object copy function...

s3_copy(aws, bucket_name, "key1";
        to_bucket = bucket_name, to_path = "key1.copy")

@test s3_get(aws, bucket_name, "key1.copy") == b"data1.v1"

url = s3_sign_url(aws, bucket_name, "key1")
curl_output = ""
@repeat 3 try
    curl_output = readstring(`curl -s -o - $url`)
catch e
    @delay_retry if true end
end
@test curl_output == "data1.v1"

fn = "/tmp/jl_qws_test_key1"
if isfile(fn)
    rm(fn)
end
@repeat 3 try
    s3_get_file(aws, bucket_name, "key1", fn)
catch e
    sleep(1)
    @retry if true end
end
@test readstring(fn) == "data1.v1"
rm(fn)

# Check exists and list objects functions...

for key in ["key1", "key2", "key3", "key1.copy"]
    @test s3_exists(aws, bucket_name, key)
    @test key in [o["Key"] for o in s3_list_objects(aws, bucket_name)]
end

# Check delete...

s3_delete(aws, bucket_name, "key1.copy")

@test !("key1.copy" in [o["Key"] for o in s3_list_objects(aws, bucket_name)])

# Check metadata...

meta = s3_get_meta(aws, bucket_name, "key1")
@test meta["ETag"] == "\"68bc8898af64159b72f349b391a7ae35\""

# Check versioned object content...

versions = s3_list_versions(aws, bucket_name, "key3")
@test length(versions) == 3
@test (s3_get(aws, bucket_name, "key3"; version = versions[3]["VersionId"])
      == b"data3.v1")
@test (s3_get(aws, bucket_name, "key3"; version = versions[2]["VersionId"])
      == b"data3.v2")
@test (s3_get(aws, bucket_name, "key3"; version = versions[1]["VersionId"])
      == b"data3.v3")

# Check pruning of old versions...

s3_purge_versions(aws, bucket_name, "key3")
versions = s3_list_versions(aws, bucket_name, "key3")
@test length(versions) == 1
@test s3_get(aws, bucket_name, "key3") == b"data3.v3"

#==============================================================================#
# End of file.
#==============================================================================#


text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.VariableTemplates: flattenednames
using CLIMA.ShallowWater
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.VTK
import CLIMA.ShallowWater:
    shallow_init_state!,
    shallow_init_aux!,
    vars_state,
    vars_aux,
    shallow_boundary_state!,
    TurbulenceClosure,
    LinearDrag,
    ConstantViscosity,
    AdvectionTerm,
    NonLinearAdvection

using CLIMAParameters
using CLIMAParameters.Planet: grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

struct GyreInABox{T} <: SWProblem
    τₒ::T
    fₒ::T # value includes τₒ, g, and ρ
    β::T
    Lˣ::T
    Lʸ::T
    H::T
end

function null_init_state!(
    ::GyreInABox,
    ::TurbulenceClosure,
    state,
    aux,
    coord,
    t,
)
    T = eltype(state.U)
    state.U = @SVector zeros(T, 3)
    state.η = 0
    return nothing
end

η_lsw(x, y, t) = cos(π * x) * cos(π * y) * cos(√2 * π * t)
u_lsw(x, y, t) = 2^(-0.5) * sin(π * x) * cos(π * y) * sin(√2 * π * t)
v_lsw(x, y, t) = 2^(-0.5) * cos(π * x) * sin(π * y) * sin(√2 * π * t)

function lsw_init_state!(
    p::GyreInABox,
    ::TurbulenceClosure,
    state,
    aux,
    coords,
    t,
)
    state.U = @SVector [
        u_lsw(coords[1], coords[2], t),
        v_lsw(coords[1], coords[2], t),
        0,
    ]

    state.η = η_lsw(coords[1], coords[2], t)

    return nothing
end

v_lkw(x, y, t) = 0
u_lkw(x, y, t) = exp(-0.5 * y^2) * exp(-0.5 * (x - t + 5)^2)
η_lkw(x, y, t) = 1 + u_lkw(x, y, t)

function lkw_init_state!(
    p::GyreInABox,
    ::TurbulenceClosure,
    state,
    aux,
    coords,
    t,
)
    state.U = @SVector [
        u_lkw(coords[1], coords[2], t),
        v_lkw(coords[1], coords[2], t),
        0,
    ]

    state.η = η_lkw(coords[1], coords[2], t)

    return nothing
end

R₋(ϵ) = (-1 - sqrt(1 + (2 * π * ϵ)^2)) / (2ϵ)
R₊(ϵ) = (-1 + sqrt(1 + (2 * π * ϵ)^2)) / (2ϵ)
D(ϵ) =
    (R₊(ϵ) * (exp(R₋(ϵ)) - 1) + R₋(ϵ) * (1 - exp(R₊(ϵ)))) /
    (exp(R₊(ϵ)) - exp(R₋(ϵ)))
R₂(x₁, ϵ) =
    (1 / D(ϵ)) * (
        (
            (R₊(ϵ) * (exp(R₋(ϵ)) - 1)) * exp(R₊(ϵ) * x₁) +
            (R₋(ϵ) * (1 - exp(R₊(ϵ)))) * exp(R₋(ϵ) * x₁)
        ) / (exp(R₊(ϵ)) - exp(R₋(ϵ)))
    )
R₁(x₁, ϵ) =
    (π / D(ϵ)) * (
        1 .+
        (
            (exp(R₋(ϵ)) - 1) * exp(R₊(ϵ) * x₁) .+
            (1 - exp(R₊(ϵ))) * exp(R₋(ϵ) * x₁)
        ) / (exp(R₊(ϵ)) - exp(R₋(ϵ)))
    )

𝒱(x₁, y₁, ϵ) = R₂(x₁, ϵ) * sin.(π * y₁)
𝒰(x₁, y₁, ϵ) = -R₁(x₁, ϵ) * cos.(π * y₁)
ℋ(x₁, y₁, ϵ, βᵖ, fₒ, γ) =
    (R₂(x₁, ϵ) / (π * fₒ)) * γ * cos(π * y₁) +
    (R₁(x₁, ϵ) / π) *
    (sin(π * y₁) * (1.0 + βᵖ * (y₁ - 0.5)) + (βᵖ / π) * cos(π * y₁))

function gyre_init_state!(p::GyreInABox, T::LinearDrag, state, aux, coords, t)
    FT = eltype(state)
    τₒ = p.τₒ
    fₒ = p.fₒ
    β = p.β
    Lˣ = p.Lˣ
    Lʸ = p.Lʸ
    H = p.H

    γ = T.λ

    βᵖ = β * Lʸ / fₒ
    ϵ = γ / (Lˣ * β)

    _grav::FT = grav(param_set)

    uˢ(ϵ) = (τₒ * D(ϵ)) / (H * γ * π)
    hˢ(ϵ) = (fₒ * Lˣ * uˢ(ϵ)) / _grav

    u = uˢ(ϵ) * 𝒰(coords[1] / Lˣ, coords[2] / Lʸ, ϵ)
    v = uˢ(ϵ) * 𝒱(coords[1] / Lˣ, coords[2] / Lʸ, ϵ)
    h = hˢ(ϵ) * ℋ(coords[1] / Lˣ, coords[2] / Lʸ, ϵ, βᵖ, fₒ, γ)

    state.U = @SVector [H * u, H * v, 0]

    state.η = h

    return nothing
end

t1(x, δᵐ) = cos((√3 * x) / (2 * δᵐ)) + (√3^-1) * sin((√3 * x) / (2 * δᵐ))
t2(x, δᵐ) = 1 - exp((-x) / (2 * δᵐ)) * t1(x, δᵐ)
t3(y, Lʸ) = π * sin(π * y / Lʸ)
t4(x, Lˣ, C) = C * (1 - x / Lˣ)

η_munk(x, y, Lˣ, Lʸ, δᵐ, C) = t4(x, Lˣ, C) * t3(y, Lʸ) * t2(x, δᵐ)

function gyre_init_state!(
    p::GyreInABox,
    V::ConstantViscosity,
    state,
    aux,
    coords,
    t,
)
    T = eltype(state.U)
    _grav::FT = grav(param_set)

    τₒ = p.τₒ
    fₒ = p.fₒ
    β = p.β
    Lˣ = p.Lˣ
    Lʸ = p.Lʸ
    H = p.H

    ν = V.ν

    δᵐ = (ν / β)^(1 / 3)
    C = τₒ / (_grav * H) * (fₒ / β)

    state.η = η_munk(coords[1], coords[2], Lˣ, Lʸ, δᵐ, C)
    state.U = @SVector zeros(T, 3)

    return nothing
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
import Distributions
using Statistics

include(joinpath(@__DIR__, "simulation_manipulation_with_unknown_parameters.jl"))

function state_res_norm(x::PhysManipulate2D)
    xVec = vec(x)[1:6];
    xVec_target = [0.,0.,pi,0.,0.,0.];
    return norm(xVec - xVec_target)
end;

function getMetricHistory(x_history)
    CostArray = Float64[];
    for t = 1:length(x_history)-1
        push!(CostArray,state_res_norm(x_history[t]))
    end
    #push!(CostArray,termCost_true(CostManipulate2D(),x_history[end],Cs));
    return CostArray
end

function evaluate()
    safe_rngnums = [];
    rngnum = 191;
    num_successful_cases = 0;

    cost_sacbp = Array{Float64}[];
    tcalc_sacbp = Array{Float64}[];
    cost_mcts = Array{Float64}[];
    tcalc_mcts = Array{Float64}[];
    cost_ilqg = Array{Float64}[];
    tcalc_ilqg = Array{Float64}[];

    time_idx_array = [];

    while num_successful_cases < 20
        println("------- RNG Seed: $(rngnum) & #Cases Tested: $(num_successful_cases) -------")
        println("1. Simulation with SACBP Controller")
        try
            y_history,tcalc_true_history,b_history,x_history,u_history,U_pool,Cs,Cu =
            simulate_main("sacbp",rng_num=rngnum, animate=false, plotfig=false, verbose=false);
            push!(cost_sacbp,getMetricHistory(x_history));
            push!(tcalc_sacbp,tcalc_true_history);
            if isempty(time_idx_array)
                push!(time_idx_array, unique([x.t for x in x_history]))
            end
        catch
            rngnum += 1
            continue
        end
        println("2. Simulation with MCTS Controller")
        try
            y_history,tcalc_true_history,b_history,x_history,u_history,U_pool,Cs,Cu =
            simulate_main("mcts",rng_num=rngnum, animate=false, plotfig=false, verbose=false);
            push!(cost_mcts,getMetricHistory(x_history));
            push!(tcalc_mcts,tcalc_true_history);
            if isempty(time_idx_array)
                push!(time_idx_array, unique([x.t for x in x_history]))
            end
        catch
            pop!(cost_sacbp)
            pop!(tcalc_sacbp)
            rngnum += 1
            continue
        end
        println("3. Simulation with Belief iLQG Controller")
        try
            y_history,tcalc_true_history,b_history,x_history,u_history,U_pool,Cs,Cu =
            simulate_main("ilqg",rng_num=rngnum, animate=false, plotfig=false, verbose=false);
            push!(cost_ilqg,getMetricHistory(x_history));
            push!(tcalc_ilqg,tcalc_true_history);
            if isempty(time_idx_array)
                push!(time_idx_array, unique([x.t for x in x_history]))
            end
        catch
            pop!(cost_sacbp)
            pop!(tcalc_sacbp)
            pop!(cost_mcts)
            pop!(tcalc_mcts)
            rngnum += 1
            continue
        end
        append!(safe_rngnums,rngnum)
        rngnum += 1;
        num_successful_cases += 1;
    end

    tcalc_sacbp = vcat(tcalc_sacbp...);
    tcalc_mcts = vcat(tcalc_mcts...);
    tcalc_ilqg = vcat(tcalc_ilqg...);

    data_dir = normpath(joinpath(@__DIR__, "../data",
                                 "manipulation_with_unknown_parameters"));

    FileIO.save(joinpath(data_dir, "statistical_data_manipulation.jld2"),
                "cost_sacbp",cost_sacbp,"tcalc_sacbp",tcalc_sacbp,
                "cost_mcts",cost_mcts,"tcalc_mcts",tcalc_mcts,
                "cost_ilqg",cost_ilqg,"tcalc_ilqg",tcalc_ilqg,
                "safe_rngnums",safe_rngnums)

    # Plot Stats Data
    cost_std_sacbp = [std([cost[ii] for cost in cost_sacbp]) for ii = 1:length(cost_sacbp[1])]
    cost_std_mcts = [std([cost[ii] for cost in cost_mcts]) for ii = 1:length(cost_mcts[1])]
    cost_std_ilqg = [std([cost[ii] for cost in cost_ilqg]) for ii = 1:length(cost_ilqg[1])];
    # # In this problem, cost is w.r.t. the true state history.
    CostHistoryMCTS = mean(cost_mcts)
    plot(time_idx_array[1][1:end-1],CostHistoryMCTS,xlabel="Time [s]",ribbon=cost_std_mcts,fillalpha=0.3,ylabel="Residual Norm",label="MCTS-DPW",color=:lightcoral,linewidth=2.,linestyle=:dash);
    CostHistoryiLQG = mean(cost_ilqg)
    plot!(time_idx_array[1][1:end-1],CostHistoryiLQG,label="Belief iLQG",color=:darkorange,ribbon=cost_std_ilqg,fillalpha=0.3,linewidth=2.,linestyle=:dashdot)
    CostHistorySACBP = mean(cost_sacbp)
    plot!(time_idx_array[1][1:end-1],CostHistorySACBP,label="SACBP",color=:darkblue,ribbon=cost_std_sacbp,fillalpha=0.3,linewidth=2.,size=(400,400),ylim=(0.,5.))
    savefig(joinpath(data_dir, "manipulation_costs.pdf"))
    # Plot Computation Time Comparison
    tcalcTrueArray = [tcalc_ilqg,tcalc_mcts,tcalc_sacbp]
    meanTcalcTrueArray = mean.(tcalcTrueArray);
    stdTcalcTrueArray = std.(tcalcTrueArray);
    bar([1,2,3],meanTcalcTrueArray,yerror=stdTcalcTrueArray,linecolor=:black,marker=stroke(4.0,:black, :dash),c=[:darkorange,:lightcoral,:darkblue],label="");
    xticks!([1,2,3],["Belief iLQG","MCTS-DPW","SACBP"]);
    plot!([0.5,3.5],[3/4*0.2,3/4*0.2],width=2.,label="Targeted",size=(300,400),color=:darkred,ylim=(-0.03,1.1))
    ylabel!("Control Computation Time [s]")
    savefig(joinpath(data_dir, "manipulation_time.pdf"))
end

if abspath(PROGRAM_FILE) == @__FILE__
    evaluate();
end


text:
@testset "Sanity Matrix/String" begin
    @test distance_matrix("", "")    == 0

    @testset "Sanity Matrix/String/Left" begin
        @test distance_matrix("a", "")   == 1
        @test distance_matrix("ab", "")  == 2
        @test distance_matrix("abc", "") == 3
    end

    @testset "Sanity Matrix/String/Right" begin
        @test distance_matrix("", "a")   == 1
        @test distance_matrix("", "ab")  == 2
        @test distance_matrix("", "abc") == 3
    end

    @testset "Sanity Matrix/String/Both" begin
        @test distance_matrix("a", "a")     == 0
        @test distance_matrix("a", "b")     == 1
        @test distance_matrix("ab", "bb")   == 1
        @test distance_matrix("abc", "bbc") == 1
        @test distance_matrix("abc", "abb") == 1
    end
end

max_stars_repo_name: jlegare/LevenshteinToolkit

max_issues_repo_name: jlegare/LevenshteinToolkit

max_forks_repo_name: jlegare/LevenshteinToolkit


text:
using Base.Cartesian
using LinearAlgebra

function non_optimized(xyz, forces)
    natoms = size(xyz, 2)
    e      = 0.0
    s      = 1.0

    for i = 1:natoms
        for j = 1:natoms
            if i == j
                continue
            end
            r = xyz[:,j] - xyz[:,i]
            d = norm(r)
            e += (s/d)^12 - (s/d)^6
            f = 12*s^12/d^13-6*s^6/d^7
            forces[:, i] += f
        end
    end

    return e
end

"""
    optimized(state::State{T}, cutoff::T, ::do_forces::Bool) where {DoF, T <: AbstractFloat}

        Calculates Lennard Jones energy of a State in a basic and inneficient
        way. If 'do_forces' is set to True, also calculate and update forces.
        Non bonded interactions are only considered for neighbouring atoms
        bellow the defined 'cutoff' (in nm). Returns the energy value.

	Example:

        naive_optimized(state, 1.2, true)
"""
function optimized(state::State{F, T}, cutoff::T, do_forces::Bool) where {F <: AbstractMatrix, T <: AbstractFloat}
    coords    = state.coords
    forces    = state.forces
    energy    = T(0)
    cutoff_sq = cutoff^2
    for i in 1:(length(coords) - 1)
        if do_forces
            forces[i] = zeros(3, 1)
            fi = zeros(3, 1)
        end
        for j in (i+1):length(coords)
            rij = coords[i] - coords[j]
            dij_sq = sum(rij.^2)
            if dij_sq < cutoff_sq
                lj2 = T(1) / dij_sq
                lj6 = lj2*lj2*lj2
                energy += T(1)*(lj6*lj6 - lj6)

                if do_forces
                    fc = T(24) * (lj6 - T(2) * lj6 * lj6) / dij_sq
                    fi .+= fc .* rij
                    forces[j:j+2] .-= fc .* rij
                end
            end
        end

        if do_forces
            forces[i:i+2] .+= fi
        end
    end

    return energy
end

"""
    serial_no_verlet(state::State{T}, cutoff::T, ::Type{Val{DoF}}) where {DoF, T <: AbstractFloat}

        Calculates Lennard Jones energy of a State. If Val{DoF} is set to True,
        also calculate and update forces. Non bonded interactions are only
        considered for neighbouring atoms bellow the defined 'cutoff' (in nm).
        Returns the energy value.

	Example:

        serial_no_verlet(state, 1.2, Val{true})
"""
@generated function serial(state::State{F, T}, cutoff::T, ::Type{Val{DoF}}) where {DoF, F <: AbstractMatrix, T <: AbstractFloat}
    quote
        coords = state.coords
        forces = state.forces
        natoms = state.n

        cutsq = convert(T, cutoff*cutoff)   # squared cutoff
        energy = T(0)                       # total energy
        σ = T(1)                            # LJ sigma

        @inbounds for i = 1:natoms-1

            # load coordinates for the i-th atom
            @nexprs 3 u -> ri_u = coords[i, u]

            #region FORCE_SECTION
            # zero the force accummulator for the i-th atom
            $(if DoF === true
                :(@nexprs 3 u -> fi_u = zero(T))
            end)
            #endregion

            for j = i+1:natoms

                # load coordinates for the j-th atom
                # and calculate the ij vector
                @nexprs 3 u -> rij_u = coords[j, u] - ri_u

                # calculate the squared distance. Skip
                # if greater than cutoff
                dij_sq = @reduce 3 (+) u -> rij_u*rij_u
                (dij_sq > cutsq) && continue

                # LJ potential
                lj2 = σ/dij_sq
                lj6 = lj2*lj2*lj2
                energy += lj6*lj6 - lj6

                #region FORCE_SECTION
                $(if DoF === true
                    quote
                        fc = T(24) * (lj6 - T(2) * lj6 * lj6) / dij_sq
                        @nexprs 3 u -> begin
                            # accumulate forces for atom i (in registers)
                            fi_u += fc * rij_u
                            # update forces for atom j
                            forces[j, u] -= fc * rij_u
                        end
                    end
                end)
                #endregion

            end
            #region FORCE_SECTION
            # update forces for the i-th atom
            $(if DoF === true
                :(@nexprs 3 u -> forces[i, u] += fi_u)
            end)
            #endregion
        end
        energy
    end
end

@generated function serial(state::State{F, T}, vlist::VerletList, ::Type{Val{DoF}}) where {DoF, F <: AbstractMatrix, T <: AbstractFloat}
    quote
        coords = state.coords
        forces = state.forces
        natoms = state.n

        cutsq = convert(T, vlist.cutoff*vlist.cutoff)   # squared cutoff
        energy = T(0)                                   # total energy
        σ = T(1)                                        # LJ sigma

        @inbounds for i = 1:natoms-1

            # ptr -> location of the first neighbor of atom i
            # ptr_stop -> location of the last neighbor of atom i
            ptr = vlist.offset[i]
            if vlist.list[ptr] < 1
                continue
            end
            ptr_stop = vlist.offset[i+1]-2

            # load coordinates for the i-th atom
            @nexprs 3 u -> ri_u = coords[i, u]

            #region FORCE_SECTION
            # zero the force accummulator for the i-th atom
            $(if DoF === true
                :(@nexprs 3 u -> fi_u = zero(T))
            end)
            #endregion

            # for j = i+1:natoms
            while ptr <= ptr_stop
                j = vlist.list[ptr]

                # load coordinates for the j-th atom
                # and calculate the ij vector
                @nexprs 3 u -> rij_u = coords[j, u] - ri_u

                # calculate the squared distance. Skip
                # if greater than cutoff
                dij_sq = @reduce 3 (+) u -> rij_u*rij_u
                # (dij_sq > cutsq) && continue

                # LJ potential
                lj2 = σ/dij_sq
                lj6 = lj2*lj2*lj2
                energy += lj6*lj6 - lj6

                #region FORCE_SECTION
                $(if DoF === true
                    quote
                        fc = T(24) * (lj6 - T(2) * lj6 * lj6) / dij_sq
                        @nexprs 3 u -> begin
                            # accumulate forces for atom i (in registers)
                            fi_u += fc * rij_u
                            # update forces for atom j
                            forces[j, u] -= fc * rij_u
                        end
                    end
                end)

                ptr += 1
                #endregion

            end
            #region FORCE_SECTION
            # update forces for the i-th atom
            $(if DoF === true
                :(@nexprs 3 u -> forces[i, u] += fi_u)
            end)
            #endregion
        end
        energy
    end
end


text:
using SparseArrays: sparse

export constant_constraint
"""
    constant_constraint(comp::CellComplex{N},
        cells::AbstractVector{Cell{N}}, k::Int) where N

Find a sparse representation of the constraint that a `k-1` form is constant on given
`k-1` dimensional cells.
"""
function constant_constraint(comp::CellComplex{N},
    cells::AbstractVector{Cell{N}}, k::Int) where N
    row_inds, col_inds, vals = Int[], Int[], Float64[]
    for (row_ind, (c1, c2)) in enumerate(zip(cells[1:end-1], cells[2:end]))
        col1 = findfirst(isequal(c1), comp.cells[k])
        col2 = findfirst(isequal(c2), comp.cells[k])
        push!(row_inds, row_ind); push!(col_inds, col1); push!(vals, 1)
        push!(row_inds, row_ind); push!(col_inds, col2); push!(vals, -1)
    end
    return sparse(row_inds, col_inds, vals, length(cells)-1, length(comp.cells[k]))
end

export zero_constraint
"""
    zero_constraint(comp::CellComplex{N}, cells::AbstractVector{Cell{N}}, k::Int) where N

Find a sparse representation of the constraint that a `k-1` form is zero on given `k-1`
dimensional cells.
"""
function zero_constraint(comp::CellComplex{N}, cells::AbstractVector{Cell{N}},
    k::Int) where N
    col_inds = [findfirst(isequal(c), comp.cells[k]) for c in cells]
    num_rows = length(col_inds)
    return sparse(1:num_rows, col_inds, ones(num_rows), num_rows, length(comp.cells[k]))
end


text:
module TimeSeriesA4P2

#Use this to print
#=
weave(Pkg.dir("$(pwd())\\FamaMacbeth.jl"),
  informat="script",
  out_path = "$(pwd())\\FamaMacbethOUT.html",
  doctype = "md2html")
=#

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD, ForwardDiff,
  Gadfly, JuMP, NLopt, HCubature, StaticArrays, CTNumerical,
  ParallelAccelerator, Measures, Formatting

importall CTIO, CTReg, CTStat

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const YAHOO_DATE_FORMAT = "m/d/yyyy"
const FRED_DATE_FORMAT = "m/d/yyyy"
const PPI_DATE_FORMAT = "m/d/yyyy"
const SP500_NAME = "GSPC"
const DBV_NAME = "DBV"
const FRED_NAME = "CPI"
const PPI_NAME = "PPI"

const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const DECIMALS = 2
const NUM_FORMAT = "{:.$(DECIMALS)f}"
const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}} #holds graphs
const OVERRIDE_STAR_STRINGS = ["\\ostar","\\dstar","\\tstar"]
const COLORS = ["DodgerBlue", "OrangeRed", "Green",
  "Brown", "Black", "BlueViolet"]

n2s(x::T where T<:Real) = num2Str(x, 3, Ints=true)

#preprocess PPI data
function preProcessPPI()::Void

  ppiDF::DataFrame = readtable("$DATA_PATH\\$PPI_NAME.csv")

  #read and pre-process the data
  dateFormat::String = PPI_DATE_FORMAT
  if :Date ∈ names(ppiDF)
    rename!(ppiDF, :Date, :DATE)
  elseif :_Date ∈ names(ppiDF)
    rename!(ppiDF, :_Date, :DATE)
  elseif :_DATE ∈ names(ppiDF)
    rename!(ppiDF, :_DATE, :DATE)
  end

  ppiDF[:,:DATE] = ((s::String)->Date(s, dateFormat)::Date).(ppiDF[:,:DATE])::DataVector{Date}
  sort!(ppiDF, cols = [:DATE])

  #difference

  ppiDF[:,:Δppi] = similar(ppiDF[:,:ppi])
  ppiDF[:,:Δppi] .= NA
  ppiDF[2:end,:Δppi] = ppiDF[2:end,:ppi] .- ppiDF[1:(end-1),:ppi]

  #log and log difference
  ppiDF[:,:lppi] = log.(ppiDF[:,:ppi])
  ppiDF[:,:Δlppi] = similar(ppiDF[:,:lppi])
  ppiDF[:,:Δlppi] .= NA
  ppiDF[2:end,:Δlppi] = ppiDF[2:end,:lppi] .- ppiDF[1:(end-1),:lppi]

  stream::IOStream = open("$DATA_PATH\\$PPI_NAME.jls", "w")
  serialize(stream, ppiDF)
  close(stream)

  return nothing
end

#to get the df
function getDF(pathStr::String, refreshData::Bool, preProcessFunc::Function)::DataFrame

  #re-process the data if needed or desired
  if !isfile(pathStr) || refreshData
    preProcessFunc()
  end

  #read the binary file
  stream::IOStream = open(pathStr)
  DF::DataFrame = deserialize(stream)
  close(stream)

  return DF
end

#problem specific graphs
function visualizeA4P2(Θ::Vector{ARMA}, ppiDF::DataFrame)::Void

  ACFRange::Range = 1:24
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  Gadfly.push_theme(:default)
  outDF = deepcopy(ppiDF)
  completecases!(outDF)

  push!(plotNames, "Summary Graphs")
  push!(plots, Vector{Plot}())
  for series::Symbol ∈ [:ppi, :Δppi, :lppi, :Δlppi]
    push!(plots[end],
      plot(outDF, x=:DATE, y=series,  Geom.line, Geom.point,
        Guide.xlabel("Date"),Guide.ylabel("$series"),
        Guide.title("Summary $series")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  push!(plotNames, "ACF Graphs")
  push!(plots, Vector{Plot}())
  for series::Symbol ∈ [:ppi, :Δppi, :lppi, :Δlppi]
    push!(plots[end],
      plot(x=collect(ACFRange), y=autoρ(Vector{Float64}(outDF[:,series]),ACFRange),
        Geom.line, Geom.point,
        Guide.xlabel("Date"),Guide.ylabel("$series ρ"),
        Guide.title("ACF $series")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  push!(plotNames, "PACF Graphs")
  push!(plots, Vector{Plot}())
  for series::Symbol ∈ [:ppi, :Δppi, :lppi, :Δlppi]
    push!(plots[end],
      plot(x=collect(ACFRange), y=pAutoρ(Vector{Float64}(outDF[:,series]),ACFRange),
        Geom.line, Geom.point,
        Guide.xlabel("Date"),Guide.ylabel("$series ∂ρ"),
        Guide.title("PACF $series")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  #write the graphs
  for i::Int ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 9inch, 7inch), plots[i])
    #display(plots[i])
  end

  return nothing
end

#script for the problem
function A4P2Script(;refreshData::Bool = true)::Void

  targetSeries::Symbol = :Δlppi

  #select the models
  testModel::ARMA =  ARMA(4,4)#ARMA(c=0.1,ϕ=[0.1,-0.2], ϑ=[0.1,0.2])
  models::Vector{ARMA} = [ARMA(1,1), ARMA(4,1), ARMA(4,4), ARMA(8,4), ARMA(8,8)]
  ppiDF::DataFrame = getDF("$DATA_PATH\\$PPI_NAME.jls", refreshData, preProcessPPI)
  #visualizeA4P2(models, ppiDF)

  #get the data in an array
  Y::Vector{Float64} = Vector{Float64}(dropna(ppiDF[:,:Δlppi]))
  generalizeARMA!(testModel)
  narrowARMA!(testModel)
  println("test log like: $(ARMALogLike(testModel, Y))")diag
  ARMAEstimate(testModel, Y)

  return nothing
end

@time begin
  A4P2Script(refreshData = true)
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module CovarianceEstimation

using Statistics
using StatsBase
using LinearAlgebra
import StatsBase: cov

export cov
export CovarianceEstimator, SimpleCovariance,
    LinearShrinkage,
    # Targets for linear shrinkage
    DiagonalUnitVariance, DiagonalCommonVariance, DiagonalUnequalVariance,
    CommonCovariance, PerfectPositiveCorrelation, ConstantCorrelation,
    # Eigendecomposition-based methods
    AnalyticalNonlinearShrinkage,
    # Biweight midcovariance
    BiweightMidcovariance

include("utils.jl")
include("biweight.jl")
include("linearshrinkage.jl")
include("nonlinearshrinkage.jl")

end # module


text:
#!/usr/bin/env julia
# Rendering slices and an isosurface of MRI data.

import GR
const gr3 = GR.gr3

function draw(mesh::Cint; x::Union{Real, Nothing}=nothing, y::Union{Real, Nothing}=nothing, z::Union{Real, Nothing}=nothing)
    gr3.clear()
    gr3.drawmesh(mesh, 1, (0,0,0), (0,0,1), (0,1,0), (1,1,1), (1,1,1))
    gr3.drawslicemeshes(data, x=x, y=y, z=z)
    GR.clearws()
    gr3.drawimage(0, 1, 0, 1, 500, 500, gr3.DRAWABLE_GKS)
    GR.updatews()
end

data = open(stream -> read!(stream, Array{UInt16}(undef, 93, 64, 64)), "mri.raw")
data = min.(data, 2000) / 2000.0 * typemax(UInt16)
data = convert(Array{UInt16, 3}, floor.(data))
data = permutedims(data, [3, 2, 1])

GR.setviewport(0, 1, 0, 1)
gr3.cameralookat(-3, 2, -2, 0, 0, 0, 0, 0, -1)

mesh = gr3.createisosurfacemesh(data, (2.0/63, 2.0/63, 2.0/92), (-1.0, -1.0, -1.0), 40000)

GR.setcolormap(1)
for z in 0:0.005:1
    draw(mesh, x=0.9, z=z)
end
for y in 1:-0.0025:0.5
    draw(mesh, x=0.9, y=y, z=1)
end
GR.setcolormap(19)
for x in 0.9:-0.003:0
    draw(mesh, x=x, y=0.5, z=1)
end
for x in 0:0.003:0.9
    draw(mesh, x=x, z=1)
end
gr3.terminate()


text:
using LinearAlgebra, JuMP, Mosek, MosekTools, SCS, LowRankApprox, Random, Dates
using TSVD, Distributions, StatsBase

include("src/BnB.jl")
include("src/SPCP.jl")


text:
# Implement kernels for polar formulation

"""
    power_balance(V, Ybus, Sbus, pv, pq)

Assembly residual function for N-R power flow.
In complex form, the balance equations writes

``
g(V) = V * (Y_{bus} * V)^* - S_{bus}
``

# Note
Code adapted from MATPOWER.
"""
function power_balance(V, Ybus, Sbus, pv, pq)
    # form mismatch vector
    mis = V .* conj(Ybus * V) - Sbus
    # form residual vector
    F = [real(mis[pv])
         real(mis[pq])
         imag(mis[pq]) ]
    return F
end

"""
    residual_jacobian(V, Ybus, pv, pq)

Compute the Jacobian w.r.t. the state `x` of the power
balance function [`power_balance`](@ref).

# Note
Code adapted from MATPOWER.
"""
function residual_jacobian(V, Ybus, pv, pq)
    n = size(V, 1)
    Ibus = Ybus*V
    diagV       = sparse(1:n, 1:n, V, n, n)
    diagIbus    = sparse(1:n, 1:n, Ibus, n, n)
    diagVnorm   = sparse(1:n, 1:n, V./abs.(V), n, n)

    dSbus_dVm = diagV * conj(Ybus * diagVnorm) + conj(diagIbus) * diagVnorm
    dSbus_dVa = 1im * diagV * conj(diagIbus - Ybus * diagV)

    j11 = real(dSbus_dVa[[pv; pq], [pv; pq]])
    j12 = real(dSbus_dVm[[pv; pq], pq])
    j21 = imag(dSbus_dVa[pq, [pv; pq]])
    j22 = imag(dSbus_dVm[pq, pq])

    J = [j11 j12; j21 j22]
end

function _state_jacobian(polar::PolarForm)
    pf = polar.network
    ref = polar.network.ref
    pv = polar.network.pv
    pq = polar.network.pq
    n = PS.get(pf, PS.NumberOfBuses())

    Y = pf.Ybus
    # Randomized inputs
    Vre = rand(n)
    Vim = rand(n)
    V = Vre .+ 1im .* Vim
    return residual_jacobian(V, Y, pv, pq)
end
_sparsity_pattern(polar::PolarForm) = findnz(_state_jacobian(polar))

"""
    get_power_injection(fr, v_m, v_a, ybus_re, ybus_im)

Computes the power injection at node `fr`.
In polar form, the power injection at node `i` satisfies
```math
p_{i} = \\sum_{j} v_{i} v_{j} (g_{ij} \\cos(\\theta_i - \\theta_j) + b_{ij} \\sin(\\theta_i - \\theta_j))
```
"""
function get_power_injection(fr, v_m, v_a, ybus_re, ybus_im)
    P = 0.0
    for c in ybus_re.colptr[fr]:ybus_re.colptr[fr+1]-1
        to = ybus_re.rowval[c]
        aij = v_a[fr] - v_a[to]
        P += v_m[fr]*v_m[to]*(ybus_re.nzval[c]*cos(aij) + ybus_im.nzval[c]*sin(aij))
    end
    return P
end

"""
    get_react_injection(fr, v_m, v_a, ybus_re, ybus_im)

Computes the reactive power injection at node `fr`.
In polar form, the power injection at node `i` satisfies
```math
q_{i} = \\sum_{j} v_{i} v_{j} (g_{ij} \\sin(\\theta_i - \\theta_j) - b_{ij} \\cos(\\theta_i - \\theta_j))
```
"""
function get_react_injection(fr::Int, v_m, v_a, ybus_re::Spmat{VI,VT}, ybus_im::Spmat{VI,VT}) where {VT <: AbstractVector, VI<:AbstractVector}
    Q = zero(eltype(v_m))
    for c in ybus_re.colptr[fr]:ybus_re.colptr[fr+1]-1
        to = ybus_re.rowval[c]
        aij = v_a[fr] - v_a[to]
        Q += v_m[fr]*v_m[to]*(ybus_re.nzval[c]*sin(aij) - ybus_im.nzval[c]*cos(aij))
    end
    return Q
end

@kernel function residual_kernel!(F, v_m, v_a,
                                  ybus_re_nzval, ybus_re_colptr, ybus_re_rowval,
                                  ybus_im_nzval, ybus_im_colptr, ybus_im_rowval,
                                  pinj, qinj, pv, pq, nbus)

    npv = size(pv, 1)
    npq = size(pq, 1)

    i = @index(Global, Linear)
    # REAL PV: 1:npv
    # REAL PQ: (npv+1:npv+npq)
    # IMAG PQ: (npv+npq+1:npv+2npq)
    fr = (i <= npv) ? pv[i] : pq[i - npv]
    F[i] -= pinj[fr]
    if i > npv
        F[i + npq] -= qinj[fr]
    end
    @inbounds for c in ybus_re_colptr[fr]:ybus_re_colptr[fr+1]-1
        to = ybus_re_rowval[c]
        aij = v_a[fr] - v_a[to]
        # f_re = a * cos + b * sin
        # f_im = a * sin - b * cos
        coef_cos = v_m[fr]*v_m[to]*ybus_re_nzval[c]
        coef_sin = v_m[fr]*v_m[to]*ybus_im_nzval[c]
        cos_val = cos(aij)
        sin_val = sin(aij)
        F[i] += coef_cos * cos_val + coef_sin * sin_val
        if i > npv
            F[npq + i] += coef_cos * sin_val - coef_sin * cos_val
        end
    end
end

function residual_polar!(F, v_m, v_a,
                         ybus_re, ybus_im,
                         pinj, qinj, pv, pq, nbus)
    npv = length(pv)
    npq = length(pq)
    if isa(F, Array)
        kernel! = residual_kernel!(CPU(), 4)
    else
        kernel! = residual_kernel!(CUDADevice(), 256)
    end
    ev = kernel!(F, v_m, v_a,
                 ybus_re.nzval, ybus_re.colptr, ybus_re.rowval,
                 ybus_im.nzval, ybus_im.colptr, ybus_im.rowval,
                 pinj, qinj, pv, pq, nbus,
                 ndrange=npv+npq)
    wait(ev)
end

@kernel function transfer_kernel!(vmag, vang, pinj, qinj,
                        x, u, p, pv, pq, ref,
                        ybus_re_nzval, ybus_re_colptr, ybus_re_rowval,
                        ybus_im_nzval, pload)
    i = @index(Global, Linear)
    npv = length(pv)
    npq = length(pq)
    nref = length(ref)

    # PV bus
    if i <= npv
        bus = pv[i]
        qacc = 0
        @inbounds for c in ybus_re_colptr[bus]:ybus_re_colptr[bus+1]-1
            to = ybus_re_rowval[c]
            aij = vang[bus] - vang[to]
            # f_re = a * cos + b * sin
            # f_im = a * sin - b * cos
            coef_cos = vmag[bus]*vmag[to]*ybus_re_nzval[c]
            coef_sin = vmag[bus]*vmag[to]*ybus_im_nzval[c]
            cos_val = cos(aij)
            sin_val = sin(aij)
            qacc += coef_cos * sin_val - coef_sin * cos_val
        end
        qinj[bus] = qacc
        vang[bus] = x[i]
        pinj[bus] = u[nref + i] - pload[bus]
        vmag[bus] = u[nref + npv + i]
    # PQ bus
    elseif i <= npv + npq
        i_pq = i - npv
        bus = pq[i_pq]
        vang[bus] = x[npv+i_pq]
        vmag[bus] = x[npv+npq+i_pq]
        pinj[bus] = p[nref + i_pq]
        qinj[bus] = p[nref + npq + i_pq]
    # REF bus
    else i <= npv + npq + nref
        i_ref = i - npv - npq
        bus = ref[i_ref]
        pacc = 0
        qacc = 0
        @inbounds for c in ybus_re_colptr[bus]:ybus_re_colptr[bus+1]-1
            to = ybus_re_rowval[c]
            aij = vang[bus] - vang[to]
            # f_re = a * cos + b * sin
            # f_im = a * sin - b * cos
            coef_cos = vmag[bus]*vmag[to]*ybus_re_nzval[c]
            coef_sin = vmag[bus]*vmag[to]*ybus_im_nzval[c]
            cos_val = cos(aij)
            sin_val = sin(aij)
            pacc += coef_cos * cos_val + coef_sin * sin_val
            qacc += coef_cos * sin_val - coef_sin * cos_val
        end
        pinj[bus] = pacc
        qinj[bus] = qacc
        vmag[bus] = u[i_ref]
        vang[bus] = p[i_ref]
    end
end

# Transfer values in (x, u, p) to buffer
function transfer!(polar::PolarForm, buffer::PolarNetworkState, x, u, p)
    if isa(x, Array)
        kernel! = transfer_kernel!(CPU(), 1)
    else
        kernel! = transfer_kernel!(CUDADevice(), 256)
    end
    nbus = length(buffer.vmag)
    pv = polar.indexing.index_pv
    pq = polar.indexing.index_pq
    ref = polar.indexing.index_ref
    ev = kernel!(
        buffer.vmag, buffer.vang, buffer.pinj, buffer.qinj,
        x, u, p,
        pv, pq, ref,
        polar.ybus_re.nzval, polar.ybus_re.colptr, polar.ybus_re.rowval,
        polar.ybus_im.nzval, polar.active_load,
        ndrange=nbus
    )
    wait(ev)
end

@kernel function active_power_kernel!(
    pg, vmag, vang, pinj, qinj,
    pv, ref, pv_to_gen, ref_to_gen,
    ybus_re_nzval, ybus_re_colptr, ybus_re_rowval,
    ybus_im_nzval, pload
)
    i = @index(Global, Linear)
    npv = length(pv)
    nref = length(ref)
    # Evaluate active power at PV nodes
    if i <= npv
        bus = pv[i]
        i_gen = pv_to_gen[i]
        pg[i_gen] = pinj[bus] + pload[bus]
    # Evaluate active power at slack nodes
    elseif i <= npv + nref
        i_ = i - npv
        bus = ref[i_]
        i_gen = ref_to_gen[i_]
        inj = 0
        @inbounds for c in ybus_re_colptr[bus]:ybus_re_colptr[bus+1]-1
            to = ybus_re_rowval[c]
            aij = vang[bus] - vang[to]
            # f_re = a * cos + b * sin
            # f_im = a * sin - b * cos
            coef_cos = vmag[bus]*vmag[to]*ybus_re_nzval[c]
            coef_sin = vmag[bus]*vmag[to]*ybus_im_nzval[c]
            cos_val = cos(aij)
            sin_val = sin(aij)
            inj += coef_cos * cos_val + coef_sin * sin_val
        end
        pg[i_gen] = inj + pload[bus]
    end
end

# Refresh active power (needed to evaluate objective)
function refresh!(polar::PolarForm, ::PS.Generator, ::PS.ActivePower, buffer::PolarNetworkState)
    if isa(buffer.vmag, Array)
        kernel! = active_power_kernel!(CPU(), 1)
    else
        kernel! = active_power_kernel!(CUDADevice(), 256)
    end
    pv = polar.indexing.index_pv
    pq = polar.indexing.index_pq
    ref = polar.indexing.index_ref
    pv_to_gen = polar.indexing.index_pv_to_gen
    ref_to_gen = polar.indexing.index_ref_to_gen

    range_ = length(pv) + length(ref)

    ev = kernel!(
        buffer.pg,
        buffer.vmag, buffer.vang, buffer.pinj, buffer.qinj,
        pv, ref, pv_to_gen, ref_to_gen,
        polar.ybus_re.nzval, polar.ybus_re.colptr, polar.ybus_re.rowval,
        polar.ybus_im.nzval, polar.active_load,
        ndrange=range_
    )
    wait(ev)
end

@kernel function reactive_power_kernel!(
    qg, vmag, vang, pinj, qinj,
    pv, ref, pv_to_gen, ref_to_gen,
    ybus_re_nzval, ybus_re_colptr, ybus_re_rowval,
    ybus_im_nzval, qload
)
    i = @index(Global, Linear)
    npv = length(pv)
    nref = length(ref)
    # Evaluate reactive power at PV nodes
    if i <= npv
        bus = pv[i]
        i_gen = pv_to_gen[i]
    # Evaluate reactive power at slack nodes
    elseif i <= npv + nref
        i_ = i - npv
        bus = ref[i_]
        i_gen = ref_to_gen[i_]
    end
    inj = 0
    @inbounds for c in ybus_re_colptr[bus]:ybus_re_colptr[bus+1]-1
        to = ybus_re_rowval[c]
        aij = vang[bus] - vang[to]
        # f_re = a * cos + b * sin
        # f_im = a * sin - b * cos
        coef_cos = vmag[bus]*vmag[to]*ybus_re_nzval[c]
        coef_sin = vmag[bus]*vmag[to]*ybus_im_nzval[c]
        cos_val = cos(aij)
        sin_val = sin(aij)
        inj += coef_cos * sin_val - coef_sin * cos_val
    end
    qg[i_gen] = inj + qload[bus]
end

function refresh!(polar::PolarForm, ::PS.Generator, ::PS.ReactivePower, buffer::PolarNetworkState)
    if isa(buffer.vmag, Array)
        kernel! = reactive_power_kernel!(CPU(), 1)
    else
        kernel! = reactive_power_kernel!(CUDADevice(), 256)
    end
    pv = polar.indexing.index_pv
    pq = polar.indexing.index_pq
    ref = polar.indexing.index_ref
    pv_to_gen = polar.indexing.index_pv_to_gen
    ref_to_gen = polar.indexing.index_ref_to_gen

    range_ = length(pv) + length(ref)
    ev = kernel!(
        buffer.qg,
        buffer.vmag, buffer.vang, buffer.pinj, buffer.qinj,
        pv, ref, pv_to_gen, ref_to_gen,
        polar.ybus_re.nzval, polar.ybus_re.colptr, polar.ybus_re.rowval,
        polar.ybus_im.nzval, polar.reactive_load,
        ndrange=range_
    )
    wait(ev)
end

@kernel function load_power_constraint_kernel!(
    g, qg, ref_to_gen, pv_to_gen, nref, npv, shift
)
    i = @index(Global, Linear)
    # Evaluate reactive power at PV nodes
    if i <= npv
        ig = pv_to_gen[i]
        g[i + nref + shift] = qg[ig]
    else i <= npv + nref
        i_ = i - npv
        ig = ref_to_gen[i_]
        g[i_ + shift] = qg[ig]
    end
end


text:
## calculatoin_tools.jl
export pack
export radial_mean, Δ_phase
export moment_proj_normed

"""
    radial_mean(data; maxbin=nothing, bin_step=nothing, pixelsize=nothing)

Calculates the radial mean of a dataset `data`.
Returns a tuple of the radial_mean and the bin_centers.

 ## Arguments:
+ `data`: data to radially average
+ `maxbin`: a maximum bin value
+ `bin_step`: optionally defines the step between the bins

## Examples
```julia-repl
julia> radial_mean([5, 1, 1, 3, 15])
([1.0, 2.0, 10.0], [0.5, 1.5, 2.5])
```
"""
function radial_mean(data; nbins=nothing, bin_step=nothing, offset=CtrFT, scale=nothing)
    if isnothing(scale)
        scale=Tuple(ones(Int,ndims(data)))
    end
    if isnothing(bin_step)
        bin_step=maximum(scale)
    end

    idx = round.(Int, rr(size(data), offset=offset, scale=scale./bin_step)) .+ 1
    if isnothing(nbins)
        nbins = maximum(idx)
    else
        idx[idx.>nbins] .= nbins
    end

    bincenters = bin_step .* collect(0.5:nbins-0.5)

    myevents = zeros(Float64, nbins)
    mysum = zeros(Float64, nbins)
    idxlist = idx[:]
    v = @view mysum[idxlist]
    v .+= data[:]
    v = @view myevents[idxlist]
    v .+= 1
    myevents[myevents .== 0] .= 1
    return mysum ./ myevents, bincenters
end

"""
    Δ_phase(arr, dim)

Calculates the relative phase slope along dimension `dim` of a non-zero array `arr` without wrap-around problems.

Arguments:
+ arr: array of which to evaluate the phase slope
+ dim: dimension along which to evaluate the slope
"""
function Δ_phase(arr, dim)
    no_shift = ((n==dim) ? (1:size(arr,dim)-1) : (:) for n in 1:ndims(arr))
    shift = ((n==dim) ? (2:size(arr,dim)) : (:) for n in 1:ndims(arr))
    return angle.(arr[no_shift...] ./ arr[shift...]) # this is a complex-valued division to obtain the relative phase angle!
end

"""
    moment_proj(data, h=3; pdims=3)

Performs a projection over the `h`-th moment of the data along dimension(s) `pdims`.
See also [`moment_proj_normed`](@ref MicroscopyTools.moment_proj_normed)

## Example
```jldoctest
julia> MicroscopyTools.moment_proj([1 3; 2 4], 2, pdims=(1,))
1×2 Matrix{Float64}:
 0.25  0.25
```
"""
function moment_proj(data, h=3; pdims=3)
    mdata = mean(data, dims=pdims)
    return mean((data .- mdata).^h, dims=pdims)
end

"""
    moment_proj_normed(data, h=3; pdims=3)

Performs a projection over the `h`-th moment of the data along dimension(s) `pdims` normed by the variance.

## Example
```jldoctest
julia> MicroscopyTools.moment_proj_normed([1 3; 2 4], 2, pdims=(1,))
1×2 Matrix{Float64}:
 0.5  0
```
"""
function moment_proj_normed(data, h=3; pdims=3)
    moment_proj(data, h, pdims=pdims) ./ (moment_proj(data,2, pdims=pdims) .^((h-1)/2))
end

max_stars_repo_name: JuliaMicroscopy/MicroscopyTools

max_issues_repo_name: JuliaMicroscopy/MicroscopyTools


text:
module BandedMatrices
using Base, FillArrays, ArrayLayouts, LinearAlgebra, SparseArrays, Random
using LinearAlgebra.LAPACK
import Base: axes, axes1, getproperty, iterate, tail
import LinearAlgebra: BlasInt, BlasReal, BlasFloat, BlasComplex, axpy!,
                        copy_oftype, checksquare, adjoint, transpose, AdjOrTrans, HermOrSym,
                        _chol!, rot180
import LinearAlgebra.BLAS: libblas
import LinearAlgebra.LAPACK: liblapack, chkuplo, chktrans
import LinearAlgebra: cholesky, cholesky!, cholcopy, norm, diag, eigvals!, eigvals, eigen!, eigen,
            qr, qr!, axpy!, ldiv!, mul!, lu, lu!, ldlt, ldlt!, AbstractTriangular,
            chkstride1, kron, lmul!, rmul!, factorize, StructuredMatrixStyle, logabsdet,
            svdvals, svdvals!, QRPackedQ, checknonsingular, ipiv2perm, tril!,
            triu!, Givens, diagzero
import SparseArrays: sparse

import Base: getindex, setindex!, *, +, -, ==, <, <=, >, isassigned,
                >=, /, ^, \, transpose, showerror, reindex, checkbounds, @propagate_inbounds

import Base: convert, size, view, unsafe_indices,
                first, last, size, length, unsafe_length, step,
                to_indices, to_index, show, fill!, promote_op,
                MultiplicativeInverses, OneTo, ReshapedArray,
                               similar, copy, convert, promote_rule, rand,
                            IndexStyle, real, imag, Slice, pointer, unsafe_convert, copyto!,
                            hcat, vcat, hvcat

import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Broadcasted, broadcasted,
                        materialize, materialize!

import ArrayLayouts: MemoryLayout, transposelayout, triangulardata,
                    conjlayout, symmetriclayout, symmetricdata,
                    triangularlayout, MatLdivVec, hermitianlayout, hermitiandata,
                    materialize!, BlasMatMulMatAdd, BlasMatMulVecAdd, BlasMatLmulVec, BlasMatLdivVec,
                    colsupport, rowsupport, symmetricuplo, MatMulMatAdd, MatMulVecAdd,
                    sublayout, sub_materialize, _fill_lmul!,
                    reflector!, reflectorApply!, _copyto!,
                    _qr!, _qr, _lu!, _lu, _factorize, TridiagonalLayout

import FillArrays: AbstractFill, getindex_value

export BandedMatrix,
       bandrange,
       brand,
       bandwidth,
       BandError,
       band,
       Band,
       BandRange,
       bandwidths,
       colrange,
       rowrange,
       isbanded,
       Zeros,
       Fill,
       Ones,
       Eye

if VERSION < v"1.2-"
    import Base: has_offset_axes
    require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
else
    import Base: require_one_based_indexing
end

if VERSION < v"1.3-"
    const _apply_ipiv_rows! = LinearAlgebra._apply_ipiv!
else
    import LinearAlgebra: _apply_ipiv_rows!
end

include("blas.jl")
include("lapack.jl")

include("generic/AbstractBandedMatrix.jl")
include("generic/broadcast.jl")
include("generic/matmul.jl")
include("generic/Band.jl")
include("generic/utils.jl")
include("generic/indexing.jl")

include("banded/BandedMatrix.jl")
include("banded/BandedLU.jl")
include("banded/bandedqr.jl")
include("banded/gbmm.jl")
include("banded/linalg.jl")

include("symbanded/symbanded.jl")
include("symbanded/ldlt.jl")
include("symbanded/BandedCholesky.jl")
include("symbanded/SplitCholesky.jl")
include("symbanded/bandedeigen.jl")

include("tribanded.jl")

include("interfaceimpl.jl")

# function _precompile_()
#     precompile(Tuple{typeof(gbmm!), Char, Char, Float64, BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}, BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}, Float64, BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}})
# end

# _precompile_()

end #module


text:
@testset "Initialization" begin
    dists = CovidPOMDPs.FitInfectionDistributions(
        CovidPOMDPs.INF_DATAFRAME,
        CovidPOMDPs.INFECTION_HORIZON,
        50
    )
    @test @isinferred first(dists)

    test_probs = [
        CovidPOMDPs.prop_above_LOD(CovidPOMDPs.VIRAL_LOADS, day, CovidPOMDPs.LIMIT_OF_DETECTION)
        for day in 1:CovidPOMDPs.INFECTION_HORIZON
    ]

    @test @isinferred first(test_probs)
    @test all(≥(0.), test_probs)
end

@testset "shifting" begin
    T = rand(Int, 4, 3)
    T′ = circshift(T, (-1,1))
    T′[:,1] .= 0
    T′[end,:] .= 0
    CovidPOMDPs.shift_test!(T)
    @test all(T .== T′)

    I = rand(Int, 14)
    I′ = circshift(I, 1)
    I′[1] = 0
    CovidPOMDPs.shift_inf!(I)
    @test all(I .== I′)
end


text:
#=

    particles.jl

    Implements the type `Particle` used in ABC SMC.

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

=#

"""
    Particle(θ::Vector,d::Number,w::Number)

SMC particle with location `θ`, distance `d`, and weight `w`.
"""
mutable struct Particle
    θ::Vector   # Location
    d::Number   # Discrepancy
    w::Number   # Weight
end
Particle(θ,d) = Particle(θ,d,1.0)

##############################################################
## Extend methods to particles
##############################################################

import Base: sort!
import Distributions: insupport
import StatsBase: sample, weights
import Statistics: cov, mean

weights(P::Vector{Particle})  = weights([p.w for p in P])
distances(P::Vector{Particle}) = [p.d for p in P]
locations(P::Vector{Particle}) = hcat([p.θ for p in P]...)
insupport(d::Distribution,p::Particle) = insupport(d,p.θ)

"""
    sort!(p::Vector{Particle})

Sort a vector of `Particle`s, `P` by distance `d`.
"""
sort!(P::Vector{Particle}) = permute!(P,sortperm(distances(P)))

"""
    cov(p::Vector{Particle})

Weighted covariance matrix of particles `P`.
"""
cov(P::Vector{Particle}) = cov(permutedims(locations(P)),weights(P))

mean(P::Vector{Particle}) = mean(locations(P),weights(P),dims=2)[:]

"""
    sample(prior,dist,N=1)

Sample particles from the `prior` distribution and attach distance using function `dist`.
"""
function StatsBase.sample(prior::Distribution,dist::Function,w::Number=1)
    θ = rand(prior)
    d = dist(θ)
    Particle(θ,d,w)
end
StatsBase.sample(prior::Distribution,dist::Function,n::Int) = [sample(prior,dist,1/n) for i in 1:n]

"""
    copy(p::Particle)

Copies a particle.
"""
Base.copy(p::Particle) = Particle(p.θ,p.d,p.w)

"""
    findmin(p::Particle)

Returns (pᵢ,i) where pᵢ is the particle with the smallest discrepancy.
"""
function Base.findmin(p::Vector{Particle})
    D = [p[i].d for i = 1:length(p)]
    idx = findmin(D)[2]
    p[idx],idx
end

"""
    minimum(p::Particle)

Finds the particle with the smallest discrepancy.
"""
Base.minimum(p::Vector{Particle}) = findmin(p)[1]

"""
    findmax(p::Particle)

Returns (pᵢ,i) where pᵢ is the particle with the largest discrepancy.
"""
function findmax(p::Vector{Particle})
    D = [p[i].d for i = 1:length(p)]
    idx = findmin(-D)[2]
    p[idx],idx
end

"""
    maximum(p::Particle)

Finds the particle with the largest discrepancy.
"""
Base.maximum(p::Vector{Particle}) = findmax(p)[1]

Base.ndims(P::Vector{Particle}) = ndims(P[1])
Base.ndims(p::Particle) = length(p.θ)

ess(P::Vector{Particle}) = 1 / sum(weights(P).^2)

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation