text:
module v2
#Tiny Yolo V2 model configuration

using Knet: Knet, progress, progress!, gpu, KnetArray, relu, minibatch, conv4, pool, softmax

import ..sigmoid, ..Settings

mutable struct v2Chain
    layers
    v2Chain(layers...) = new(layers)
end

mutable struct Conv #Define convolutional layer
    w
    b
    stride
    padding
    f
end
mutable struct YoloPad #Define Yolo padding layer (assymetric padding).
    w
    minibatch_size::Int
end
struct Pool # Define pool layer
    size
    stride
    pad
end

YoloPad(w1::Int, w2::Int, cx::Int, cy::Int, minibatch_size::Int) =
    YoloPad(zeros(Float32, w1, w2, cx, cy), minibatch_size)#Constructor for Yolopad

Conv(w1::Int, w2::Int, cx::Int, cy::Int, st, pd, f) =
    Conv(randn(Float32, w1, w2, cx, cy), randn(Float32, 1, 1, cy, 1), st, pd, f)#Constructor for convolutional layer

#Assymetric padding function
function (y::YoloPad)(x)
    x = reshape(x, y.grid_x + 1, y.grid_y + 1, 1, 512 * y.minibatch_size)
    return reshape(conv4(y.w, x; stride = 1), y.grid_x, y.grid_y, 512, y.minibatch_size)
end

(p::Pool)(x) = pool(x; window = p.size, stride = p.stride, padding = p.pad) #pool function

(c::v2Chain)(x) = (for l in c.layers #chain function
    x = l(x)
end; x)

(c::Conv)(x) =
    c.f.(conv4(c.w, x; stride = c.stride, padding = c.padding) .+ c.b) #convolutional layer function

# Define a custom chain of layers:
struct Chain_layers; layers; Chain_layers(args...)=new(args); end
function (c::Chain_layers)(img)
    res_1_17   = c.layers[1](img)
    res_18_25  = c.layers[2](res_1_17)
    res_26_27  = c.layers[3](res_1_17)
    res_28     = reorg(res_26_27, 2)
    res_29     = concat(res_28, res_18_25)
    res_30_31  = c.layers[4](res_29)
end
(c::Chain_layers)(x,y) = nll(c(x),y)

concat(x,y) = (res = cat(x,y,dims=3))
function reorg(x,stride)
    d1,d2,nr_filt,nr_batch = size(x)
    reshape(x,div(d1,stride),div(d2,stride),nr_filt*(stride^2),nr_batch)
end

#leaky function
function leaky(x)
    return max(Float32(0.1) * x, x)
end

"""
    load(sets::Settings)

Loads the YOLOv2 model.

Note on Asymmetric Padding (see https://github.com/Ybakman/YoloV2#asymmetric-padding)
An alternative solution for asymmetric padding on pooling. We solved the asymmetric padding process in the following 4 steps: Let's assume we have a matrix as (d1,d2,depth,minibatch_size)

1- Apply symmetric padding = 1 with the pooling function

2- reshape it into => (d1,d2,1,minibatch_size*depth)

3- apply the following convolutional layer => 1 0  with stride = 1, padding = 0
                                              0 0

4- Lastly, reshape again into => (d1,d2,depth,minibatch_size)

Now, we made pooling with asymmetric padding on the right and bottom side. If you want to
make it on left and top side, apply this conv layer =>  0 0
                                                        0 1

The last 3 steps are implemented in YoloPad(x) function
The symmetric padding part is implemented in the pooling function just before Yolopad(x)
"""
function load(sets::Settings)
    # comments are line numbers corresponding
    # with https://github.com/pjreddie/darknet/blob/master/cfg/yolov2-voc.cfg
c1_17 = v2Chain(
    Conv(3, 3, 3, 32, 1, 1, leaky), #25 l1
    Pool(2, 2, 0), #33 l2
    Conv(3, 3, 32, 64, 1, 1, leaky), #37 l3
    Pool(2, 2, 0), #45 l4
    Conv(3, 3, 64, 128, 1, 1, leaky), #49 l5
    Conv(1, 1, 128, 64, 1, 0, leaky), #57 l6
    Conv(3, 3, 64, 128, 1, 1, leaky), #65 l7
    Pool(2, 2, 0), #73 l8
    Conv(3, 3, 128, 256, 1, 1, leaky), #77 l9
    Conv(1, 1, 256, 128, 1, 0, leaky), #85 l10
    Conv(3, 3, 128, 256, 1, 1, leaky), #93 l11
    Pool(2, 2, 0), #101 l12
    Conv(3, 3, 256, 512, 1, 1, leaky), #105 l13
    Conv(1, 1, 512, 256, 1, 0, leaky), #113 l14
    Conv(3, 3, 256, 512, 1, 1, leaky), #121 l15
    Conv(1, 1, 512, 256, 1, 0, leaky), #129 l16
    Conv(3, 3, 256, 512, 1, 1, leaky) #137 l17
    )
c18_24 = v2Chain(
	Pool(2, 2, 0), #145 l18
	Conv(3, 3, 512, 1024, 1, 1, leaky), #149 l19
	Conv(1, 1, 1024, 512, 1, 0, leaky), #157 l20
	Conv(3, 3, 512, 1024, 1, 1, leaky), #165 l21
	Conv(1, 1, 1024, 512, 1, 0, leaky), #173 l22
	Conv(3, 3, 512, 1024, 1, 1, leaky), #181 l23
	######
	Conv(3, 3, 1024, 1024, 1, 1, leaky), #192 l24
	Conv(3, 3, 1024, 1024, 1, 1, leaky) #200 l25
	)
c26_27 = v2Chain(
	Conv(1, 1, 512, 64, 1, 0, leaky)
	)

c30_31 = v2Chain(
	Conv(3, 3,1280 , 1024, 1, 1, leaky), #225 l39
	Conv(1, 1, 1024, sets.cell_bboxes * (5 + sets.num_classes), 1, 0, identity), #233 l31
	)

return Chain_layers(c1_17, c18_24, c26_27, c30_31)

end

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

end #module
# ## initial
# return v2Chain(
#     Conv(3, 3, 3, 32, 1, 1, leaky), #25 l1
#     Pool(2, 2, 0), #33 l2
#     Conv(3, 3, 32, 64, 1, 1, leaky), #37 l3
#     Pool(2, 2, 0), #45 l4
#     Conv(3, 3, 64, 128, 1, 1, leaky), #49 l5
#     Conv(1, 1, 128, 64, 1, 1, leaky), #57 l6
#     Conv(3, 3, 64, 128, 1, 1, leaky), #65 l7
#     Pool(2, 2, 0), #73 l8
#     Conv(3, 3, 128, 256, 1, 1, leaky), #77 l9
#     Conv(1, 1, 256, 128, 1, 1, leaky), #85 l10
#     Conv(3, 3, 128, 256, 1, 1, leaky), #93 l11
#     Pool(2, 2, 0), #101 l12
#     Conv(3, 3, 256, 512, 1, 1, leaky), #105 l13
#     Conv(1, 1, 512, 256, 1, 1, leaky), #113 l14
#     Conv(3, 3, 256, 512, 1, 1, leaky), #121 l15
#     Conv(1, 1, 512, 256, 1, 1, leaky), #129 l16
#     Conv(3, 3, 256, 512, 1, 1, leaky), #137 l17
#     Pool(2, 2, 0), #145 l18
#     Conv(3, 3, 512, 1024, 1, 1, leaky), #149 l19
#     Conv(1, 1, 1024, 512, 1, 1, leaky), #157 l20
#     Conv(3, 3, 512, 1024, 1, 1, leaky), #165 l21
#     Conv(1, 1, 1024, 512, 1, 1, leaky), #173 l22
#     Conv(3, 3, 512, 1024, 1, 1, leaky), #181 l23
#     ######
#     Conv(3, 3, 1024, 1024, 1, 1, leaky), #192 l24
#     Conv(3, 3, 1024, 1024, 1, 1, leaky), #200 l25
#     #TODO: [route] #208 https://github.com/pjreddie/darknet/blob/61c9d02ec461e30d55762ec7669d6a1d3c356fb2/cfg/yolov2.cfg#L208 l26
#     Conv(1, 1, 512, 64, 1, 1, leaky), #211 l27
#     #TODO: [reorg] #219 l28
#     #TODO: [route] #222 l29
#     Conv(3, 3,1280 , 1024, 1, 1, leaky), #225 l39
#     Conv(1, 1, 1024, sets.cell_bboxes * (5 + sets.num_classes), 1, 0, identity), #233 l31
# )


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@testset "basic path functions for string type $S" for S in (String, GenericString)
    dir = pwd()
    sep = Base.Filesystem.path_separator
    @testset "isabspath and abspath" begin
        @test abspath(S("foo")) == "$dir$(sep)foo"
        @test abspath(S("foo"), S("bar")) == "$dir$(sep)foo$(sep)bar"
        @test isabspath(S(homedir()))
        @test !isabspath(S("foo"))
    end
    @test basename(S("foo$(sep)bar")) == "bar"
    @test dirname(S("foo$(sep)bar")) == "foo"

    @testset "expanduser" begin
        @test expanduser(S("")) == ""
        @test expanduser(S("x")) == "x"
        @test expanduser(S("~")) == (Sys.iswindows() ? "~" : homedir())
    end
    @testset "Base.contractuser" begin
        @test Base.contractuser(S(homedir())) == (Sys.iswindows() ? homedir() : "~")
        @test Base.contractuser(S(joinpath(homedir(), "x"))) ==
              (Sys.iswindows() ? joinpath(homedir(), "x") : "~$(sep)x")
        @test Base.contractuser(S("/foo/bar")) == "/foo/bar"
    end
    @testset "isdirpath" begin
        @test !isdirpath(S("foo"))
        @test isdirpath(S("foo$sep"))
        @test isdirpath(S(""))
        @test isdirpath(S("."))
        @test isdirpath(S(".."))
    end
    @testset "joinpath" begin
        @test joinpath(S("")) == ""
        @test joinpath(S("foo")) == "foo"
        @test joinpath(S("foo"), S("bar")) == "foo$(sep)bar"
        @test joinpath(S("foo"), S("bar"), S("baz")) == "foo$(sep)bar$(sep)baz"
        @test joinpath(S("foo"), S(""), S("baz")) == "foo$(sep)baz"
        @test joinpath(S("foo"), S(""), S("")) == "foo$(sep)"
        @test joinpath(S("foo"), S(""), S(""), S("bar")) == "foo$(sep)bar"

        @test joinpath(S("foo"), S(homedir())) == homedir()
        @test joinpath(S(abspath("foo")), S(homedir())) == homedir()

        if Sys.iswindows()
            @test joinpath(S("foo"),S("bar:baz")) == "bar:baz"
            @test joinpath(S("C:"),S("foo"),S("D:"),S("bar")) == "D:bar"
            @test joinpath(S("C:"),S("foo"),S("D:bar"),S("baz")) == "D:bar$(sep)baz"

            # relative folders and case-insensitive drive letters
            @test joinpath(S("C:\\a\\b"), S("c:c\\e")) == "C:\\a\\b\\c\\e"

            # UNC paths
            @test joinpath(S("\\\\server"), S("share")) == "\\\\server\\share"
            @test joinpath(S("\\\\server"), S("share"), S("a")) == "\\\\server\\share\\a"
            @test joinpath(S("\\\\server\\"), S("share"), S("a")) == "\\\\server\\share\\a"
            @test joinpath(S("\\\\server"), S("share"), S("a"), S("b")) == "\\\\server\\share\\a\\b"
            @test joinpath(S("\\\\server\\share"),S("a")) == "\\\\server\\share\\a"
            @test joinpath(S("\\\\server\\share\\"), S("a")) == "\\\\server\\share\\a"

        elseif Sys.isunix()
            @test joinpath(S("foo"),S("bar:baz")) == "foo$(sep)bar:baz"
            @test joinpath(S("C:"),S("foo"),S("D:"),S("bar")) == "C:$(sep)foo$(sep)D:$(sep)bar"
            @test joinpath(S("C:"),S("foo"),S("D:"),S("bar"),S("baz")) == "C:$(sep)foo$(sep)D:$(sep)bar$(sep)baz"
        end
    end
    @testset "normpath" begin
        @test normpath(S(joinpath("."))) == "."
        @test normpath(S(joinpath(".."))) == ".."
        @test normpath(S(joinpath("..","."))) == ".."
        @test normpath(S(joinpath(".",".."))) == ".."
        @test normpath(S(joinpath("..",".."))) == "..$(sep).."
        @test normpath(S(joinpath(".","..",".."))) == "..$(sep).."
        @test normpath(S(joinpath("..",".",".."))) == "..$(sep).."
        @test normpath(S(joinpath("..","..","."))) == "..$(sep).."

        @test normpath(S(joinpath("foo","."))) == "foo$sep"
        @test normpath(S(joinpath("foo",".."))) == "."
        @test normpath(S(joinpath("foo","..","."))) == "."
        @test normpath(S(joinpath("foo",".",".."))) == "."
        @test normpath(S(joinpath("foo","..",".."))) == ".."
        @test normpath(S(joinpath("foo",".","..",".."))) == ".."
        @test normpath(S(joinpath("foo","..",".",".."))) == ".."
        @test normpath(S(joinpath("foo","..","..","."))) == ".."

        @test normpath(S(joinpath(".","bar"))) == "bar"
        @test normpath(S(joinpath("..","bar"))) == "..$(sep)bar"
        @test normpath(S(joinpath("..",".","bar"))) == "..$(sep)bar"
        @test normpath(S(joinpath(".","..","bar"))) == "..$(sep)bar"
        @test normpath(S(joinpath("..","..","bar"))) == "..$(sep)..$(sep)bar"
        @test normpath(S(joinpath(".","..","..","bar"))) == "..$(sep)..$(sep)bar"
        @test normpath(S(joinpath("..",".","..","bar"))) == "..$(sep)..$(sep)bar"
        @test normpath(S(joinpath("..","..",".","bar"))) == "..$(sep)..$(sep)bar"

        @test normpath(S(joinpath("foo",".","bar"))) == "foo$(sep)bar"
        @test normpath(S(joinpath("foo","..","bar"))) == "bar"
        @test normpath(S(joinpath("foo","..",".","bar"))) == "bar"
        @test normpath(S(joinpath("foo",".","..","bar"))) == "bar"
        @test normpath(S(joinpath("foo","..","..","bar"))) == "..$(sep)bar"
        @test normpath(S(joinpath("foo",".","..","..","bar"))) == "..$(sep)bar"
        @test normpath(S(joinpath("foo","..",".","..","bar"))) == "..$(sep)bar"
        @test normpath(S(joinpath("foo","..","..",".","bar"))) == "..$(sep)bar"
    end
    @test relpath(S(joinpath("foo","bar")), S("foo")) == "bar"

    @testset "splitpath" begin
        @test splitpath(S(joinpath("a","b","c"))) == ["a", "b", "c"]
        @test splitpath("") == [""]

        @test splitpath(S(joinpath("cats are", "gr8t"))) == ["cats are", "gr8t"]
        @test splitpath(S(joinpath("  ", " "))) == ["  ", " "]

        # Unix-style paths are understood by all systems.
        @test splitpath(S("/a/b")) == ["/", "a", "b"]
        @test splitpath(S("/")) == ["/"]
        @test splitpath(S("a/")) == ["a"]
        @test splitpath(S("a/b/")) == ["a", "b"]
        @test splitpath(S("a.dir/b.txt")) == ["a.dir", "b.txt"]
        @test splitpath(S("///")) == ["/"]
        @test splitpath(S("///a///b///")) == ["/", "a", "b"]

        if Sys.iswindows()
            @test splitpath(S("C:\\\\a\\b\\c")) == ["C:\\", "a", "b", "c"]
            @test splitpath(S("C:\\\\")) == ["C:\\"]
            @test splitpath(S("J:\\")) == ["J:\\"]
            @test splitpath(S("C:")) == ["C:"]
            @test splitpath(S("C:a")) == ["C:a"]
            @test splitpath(S("C:a\\b")) == ["C:a", "b"]

            @test splitpath(S("a\\")) == ["a"]
            @test splitpath(S("a\\\\b\\\\")) == ["a","b"]
            @test splitpath(S("a.dir\\b.txt")) == ["a.dir", "b.txt"]
            @test splitpath(S("\\a\\b\\")) == ["\\", "a","b"]
            @test splitpath(S("\\\\a\\b")) == ["\\\\a\\b"]  # This is actually a valid drive name in windows.

            @test splitpath(S("/a/b\\c/d\\\\e")) == ["/", "a", "b", "c", "d", "e"]
            @test splitpath(S("/\\/\\")) == ["/"]
            @test splitpath(S("\\/\\a/\\//b")) == ["\\","a","b"]
        end
    end

    @testset "splitdir, splitdrive" begin
        @test joinpath(splitdir(S(homedir()))...) == homedir()
        @test string(splitdrive(S(homedir()))...) == homedir()
        @test splitdrive("a\nb") == ("", "a\nb")

        if Sys.iswindows()
            @test splitdrive(S("\\\\servername\\hello.world\\filename.ext")) ==
                ("\\\\servername\\hello.world","\\filename.ext")
            @test splitdrive(S("\\\\servername.com\\hello.world\\filename.ext")) ==
                ("\\\\servername.com\\hello.world","\\filename.ext")
            @test splitdrive(S("C:\\foo\\bar")) ==
                ("C:","\\foo\\bar")
        end

        @test splitext(S("")) == ("", "")
        @test splitext(S(".")) == (".", "")
        @test_broken splitext(S("..")) == ("..", "")
        @test_broken splitext(S("...")) == ("...", "")
        @test splitext(S("foo")) == ("foo", "")
        @test splitext(S("foo.")) == ("foo", ".")
        @test_broken splitext(S("foo..")) == ("foo", "..")
        @test_broken splitext(S("foo...")) == ("foo", "...")
        @test splitext(S("foo.bar")) == ("foo", ".bar")
        @test splitext(S(".foo")) == (".foo", "")
        @test splitext(S(".foo.")) == (".foo", ".")
        @test_broken splitext(S(".foo..")) == (".foo", "..")
        @test_broken splitext(S(".foo...")) == (".foo", "...")
        @test splitext(S(".foo.bar")) == (".foo", ".bar")
    end

    @testset "isabspath" begin
        @test isabspath(S("~")) == false
        @test isabspath(S("/")) == true # on windows, this is relatively absolute
        @test isabspath(S("A:/")) == Sys.iswindows()
        @test isabspath(S("B:\\")) == Sys.iswindows()
        @test isabspath(S("./")) == false
        @test isabspath(S("C:")) == false
        @test isabspath(S("C:.")) == false
        @test isabspath(S("α:/")) == false
        @test isabspath(S(".:/")) == false
        #@test isabspath(S("_:/")) == false # FIXME?
        #@test isabspath(S("AB:/")) == false # FIXME?
        @test isabspath(S("\\\\")) == Sys.iswindows()
        if Sys.isunix()
            @test isabspath(S(expanduser("~"))) == true
            @test startswith(expanduser(S("~")), homedir())
        else
            @test expanduser(S("~")) == "~"
        end
    end

    @testset "relpath" begin
        function test_relpath()
            sep = Base.Filesystem.path_separator
            filepaths = [
                "$(sep)home$(sep)user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)Test1.md",
                "$(sep)home$(sep)user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "$(sep)home$(sep)user$(sep).julia$(sep)测试2$(sep)docs$(sep)api$(sep)测试2.md",
                "$(sep)home$(sep)user$(sep)dir_withendsep$(sep)",
                "$(sep)home$(sep)dir2_withendsep$(sep)",
                "$(sep)home$(sep)test.md",
                "$(sep)home",
                # Special cases
                "$(sep)",
                "$(sep)home$(sep)$(sep)$(sep)"
            ]
            startpaths = [
                "$(sep)home$(sep)user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)genindex.md",
                "$(sep)multi_docs$(sep)genindex.md",
                "$(sep)home$(sep)user$(sep)dir_withendsep$(sep)",
                "$(sep)home$(sep)dir2_withendsep$(sep)",
                "$(sep)home$(sep)test.md",
                "$(sep)home",
                # Special cases
                "$(sep)",
                "$(sep)home$(sep)$(sep)$(sep)"
            ]
            relpath_expected_results = [
                "..$(sep)Test1.md",
                "..$(sep)..$(sep)home$(sep)user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)Test1.md",
                "..$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)Test1.md",
                "..$(sep)user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)Test1.md",
                "..$(sep)user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)Test1.md",
                "user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)Test1.md",
                "home$(sep)user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)Test1.md",
                "user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)Test1.md",
                "..$(sep)lib$(sep)file1.md",
                "..$(sep)..$(sep)home$(sep)user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "..$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "..$(sep)user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "..$(sep)user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "home$(sep)user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "user$(sep).julia$(sep)Test1$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "..$(sep)..$(sep)..$(sep)..$(sep)测试2$(sep)docs$(sep)api$(sep)测试2.md",
                "..$(sep)..$(sep)home$(sep)user$(sep).julia$(sep)测试2$(sep)docs$(sep)api$(sep)测试2.md",
                "..$(sep).julia$(sep)测试2$(sep)docs$(sep)api$(sep)测试2.md",
                "..$(sep)user$(sep).julia$(sep)测试2$(sep)docs$(sep)api$(sep)测试2.md",
                "..$(sep)user$(sep).julia$(sep)测试2$(sep)docs$(sep)api$(sep)测试2.md",
                "user$(sep).julia$(sep)测试2$(sep)docs$(sep)api$(sep)测试2.md",
                "home$(sep)user$(sep).julia$(sep)测试2$(sep)docs$(sep)api$(sep)测试2.md",
                "user$(sep).julia$(sep)测试2$(sep)docs$(sep)api$(sep)测试2.md",
                "..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)dir_withendsep",
                "..$(sep)..$(sep)home$(sep)user$(sep)dir_withendsep",".","..$(sep)user$(sep)dir_withendsep",
                "..$(sep)user$(sep)dir_withendsep","user$(sep)dir_withendsep",
                "home$(sep)user$(sep)dir_withendsep","user$(sep)dir_withendsep",
                "..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)dir2_withendsep",
                "..$(sep)..$(sep)home$(sep)dir2_withendsep","..$(sep)..$(sep)dir2_withendsep",".",
                "..$(sep)dir2_withendsep","dir2_withendsep","home$(sep)dir2_withendsep","dir2_withendsep",
                "..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)test.md","..$(sep)..$(sep)home$(sep)test.md",
                "..$(sep)..$(sep)test.md","..$(sep)test.md",".","test.md","home$(sep)test.md","test.md",
                "..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..","..$(sep)..$(sep)home","..$(sep)..",
                "..","..",".","home",".","..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..","..$(sep)..",
                "..$(sep)..$(sep)..","..$(sep)..","..$(sep)..","..",".","..",
                "..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..","..$(sep)..$(sep)home","..$(sep)..",
                "..","..",".","home","."
            ]
            idx = 0
            for filep in filepaths
                for startp in startpaths
                    res = relpath(filep, startp)
                    idx += 1
                    @test res == relpath_expected_results[idx]
                    if Sys.iswindows()
                        @test relpath("e:$filep", "e:$startp") == relpath_expected_results[idx]
                        @test relpath("e:$filep", "E:$startp") == relpath_expected_results[idx]
                        @test relpath("E:$filep", "e:$startp") == relpath_expected_results[idx]
                        @test relpath("E:$filep", "E:$startp") == relpath_expected_results[idx]
                    end
                end
            end
            # Additional cases
            @test_throws ArgumentError relpath(S("$(sep)home$(sep)user$(sep)dir_withendsep$(sep)"), "")
            @test_throws ArgumentError relpath(S(""), S("$(sep)home$(sep)user$(sep)dir_withendsep$(sep)"))
        end
        test_relpath()
    end

    if Sys.iswindows()
        @testset "issue #23646" begin
            @test lowercase(relpath("E:\\a\\b", "C:\\c")) == "e:\\a\\b"
            @test lowercase(relpath("E:\\a\\b", "c:\\c")) == "e:\\a\\b"
            @test lowercase(relpath("e:\\a\\b", "C:\\c")) == "e:\\a\\b"
            @test lowercase(relpath("e:\\a\\b", "c:\\c")) == "e:\\a\\b"

            @test relpath("C:\\a\\b", "c:\\a\\b") == "."
            @test relpath("c:\\a\\b", "C:\\a\\b") == "."
            @test lowercase(relpath("C:\\a\\b", "c:\\c\\d")) == "..\\..\\a\\b"
            @test lowercase(relpath("c:\\a\\b", "C:\\c\\d")) == "..\\..\\a\\b"
        end
    end

    @testset "type stability" begin
        @test isa(joinpath(S("a"), S("b")), String)
        @test isa(joinpath(S(abspath("a")), S("b")), String)
    end
end

@testset "homedir" begin
    var = Sys.iswindows() ? "USERPROFILE" : "HOME"
    AVG_PATH = Base.Filesystem.AVG_PATH - 1 # null-termination character
    for i = 0:9
        local home = " "^AVG_PATH * "123456789"[1:i]
        @test withenv(var => home) do
            homedir()
        end == home
    end
    @test isabspath(withenv(homedir, var => nothing))
end

ext: jl

max_stars_repo_name: udesou/julia

max_issues_repo_name: udesou/julia

max_forks_repo_name: udesou/julia


text:
#---------------------------------------------------------------------#
#This code computes the Chebyshev points & weights
#Written by F.X. Giraldo on 7/2014
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function chebyshev_gauss(P::Integer)

    #Initialize arrays
    xgl=zeros(P)
    wgl=zeros(P)

    #Construct arrays
    for i=1:P
        xgl[P+1-i]=cos( (2*i-1)*pi/(2*P) );
        wgl[P+1-i]=pi/P
    end
    return (xgl,wgl)
end #function


text:
## Ray Tracing Demonstrations

## Preamble
using Plots

include("../AcousticPropagation.jl")

## Demonstrative Scenario
# Altimetry
zAtiMin = -10
zAtiMax = 50
zAti(r) = zAtiMin + (zAtiMax - zAtiMin)*(sin(r/1e3) + 1.)/2
ati = AcousticPropagation.Boundary(zAti)

# Bathymetry
rPeak = 5e3
rMax = 10e3
zMax = 1e3
zMin = 8e2
Aᵣ = (2rPeak/3)^2/log((9zMax - 11zMin)/(10(zMax - zMin)))
zBty(r) = zMax - (zMax - zMin)*exp(-(r - rPeak)^2/4e5)
bty = AcousticPropagation.Boundary(zBty)

# Ocean
cMin = 1500
cMax = 1600
C(r) = [1 zAti(r) zAti(r)^2
	1 (zAti(r) + zBty(r))/2 ((zAti(r) + zBty(r))/2)^2
	1 zBty(r) zBty(r)^2]
c_(r) = C(r)\[cMax, cMin, cMax]
c₀(r) = c_(r)[1]
c₁(r) = c_(r)[2]
c₂(r) = c_(r)[3]
c(r, z) = c₀(r) + c₁(r)*z + c₂(r)*z^2
ocn = AcousticPropagation.Medium(c, rMax)

# Initial Conditions
r₀ = 0.
z₀ = (zBty(r₀) + zAti(r₀))/2
src = AcousticPropagation.Entity.(r₀, z₀)

θ₀ = acos(c(r₀, z₀)/cMax).*(-1.5:0.125:1.5)
# θ₀ = acos(c(r₀, z₀)/cMax).*[-1.5, 1.5]
rays = AcousticPropagation.Ray.(θ₀, src, ocn, bty, ati)

println("Now plotting.")
pt = plot(
	xaxis = "Range (m)",
	yaxis = ("Depth (m)", :flip))
r = range(0, rMax, length = 100)
# z = range(zAtiMin, zMax, length = 51)
plot!(r, zAti, label = "Altimetry")
plot!(r, zBty, label = "Bathymetry")
for nRay = 1:length(rays)
	plot!(rays[nRay].Sol, vars = (1, 2), label = "")
end
display(pt)

savefig(pt, "img/RayTrace_FirstExample.png")

## Parabolic Bathymetry
cVal = 250
R = 20e3
ocn = AcousticPropagation.Medium((r, z) -> cVal, R)
ati = AcousticPropagation.Boundary(r -> 0)
bty = AcousticPropagation.Boundary(r -> 2e-3*2.5e5sqrt(1 + r/cVal))
src = AcousticPropagation.Entity(0, 0)
θ₀ = range(atan(5e3/2e3), atan(5e3/20e3), length = 30)

rays = AcousticPropagation.Ray.(θ₀, src, ocn, bty, ati)

pt = plot(
	xaxis = "Range (m)",
	yaxis = ("Depth(m)", :flip))
plot!(range(0, R, length = 101), bty.z, label = "Bathymetry")
# plot!(ray.Sol, vars = (1, 2))
for nRay = 1:length(rays)
	plot!(rays[nRay].Sol, vars = (1, 2), label = "")
end
display(pt)

savefig(pt, "img/RayTrace_ParabolicBoundary.png")

## Uniformly Increasing Celerity
ati = AcousticPropagation.Boundary(0)
bty = AcousticPropagation.Boundary(5e3)
ocn = AcousticPropagation.Medium((r, z) -> 1500 + 100z/5e3, 1e5)
src = AcousticPropagation.Entity(0, 0)
θ_crit = acos(ocn.c(0, 0)/ocn.c(0, 5e3))
# ocn.c(0, 0)
# ocn.c(0, 5e3)
# θ_crit = acos(1500. /1600.)
θ₀ = range(0.1θ_crit, θ_crit, length = 10)
# θ₀ = range(0.5θ_crit, θ_crit, length = 2)
rays = AcousticPropagation.Ray.(θ₀, src, ocn, bty, ati)

pt = plot(
	xaxis = "Range (m)",
	yaxis = ("Depth (m)", :flip))
for nRay = 1:length(rays)
	plot!(rays[nRay].Sol, vars = (1, 2), label = "")
end
display(pt)

savefig(pt, "img/RayTrace_UpwardRefracting.png")

## Convergence Zones
cVec = [1520, 1500, 1515, 1495, 1545.]
zVec = [0., 300., 1200., 2e3, 5000.]

# cMat = vcat([0 0 250e3], hcat(zVec, cVec, cVec))

ati = AcousticPropagation.Boundary(0)
bty = AcousticPropagation.Boundary(5e3)
# ocn = AcousticPropagation.Medium(cMat, 250e3)
ocn = AcousticPropagation.Medium(zVec, cVec, 250e3)
src = AcousticPropagation.Entity(0, 0)
θ_crit = acos(ocn.c(0, 0)/ocn.c(0, 5e3))
θ₀ = range(0.01θ_crit, θ_crit, length = 10)

rays = AcousticPropagation.Ray.(θ₀, src, ocn, bty, ati)

pt = plot(
	xaxis = "Range (m)",
	yaxis = ("Depth (m)", :flip))
for nRay = 1:length(rays)
	plot!(rays[nRay].Sol, vars = (1, 2), label = "")
end
display(pt)

savefig(pt, "img/RayTrace_ConvergenceZones.png")


text:
# https://github.com/stan-dev/blob/master/basic_estimators/bernoulli.stan

const berstanmodel = "
data {
  int<lower=0> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  for (n in 1:N)
    y[n] ~ bernoulli(theta);
}
"


text:
function QUBO_linear_sum(n, σ)
    N = n^2
    Q = zeros(N,N)
    for i in 1:n, k in 1:n
        a = (i - 1) * n + k
        Q[a, a] = -2 * σ * (k - 1)
        for j in 1:n, l in 1:n
            b = (j - 1) * n + l
            Q[a, b] += (k - 1) * (l - 1)
        end
    end
    return Q
end


text:
# Plots a 1D or 2D GPE object.
#
# Keyword arguments:
#   - β: level of confidence interval for ribbon (1D only)
#   - obsv: plot observation points (1D only)
#   - var: plot predicted variance rather than mean (2D only)
@recipe function f(gp::GPE; β=0.95, obsv=true, var=false)
    @assert gp.dim ∈ (1,2)
    if gp.dim == 1
        xlims --> (minimum(gp.x), maximum(gp.x))
        xmin, xmax = plotattributes[:xlims]
        x = range(xmin, stop=xmax, length=100)
        μ, Σ = predict_f(gp, x)
        y = μ
        err = invΦ((1+β)/2)*sqrt.(Σ)

        @series begin
            seriestype := :path
            ribbon := err
            fillcolor --> :lightblue
            color --> :black
            x,y
        end
        if obsv
            @series begin
                seriestype := :scatter
                markershape := :circle
                markercolor := :black
                gp.x', gp.y
            end
        end
    else
        xlims --> (minimum(gp.x[1,:]), maximum(gp.x[1,:]))
        ylims --> (minimum(gp.x[2,:]), maximum(gp.x[2,:]))
        xmin, xmax = plotattributes[:xlims]
        ymin, ymax = plotattributes[:ylims]
        x = range(xmin, stop=xmax, length=50)
        y = range(ymin, stop=ymax, length=50)
        xgrid = repeat(x', 50, 1)
        ygrid = repeat(y, 1, 50)
        μ, Σ = predict_f(gp,[vec(xgrid)';vec(ygrid)'])
        if var
            zgrid  = reshape(Σ,50,50)
        else
            zgrid  = reshape(μ,50,50)
        end
        x, y, zgrid
    end
end


text:
using Plots, Test

@testset "Shorthands" begin
    @testset "Set Lims" begin
        p = plot(rand(10))

        xlims!((1, 20))
        @test xlims(p) == (1, 20)

        ylims!((-1, 1))
        @test ylims(p) == (-1, 1)

        zlims!((-1, 1))
        @test zlims(p) == (-1, 1)

        xlims!(-1, 11)
        @test xlims(p) == (-1, 11)

        ylims!((-10, 10))
        @test ylims(p) == (-10, 10)

        zlims!((-10, 10))
        @test zlims(p) == (-10, 10)
    end

    @testset "Set Ticks" begin
        p = plot([0, 2, 3, 4, 5, 6, 7, 8, 9, 10])

        xticks = 2:6
        xticks!(xticks)
        @test Plots.get_subplot(current(), 1).attr[:xaxis][:ticks] == xticks

        yticks = 0.2:0.1:0.7
        yticks!(yticks)
        @test Plots.get_subplot(current(), 1).attr[:yaxis][:ticks] == yticks

        xticks = [5, 6, 7.5]
        xlabels = ["a", "b", "c"]

        xticks!(xticks, xlabels)
        @test Plots.get_subplot(current(), 1).attr[:xaxis][:ticks] == (xticks, xlabels)

        yticks = [0.5, 0.6, 0.75]
        ylabels = ["z", "y", "x"]
        yticks!(yticks, ylabels)
        @test Plots.get_subplot(current(), 1).attr[:yaxis][:ticks] == (yticks, ylabels)
    end
end


text:
mutable struct EpisodeLength <: LoggerKeyData
    episode_length::Array{Int64}
    step_counter::Int

    function EpisodeLength(logger_init_info)
        new([],0)
    end
end

function lg_step!(self::EpisodeLength, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    self.step_counter += 1
    if is_terminal == true
        push!(self.episode_length, self.step_counter)
        self.step_counter = 0
    # else
    end

end

function lg_episode_end!(self::EpisodeLength, cur_step_in_episode, cur_step_total)
    # Environment must have cut off trajectory for exceeding max length
    if self.step_counter > 0
        # subtract a step since the counting in step is forward
        push!(self.episode_length, self.step_counter)
        self.step_counter = 0
    end
end

function save_log(self::EpisodeLength, save_dict::Dict)
    save_dict[:episode_length] = self.episode_length
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
module DifferentialGeometry

using Match
using LinearAlgebra, Polynomials, Statistics
using MiniQhull

include("io.jl")
using .DataIO: read_obj, read_ply

export Mesh
export mesh

export Surface
export basis, pullback

export Manifold
export manifold, makefield

numdims(x) = sum(1 for d in 1:ndims(x) if size(x,d) > 1)

# ------------------------------------------------------------------------
# mesh type

"""
    struct Mesh{T <: Real}
        r   :: Array{T, 2}
        vₙ  :: Array{T, 2}
        tri :: Array{Int, 2}
    end

Store a 2 dimensional triangular mesh, embedded into arbitrary dimensions.
`r` and `vₙ` denote vertex positions and normals respectively.
`tri` denotes (non-oriented) triangular faces.
"""
struct Mesh{T <: Real}
    r   :: Array{T, 2}     # vertex positions
    vₙ  :: Array{T, 2}     # vertex normals
    tri :: Array{Int, 2}   # triangular faces
end

"""
    mesh(io::IO, type::Symbol=:obj)

Load a `Mesh` object from stream `io` formatted with `type`.
As of now, only `.obj` files are supported.
"""
function mesh(io::IO, type::Symbol=:obj)
    if type != :obj
        error("not implemented")
    end

    data = read_obj(io)

    return Mesh{eltype(data.verts[1])}(data.verts, data.norms, data.faces)
end

# ------------------------------------------------------------------------
# surface type

# internal functions
"""
    ellipse(r)

Fit an ellipse to 2D point cloud `r`.
`r` is assumed to be sized ``N \\times 2``.
"""
function ellipse(r)
    # constrained least squares (constraint: 4AC-B² = 1)
    D = hcat(r[:,1].*r[:,1], r[:,1].*r[:,2], r[:,2].*r[:,2], r[:,1], r[:,2], ones(size(r, 1),1))
    S = D'*D

    C      = zeros(6,6)
    C[1,3] = C[3,1] = 2
    C[2,2] = -1

    λ, ν = eigen(inv(S) * C)
    n    = argmax(λ)
    p    = real(ν[:,n])

    # equations taken from wikipedia
    A, B, C, D, E, F = p[1], p[2]/2, p[3], p[4]/2, p[5]/2, p[6]
    num = B^2 - A*C

    x₀ = (C*D-B*E)/num
    y₀ = (A*E-B*D)/num
    ϕ  = 0.5*atan(2*B/(A-C))

    up    = 2*(A*E^2+C*D^2+F*B^2-2*B*D*E-A*C*F)
    down1 = (B^2-A*C)*( (C-A)*√(1+4*B^2/((A-C)*(A-C)))-(C+A))
    down2 = (B^2-A*C)*( (A-C)*√(1+4*B^2/((A-C)*(A-C)))-(C+A))
    a, b  = √(abs(up/down1)), √(abs(up/down2))

    return [x₀, y₀, a, b, ϕ]
end

"""
    struct Surface{T <: Real}
        Θ  :: Array{T, 2}
        x  :: Array{T}
        L  :: T
        Λ  :: Array{Polynomial{T}}
        ∂Λ :: Array{Polynomial{T}}
    end

Store a representation of a 2D surface embedded into higher dimensional Euclidean space.
Fits the surface by:
  1. Partition the x axis so that each resultant interval has, on average, 20 points.
  2. Fit points within each partition to a 2D ellipse.
  3. Fit a polynomial function to each elliptical parameter over all partitions.
  4. Use polynomials to estimate tangent vectors.
`Θ` denotes the elliptical parameters fit per partition.
`x` denotes the input data.
`L` denotes the length along the x axis.
`Λ` denotes the polynomials for each elliptical parameter.
`∂Λ` denotes the derivative of polynomials for each elliptical parameter.
"""
struct Surface{T <: Real}
    Θ  :: Array{T, 2}
    x  :: Array{T}
    L  :: T
    Λ  :: Array{Polynomial{T}}
    ∂Λ :: Array{Polynomial{T}}

    # constructors
    Surface(Θ, x, L, Λ, ∂Λ) = new{eltype(x)}(Θ, x, L, Λ, ∂Λ)

    function Surface(x; order=12, debug=false)
        nbin = size(x,1)÷20
        xbin = LinRange(minimum(x[:,1]), maximum(x[:,1]), nbin+1)

        select(x, index) = hcat([x[index,d] for d in 2:size(x,2)]...)

        T  = eltype(x)
        x₀ = T[]
        Θ  = Array{T}[]
        for (lb, ub) in zip(xbin[1:end-1],xbin[2:end])
            ι = (x[:,1] .>= lb) .& (x[:,1] .<= ub)
            if sum(ι) < 6
                continue
            end

            push!(x₀, mean(x[ι,1]))
            push!(Θ,  ellipse(select(x, ι)))
        end

        Θ = hcat(Θ...)' # [x₀, y₀, a, b, ϕ]
        Λ = Array{Polynomial}(undef, size(Θ,2))
        L = (maximum(x₀) - minimum(x₀))./2
        for p in 1:size(Θ,2)
            M    = hcat([(x₀./L).^n for n in 0:order]...)
            α    = M \ Θ[:,p]
            Λ[p] = Polynomial(α)
            if debug
                χ = range(minimum(x₀), maximum(x₀), length=1000)
                Plots.scatter(x₀, Θ[:,p])
                Plots.plot!(χ, Λ[p].(χ./L), ylims=(minimum(Θ[:,p]),maximum(Θ[:,p])))
                gui()
                readline()
            end
        end

        ∂Λ = [derivative(Λ[i]) for i in 1:length(Λ)]

        new{T}(Θ, x₀, L, Λ, ∂Λ)
    end
end

"""
    function pullback(s::Surface, r)

Transform Euclidean point cloud `r` into a 2D cylindrical projection defined by `s`.
Cylindrical coordinates are [x,φ].
Assumes all points within `r` are distributed over the surface.
Assumes `r` is sized ``N \\times 3``.
"""
function pullback(s::Surface, r)
    z₀ = r[:,1] ./ s.L
    x₀ = s.Λ[1].(z₀)
    y₀ = s.Λ[2].(z₀)
    a  = s.Λ[3].(z₀)
    b  = s.Λ[4].(z₀)
    ϕ  = s.Λ[5].(z₀)

    v = Array{eltype(r),2}(undef, length(a), 2)
    for (i, p) in enumerate(ϕ)
        v[i,:] = [cos(p) +sin(p); -sin(p) cos(p)]*(r[i,2:3] .- [x₀[i], y₀[i]])
    end
    θ = mod2pi.(atan.(v[:,2]./b, v[:,1]./a) .+ π .- 1.65) .- π

    v[:,1] = z₀
    v[:,2] = θ

    return v
end

"""
    function basis(s::Surface, r)

Compute the tangent vectors ``\\hat{\\bm{e}}_\\phi, \\hat{\\bm{e}}_x`` associated to each point within `r`.
Assumes all points within `r` are distributed over the surface.
Assumes `r` is sized ``N \\times 3``.
"""
function basis(s::Surface, r; surface=false, normalized=true)
    x  = surface ? r : pullback(s, r)

    # function
    x₀ = s.Λ[1].(x[:,1])
    y₀ = s.Λ[2].(x[:,1])
    a  = s.Λ[3].(x[:,1])
    b  = s.Λ[4].(x[:,1])
    ϕ  = s.Λ[5].(x[:,1])

    # derivatives
    ∂x₀ = s.∂Λ[1].(x[:,1])./s.L
    ∂y₀ = s.∂Λ[2].(x[:,1])./s.L
    ∂a  = s.∂Λ[3].(x[:,1])./s.L
    ∂b  = s.∂Λ[4].(x[:,1])./s.L
    ∂ϕ  = s.∂Λ[5].(x[:,1])./s.L

    eθ = hcat(zeros(length(a)),
              -a.*cos.(ϕ).*sin.(x[:,2]) .- b.*sin.(ϕ).*cos.(x[:,2]) .+ x₀,
              -a.*sin.(ϕ).*sin.(x[:,2]) .+ b.*cos.(ϕ).*cos.(x[:,2]) .+ y₀)

    ez = hcat(ones(length(a)),
              +(∂a.*cos.(ϕ).-a.*sin.(ϕ).*∂ϕ).*cos.(x[:,2]) .- (∂b.*sin.(ϕ).+b.*cos.(ϕ).*∂ϕ).*sin.(x[:,2]) .+ ∂x₀,
              +(∂a.*sin.(ϕ).+a.*cos.(ϕ).*∂ϕ).*cos.(x[:,2]) .+ (∂b.*cos.(ϕ).-b.*sin.(ϕ).*∂ϕ).*sin.(x[:,2]) .+ ∂y₀)

    if normalized
        ez = ez ./ sqrt.(ez⋅ez)
        eθ = eθ ./ sqrt.(eθ⋅eθ)
    end

    return ez, eθ
end

# ------------------------------------------------------------------------
# manifold type

"""
    struct Manifold{T <: Real}
        mesh :: Mesh{T}
        surf :: Surface{T}
    end

Store the representation of a differential geometry object.
`mesh` is the empirical point cloud.
`surf` is the estimated differentiable surface.
"""
struct Manifold{T <: Real}
    mesh :: Mesh{T}
    surf :: Surface{T}
end

"""
    function pullback(ℳ::Manifold)

Returns the 2D cylindrical projection of the mesh, as estimated by the internal surface.
"""
function pullback(ℳ::Manifold)
    x    = pullback(ℳ.surf, ℳ.mesh.r)
    maxₗ = maximum(hcat(abs.(x[ℳ.mesh.tri[:,1],2]-x[ℳ.mesh.tri[:,2],2]),
                        abs.(x[ℳ.mesh.tri[:,2],2]-x[ℳ.mesh.tri[:,3],2]),
                        abs.(x[ℳ.mesh.tri[:,3],2]-x[ℳ.mesh.tri[:,1],2])),
                   dims=2)
    return x, ℳ.mesh.tri[getindex.(findall(maxₗ .> π/2),1), :]
end

pullback(ℳ::Manifold, r) = pullback(ℳ.surf, r)

"""
    function order(tri, r)

Reorder the triangulation so that all labels are counterclockwise.
"""
function order(tri, r)
    r21 = r[tri[:,2],:] .- r[tri[:,1],:]
    r32 = r[tri[:,3],:] .- r[tri[:,2],:]

    det = r21[:,1].*r32[:,2] - r21[:,2].*r32[:,1]
    idx = findall(det .< 0)

    tmp = tri[idx, 2]
    tri[idx, 2] .= tri[idx,3]
    tri[idx, 3] .= tmp

    return tri
end

"""
    function triangulation(r₀, rᵢ)

Return the triangle containing each point given in `rᵢ`.
The triangulation used is delaunay triangulation defined by point cloud `r₀`.
"""
function triangulation(r₀, rᵢ)
    tri = order(delaunay(Array{eltype(r₀),2}(r₀'))', r₀)

    # compute containing faces
    t1  = r₀[tri[:,2],:] - r₀[tri[:,1],:]
    t2  = r₀[tri[:,3],:] - r₀[tri[:,2],:]
    t3  = r₀[tri[:,1],:] - r₀[tri[:,3],:]

    t1 /= sqrt.(t1⋅t1)
    t2 /= sqrt.(t2⋅t2)
    t3 /= sqrt.(t3⋅t3)

    t1  = t1*[0 +1; -1 0]
    t2  = t2*[0 +1; -1 0]
    t3  = t3*[0 +1; -1 0]

    e1  = .5*(r₀[tri[:,2],:] + r₀[tri[:,1],:])
    e2  = .5*(r₀[tri[:,3],:] + r₀[tri[:,2],:])
    e3  = .5*(r₀[tri[:,1],:] + r₀[tri[:,3],:])

    θ1  = t1*rᵢ' .- sum(t1.*e1,dims=2)
    θ2  = t2*rᵢ' .- sum(t2.*e2,dims=2)
    θ3  = t3*rᵢ' .- sum(t3.*e3,dims=2)

    τ(i)  = findfirst((θ1[:,i] .> 0) .& (θ2[:,i] .> 0) .& (θ3[:,i] .> 0))
    faces = [τ(i) for i in 1:size(rᵢ,1)]
    return tri, faces
end

function continuum(tri, ϕ, r)
    function fit(ϕ, r)
        M = hcat(r[:,1], r[:,2], ones(size(r,1)))
        ρ = M \ ϕ

        f(x) = ρ[1]*x[1] + ρ[2]*x[2] + ρ[3]
        return f
    end

    return [fit(ϕ[tri[t,:]], r[tri[t,:],:]) for t in 1:size(tri,1)]
end

"""
    function interpolate(M, ϕ, r)

Interpolate the tensor field `ϕ`, defined at points `r` onto the vertices of manifold `M`.
Interpolation is computed by finding the containing triangle within the mesh of `M` for each point `r`.
Linear interpolation is performed per face.
"""
function interpolate(M, ϕ, r)
    x, _ = pullback(M)
    y    = pullback(M.surf, r)

    tri, face = triangulation(y, x)
    dist = sum( (x[:,d] .- y[:,d]').^2 for d in 1:size(x,2) )

    @match numdims(ϕ) begin
        1 => begin
            field  = continuum(tri, ϕ, y)
            result = Array{eltype(ϕ)}(undef, size(x,1))
            for i in 1:size(x,1)
                if isnothing(face[i])
                    result[i] = ϕ[argmin(dist[i,:])]
                else
                    result[i] = field[face[i]](x[i,:])
                end
            end

            return result
        end
        2 => begin
            ψ = Array{eltype(ϕ)}(undef, size(x,1), 3)
            for d in 1:3
                field  = continuum(tri, ϕ[:,d], y)
                for i in 1:size(x,1)
                    if isnothing(face[i])
                        ψ[i,d] = ϕ[argmin(dist[i,:]), d]
                    else
                        ψ[i,d] = [field[face[i]](x[i,:]) for i in 1:size(x,1)]
                    end
                end
            end
            return ψ
        end
        _ => return error("unsupported field dimension = $(size(ϕ))")
    end
end

"""
    rescale(x)

Rescale array `x` to run between ``[0,1]``.
"""
function rescale(x)
    min = minimum(x)
    max = maximum(x)
    return (x.-min)./(max.-min)
end

"""
    scalar(M::Manifold, ϕ, field::Symbol)

Return a scalar field `ϕ` interpolated onto either the embedded space of manifold `M` or the pullback.
If `field` is :ℝ² the pullback is computed.
If `field` is :ℝ³ the embedding space is returned.
"""
function scalar(M::Manifold, ϕ, field::Symbol)
    @match field begin
        :ℝ² => begin
            q, tri = pullback(ℳ)
            return (
                x1  = q[:,1],
                x2  = q[:,2],
                tri = tri',
                ϕ   = ϕ,
            )
        end
        :ℝ³ => begin
            return (
                x1  = M.mesh.r[:,1],
                x2  = M.mesh.r[:,2],
                x3  = M.mesh.r[:,3],
                tri = M.mesh.tri',
                ϕ   = ϕ,
            )
        end
        _   => return error("unsupported field $(field)")
    end
end

"""
    scalar(M::Manifold, ϕ, field::Symbol)

Return a vector field field `ϕ` interpolated onto either the embedded space of manifold `M` or the pullback.
If `field` is :ℝ² the pullback is computed. Only the tangent space component is kept.
If `field` is :ℝ³ the embedding space is returned.
"""
function vector(ℳ::Manifold, ϕ, field::Symbol)
    @match field begin
        :ℝ² => begin
                q, tri = pullback(ℳ)
                ez, eθ = basis(ℳ.surf, ℳ.mesh.r)
                ϕz, ϕθ = sum(ϕ.*ez,dims=2), sum(ϕ.*eθ,dims=2)

                return (
                    x1  = q[:,1],
                    x2  = q[:,2],
                    tri = tri',
                    v1  = ϕz,
                    v2  = ϕθ,
                )
            end
        :ℝ³ => begin
                return (
                    x1  = M.mesh.r[:,1],
                    x2  = M.mesh.r[:,2],
                    x3  = M.mesh.r[:,3],
                    tri = M.mesh.tri',
                    v1  = ϕ[:,1],
                    v2  = ϕ[:,2],
                    v3  = ϕ[:,3],
                )
            end
        _   => return error("unsupported field $field")
    end
end

"""
    makefield(ℳ::Manifold, ϕ::AbstractArray; field::Symbol = :ℝ³)

Return an arbitrary tensor field defined over _either_ the embedding space or the cylindrical pullback of manifold `M`.
Higher level function than either [`scalar`](@ref) or [`vector`](@ref).
"""
function makefield(ℳ::Manifold, ϕ::AbstractArray; field::Symbol = :ℝ³)
    @match numdims(ϕ) begin
        1                         => return scalar(ℳ, ϕ, field)
        2, if size(ϕ,2) == 3 end  => return vector(ℳ, ϕ, field)
        _                         => return error("unsupported field dimension = $(size(ϕ))")
    end
end

function manifold(io::IO, type::Symbol=:obj)
    m = mesh(io, type)
    s = Surface(m.r)

    return Manifold(m, s)
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
#=
This file is auto-generated. Do not edit.
=#

mutable struct TechHydro <: TechnicalParams
    rating::Float64  # Thermal limited MVA Power Output of the unit. <= Capacity
    primemover::PrimeMovers  # PrimeMover Technology according to EIA 923
    activepowerlimits::Min_Max
    reactivepowerlimits::Union{Nothing, Min_Max}
    ramplimits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    timelimits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    internal::InfrastructureSystemsInternal
end

function TechHydro(rating, primemover, activepowerlimits, reactivepowerlimits, ramplimits, timelimits, )
    TechHydro(rating, primemover, activepowerlimits, reactivepowerlimits, ramplimits, timelimits, InfrastructureSystemsInternal())
end

function TechHydro(; rating, primemover, activepowerlimits, reactivepowerlimits, ramplimits, timelimits, )
    TechHydro(rating, primemover, activepowerlimits, reactivepowerlimits, ramplimits, timelimits, )
end

# Constructor for demo purposes; non-functional.

function TechHydro(::Nothing)
    TechHydro(;
        rating=0.0,
        primemover=HY::PrimeMovers,
        activepowerlimits=(min=0.0, max=0.0),
        reactivepowerlimits=nothing,
        ramplimits=nothing,
        timelimits=nothing,
    )
end

"""Get TechHydro rating."""
get_rating(value::TechHydro) = value.rating
"""Get TechHydro primemover."""
get_primemover(value::TechHydro) = value.primemover
"""Get TechHydro activepowerlimits."""
get_activepowerlimits(value::TechHydro) = value.activepowerlimits
"""Get TechHydro reactivepowerlimits."""
get_reactivepowerlimits(value::TechHydro) = value.reactivepowerlimits
"""Get TechHydro ramplimits."""
get_ramplimits(value::TechHydro) = value.ramplimits
"""Get TechHydro timelimits."""
get_timelimits(value::TechHydro) = value.timelimits
"""Get TechHydro internal."""
get_internal(value::TechHydro) = value.internal


text:
using Parameters

####### convergence constants #######
@inline function chargecon()
    return 18.22261720426243437986
end

@inline function debye()
    return 4.8033324 # eÅ
end

abstract type TTM_Constants end

@with_kw struct TTM21_Constants <: TTM_Constants
    # vdw
    vdwA::Float64 = -1.329565985e+6
    vdwB::Float64 =  3.632560798e+5
    vdwC::Float64 = -2.147141323e+3
    vdwD::Float64 =  1.0e+13
    vdwE::Float64 =  13.2

    # M-site positioning
    γ_M::Float64 = 0.426706882
    γ_1::Float64 = 1.0 - γ_M
    γ_2::Float64 = γ_M / 2

    # polarizability
    α_O::Float64 = 0.837
    α_H::Float64 = 0.496
    α_M::Float64 = 0.0

    # Thole damping factors
    damping_factor_O::Float64 = α_O
    damping_factor_H::Float64 = α_H
    damping_factor_M::Float64 = α_O

    name::Symbol = :ttm21
end

@with_kw struct TTM22_Constants <: TTM_Constants
    # vdw
    vdwA::Float64 = -1.329565985e+6
    vdwB::Float64 =  3.632560798e+5
    vdwC::Float64 = -2.147141323e+3
    vdwD::Float64 =  1.0e+13
    vdwE::Float64 =  13.2

    # dms (start with the ttm3 modified dms)
    dms_param1::Float64 = 0.5
    dms_param2::Float64 = 0.9578
    dms_param3::Float64 = 0.012

    # M-site positioning
    γ_M::Float64 = 0.426706882
    γ_1::Float64 = 1.0 - γ_M
    γ_2::Float64 = γ_M / 2

    # polarizability
    α_O::Float64 = 0.837
    α_H::Float64 = 0.496
    α_M::Float64 = 0.0

    # Thole damping factors
    damping_factor_O::Float64 = α_O
    damping_factor_H::Float64 = α_H
    damping_factor_M::Float64 = α_O

    name::Symbol = :ttm22
end

@with_kw struct TTM3_Constants <: TTM_Constants
    # vdw
    vdwC::Float64 = -0.72298855E+03
    vdwD::Float64 =  0.10211829E+06
    vdwE::Float64 =  0.37170376E+01

    # dms
    dms_param1::Float64 = 0.5
    dms_param2::Float64 = 0.9578
    dms_param3::Float64 = 0.012

    # M-site positioning
    γ_M::Float64 = 0.46
    γ_1::Float64 = 1.0 - γ_M
    γ_2::Float64 = γ_M / 2

    # polarizability
    α_O::Float64 = 0.0
    α_H::Float64 = 0.0
    α_M::Float64 = 1.444

    # Thole damping factors
    damping_factor_O::Float64 = 0.837
    damping_factor_H::Float64 = 0.496
    damping_factor_M::Float64 = 0.837

    name::Symbol = :ttm3
end

@with_kw struct TTM4_Constants <: TTM_Constants
    # vdw
    A6::Float64 = -0.126503e+5
    A8::Float64 =  0.526347e+6
    A10::Float64 = -0.964270e+7
    A12::Float64 =  0.877792e+8
    A14::Float64 = -0.367476e+9
    A16::Float64 =  0.572395e+9

    # M-site positioning
    γ_M = 0.426706882
    γ_1 = 1.0 - gammaM
    γ_2 = gammaM / 2

    # polarizability
    α_O::Float64 = 1.310
    α_H::Float64 = 0.294
    α_M::Float64 = 0.0

    # Thole damping factors
    damping_factor_O::Float64 = α_O
    damping_factor_H::Float64 = α_H
    damping_factor_M::Float64 = α_O

    name::Symbol = :ttm4
end

# probably factor all of the qtip4pf stuff into it's own directory
@with_kw struct qtip4pf_Constants
    qM::Float64 = -1.1128 * chargecon()
    qH::Float64 = -qM/2

    rOHeq::Float64= 0.9419 # A
    aHOHeq::Float64 = 107.4 # degree

    epsilon::Float64 = 0.1852 # kcal/mol
    sigma::Float64 = 3.1589 # A

    alpha_r::Float64 = 2.287 # A-1
    D_r::Float64 = 116.09 # kcal/mol
    k_theta::Float64 = 87.85 # kcal/mol/rad^2

    gammaM::Float64 = 0.73612
    gamma1::Float64 = (1.0 - gammaM)/2
end

function get_constants_type(potential::Symbol)
    if potential == :ttm3
        return TTM3_Constants()
    elseif potential == :ttm21
        return TTM21_Constants()
    elseif potential == :ttm4
        return TTM4_Constants()
    else
        @assert false "Didn't receive a valid potential. Pass as a symbol :ttm3, :ttm21, or :ttm4."
    end
end

max_stars_repo_name: heindelj/MolecularPotentials

max_issues_repo_name: heindelj/MolecularPotentials

max_forks_repo_name: heindelj/MolecularPotentials


text:
include("filtering.jl");
gr(size=(600,600), tick_orientation = :out, grid = false,
    linecolor = :black,
    markerstrokecolor = :black,
    thickness_scaling = 2,
    markersize = 6)
check = combine(groupby(pokes,[:Phase,:Day]),:Treatment => t -> [union(t)],:Stim => t -> [union(t)])
    open_html_table(sort!(check,:Day))
    check = combine(groupby(streaks,[:Phase,:Day]),:Treatment => t -> [union(t)],:Stim => t -> [union(t)])
    open_html_table(sort!(check,:Day))
################################ Adjust Streaks table ##################################
list = ["PreVehicle",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Methysergide",
    "Citalopram",
    "SB242084_opt",
    "Saline"]
for df in [pokes, streaks]
    filter!(r->r.Treatment in list &&
        r.Trial < 51 &&
        r.MouseID != "pc7",
        df)
    (df.Treatment .== "Saline") .& (df.Phase .== "Optogenetic")
    df[df.Treatment .== "PreVehicle",:Treatment] .= "Control"
    (df.Treatment .== "Saline") .& (df.Phase .== "Optogenetic")
    df[(df.Treatment .== "Saline") .& (df.Phase .== "Optogenetic"),:Treatment] =
        [o ? "Optogenetic" : "Control" for o in df[(df.Treatment .== "Saline") .& (df.Phase .== "Optogenetic"),:Stim]]
    df[df.Treatment .== "SB242084_opt",:Phase] .=  "SB242084_opt"
    df[df.Treatment .== "SB242084_opt",:Treatment] = [o ? "SB242084_opt" : "Control" for o in df[df.Treatment .== "SB242084_opt",:Stim]]
    df[!,:Treatment] = categorical(df.Treatment, ordered = false)
end
streaks[!,:ROILeavingTime] = streaks.Stop_trial .- streaks.Stop_poking
for df in [pokes, streaks]
    levels!(df.Treatment,["Control",
        "Altanserin",
        "SB242084",
        "Way_100135",
        "Citalopram",
        "Optogenetic",
        "Methysergide",
        "SB242084_opt",
        "Saline"
        ])
end
levels(streaks.Treatment)

check = combine(groupby(pokes,[:Phase,:Day]),:Treatment => t -> [union(t)],:Stim => t -> [union(t)])
    open_html_table(sort!(check,:Day))
    check = combine(groupby(streaks,[:Phase,:Day]),:Treatment => t -> [union(t)],:Stim => t -> [union(t)])
    open_html_table(sort!(check,:Day))
pcheck = filter(r-> r.Treatment == "Citalopram", pokes)
union(pcheck.MouseID)
union(pcheck.Treatment)

################################ Adjust Pokes table ##################################
filt_1 = filter(r->r.Treatment == "Control",pokes)
gd1 = groupby(filt_1,[:MouseID,:TimeFromLeaving,:Protocol])
df2 = combine(gd1, :CumRewTrial => mean => :CumRewTrial)
gd2 = groupby(df2,[:Protocol,:TimeFromLeaving])
df3 = combine(gd2, :CumRewTrial => mean, :CumRewTrial => sem)
Protocol_colors!(df3)

@df df3 scatter(:TimeFromLeaving,:CumRewTrial_mean, group = :Protocol, xlims = (0,30), markersize = 4, color = :color)
##
gd1 = groupby(streaks,[:Protocol,:MouseID,:Phase,:Treatment])
df1 = combine(gd1,
    :Leaving_NextPrew => mean => :Leaving_NextPrew,
    :AverageRewRate => mean => :AverageRewRate)

gd2 = groupby(df1,:Treatment)
df2 = combine(gd2,
    :AverageRewRate .=> [mean, median], :AverageRewRate .=> [sem,CIq],
    :Leaving_NextPrew .=> [mean, median], :Leaving_NextPrew .=> [sem, CIq])

filt_2 = filter(r->r.Treatment in ["Altanserin", "SB242084", "Control", "Citalopram", "Optogenetic"],df2)
res = combine(groupby(filt_2, :Treatment)) do dd
    DataFrame(Condition = ["Average", "At Leaving"],
        Mean = [dd[1,:AverageRewRate_mean], dd[1,:Leaving_NextPrew_mean]],
        Sem = [dd[1,:AverageRewRate_sem], dd[1,:Leaving_NextPrew_sem]],
        Median = [dd[1,:AverageRewRate_median], dd[1,:Leaving_NextPrew_median]],
        CIq = [dd[1,:AverageRewRate_CIq], dd[1,:Leaving_NextPrew_CIq]])
    end
Drug_colors!(res)
@df res scatter(:Condition,:Mean, yerror = :Sem, color = :color, xlims = (0.25,1.75), markeralpha=0.7, group = :Treatment, legend = false)
@df res scatter(:Condition,:Median, yerror = :Sem, color = :color, xlims = (0.25,1.75), markeralpha=0.7, group = :Treatment, legend = false)
##
#=
1- bin cumulative rewards over time from trial start
1b- skip interpoke interval cost
    1- accumulate poke duration over trial
2- fit the cumulative with polyfit LsqFit
3- find average travel speed
4- find a tangent
=#

filt_3  = filter(r->r.Treatment in ["Altanserin", "SB242084", "Control", "Citalopram", "Optogenetic"],pokes)
gd = groupby(filt_3, [:Day,:MouseID, :Trial])
transform!(gd,[:PokeIn,:PokeOut] => ((i,o) -> (round.(o .- i[1], digits = 1))) => :BinnedTrialTime,
    [:PokeIn,:PokeOut] => ((i,o) -> round.(cumsum(o .- i), digits = 1)) => :ForagingTime,
    :Reward => cumsum => :CumTrialReward)

timedf = combine(groupby(filt_3, [:Day,:MouseID, :Trial, :Treatment])) do dd
    timerange = 0:0.1:maximum(dd.BinnedTrialTime)
    # timerange = 0:0.1:maximum(dd.ForagingTime)
    df = DataFrame(Time = collect(timerange), Rew = 0)
    for r in eachrow(dd)
        df[findfirst(timerange .== r.BinnedTrialTime),:Rew] = r.Reward
        # df[findfirst(timerange2 .== r.ForagingTime),:Rew] = r.Reward
    end
    df.CumRew = cumsum(df.Rew)
    return df
end
timedf2 = combine(groupby(filt_3, [:Day,:MouseID, :Trial, :Treatment])) do dd
    # timerange = 0:0.1:maximum(dd.BinnedTrialTime)
    timerange = 0:0.1:maximum(dd.ForagingTime)
    df = DataFrame(Time = collect(timerange), Rew = 0)
    for r in eachrow(dd)
        # df[findfirst(timerange .== r.BinnedTrialTime),:Rew] = r.Reward
        df[findfirst(timerange .== r.ForagingTime),:Rew] = r.Reward
    end
    df.CumRew = cumsum(df.Rew)
    return df
end
# open_html_table(timedf[2400:3300,:])
union(timedf.Treatment)
## Fit cumulative per session
using LsqFit
function fit_cumulative(cum,time)
    model(x,p) = @. (1 - exp(-p[1]*x)) *p[2]
    p0 = [0.5,1]
    fitting = curve_fit(model, cum, time, p0)
    # return DataFrame(Coef = fitting.param[1], Scale = fitting.param[2])
    # return (fitting.param...,)
    return fitting.param
end
using Roots
function MVTtangent(coef, P)
    # P is the point where the line pass through
    # cumulative to be tangential to
    c(x) = @. (1 - exp(-coef[1]*x))*coef[2]
    #derivative of the cum, equal to the slope
    d(x) = @. coef[1] * exp(-coef[1]*x) *coef[2]
    #solve using slope-intercept form of the line
    # y-y0 = m(x-x0)
    # x0 has to be a point on the curve so express it in this way
    # y-c(x0) = d(x0)(x-x0)
    # solve passing for the point P
    # yp -c(x0) = d(x0)(xp-x0)
    # yp -c(x0) - d(x0)(xp-x0) = 0
    # f(x) = P[2] - (1 - exp(-coef[1]*x))*coef[2] - (coef[1] * exp(-coef[1]*x))*coef[2]*(P[1]-x)
    f(x) = P[2] - c(x) - d(x)*(P[1]-x)
    x0 = Roots.find_zero(f,0) # search x0 in a reasonable range of trial time
    #find inercept using P and m = d(x0)
    # y = mx +q
    # P[2] = d(x0)*P[1] + q
    # q = -d(x0)*P[1] - P[2]
    q = -d(x0)*P[1] - P[2]
    # sol(x) = @. x * d(x0) + q
    # return DataFrame(OptLeave = x0, Slope = d(x0), Intercept = q)
    return (x0,d(x0),q)
end

control = filter(r->r.Treatment == "Control", timedf)
coef = fit_cumulative(control.Time, control.CumRew)
control2 = filter(r->r.Treatment == "Control", timedf2)
coef2 = fit_cumulative(control2.Time, control2.CumRew)
GainModel(x,p) = @. (1 - exp(-p[1]*x)) *p[2]
xaxis = 0:0.5:35
plot(xaxis,GainModel(xaxis,coef), label = "TrialTime", legend = :bottomright)
plot!(xaxis,GainModel(xaxis,coef2), label = "ForagingTime", linecolor = :purple)
## Different gain models per animal
coef3 = combine(groupby(control,[:MouseID]), [:CumRew,:Time] => ((cr,t) -> [fit_cumulative(cr,t)])=> :Coeff)
coef3[!,:GlobCoeff] .= [coef]
open_html_table(coef3)
##
# timedf2 = combine(groupby(timedf,[:MouseID,:Day, :Treatment]),
#     # [:Time,:CumRew] => ((t,r)-> fit_cumulative(t,r)) => [:Coef,:Scale])
#     [:Time,:CumRew] => ((t,r)-> fit_cumulative(t,r)) => :Params)
## Costs
# travel
pretravel = filter(r->r.Treatment in ["Altanserin", "SB242084", "Control", "Citalopram", "Optogenetic"],streaks)
trav1 = combine(groupby(pretravel,[:MouseID,:Treatment]),
    :Poking_Travel_to => median => :Travel,
    :Poking_duration => median => :Leaving)
# Poking

function pokecost(rvec,pokevec)
    v = pokevec[.!rvec]
    median(skipmissing(v))
end
filt_3.Pre_Interpoke
pcost = combine(groupby(filt_3, [:MouseID, :Treatment]), [:Reward,:PokeDuration] => ((r,p) -> pokecost(r,p)) => :PokeCost,
    :Pre_Interpoke => (i -> median(skipmissing(i))) => :InterpokeCost)
cost = leftjoin(trav1,pcost, on=[:MouseID, :Treatment])
params = leftjoin(cost, coef3; on=:MouseID)
## Optimal leaving
timedf3 = combine(groupby(params,[:MouseID,:Treatment]),
    [:Travel, :PokeCost, :Coeff] => ByRow(((trav,pcost, coeff) -> MVTtangent(coeff, [-trav+pcost,0.0]))) => [:OptLeave, :Slope, :Intercept],
    :PokeCost,:Travel,:Leaving)
timedf3.Diff = timedf3.Leaving .- timedf3.OptLeave
sort!(timedf3,[:Treatment,:MouseID])
open_html_table(timedf3)
timedf4 = combine(groupby(timedf3,:Treatment)) do dd
    x = Vector(dd.Diff)
    wilcoxon(x)
end
open_html_table(timedf4)
Drug_colors!(timedf4)
@df timedf4 scatter(string.(:Treatment), :Median, yerror = :CI, color = :color, legend = false)
Plots.abline!(0,0,color = :black, linestyle = :dash)
##GLM
fm = @formula(Diff ~ 1 + Treatment + (1|MouseID+Treatment))
fm1 = fit(MixedModel, fm, timedf3)
## Leaving by Poke_within_Trial
#= using median interpoke and median unrewarded pokes duration we can normalize median travel time on the same
abstract unit measure =#
#1 filter pokes and streaks
fpokes = filter(r->r.Treatment in ["Altanserin", "SB242084", "Control", "Citalopram", "Optogenetic"],pokes)
fstreaks = filter(r->r.Treatment in ["Altanserin", "SB242084", "Control", "Citalopram", "Optogenetic"],streaks)
#2 trial cumulative reward in pokes df
gd = groupby(fpokes, [:Day,:MouseID, :Trial])
transform!(gd,:Reward => cumsum => :CumTrialReward)
#3 calculate median travel, per session
travcost = combine(groupby(fstreaks,[:MouseID,:Treatment,:Day]),
    # :Poking_Travel_to => (t->median(t[t.>=1])) => :Travel)
    :Poking_Travel_to => median => :Travel,
    :Num_pokes => mean => :Leave
    )
#4 calculate median poke duration and median interpoke duration
pcost = combine(groupby(fpokes,[:MouseID,:Treatment,:Day]),
    [:Reward,:PokeDuration] => ((r,p) -> pokecost(r,p)) => :PokeCost,
    :Pre_Interpoke => (i -> median(skipmissing(i))) => :InterpokeCost)
#5 join costs datasets
costs = leftjoin(travcost,pcost, on = [:MouseID,:Treatment,:Day])
costs[!,:Pcost] = costs.InterpokeCost .+ costs.PokeCost
costs[!,:NormTravel] = -costs.Travel ./ costs.Pcost
#6 calculate cumulative reward per ...
GainModel(x,p) = @. (1 - exp(-p[1]*x)) *p[2]
# coef = combine(groupby(fpokes,[:Treatment,:Day]),
#     [:CumTrialReward,:Poke_within_Trial] => ((cr,t) -> [fit_cumulative(cr,t)])=> :Coeff)
coef = combine(groupby(filter(r-> r.Treatment == "Control", fpokes),[:MouseID]),
    [:CumTrialReward,:Poke_within_Trial] => ((cr,t) -> [fit_cumulative(cr,t)])=> :Coeff)
res = leftjoin(costs, coef, on = [:MouseID])
#7 find optimal leaving time per session
groups = groupby(res, [:MouseID,:Treatment,:Day])
transform!(groups,
    [:NormTravel, :Coeff] => ByRow(((ntrav, coeff) -> MVTtangent(coeff, [ntrav,0.0]))) => [:OptLeave, :Slope, :Intercept]
    )

open_html_table(sort!(res,:Treatment))
#8 calculate difference between otpimal and actual leaving time
res[!,:Diff] = res.Leave .- res.OptLeave
fm = @formula(Diff ~ 1 + Treatment + (1|MouseID+Treatment))
fm1 = fit(MixedModel, fm, res)
##
reswilc = combine(groupby(res,:Treatment)) do dd
    x = Vector(dd.Diff)
    wilcoxon(x)
end
Drug_colors!(reswilc)
@df reswilc scatter(string.(:Treatment), :Median, yerror = :CI, color = :color, legend = false)
Plots.abline!(0,0,color = :black, linestyle = :dash)


text:
function value_of(::Type{Keyword}, arg0::JString)
    return jcall(Keyword, "valueOf", Keyword, (JString,), arg0)
end

function values(::Type{Keyword})
    return jcall(Keyword, "values", Vector{Keyword}, ())
end


text:
println("=======================================================================")
println("Beginning GMHExamples.jl/ode/fitzhughnagumo/scripts/FitzHughNagumo1.jl")
println("=======================================================================")

import GeneralizedMetropolisHastings #pre-include to avoid race condition when running in parallel
import GMHExamples

println("Number of parallel processes running: ",nprocs())
println("Use addprocs() in the REPL if you want to run on more processes")

@everywhere using GeneralizedMetropolisHastings
@everywhere using GMHExamples
using PyPlot

#Standard M-H for nproposals == 1
#Generalized M-H for nproposals > 1
nproposals = 30

#MCMC iteration specifications
nburnin = 1000
niterations = 1000
ntunerperiod = 50

###Initial conditions for the ODE (membrane potential and refractory variable)
initial = [-1.0,1.0]

###Default values of the parameters (a,b,c) and prior boundaries
defaults = [0.3,0.3,2.0]
lows = zeros(3)
highs = [5.0,5.0,5.0]

###The variance of the noise on the input data
variance = [0.02,0.005]

println("==========================================")
println("Simulation parameters defined successfully")
println("==========================================")

###Create a FitzHugh-Nagumo model with measurement data, ODE function and parameters with default values and priors
m = fitzhughnagumomodel(initial,variance,lows,highs,defaults)

###Show the model
println("==========================")
println("Model defined successfully")
println("==========================")
show(m)

###Plot the measurement data (simmulated data + noise)
figure(string("FitzHughNagumo1-",nprocs())) ; clf()
subplot(231)
plot(dataindex(m),measurements(m)[:,1];label="membrane potential")
plot(dataindex(m),measurements(m)[:,2];label="refractory variable")
xlabel("Time")
ylabel("Amplitude")
title("FitzHugh-Nagumo measurement data")
grid("on")
legend(loc="upper right",fancybox="true")

###Create a Metropolis sampler with a Normal proposal density
s = sampler(:mh,:normal,0.5,ones(3))
println("============================")
println("Sampler defined successfully")
println("============================")
show(s)

###Create a tuner that scales the proposal density
t = tuner(:scale,ntunerperiod,0.5,:erf)
println("==========================")
println("Tuner defined successfully")
println("==========================")
show(t)

###Create a Generalized Metropolis-Hastings runner (which will default to Standard MH when nproposals=1)
p = policy(:mh,nproposals;initialize=:prior)
r = runner(p,niterations,nproposals;numburnin = nburnin)
println("===========================")
println("Runner defined successfully")
println("===========================")
show(r)

###Run the MCMC (can take quite a bit of time)
println("=======================")
println("Run the MCMC simulation")
println("=======================")
@time c = run!(r,m,s,t)
println("==========================")
println("Completed MCMC simulation")
println("=========================")

###Show the results of the simulations
show(c)

nparas = numparas(m)
meanparamvals = mean(samples(c),2)
stdparamvals = std(samples(c),2)

println("Results of the MCMC simulation:")
for i=1:nparas
    println(" parameter $(parameters(m)[i].key):  mean = ",meanparamvals[i]," std = ",stdparamvals[i])
end

println("================")
println("Plotting results")
println("================")

###Plot the loglikelihood values across samples
###After an initial few low values, this should remain relatively high
subplot(232)
plot(1:numsamples(c),logposterior(c,m))
title("Log-Posterior values across samples")
xlabel("Samples")
ylabel("Log-Posterior")

###Plot the histograms of a,b,c values
for i=1:nparas
  subplot(233 + i)
    h = PyPlot.plt[:hist](vec(getindex(samples(c),i,:)),20)
  grid("on")
  title("Parameter $(parameters(m)[i].key)")
  xlabel("Values")
  ylabel("Samples")
end

subplot(233)
modeldata = evaluate!(m,vec(meanparamvals))
plot(dataindex(m),measurements(m)[:,1];label="membrane potential")
plot(dataindex(m),measurements(m)[:,2];label="refractory variable")
plot(dataindex(m),modeldata[:,1])
plot(dataindex(m),modeldata[:,2])
xlabel("Time")
ylabel("Amplitude")
title("FitzHugh-Nagumo model data")
grid("on")
legend(loc="upper right",fancybox="true")


text:
module CvxCompress

using CvxCompress_jll

import Base.copy

"""
    CvxCompressor{N}

Struct used for managing the compression target
"""
mutable struct CvxCompressor{N}
    br::NTuple{N,Int}
    scale::Float32
    compressed_length::Base.RefValue{Clong}
end

"""
    CvxCompressor(b::NTuple{N,Int}, scale=1e-2)

Create a compressor struc instance
"""
CvxCompressor(b::NTuple{N,Int}, scale=1e-2) where {N} = CvxCompressor{N}(b, scale, Ref{Clong}(0))

"""
    CvxCompressor(b::Vararg{Int,N})

Create a compressor struc instance
"""
CvxCompressor(b::Vararg{Int,N}) where {N} = CvxCompressor(b)

"""
    copy(c::CvxCompressor{N})

Create a copy of a CvxCompressor struct
"""
copy(c::CvxCompressor{N}) where {N} = CvxCompressor{N}(c.br, c.scale, Ref{Clong}(0))

"""
    compress!(compressed_volume::Array{UInt32,1}, c::CvxCompressor{3}, volume::Array{Float32,3})

3D compression of a volume
"""
function compress!(compressed_volume::Array{UInt32,1}, c::CvxCompressor{3}, volume::Array{Float32,3})
    nz, ny, nx = size(volume)
    ccall((:cvx_compress, libcvxcompress),
          Float32,
          (Cfloat, Ptr{Cfloat}, Cint, Cint, Cint, Cint,    Cint,    Cint,    Ptr{Cuint},        Ref{Clong}),
          c.scale, volume,      nz,   ny,   nx,   c.br[1], c.br[2], c.br[3], compressed_volume, c.compressed_length)
    c.compressed_length[]
end

"""
    compress!(compressed_volume::Array{UInt32,1}, c::CvxCompressor{3}, volume::Array{Float64,3})

3D compression of a volume
"""
compress!(compressed_volume::Array{UInt32,1}, c::CvxCompressor{3}, volume::Array{Float64,3}) = compress!(compressed_volume, c, convert(Array{Float32,3}, volume))

"""
    decompress!(volume::Array{Float32,3}, c::CvxCompressor{3}, compressed_volume::Array{UInt32,1}, compressed_length::Integer)

3D de-compression of a volume
"""
function decompress!(volume::Array{Float32,3}, c::CvxCompressor{3}, compressed_volume::Array{UInt32,1}, compressed_length::Integer)
    nz, ny, nx = size(volume)
    ccall((:cvx_decompress_inplace, libcvxcompress),
          Cvoid,
          (Ptr{Cfloat}, Cint, Cint, Cint, Ptr{Cuint},        Clong),
          volume,       nz,   ny,   nx,   compressed_volume, compressed_length)
end

"""
    decompress!(volume::Array{Float64,3}, c::CvxCompressor{3}, compressed_volume::Array{UInt32,1}, compressed_length::Integer)

3D de-compression of a volume
"""
function decompress!(volume::Array{Float64,3}, c::CvxCompressor{3}, compressed_volume::Array{UInt32,1}, compressed_length::Integer)
    v = Array{Float32}(undef, size(volume))
    decompress!(v, c, compressed_volume, compressed_length)
    volume[:] = v[:]
end

"""
    compress!(compressed_volume::Array{UInt32,1}, c::CvxCompressor{2}, volume::Array{Float32,2})

2D compression of a volume
"""
function compress!(compressed_volume::Array{UInt32,1}, c::CvxCompressor{2}, volume::Array{Float32,2})
    nz, nx = size(volume)
    ccall((:cvx_compress, libcvxcompress),
          Float32,
          (Cfloat, Ptr{Cfloat}, Cint, Cint, Cint, Cint,    Cint,    Cint,    Ptr{Cuint},        Ref{Clong}),
          c.scale, volume,      nz,   nx,   1,    c.br[1], c.br[2], 1,       compressed_volume, c.compressed_length)
    c.compressed_length[]
end

"""
    compress!(compressed_volume::Array{UInt32,1}, c::CvxCompressor{2}, volume::Array{Float64,2})

2D compression of a volume
"""
compress!(compressed_volume::Array{UInt32,1}, c::CvxCompressor{2}, volume::Array{Float64,2}) = compress!(compressed_volume, c, convert(Array{Float32,2}, volume))

"""
    decompress!(volume::Array{Float32,2}, c::CvxCompressor{2}, compressed_volume::Array{UInt32,1}, compressed_length::Integer)

2D de-compression of a volume
"""
function decompress!(volume::Array{Float32,2}, c::CvxCompressor{2}, compressed_volume::Array{UInt32,1}, compressed_length::Integer)
    nz, nx = size(volume)
    ccall((:cvx_decompress_inplace, libcvxcompress),
          Cvoid,
          (Ptr{Cfloat}, Cint, Cint, Cint, Ptr{Cuint},        Clong),
          volume,       nz,   nx,   1,    compressed_volume, compressed_length)
end

"""
    decompress!(volume::Array{Float64,2}, c::CvxCompressor{2}, compressed_volume::Array{UInt32,1}, compressed_length::Integer)

2D de-compression of a volume
"""
function decompress!(volume::Array{Float64,2}, c::CvxCompressor{2}, compressed_volume::Array{UInt32,1}, compressed_length::Integer)
    v = Array{Float32}(size(volume))
    decompress!(v, c, compressed_volume, compressed_length)
    volume[:] = v[:]
end

export CvxCompressor
export CvxCompressor2D
export CvxCompressor3D
export compress!
export decompress!

end