text:
using FactCheck, Images, Colors, FixedPointNumbers

testing_units = Int == Int64
if testing_units
    using SIUnits, SIUnits.ShortUnits
end

const scalar_getindex_new = VERSION >= v"0.5.0-dev+1195"

facts("Core") do
    a = rand(3,3)
    @inferred(Image(a))
    # support integer-valued types, but these are NOT recommended (use UFixed)
    B = rand(convert(UInt16, 1):convert(UInt16, 20), 3, 5)
    # img, imgd, and imgds will be used in many more tests
    # Thus, these must be defined as local if reassigned in any context() blocks
    cmap = reinterpret(RGB, repmat(reinterpret(UFixed8, round(UInt8, linspace(12, 255, 20)))', 3, 1))
    img = ImageCmap(copy(B), cmap, Dict{ASCIIString, Any}([("pixelspacing", [2.0, 3.0]), ("spatialorder", Images.yx)]))
    imgd = convert(Image, img)
    if testing_units
        imgd["pixelspacing"] = [2.0mm, 3.0mm]
    end
    imgds = separate(imgd)

    context("Constructors of Image types") do
        @fact colorspace(B) --> "Gray"
        @fact colordim(B) --> 0
        local img = Image(B, colorspace="RGB", colordim=1)  # keyword constructor
        @fact colorspace(img) --> "RGB"
        @fact colordim(img) --> 1
        img = grayim(B)
        @fact colorspace(img) --> "Gray"
        @fact colordim(B) --> 0
        @fact grayim(img) --> img
        # this is recommended for "integer-valued" images (or even better, directly as a UFixed type)
        Bf = grayim(round(UInt8, B))
        @fact eltype(Bf) --> UFixed8
        @fact colorspace(Bf) --> "Gray"
        @fact colordim(Bf) --> 0
        Bf = grayim(B)
        @fact eltype(Bf) --> UFixed16
        # colorspace encoded as a Color (enables multiple dispatch)
        BfCV = reinterpret(Gray{UFixed8}, round(UInt8, B))
        @fact colorspace(BfCV) --> "Gray"
        @fact colordim(BfCV) --> 0
        Bf3 = grayim(reshape(convert(UInt8,1):convert(UInt8,36), 3,4,3))
        @fact eltype(Bf3) --> UFixed8
        Bf3 = grayim(reshape(convert(UInt16,1):convert(UInt16,36), 3,4,3))
        @fact eltype(Bf3) --> UFixed16
        Bf3 = grayim(reshape(1.0f0:36.0f0, 3,4,3))
        @fact eltype(Bf3) --> Float32
    end

    context("Colorim") do
        C = colorim(rand(UInt8, 3, 5, 5))
        @fact eltype(C) --> RGB{UFixed8}
        @fact colordim(C) --> 0
        @fact colorim(C) --> C
        C = colorim(rand(UInt16, 4, 5, 5), "ARGB")
        @fact eltype(C) --> ARGB{UFixed16}
        C = colorim(rand(1:20, 3, 5, 5))
        @fact eltype(C) --> Int
        @fact colordim(C) --> 1
        @fact colorspace(C) --> "RGB"
        @fact eltype(colorim(rand(UInt16, 3, 5, 5))) --> RGB{UFixed16}
        @fact eltype(colorim(rand(3, 5, 5))) --> RGB{Float64}
        @fact colordim(colorim(rand(UInt8, 5, 5, 3))) --> 3
        @fact spatialorder(colorim(rand(UInt8, 3, 5, 5))) --> ["x", "y"]
        @fact spatialorder(colorim(rand(UInt8, 5, 5, 3))) --> ["y", "x"]
        @fact eltype(colorim(rand(UInt8, 4, 5, 5), "RGBA")) --> RGBA{UFixed8}
        @fact eltype(colorim(rand(UInt8, 4, 5, 5), "ARGB")) --> ARGB{UFixed8}
        @fact colordim(colorim(rand(UInt8, 5, 5, 4), "RGBA")) --> 3
        @fact colordim(colorim(rand(UInt8, 5, 5, 4), "ARGB")) --> 3
        @fact spatialorder(colorim(rand(UInt8, 4, 5, 5), "ARGB")) --> ["x", "y"]
        @fact spatialorder(colorim(rand(UInt8, 5, 5, 4), "ARGB")) --> ["y", "x"]
        @fact_throws ErrorException colorim(rand(UInt8, 3, 5, 3))
        @fact_throws ErrorException colorim(rand(UInt8, 4, 5, 5))
        @fact_throws ErrorException colorim(rand(UInt8, 5, 5, 4))
        @fact_throws ErrorException colorim(rand(UInt8, 4, 5, 4), "ARGB")
        @fact_throws ErrorException colorim(rand(UInt8, 5, 5, 5), "foo")
        @fact_throws ErrorException colorim(rand(UInt8, 2, 2, 2), "bar")
    end

    context("Indexed color") do
        local cmap = linspace(RGB(0x0cuf8, 0x00uf8, 0x00uf8), RGB(0xffuf8, 0x00uf8, 0x00uf8), 20)
        local img_ = ImageCmap(copy(B), cmap, Dict{ASCIIString, Any}([("spatialorder", Images.yx)]))
        @fact colorspace(img_) --> "RGB"
        img_ = ImageCmap(copy(B), cmap, spatialorder=Images.yx)
        @fact colorspace(img_) --> "RGB"
        # Note: img from opening of facts() block
        # TODO: refactor whole block
        @fact eltype(img) --> RGB{UFixed8}
        @fact eltype(imgd) --> RGB{UFixed8}
    end

    context("Basic information") do
        @fact size(img) --> (3,5)
        @fact size(imgd) --> (3,5)
        @fact ndims(img) --> 2
        @fact ndims(imgd) --> 2
        @fact size(img,"y") --> 3
        @fact size(img,"x") --> 5
        @fact strides(img) --> (1,3)
        @fact strides(imgd) --> (1,3)
        @fact strides(imgds) --> (1,3,15)
        @fact nimages(img) --> 1
        @fact ncolorelem(img) --> 1
        @fact ncolorelem(imgd) --> 1
        @fact ncolorelem(imgds) --> 3
        @fact sort(vec(img)).data[1] --> minimum(img)  # resize! necessary
        @fact sort(vec(imgds)).data[1] --> minimum(imgds)
    end

    context("1-dimensional images") do
        @fact colordim(1:5) --> 0
        @fact nimages(1:5) --> 1
        @fact colorspace(1:5) --> "Gray"
        img1 = Image(map(Gray, 0.1:0.1:0.5), spatialorder=["z"])
        @fact colordim(img1) --> 0
        @fact colorspace(img1) --> "Gray"
        @fact sdims(img1) --> 1
        @fact convert(Vector{Gray{Float64}}, img1) --> map(Gray, 0.1:0.1:0.5)
        @fact size(img1, "z") --> 5
        @fact_throws ErrorException size(img1, "x")
    end

    context("Printing") do
        iob = IOBuffer()
        show(iob, img)
        show(iob, imgd)
    end

    context("Copy / similar") do
        A = randn(3,5,3)
        imgc = @inferred(copy(img))
        @fact imgc.data --> img.data
        imgc = copyproperties(imgd, A)
        @fact imgc.data --> A
        img2 = @inferred(similar(img))
        @fact isa(img2, ImageCmap) --> true
        @fact (img2.data == img.data) --> false
        img2 = @inferred(similar(imgd))
        @fact isa(img2, Image) --> true
        img2 = similar(img, (4,4))
        @fact isa(img2, ImageCmap) --> true
        @fact size(img2) --> (4,4)
        img2 = similar(imgd, (3,4,4))
        @fact isa(img2, Image) --> true
        @fact size(img2) --> (3,4,4)
        @fact copyproperties(B, A) --> A
        @fact shareproperties(A, B) --> B
        @fact shareproperties(img, B) --> B
    end

    context("Getindex / setindex!") do
        prev = img[4]
        @fact prev --> B[4]
        img[4] = prev+1
        @fact img.data[4] --> prev+1
        @fact img[4] --> prev+1
        @fact img[1,2] --> prev+1
        img[1,2] = prev
        @fact img[4] --> prev
    end

    context("Properties") do
        @fact colorspace(img) --> "RGB"
        @fact colordim(img) --> 0
        @fact colordim(imgds) --> 3
        @fact timedim(img) --> 0
        @fact pixelspacing(img) --> [2.0, 3.0]
        @fact spacedirections(img) --> Vector{Float64}[[2.0, 0], [0, 3.0]]
        if testing_units
            @fact pixelspacing(imgd) --> [2.0mm, 3.0mm]
            @fact spacedirections(imgd) -->
                Vector{SIUnits.SIQuantity{Float64,1,0,0,0,0,0,0}}[[2.0mm, 0.0mm], [0.0mm, 3.0mm]]
        end
    end

    context("Dims and ordering") do
        @fact sdims(img) --> sdims(imgd)
        @fact coords_spatial(img) --> coords_spatial(imgd)
        @fact size_spatial(img) --> size_spatial(imgd)
        A = randn(3,5,3)
        tmp = Image(A, Dict{ASCIIString,Any}())
        copy!(tmp, imgd, "spatialorder")
        @fact properties(tmp) --> Dict{ASCIIString,Any}([("spatialorder",Images.yx)])
        copy!(tmp, imgd, "spatialorder", "pixelspacing")
        if testing_units
            @fact tmp["pixelspacing"] --> [2.0mm, 3.0mm]
        end

        @fact storageorder(img) --> Images.yx
        @fact storageorder(imgds) --> [Images.yx; "color"]

        A = rand(4,4,3)
        @fact colordim(A) --> 3
        Aimg = permutedims(convert(Image, A), [3,1,2])
        @fact colordim(Aimg) --> 1
    end

    context("Sub / slice") do
        s = sub(img, 2, 1:4)
        @fact ndims(s) --> 2
        @fact sdims(s) --> 2
        @fact size(s) --> (1,4)
        @fact data(s) --> sub(B, 2, 1:4)
        if scalar_getindex_new
            s = getindexim(img, 2, 1:4)
            @fact ndims(s) --> 1
            @fact sdims(s) --> 1
            @fact size(s) --> (4,)
            @fact data(s) --> B[2, 1:4]
            s = getindexim(img, 2:2, 1:4)
            @fact ndims(s) --> 2
            @fact sdims(s) --> 2
            @fact size(s) --> (1,4)
            @fact data(s) --> B[2:2, 1:4]
        else
            s = getindexim(img, 2, 1:4)
            @fact ndims(s) --> 2
            @fact sdims(s) --> 2
            @fact size(s) --> (1,4)
            @fact data(s) --> B[2, 1:4]
        end
        s = subim(img, 2, 1:4)
        @fact ndims(s) --> 2
        @fact sdims(s) --> 2
        @fact size(s) --> (1,4)
        s = subim(img, 2, [1,2,4])
        @fact ndims(s) --> 2
        @fact sdims(s) --> 2
        @fact size(s) --> (1,3)
        s = sliceim(img, 2, 1:4)
        @fact ndims(s) --> 1
        @fact sdims(s) --> 1
        @fact size(s) --> (4,)
        s = sliceim(img, 2, [1,2,4])
        @fact ndims(s) --> 1
        @fact sdims(s) --> 1
        @fact size(s) --> (3,)
        s = sliceim(imgds, 2, 1:4, 1:3)
        @fact ndims(s) --> 2
        @fact sdims(s) --> 1
        @fact colordim(s) --> 2
        @fact spatialorder(s) --> ["x"]
        s = sliceim(imgds, 2, :, 1:3)
        @fact ndims(s) --> 2
        @fact sdims(s) --> 1
        @fact colordim(s) --> 2
        @fact spatialorder(s) --> ["x"]
        s = sliceim(imgds, 2:2, 1:4, 1:3)
        @fact ndims(s) --> 3
        @fact sdims(s) --> 2
        @fact colordim(s) --> 3
        @fact colorspace(s) --> "RGB"
        s = getindexim(imgds, 2:2, 1:4, 2)
        @fact ndims(s) --> 2
        @fact sdims(s) --> 2
        @fact colordim(s) --> 0
        @fact colorspace(s) --> "Unknown"
        s = sliceim(imgds, 2:2, 1:4, 2)
        @fact ndims(s) --> 2
        @fact sdims(s) --> 2
        @fact colordim(s) --> 0
        @fact colorspace(s) --> "Unknown"
        s = sliceim(imgds, 2:2, 1:4, 1:2)
        @fact ndims(s) --> 3
        @fact sdims(s) --> 2
        @fact colordim(s) --> 3
        @fact colorspace(s) --> "Unknown"
        @fact spatialorder(s) --> ["y","x"]
        s = sub(img, "y", 2)
        @fact ndims(s) --> 2
        @fact sdims(s) --> 2
        @fact size(s) --> (1,5)
        s = slice(img, "y", 2)
        @fact ndims(s) --> 1
        @fact size(s) --> (5,)
        @fact size(getindexim(imgds, :, 1:2, :)) --> (size(imgds,1), 2, 3)

        s = permutedims(imgds, (3,1,2))
        @fact colordim(s) --> 1
        ss = getindexim(s, 2:2, :, :)
        @fact colorspace(ss) --> "Unknown"
        @fact colordim(ss) --> 1
        sss = squeeze(ss, 1)
        @fact colorspace(ss) --> "Unknown"
        @fact colordim(sss) --> 0
        ss = getindexim(imgds, 2:2, :, :)
        @fact colordim(ss) --> 3
        @fact spatialorder(ss) --> ["y", "x"]
        sss = squeeze(ss, 1)
        @fact colordim(sss) --> 2
        @fact spatialorder(sss) --> ["x"]
    end

# # reslicing
# D = randn(3,5,4)
# sd = SliceData(D, 2)
# C = slice(D, sd, 2)
# @fact C --> reshape(D[1:end, 2, 1:end], size(C))
# reslice!(C, sd, 3)
# @fact C --> reshape(D[1:end, 3, 1:end], size(C))
# sd = SliceData(D, 3)
# C = slice(D, sd, 2)
# @fact C --> reshape(D[1:end, 1:end, 2], size(C))
#
# sd = SliceData(imgds, 2)
# s = sliceim(imgds, sd, 2)
# @fact colordim(s) --> 2
# @fact colorspace(s) --> "RGB"
# @fact spatialorder(s) --> ["y"]
# @fact s.data --> reshape(imgds[:,2,:], size(s))
# sd = SliceData(imgds, 3)
# s = sliceim(imgds, sd, 2)
# @fact colordim(s) --> 0
# @fact colorspace(s) --> "Unknown"
# @fact spatialorder(s) --> Images.yx
# @fact s.data --> imgds[:,:,2]
# reslice!(s, sd, 3)
# @fact s.data --> imgds[:,:,3]

    context("Named indexing") do
        @fact dimindex(imgds, "color") --> 3
        @fact dimindex(imgds, "y") --> 1
        @fact dimindex(imgds, "z") --> 0
        imgdp = permutedims(imgds, [3,1,2])
        @fact dimindex(imgdp, "y") --> 2
        @fact coords(imgds, "x", 2:4) --> (1:3, 2:4, 1:3)
        @fact coords(imgds, x=2:4, y=2:3) --> (2:3, 2:4, 1:3)
        @fact img["y", 2, "x", 4] --> B[2,4]
        @fact img["x", 4, "y", 2] --> B[2,4]
        chan = imgds["color", 2]
        Blookup = reshape(green(cmap[B[:]]), size(B))
        @fact chan --> Blookup

        sd = SliceData(imgds, "x")
        s = sliceim(imgds, sd, 2)
        @fact spatialorder(s) --> ["y"]
        @fact s.data --> reshape(imgds[:,2,:], size(s))
        sd = SliceData(imgds, "y")
        s = sliceim(imgds, sd, 2)
        @fact spatialorder(s) --> ["x"]
        @fact s.data --> reshape(imgds[2,:,:], size(s))
        sd = SliceData(imgds, "x", "y")
        s = sliceim(imgds, sd, 2, 1)
        @fact s.data --> reshape(imgds[1,2,:], 3)
    end

    context("Spatial order, width/ height, and permutations") do
        @fact spatialpermutation(Images.yx, imgds) --> [1,2]
        @fact widthheight(imgds) --> (5,3)
        C = convert(Array, imgds)
        @fact C --> imgds.data
        imgds["spatialorder"] = ["x", "y"]
        @fact spatialpermutation(Images.xy, imgds) --> [1,2]
        @fact widthheight(imgds) --> (3,5)
        C = convert(Array, imgds)
        @fact C --> permutedims(imgds.data, [2,1,3])
        imgds.properties["spatialorder"] = ["y", "x"]
        @fact spatialpermutation(Images.xy, imgds) --> [2,1]
        imgds.properties["spatialorder"] = ["x", "L"]
        @fact spatialpermutation(Images.xy, imgds) --> [1,2]
        imgds.properties["spatialorder"] = ["L", "x"]
        @fact spatialpermutation(Images.xy, imgds) --> [2,1]
        A = randn(3,5,3)
        @fact spatialpermutation(Images.xy, A) --> [2,1]
        @fact spatialpermutation(Images.yx, A) --> [1,2]

        imgds.properties["spatialorder"] = Images.yx
        imgp = permutedims(imgds, ["x", "y", "color"])
        @fact imgp.data --> permutedims(imgds.data, [2,1,3])
        imgp = permutedims(imgds, ("color", "x", "y"))
        @fact imgp.data --> permutedims(imgds.data, [3,2,1])
        if testing_units
            @fact pixelspacing(imgp) --> [3.0mm, 2.0mm]
        end
        imgc = copy(imgds)
        imgc["spacedirections"] = spacedirections(imgc)
        delete!(imgc, "pixelspacing")
        imgp = permutedims(imgc, ["x", "y", "color"])
        if testing_units
            @fact spacedirections(imgp) --> Vector{SIUnits.SIQuantity{Float64,1,0,0,0,0,0,0}}[[0.0mm, 3.0mm],[2.0mm, 0.0mm]]
            @fact pixelspacing(imgp) --> [3.0mm, 2.0mm]
        end
    end

    context("Reinterpret, separate, more convert") do
        a = RGB{Float64}[RGB(1,1,0)]
        af = reinterpret(Float64, a)
        @fact vec(af) --> [1.0,1.0,0.0]
        @fact size(af) --> (3,1)
        @fact_throws ErrorException reinterpret(Float32, a)
        anew = reinterpret(RGB, af)
        @fact anew --> a
        anew = reinterpret(RGB, vec(af))
        @fact anew[1] --> a[1]
        @fact ndims(anew) --> 0
        anew = reinterpret(RGB{Float64}, af)
        @fact anew --> a
        @fact_throws ErrorException reinterpret(RGB{Float32}, af)
        Au8 = rand(0x00:0xff, 3, 5, 4)
        A8 = reinterpret(UFixed8, Au8)
        rawrgb8 = reinterpret(RGB, A8)
        @fact eltype(rawrgb8) --> RGB{UFixed8}
        @fact reinterpret(UFixed8, rawrgb8) --> A8
        @fact reinterpret(UInt8, rawrgb8) --> Au8
        rawrgb32 = convert(Array{RGB{Float32}}, rawrgb8)
        @fact eltype(rawrgb32) --> RGB{Float32}
        @fact ufixed8(rawrgb32) --> rawrgb8
        @fact reinterpret(UFixed8, rawrgb8) --> A8
        imrgb8 = convert(Image, rawrgb8)
        @fact spatialorder(imrgb8) --> Images.yx
        @fact convert(Image, imrgb8) --> exactly(imrgb8)
        @fact convert(Image{RGB{UFixed8}}, imrgb8) --> exactly(imrgb8)
        im8 = reinterpret(UFixed8, imrgb8)
        @fact data(im8) --> A8
        @fact permutedims(reinterpret(UFixed8, separate(imrgb8)), (3, 1, 2)) --> im8
        @fact reinterpret(UInt8, imrgb8) --> Au8
        @fact reinterpret(RGB, im8) --> imrgb8
        ims8 = separate(imrgb8)
        @fact colordim(ims8) --> 3
        @fact colorspace(ims8) --> "RGB"
        @fact convert(Image, ims8) --> exactly(ims8)
        @fact convert(Image{UFixed8}, ims8) --> exactly(ims8)
        @fact separate(ims8) --> exactly(ims8)
        imrgb8_2 = convert(Image{RGB}, ims8)
        @fact isa(imrgb8_2, Image{RGB{UFixed8}}) --> true
        @fact imrgb8_2 --> imrgb8
        A = reinterpret(UFixed8, UInt8[1 2; 3 4])
        imgray = convert(Image{Gray{UFixed8}}, A)
        @fact spatialorder(imgray) --> Images.yx
        @fact data(imgray) --> reinterpret(Gray{UFixed8}, [0x01 0x02; 0x03 0x04])
        @fact eltype(convert(Image{HSV{Float32}}, imrgb8)) --> HSV{Float32}
        @fact eltype(convert(Image{HSV}, float32(imrgb8))) --> HSV{Float32}

        @fact eltype(convert(Array{Gray}, imrgb8)) --> Gray{U8}
        @fact eltype(convert(Image{Gray}, imrgb8)) --> Gray{U8}
        @fact eltype(convert(Array{Gray}, data(imrgb8))) --> Gray{U8}
        @fact eltype(convert(Image{Gray}, data(imrgb8))) --> Gray{U8}
        # Issue 232
        local img = Image(reinterpret(Gray{UFixed16}, rand(UInt16, 5, 5)))
        imgs = subim(img, :, :)
        @fact isa(minfinite(imgs), UFixed16) --> true
        # Raw
        imgdata = rand(UInt16, 5, 5)
        img = Image(reinterpret(Gray{UFixed16}, imgdata))
        @fact all(raw(img) .== imgdata) --> true
        @fact typeof(raw(img)) --> Array{UInt16,2}
        @fact typeof(raw(Image(rawrgb8))) --> Array{UInt8,3}  # check color images
        @fact size(raw(Image(rawrgb8))) --> (3,5,4)
        @fact typeof(raw(imgdata)) --> typeof(imgdata)  # check array fallback
        @fact all(raw(imgdata) .== imgdata) --> true
    end
end


text:
# SECTION 4.1

using NetworkDesignModels, Gurobi, JuMP
using DataFrames, JLD2

# load data
env = Gurobi.Env()
results = DataFrame(
    budget = Float64[], niterations = Int[], runtime = Float64[],
    line = Int[], order = Int[], stop_id = String[], freq = Float64[],
    pctxfr = Float64[], obj = Float64[])
JLD2.@load "data/transit-network-toy.jld2" stns np

# create master problem
options = NetworkDesignModels.MasterOptions(
    nlegs = 2,
    nfreqs = 2, costwts = [1.0, 0.5], freqwts = [1.0, 0.75], xfrwts = [1.0, 0.25])
rmp = NetworkDesignModels.MasterProblem(
    np,
    linelist = Vector{Int}[],
    options = options,
    solver = GurobiSolver(env, OutputFlag = 0)
);
scale = 1.0
budget = ceil(scale*sum(NetworkDesignModels.linecost(np,line) for line in np.lines))
NetworkDesignModels.optimize(rmp, budget)

# create subproblem objects
sp1 = NetworkDesignModels.SubProblemCP(rmp, maxdist = 1.45, traveltimes = true,
    solver = GurobiSolver(env, OutputFlag = 0))
sp2 = NetworkDesignModels.SubProblemCP(rmp, nlegs = 1, maxdist = 1.45, traveltimes = true,
    solver = GurobiSolver(env, OutputFlag = 0))

# generate columns
t0 = time()
nit = 0
for i in 1:30

    global nit += 1

    # solve master
    NetworkDesignModels.optimize(rmp, budget)
    println(getobjectivevalue(rmp.model))

    # create subproblem
    global sp1 = NetworkDesignModels.SubProblemCP(rmp, maxdist = 1.45, traveltimes = true,
        solver = GurobiSolver(env, OutputFlag = 0))
    path = NetworkDesignModels.generatecolumn(sp1, rmp)

    if length(path) > 0 && !NetworkDesignModels.linein(path, rmp.linelist)
        println("$(i) adding high frequency column")
        for k1 in 1:(length(path)-1), k2 in (k1+1):length(path)
            if !NetworkDesignModels.linein(path[k1:k2], rmp.linelist)
                NetworkDesignModels.addcolumn!(rmp, path[k1:k2])
            end
        end
    elseif length(path) > 0
        error("Bad column")
        break
    else
        global sp2 = NetworkDesignModels.SubProblemCP(rmp, nlegs = 1, maxdist = 1.45, traveltimes = true,
            solver = GurobiSolver(env, OutputFlag = 0))
        path = NetworkDesignModels.generatecolumn(sp2, rmp, f=2)
        if length(path) > 0 && !NetworkDesignModels.linein(path, rmp.linelist)
            println("$(i) adding low frequency column")
            for k1 in 1:(length(path)-1), k2 in (k1+1):length(path)
                if !NetworkDesignModels.linein(path[k1:k2], rmp.linelist)
                    NetworkDesignModels.addcolumn!(rmp, path[k1:k2])
                end
            end
        elseif length(path) > 0
            error("bad column ", path)
            break
        else
            break
        end
    end
end
t0 = time() - t0

# solve IP
options.modeltype = :ip
rmpip = NetworkDesignModels.MasterProblem(
    np,
    linelist = rmp.linelist,
    options = options,
    solver = GurobiSolver(env, OutputFlag = 0)
);
NetworkDesignModels.optimize(rmpip, budget)
println(getobjectivevalue(rmpip.model))

# save data
twolines = 0
total = 0
for (u,v) in NetworkDesignModels.commutes(np)
    uv1 = getvalue(rmpip.θ[(u,v)])
    if length(np.xfrstns[(u,v)]) > 0
        uv2 = getvalue(sum(rmpip.model[:twoline][(u,v),w] for w in np.xfrstns[(u,v)]))
        global twolines += uv2
    end
    global total += uv1
end
println(twolines/total)

finallist = Vector{Vector{Int}}()
for l in 1:length(rmpip.linelist), f in 1:options.nfreqs
    if getvalue(rmpip.x[l,f]) > 1e-3
        println(f, "\t", rmpip.linelist[l])
        push!(finallist, rmpip.linelist[l])
        for j in 1:length(rmpip.linelist[l])
            push!(results, [budget, nit, t0,
                l, j, stns[rmpip.linelist[l][j]], f,
                twolines/total,
                getobjectivevalue(rmpip.model)])
        end
    end
end

ext: jl


text:
include("Models.jl")
include("Scores.jl")

import Random.randperm
import Random.shuffle!
using Statistics
using Optim

function splitn(x::AbstractArray, n::Int)
    shuffle!(x)
    s = length(x) / n
    [x[round(Int64, (i-1)*s)+1:min(length(x),round(Int64, i*s))] for i=1:n]
end

function Kfold(cartesianinds::Array, setting::Tuple, nfolds)
    if setting == Tuple([false for i in 1:length(setting)]) #setting A
        buckets = splitn(cartesianinds, nfolds)
    elseif setting == Tuple([true for i in 1:length(setting)])
        throw("kfolding for this setting not yet implemented")
    else
        d = Dict() # keys: contain indices where to freely select, values: all observed CartesianIndices according to a key
        for I in cartesianinds
            key = [I[i] for i in 1:length(setting) if setting[i]==false]
            if key in keys(d)
                push!(d[key], I)
            else
                d[key]=[I]
            end
        end
        buckets_keys = splitn(collect(keys(d)), nfolds)
        buckets = [[] for i in 1:length(buckets_keys)]
        for i in 1:length(buckets)
            for j in 1:length(buckets_keys[i])
                buckets[i] = cat(buckets[i], d[buckets_keys[i][j]], dims=1)
            end
        end
    end
    folds = []
    for i in 1:nfolds
        a = (train, test) = vcat(buckets[1:end .!= i]...), buckets[i]
        push!(folds,a)
    end
    return folds

end
function CVestimate(model::MultilinearKroneckerModel, kernels, Y, inds, folding, setting::Tuple)
    #determine the folds
    if folding == "LOO"
        return predict_LOO(model, kernels, Y, inds, setting)
    elseif typeof(folding) <: Int
        folds = Kfold(inds, setting, folding)
    elseif typeof(folding) <: AbstractArray
        folds = folding
    else
        throw("Input for nfolds not valid")
    end
    Ycvest = Array{Float64}(undef, size(Y))
    for fold in folds
        traininds, testinds = fold
        model.A_[:]=0.001*randn(size(model.A_))[:]
        fit(model, kernels, Y, traininds)
        Yest = predict(model, kernels)
        for I in testinds
            Ycvest[I] = Yest[I]
        end
    end
    return Ycvest
end

function CVscore(model::MultilinearKroneckerModel, kernels, Y, inds, folding, setting::Tuple, score)
    Ycvest = CVestimate(model, kernels, Y, inds, folding, setting)
    Y = [Y[inds[i]] for i in 1:length(inds)]
    Ycvest = [Ycvest[inds[i]] for i in 1:length(inds)]
    return score(Ycvest, Y)
end

function CVscore(model::MultilinearKroneckerModel, kernels, Y, inds, folding, setting::Tuple, score, times)
    scores = []
    for i in 1:times
        sc = CVscore(model, kernels, Y, inds, folding, setting, score)
        push!(scores, sc)
    end
    return mean(scores),3*std(scores)
end

# when all hyperparameters the same:
function optimizeHyperParameters(model::MultilinearKroneckerModel, kernels, Y, inds, folding, setting::Tuple, score)
    lower = 0.001
    upper = 10000
    optimizer = Brent()
    optimizeHyperParameters(model, kernels, Y, inds, folding, setting, score, lower, upper, optimizer)
end
function optimizeHyperParameters(model::MultilinearKroneckerModel, kernels, Y, inds, folding, setting::Tuple, score, lower, upper, optimizer)
    function objective(λ, model::MultilinearKroneckerModel, kernels, Y, inds, folding, setting::Tuple, score)
        model.λ_[:]=λ[:]
        return -CVscore(model, kernels, Y, inds, folding, setting, score)
    end
    opt = optimize(λ -> objective([λ for i in 1:length(model.λ_)],model, kernels,Y, inds,folding,setting,score), lower, upper, optimizer, store_trace=true)
    model.λ_[:]=[opt.minimizer for i in 1:length(model.λ_)]
    return opt
end

#when differetn hyperparemeters
function optimizeHyperParameters(model::MultilinearKroneckerModel, kernels, Y, inds, folding, setting::Tuple, score, optimizer, init)
    function objective(λ, model::MultilinearKroneckerModel, kernels, Y, inds, folding, setting::Tuple, score)
        model.λ_[:]=λ[:]
        #println(model.λ_)
        return -CVscore(model, kernels, Y, inds, folding, setting, score)
    end
    lower = zeros(length(model.λ_))
    upper = [Inf for i in 1:length(model.λ_)]
    inner_optimizer = optimizer
    opt = optimize(λ -> objective(λ,model, kernels,Y, inds,folding,setting,score), lower, upper, init, Fminbox(inner_optimizer), Optim.Options(store_trace=true))
    model.λ_[:]=opt.minimizer[:]
    return opt
end

max_stars_repo_name: padwulf/MyExample

max_issues_repo_name: padwulf/MyExample

max_forks_repo_name: padwulf/MyExample


text:
#!/usr/bin/env julia
using Plots
#=
This is to plot the results about the imputation results with HD genotypes
I also take this chance to practise using Julia instead of bash for
scripting.
=#
only = "828-only"
vran = "inc8d-vs-random"
ext8 = "extra-8d-vs-random"

#= Figure 1:
Randomly sample 50 ID to be masked
Impute them back with:
    1. 828 data
    2. random 828 from 17k data
and compare
=#
rst = read(pipeline(`grep error $only/rates`,
                    `gawk '{print $NF}'`), String)
dat = reshape([parse(Float64, x) for x in split(rst)], 2, :)'

plot(dat,
     label=["Random", "data 828"],
     dpi=300,
     xlabel="Repeat",
     ylabel="Error rate",
     legend=:left)

savefig("828-vs-ran.png")

#= Figure 2:
Randomly sample 50 ID to be masked
Impute them back with

1. 828+extr random ID as reference
2. the above number of random ID from 17k data as reference

Impute and compare
=#
rst = read(pipeline(`grep error $vran/rates`,
                    `gawk '{print $NF}'`), String)
dat = reshape([parse(Float64, x) for x in split(rst)], 2, :)'
nid = 50+828:200:1500+828

plot(nid, dat,
     label=["Random", "828+(nran-828)"],
     dpi=300,
     xlabel="Reference set size",
     ylabel="Error rate")

savefig("inc828-vs-ran.png")

#= Figure 3:
Randomly sample 50 ID to be masked
Impute them back with

1. random select ID from 17k data
2. using above, plus 828 extra

Imputed and compare
=#
rst = read(pipeline(`grep error $ext8/rates`,
                    `gawk '{print $NF}'`), String)
dat = reshape([parse(Float64, x) for x in split(rst)], 2, :)'
nid = 1000:100:1500

plot(nid, dat,
     label=["Random", "Random + extr 828"],
     dpi=300,
     xlabel="Random set size, (the other has extra 828ID)",
     ylabel="Error rate",
     legend=:bottom)

savefig("inc828-vs-ran.png")

max_stars_repo_name: xijiang/Refactered-NSG

max_issues_repo_name: xijiang/Refactered-NSG

max_forks_repo_name: xijiang/Refactered-NSG


text:
@testset "dfvstack.jl" begin
    using DataFrames
    df1 = DataFrame("t" => 1:10, "f(t)" => randn(10))
    df2 = DataFrame("t" => 5:0.5:13, "f(t)" => randn(17))
    df3 = DataFrame("t" => 14:23, "f(t)" => randn(10))
    @test try dfvstack(["s1", "s2", "s3"] .=> [df1, df2, df3], "t", "f(t)"); true; catch; false; end
    dflong = dfvstack(["s1", "s2", "s3"] .=> [df1, df2, df3], "t", "f(t)")
    gd = groupby(dflong, "sheet")
    @test gd[1][!, :t] == df1[!, :t]
    @test gd[1][!, "f(t)"] == df1[!, "f(t)"]
    @test gd[2][!, :t] == df2[!, :t]
    @test gd[2][!, "f(t)"] == df2[!, "f(t)"]
    @test gd[3][!, :t] ==  df3[!, :t]
    @test gd[3][!, "f(t)"] == df3[!, "f(t)"]
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "xdm"
version = v"2.4.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/xyce/xdm.git", "c87548b0bdd4d696ea103008d452082907951fc3"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/xdm
atomic_patch -p1 ../patches/pyinterp.patch
atomic_patch -p1 ../patches/musl.patch
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DBoost_NO_BOOST_CMAKE=ON
make -j$nproc
make install
# xdm libraries are build with .so extension, because python (for some reason).
# Symlink them to .dylib, so that oour audit system can find and verify them.
if [[ $target == *apple* ]]; then
ln -s SpiritExprCommon.so $prefix/xdm_bundle/SpiritExprCommon.dylib
ln -s XdmRapidXmlReader.so $prefix/xdm_bundle/XdmRapidXmlReader.dylib
ln -s SpiritCommon.so $prefix/xdm_bundle/SpiritCommon.dylib
fi
"""

# No windows python support in Yggdrasil at the moment
platforms = filter(!Sys.iswindows, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("SpiritExprCommon", :SpiritExprCommon, "xdm_bundle"),
    LibraryProduct("XdmRapidXmlReader", :XdmRapidXmlReader, "xdm_bundle"),
    LibraryProduct("SpiritCommon", :SpiritCommon, "xdm_bundle")
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="boostpython_jll", uuid="398de629-0a17-50a6-9837-8b3a70a53854"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using DiffEqJump, DiffEqBase
using Test

rate = (u, p, t) -> u[1]
affect! = function (integrator)
    integrator.u[1] -= 1
    integrator.u[2] += 1
end
jump = ConstantRateJump(rate, affect!)

prob = DiscreteProblem([0.0, 0.0], (0.0, 10.0))
jump_prob = JumpProblem(prob, Direct(), jump)

sol = solve(jump_prob, SSAStepper())

@test sol.t == [0.0, 10.0]
@test sol.u == [[0.0, 0.0], [0.0, 0.0]]

condition(u,t,integrator) = t == 5
function fuel_affect!(integrator)
  integrator.u[1] += 100
  reset_aggregated_jumps!(integrator)
end
cb = DiscreteCallback(condition, fuel_affect!, save_positions=(false, true))

sol = solve(jump_prob, SSAStepper(), callback=cb, tstops=[5])

@test sol.t[1:2] == [0.0, 5.0] # no jumps between t=0 and t=5
@test sol(5 + 1e-10) == [100, 0] # state just after fueling before any decays can happen


text:
module MultivariateStats_

export RidgeRegressor, PCA, KernelPCA, ICA, LDA, BayesianLDA

import MLJModelInterface
import MLJModelInterface: @mlj_model, metadata_pkg, metadata_model,
    Table, Continuous, Finite

const MMI = MLJModelInterface

import StatsBase: CovarianceEstimator
using Distances
using LinearAlgebra

import MultivariateStats

const MS = MultivariateStats

struct LinearFitresult{F} <: MMI.MLJType
    coefficients::Vector{F}
    intercept::F
end

const RIDGE_DESCR = "Ridge regressor with regularization parameter lambda. Learns a linear regression with a penalty on the l2 norm of the coefficients."

const PCA_DESCR = "Principal component analysis. Learns a linear transformation to project the data  on a lower dimensional space while preserving most of the "*
                 "initial variance."

const KPCA_DESCR = "Kernel principal component analysis."

const ICA_DESCR = "Independent component analysis."

const LDA_DESCR = "Multiclass linear discriminant analysis. The algorithm learns a projection matrix `P` that projects a feature matrix `Xtrain` onto a"*
                " lower dimensional space of dimension `out_dim` such that the trace of the transformed between-class scatter matrix(`Pᵀ*Sb*P`) is maximized"*
                " relative to the trace of the transformed within-class scatter matrix (`Pᵀ*Sw*P`).The projection matrix is scaled such that `Pᵀ*Sw*P=I` or "*
                "`Pᵀ*Σw*P=I`(where `Σw` is the within-class covariance matrix) . \nPredicted class posterior probability for feature matrix `Xtest` are derived "*
                "by applying a softmax transformation to a matrix `Pr`, such that  rowᵢ of `Pr` contains computed distances(based on a distance metric) in the "*
                "transformed space of rowᵢ in `Xtest` to the centroid of each class. "

const BayesianLDA_DESCR = "Bayesian Multiclass linear discriminant analysis. The algorithm learns a projection matrix `P` that projects a feature matrix `Xtrain` "*
                        "onto a lower dimensional space of dimension `out_dim` such that the trace of the transformed between-class scatter matrix(`Pᵀ*Sb*P`) is "*
                        "maximized relative to the trace of the transformed within-class scatter matrix (`Pᵀ*Sw*P`). The projection matrix is scaled such that "*
                        "`Pᵀ*Sw*P = n` or `Pᵀ*Σw*P=I` (Where `n` is the number of training samples and `Σw` is the within-class covariance matrix).\nPredicted "*
                        "class posterior probability distibution are derived by applying Bayes rule with a multivariate Gaussian class-conditional distribution."

const BayesianSubspaceLDA_DESCR = "Bayesian Multiclass linear discriminant analysis. Suitable for high dimensional data(Avoids computing scatter matrices"*
                                " `Sw` ,`Sb`). The algorithm learns a projection matrix `P = W*L` (`Sw`), that projects a feature matrix `Xtrain` onto a lower "*
                                "dimensional space of dimension `nc-1` such that the trace of the transformed between-class scatter matrix(`Pᵀ*Sb*P`) is maximized"*
                                " relative to the trace of the transformed within-class scatter matrix (`Pᵀ*Sw*P`). The projection matrix is scaled such that"*
                                " `Pᵀ*Sw*P = mult*I` or `Pᵀ*Σw*P=mult/(n-nc)*I` (where `n` is the number of training samples, `mult` is  one of `n` or `1` "*
                                "depending on whether `Sb` is normalized, `Σw` is the within-class covariance matrix, and `nc` is the number of unique classes "*
                                "in `y`) and also obeys `Wᵀ*Sb*p = λ*Wᵀ*Sw*p`, for every column `p` in `P`. \nPosterior class probability distibution are derived "*
                                "by applying Bayes rule with a multivariate Gaussian class-conditional distribution"
const SubspaceLDA_DESCR = "Multiclass linear discriminant analysis. Suitable for high dimensional data (Avoids computing scatter matrices `Sw` ,`Sb`). The "*
                          "algorithm learns a projection matrix `P = W*L` that projects a feature matrix `Xtrain` onto a lower dimensional space of dimension "*
                          "`nc - 1 ` such that the trace of the transformed between-class scatter matrix(`Pᵀ*Sb*P`) is maximized relative to the trace of the "*
                          "transformed within-class scatter matrix (`Pᵀ*Sw*P`). The projection matrix is scaled such that `Pᵀ*Sw*P = mult*I` or "*
                          "`Pᵀ*Σw*P=mult/(n-nc)*I` (where `n` is the number of training samples, mult` is  one of `n` or `1` depending on whether `Sb` is "*
                          "normalized, `Σw` is the within-class covariance matrix, and `nc` is the number of unique classes in `y`) and also obeys"*
                          " `Wᵀ*Sb*p = λ*Wᵀ*Sw*p`, for every column `p` in `P`.\nPredicted class posterior probability for feature matrix `Xtest` are derived by"*
                          " applying a softmax transformation to a matrix `Pr`, such that  rowᵢ of `Pr` contains computed distances(based on a distance metric)"*
                          " in the transformed space of rowᵢ in `Xtest` to the centroid of each class. "

####
#### RIDGE
####

"""
    RidgeRegressor(; lambda=1.0)

$RIDGE_DESCR

## Parameters

* `lambda=1.0`: non-negative parameter for the regularization strength.
"""
@mlj_model mutable struct RidgeRegressor <: MMI.Deterministic
    lambda::Real = 1.0::(_ ≥ 0)
end

function MMI.fit(model::RidgeRegressor, verbosity::Int, X, y)
    Xmatrix   = MMI.matrix(X)
    features  = MMI.schema(X).names
    θ         = MS.ridge(Xmatrix, y, model.lambda)
    coefs     = θ[1:end-1]
    intercept = θ[end]

    fitresult = LinearFitresult(coefs, intercept)
    report    = NamedTuple()
    cache     = nothing

    return fitresult, cache, report
end

MMI.fitted_params(::RidgeRegressor, fr) =
    (coefficients=fr.coefficients, intercept=fr.intercept)

function MMI.predict(::RidgeRegressor, fr, Xnew)
    Xmatrix = MMI.matrix(Xnew)
    return Xmatrix * fr.coefficients .+ fr.intercept
end

####
#### PCA
####

const PCAFitResultType = MS.PCA

"""
    PCA(; maxoutdim=0, method=:auto, pratio=0.99, mean=nothing)

$PCA_DESCR

## Parameters

* `maxoutdim=0`:       maximum number of output dimensions, unconstrained if 0
* `method=:auto`:      method to use to solve the problem, one of `:auto`,
                       `:cov` or `:svd`
* `pratio=0.99`:       ratio of variance preserved
* `mean=nothing`:      if set to nothing centering will be computed and
                       applied, if set to `0` no centering (assumed pre-
                       centered), if a vector is passed, the centering is done
                       with that vector.
"""
@mlj_model mutable struct PCA <: MMI.Unsupervised
    maxoutdim::Int = 0::(_ ≥ 0)
    method::Symbol  = :auto::(_ in (:auto, :cov, :svd))
    pratio::Float64 = 0.99::(0.0 < _ ≤ 1.0)
    mean::Union{Nothing, Real, Vector{Float64}} = nothing::(_ === nothing || (_ isa Real && iszero(_)) || true)
end

function MMI.fit(model::PCA, verbosity::Int, X)
    Xarray = MMI.matrix(X)
    mindim = minimum(size(Xarray))

    maxoutdim = model.maxoutdim == 0 ? mindim : model.maxoutdim

    # NOTE: copy/transpose
    fitresult = MS.fit(MS.PCA, permutedims(Xarray);
                       method=model.method,
                       pratio=model.pratio,
                       maxoutdim=maxoutdim,
                       mean=model.mean)

    cache = nothing
    report = (indim=MS.indim(fitresult),
              outdim=MS.outdim(fitresult),
              mean=MS.mean(fitresult),
              principalvars=MS.principalvars(fitresult),
              tprincipalvar=MS.tprincipalvar(fitresult),
              tresidualvar=MS.tresidualvar(fitresult),
              tvar=MS.tvar(fitresult))

    return fitresult, cache, report
end

MMI.fitted_params(::PCA, fr) = (projection=MS.projection(fr),)

function MMI.transform(::PCA, fr::PCAFitResultType, X)
    # X is n x d, need to transpose and copy twice...
    Xarray = MMI.matrix(X)
    Xnew   = permutedims(MS.transform(fr, permutedims(Xarray)))
    return MMI.table(Xnew, prototype=X)
end

####
#### KernelPCA
####

const KernelPCAFitResultType = MS.KernelPCA

default_kernel = (x, y) -> x'y

"""
    KernelPCA(; maxoutdim=nothing, kernel=(x,y)->x'y, solver=:auto, pratio=0.99, mean=nothing)

$KPCA_DESCR

## Parameters

* `maxoutdim=0`:       maximum number of output dimensions, unconstrained if 0
* `kernel=nothing`:    kernel function of 2 vector arguments x and y, returns a
                       scalar value, (x,y)->x'y if nothing
* `solver=:auto`:      solver to use for the eigenvalues, one of `:eig`
                       (default), `:eigs`
* `inverse=false`:     perform calculation for inverse transform
* `beta=1.0`:          strength of the ridge regression that learns the inverse
                       transform when inverse is true
* `tol=0.0`:           Convergence tolerance for eigs solver
* `maxiter=300`:       maximum number of iterations for eigs solver
"""
@mlj_model mutable struct KernelPCA <: MMI.Unsupervised
    maxoutdim::Int = 0::(_ ≥ 0)
    kernel::Function = default_kernel
    solver::Symbol   = :eig::(_ in (:eig, :eigs))
    inverse::Bool    = false
    beta::Real       = 1.0::(_ ≥ 0.0)
    tol::Real        = 1e-6::(_ ≥ 0.0)
    maxiter::Int     = 300::(_ ≥ 1)
end

function MMI.fit(model::KernelPCA, verbosity::Int, X)
    Xarray = MMI.matrix(X)
    mindim = minimum(size(Xarray))
    # default max out dim if not given
    maxoutdim = model.maxoutdim == 0 ? mindim : model.maxoutdim

    fitresult = MS.fit(MS.KernelPCA, permutedims(Xarray);
                       kernel=model.kernel,
                       maxoutdim=maxoutdim,
                       solver=model.solver,
                       inverse=model.inverse,
                       β=model.beta,
                       tol=model.tol,
                       maxiter=model.maxiter)

    cache  = nothing
    report = (indim=MS.indim(fitresult),
              outdim=MS.outdim(fitresult),
              projection=MS.projection(fitresult),
              principalvars=MS.principalvars(fitresult))

    return fitresult, cache, report
end

MMI.fitted_params(::KernelPCA, fr) = (projection=fr,)

function MMI.transform(::KernelPCA, fr::KernelPCAFitResultType, X)
    # X is n x d, need to transpose and copy twice...
    Xarray = MMI.matrix(X)
    Xnew   = permutedims(MS.transform(fr, permutedims(Xarray)))
    return MMI.table(Xnew, prototype=X)
end

####
#### ICA
####

const ICAFitResultType = MS.ICA

"""
    ICA(; k=0, alg=:fastica, fun=:tanh, do_whiten=true, maxiter=100,
      tol=1e-6, mean=nothing, winit=zeros(0,0))

$ICA_DESCR

## Parameters

* `k=0`:              number of independent components to recover, set
                      automatically if `0`
* `alg=:fastica`:     algorithm to use (only `:fastica` is supported at the
                      moment)
* `fun=:tanh`:        approximate neg-entropy functor, via the function
                      `MultivariateStats.icagfun`, one of `:tanh` and `:gaus`
* `do_whiten=true`:   whether to perform pre-whitening
* `maxiter=100`:      maximum number of iterations
* `tol=1e-6`:         convergence tolerance for change in matrix W
* `mean=nothing`:     mean to use, if nothing centering is computed and
                      applied, if zero, no centering, a vector of means can be
                      passed
* `winit=zeros(0,0)`: initial guess for matrix `W` either an empty matrix
                      (random init), a matrix of size `k × k` (if `do_whiten`
                      is true), a matrix of size `m × k` otherwise
"""
@mlj_model mutable struct ICA <: MMI.Unsupervised
    k::Int          = 0::(_ ≥ 0)
    alg::Symbol     = :fastica::(_ in (:fastica,))
    fun::Symbol     = :tanh::(_ in (:tanh, :gaus))
    do_whiten::Bool = true
    maxiter::Int    = 100
    tol::Real       = 1e-6::(_ ≥ 0.0)
    winit::Union{Nothing,Matrix{<:Real}} = nothing
    mean::Union{Nothing, Real, Vector{Float64}} = nothing::(_ === nothing || (_ isa Real && iszero(_)) || true)
end

function MMI.fit(model::ICA, verbosity::Int, X)
    Xarray = MMI.matrix(X)
    n, p   = size(Xarray)

    m = min(n, p)
    k = ifelse(model.k ≤ m, model.k, m)

    fitresult = MS.fit(MS.ICA, permutedims(Xarray), k;
                       alg=model.alg,
                       fun=MS.icagfun(model.fun),
                       do_whiten=model.do_whiten,
                       maxiter=model.maxiter,
                       tol=model.tol,
                       mean=model.mean,
                       winit=ifelse(model.winit === nothing, zeros(0,0), model.winit))

    cache = nothing
    report = (indim=MS.indim(fitresult),
              outdim=MS.outdim(fitresult),
              mean=MS.mean(fitresult))

    return fitresult, cache, report
end

MMI.fitted_params(::ICA, fr) = (projection=fr,)

function MMI.transform(::ICA, fr::ICAFitResultType, X)
    # X is n x d, need to transpose and copy twice...
    Xarray = MMI.matrix(X)
    Xnew   = permutedims(MS.transform(fr, permutedims(Xarray)))
    return MMI.table(Xnew, prototype=X)
end

####
#### MulticlassLDA
####

"""
    LDA(; kwargs...)

$LDA_DESCR

## Parameters

* `method=:gevd`:      choice of solver, one of `:gevd` or `:whiten` methods
* `cov_w=SimpleCovariance()`: an estimator for the within-class covariance (used in computing
                              within-class scatter matrix, Sw), by
                              default set to the standard
                              `MultivariateStats.CovarianceEstimator` but could
                              be set to any robust estimator from `CovarianceEstimation.jl`.
* `cov_b=SimpleCovariance()`: same as `cov_w` but for the between-class
                              covariance(used in computing between-class scatter matrix, Sb).
* `out_dim`:          the output dimension, i.e dimension of the transformed
                      space, automatically set if 0 is given (default).
* `regcoef`:          regularization coefficient (default value 1e-6). A
                      positive value `regcoef * eigmax(Sw)` where `Sw` is the
                      within-class scatter matrix, is added to the
                      diagonal of Sw to improve numerical stability. This can
                      be useful if using the standard covariance estimator.
* `dist=SqEuclidean`: the distance metric to use when performing classification
                      (to compare the distance between a new point and
                      centroids in the transformed space), an alternative
                      choice can be the `CosineDist`.

See also the [package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Li, Zhu and Ogihara, [Using Discriminant Analysis for Multi-class Classification: An Experimental Investigation](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7068&rep=rep1&type=pdf).
"""
@mlj_model mutable struct LDA <: MMI.Probabilistic
    method::Symbol   = :gevd::(_ in (:gevd, :whiten))
    cov_w::CovarianceEstimator = MS.SimpleCovariance()
    cov_b::CovarianceEstimator = MS.SimpleCovariance()
    out_dim::Int     = 0::(_ ≥ 0)
    regcoef::Float64    = 1e-6::(_ ≥ 0)
    dist::SemiMetric = SqEuclidean()
end

function MMI.fit(model::LDA, ::Int, X, y)
    class_list = MMI.classes(y[1]) #class list containing entries in pool of y
    classes_seen  = filter(in(unique(y)), class_list) #class list containing unique entries in y
    nc   = length(classes_seen) #number of classes in pool of y
    nclasses = length(class_list)
    nc == nclasses || (integers_seen = MMI.int(classes_seen))

    # NOTE: copy/transpose
    Xm_t   = MMI.matrix(X, transpose=true) # now p x n matrix
    yplain = MMI.int(y) # vector of n ints in {1,..., nclasses}
    p, n   = size(Xm_t)
    #recode yplain to be in {1,..., nc}
    nc == nclasses ||  replace!(yplain, (integers_seen .=> 1:nc)...)

    #check to make sure we have more than one class in training sample
    # This is to prevent Sb from being a zero matrix
    nc >= 2 || throw(ArgumentError("The number of unique classes in traning sample"*
                                        " `nc` has to be greater than one"))

    #check to make sure we have more samples than classes
    # This is to prevent Sw from being the zero matrix
    n > nc || throw(ArgumentError("The number of training samples `n` has to be"*
                                    " greater than the number of unique classes `nc`"))

    # check output dimension default is min(p, nc-1)
    def_outdim = min(p, nc - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ p || throw(ArgumentError("`out_dim` must not be larger than `p`"*
                                        "where `p` is the number of features in `X`"))

    core_res = MS.fit(MS.MulticlassLDA, nc, Xm_t, Int.(yplain);
                      method=model.method,
                      outdim=out_dim,
                      regcoef=model.regcoef,
                      covestimator_within=model.cov_w,
                      covestimator_between=model.cov_b)

    cache     = nothing
    report    = (classes       = classes_seen,
                 out_dim       = MS.outdim(core_res),
                 class_means   = MS.classmeans(core_res),
                 mean          = MS.mean(core_res),
                 class_weights = MS.classweights(core_res),
                 Sw            = MS.withclass_scatter(core_res),
                 Sb            = MS.betweenclass_scatter(core_res),
                 nc            = nc)
    fitresult = (core_res, classes_seen)

    return fitresult, cache, report
end

function MMI.fitted_params(::LDA, (core_res, classes_seen))
    return (projected_class_means = MS.classmeans(core_res),
            projection_matrix     = MS.projection(core_res))
end

function MMI.predict(m::LDA, (core_res, classes_seen), Xnew)
    # projection of Xnew, XWt is nt x o  where o = number of out dims
    # nt = number ot test samples
    XWt = MMI.matrix(Xnew) * core_res.proj
    # centroids in the transformed space, nc x o
    centroids = permutedims(core_res.pmeans)

    # compute the distances in the transformed space between pairs of rows
    # the probability matrix is `n x nc` and normalised accross rows
    P = pairwise(m.dist, XWt, centroids, dims=1)
    # apply a softmax transformation
    P .-= maximum(P, dims=2)
    P  .= exp.(-P)
    P ./= sum(P, dims=2)

    return MMI.UnivariateFinite(classes_seen, P)
end

####
#### BayesianLDA
####

"""
    BayesianLDA(; kwargs...)

$BayesianLDA_DESCR

## Parameters

* `method=:gevd`:    choice of solver, one of `:gevd` or `:whiten`
                    methods
* `cov_w=SimpleCovariance()`: an estimator for the within-class covariance (used in computing
                              within-class scatter matrix, Sw), by
                              default set to the standard
                              `MultivariateStats.CovarianceEstimator` but could
                              be set to any robust estimator from `CovarianceEstimation.jl`.
* `cov_b=SimpleCovariance()`: same as `cov_w` but for the between-class
                              covariance(used in computing between-class scatter matrix, Sb).
* `out_dim`:                  the output dimension, i.e dimension of the
                    transformed space, automatically set if 0 is
                    given (default).
* `regcoef`:        regularization coefficient (default value 1e-6).
                    A positive value `regcoef * eigmax(Sw)` where `Sw` is the
                    within-class covariance estimator, is added to the diagonal
                    of Sw to improve numerical stability. This can be useful if
                    using the standard covariance estimator.
* `priors=nothing`: For use in prediction with Baye's rule. If `priors = nothing` then `priors` are estimated
                    from the class proportions in the training data. Otherwise it requires a `Vector` containing
                    class probabilities with probabilities specified using the order given by `levels(y)` where y is the
                    target vector.
See also the [package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Li, Zhu and Ogihara, [Using Discriminant Analysis for Multi-class Classification: An Experimental Investigation](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7068&rep=rep1&type=pdf).
"""
@mlj_model mutable struct BayesianLDA <: MMI.Probabilistic
    method::Symbol   = :gevd::(_ in (:gevd, :whiten))
    cov_w::CovarianceEstimator = MS.SimpleCovariance()
    cov_b::CovarianceEstimator = MS.SimpleCovariance()
    out_dim::Int     = 0::(_ ≥ 0)
    regcoef::Float64    = 1e-6::(_ ≥ 0)
    priors::Union{Nothing, Vector{Float64}} = nothing
end

function MMI.fit(model::BayesianLDA, ::Int, X, y)
    class_list = MMI.classes(y[1]) #class list containing entries in pool of y
    nclasses = length(class_list)
    classes_seen  = filter(in(unique(y)), class_list) #class list containing unique entries in y
    nc   = length(classes_seen) #number of classes in pool of y
    nc == nclasses || (integers_seen = MMI.int(classes_seen))

    ## If piors are specified check if they makes sense.
    ## put this here to through errors earlier
    if isa(model.priors, Vector)
        length(model.priors) == nclasses || throw(ArgumentError("Invalid size of `priors`."))
        isapprox(sum(model.priors), 1) || throw(ArgumentError("probabilities specified in `priors` must sum to 1"))
        all(model.priors .>= 0) || throw(ArgumentError("probabilities specified in `priors` must non-negative"))
        #Select priors for unique classes in `y` (For resampling purporses)
        priors = nc == nclasses ? model.priors : @view model.priors[integers_seen]
    end

    # NOTE: copy/transpose
    Xm_t   = MMI.matrix(X, transpose=true) # now p x n matrix
    yplain = MMI.int(y) # vector of n ints in {1,..., nclasses}
    p, n   = size(Xm_t)
    #recode yplain to be in {1,..., nc}
    nc == nclasses ||  replace!(yplain, (integers_seen .=> 1:nc)...)

    #check to make sure we have more than one class in training sample
    #This is to prevent Sb from being a zero matrix
    nc >= 2 || throw(ArgumentError("The number of unique classes in "*
                                            "traning sample has to be greater than one"))
    #check to make sure we have more samples than classes
    #This is to prevent Sw from being the zero matrix
    n > nc || throw(ArgumentError("The number of training samples `n` has"*
                                    " to be greater than number of unique classes `nc`"))

    # check output dimension default is min(p, nc-1)
    def_outdim = min(p, nc - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ p || throw(ArgumentError("`out_dim` must not be larger than `p`"*
                                            "where `p` is the number of features in `X`"))

    core_res = MS.fit(MS.MulticlassLDA, nc, Xm_t, Int.(yplain);
                      method=model.method,
                      outdim=out_dim,
                      regcoef= model.regcoef,
                      covestimator_within=model.cov_w,
                      covestimator_between=model.cov_b)

    ## Estimates prior probabilities if specified by user.
    ## Put it here to avoid recomputing as the fitting process does this already
    if isa(model.priors, Nothing)
        weights = MS.classweights(core_res)
        total = core_res.stats.tweight
        priors = weights ./ total
    end

    cache     = nothing
    report    = (classes       = classes_seen,
                 out_dim       = MS.outdim(core_res),
                 class_means   = MS.classmeans(core_res),
                 mean          = MS.mean(core_res),
                 class_weights = MS.classweights(core_res),
                 Sw            = MS.withclass_scatter(core_res),
                 Sb            = MS.betweenclass_scatter(core_res),
                 nc            = nc)

    fitresult = (core_res, classes_seen, priors, n)

    return fitresult, cache, report
end

function MMI.fitted_params(::BayesianLDA, (core_res, classes_seen, priors, n))
   return ( projected_class_means = MS.classmeans(core_res),
            projection_matrix     = MS.projection(core_res),
            priors                = priors)
end

function MMI.predict(m::BayesianLDA, (core_res, classes_seen, priors, n), Xnew)
     # projection of Xnew, XWt is nt x o  where o = number of out dims
    # nt = number ot test samples
    XWt = MMI.matrix(Xnew) * core_res.proj
    # centroids in the transformed space, nc x o
    centroids = permutedims(core_res.pmeans)

    # compute the distances in the transformed space between pairs of rows
    # The discriminant matrix `Pr` is of dimension `nt x nc`
    # Pr[i,k] = -0.5*(xᵢ −  µₖ)ᵀ(Σw⁻¹)(xᵢ −  µₖ) + log(priorsₖ) where (Σw = Sw/n)
    # In the transformed space this becomes
    # Pr[i,k] = -0.5*(Pᵀxᵢ −  Pᵀµₖ)ᵀ(PᵀΣw⁻¹P)(Pᵀxᵢ −  Pᵀµₖ) + log(priorsₖ)
    # PᵀSw⁻¹P = I and PᵀΣw⁻¹P = n*I due to the nature of the projection_matrix, P
    # Giving Pr[i,k] = -0.5*n*(Pᵀxᵢ −  Pᵀµₖ)ᵀ(Pᵀxᵢ −  Pᵀµₖ) + log(priorsₖ)
    # (Pᵀxᵢ −  Pᵀµₖ)ᵀ(Pᵀxᵢ −  Pᵀµₖ) is the SquaredEquclidean distance in the transformed space
    Pr = pairwise(SqEuclidean(), XWt, centroids, dims=1)
    Pr .*= (-0.5*n)
    Pr .+= log.(priors)'

    # apply a softmax transformation to convert Pr to a probability matrix
    Pr .-= maximum(Pr, dims=2)
    Pr  .= exp.(Pr)
    Pr ./= sum(Pr, dims=2)

    return MMI.UnivariateFinite(classes_seen, Pr)
end

####
#### BayesianSubspaceLDA
####

"""
    BayesianSubspaceLDA(; kwargs...)

$BayesianSubspaceLDA_DESCR

## Parameters

* `normalize=true`: Option to normalize the between class variance for the number of
                    observations in each class, one of `true` or `false`.
* `out_dim`:        the dimension of the transformed space
                    to be used by `predict` and `transform` methods, automatically set if 0 is given (default).
* `priors=nothing`: For use in prediction with Baye's rule. If `priors = nothing` then `priors` are estimated
                    from the class proportions in the training data. Otherwise it requires a `Vector` containing
                    class probabilities with probabilities specified using the order given by `levels(y)` where y is the
                    target vector.

For more information about the algorithm, see the paper by Howland & Park (2006),
"Generalizing discriminant analysis using the generalized singular value decomposition",IEEE Trans. Patt.
 Anal. & Mach. Int., 26: 995-1006.
"""
@mlj_model mutable struct BayesianSubspaceLDA <: MMI.Probabilistic
    normalize::Bool=false
    out_dim::Int   = 0::(_ ≥ 0)
    priors::Union{Nothing, Vector{Float64}} = nothing
end

function MMI.fit(model::BayesianSubspaceLDA, ::Int, X, y)
    class_list = MMI.classes(y[1]) #class list containing entries in pool of y
    classes_seen  = filter(in(unique(y)), class_list) #class list containing unique entries in y
    nc   = length(classes_seen) #number of classes in pool of y
    nclasses = length(class_list)
    nc == nclasses || (integers_seen = MMI.int(classes_seen))

    ## If piors are specified check if they makes sense.
    ## put this here to through errors earlier
    if isa(model.priors, Vector)
        length(model.priors) == nclasses || throw(ArgumentError("Invalid size of `priors`."))
        isapprox(sum(model.priors), 1) || throw(ArgumentError("probabilities specified in `priors` must sum to 1"))
        all(model.priors .>= 0) || throw(ArgumentError("probabilities specified in `priors` must non-negative"))
        #Select priors for unique classes in `y`(For resampling purporses)
        priors = nc == nclasses ? model.priors : @view model.priors[integers_seen]
    end
    # NOTE: copy/transpose
    Xm_t   = MMI.matrix(X, transpose=true) # now p x n matrix
    yplain = MMI.int(y) # vector of n ints in {1,..., nclasses}
    p, n   = size(Xm_t)
    #recode yplain to be in {1,..., nc}
    nc == nclasses ||  replace!(yplain, (integers_seen .=> 1:nc)...)

    #check to make sure we have more than one class in training sample
    #This is to prevent Sb from being a zero matrix
    nc >= 2 || throw(ArgumentError("The number of unique classes in traning sample"*
                                    " has to be greater than one"))
    #check to make sure we have more samples than classes
    #This is to prevent Sw from being the zero matrix
    n > nc || throw(ArgumentError("The number of training samples has to be greater"*
                                    " than the number of unique classes in `y`"))

    # check output dimension default is  min(p, nc-1)
    def_outdim = min(p, nc - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ nc - 1 || throw(ArgumentError("`out_dim` must not be larger than `nc - 1`"*
                                            " where  `nc` is the number of unique classes in `y`"))

    core_res = MS.fit(MS.SubspaceLDA, Xm_t,
                        Int.(yplain),
                        nc;
                        normalize = model.normalize)

    λ = core_res.λ # λ is a (nc -1) x 1 vector containing the eigen values sorted in descending order.
    explained_variance_ratio = λ ./ sum(λ) #proportions of variance

    mult = model.normalize ? n : 1 #used in prediction

    ## Estimates prior probabilities if specified by user.
    ## Put it here to avoid recomputing as the fitting process does this already
    if isa(model.priors, Nothing)
        weights = MS.classweights(core_res)
        priors = weights ./ n
    end

    cache     = nothing
    report    = (explained_variance_ratio  = explained_variance_ratio,
                 classes       = classes_seen,
                 class_means   = MS.classmeans(core_res),
                 mean          = MS.mean(core_res),
                 class_weights = MS.classweights(core_res),
                 nc            = nc)
    fitresult = (core_res, out_dim, classes_seen, priors, n, mult)

    return fitresult, cache, report
end

function MMI.fitted_params(::BayesianSubspaceLDA, (core_res, _, _, priors,_))
    return (projected_class_means  = MS.classmeans(core_res),
            projection_matrix      = MS.projection(core_res),
            priors                 = priors)
end

function MMI.predict(m::BayesianSubspaceLDA, (core_res, out_dim, classes_seen, priors, n, mult), Xnew)
    # projection of Xnew, XWt is nt x o  where o = number of out dims
    # nt = number ot test samples
    proj = core_res.projw * view(core_res.projLDA, :, 1:out_dim) #proj is the projection_matrix
    XWt = MMI.matrix(Xnew) * proj

    # centroids in the transformed space, nc x o
    centroids = permutedims(proj' * core_res.cmeans)
    nc = length(classes_seen)

    # compute the distances in the transformed space between pairs of rows
    # The discriminant matrix `Pr` is of dimension `nt x nc`
    # Pr[i,k] = -0.5*(xᵢ −  µₖ)ᵀ(Σw⁻¹)(xᵢ −  µₖ) + log(priorsₖ) where (Σw = Sw/(n-nc)
    # In the transformed space this becomes
    # Pr[i,k] = -0.5*(Pᵀxᵢ −  Pᵀµₖ)ᵀ(PᵀΣw⁻¹P)(Pᵀxᵢ −  Pᵀµₖ) + log(priorsₖ)
    # PᵀSw⁻¹P = (1/mult)*I and PᵀΣw⁻¹P = (n-nc)/mult*I due to the nature of the projection_matrix, P
    # Giving Pr[i,k] = -0.5*n*(Pᵀxᵢ −  Pᵀµₖ)ᵀ(Pᵀxᵢ −  Pᵀµₖ) + log(priorsₖ)
    # (Pᵀxᵢ −  Pᵀµₖ)ᵀ(Pᵀxᵢ −  Pᵀµₖ) is the SquaredEquclidean distance in the transformed space
    Pr = pairwise(SqEuclidean(), XWt, centroids, dims=1)
    Pr .*= (-0.5 * (n-nc)/mult)
    Pr .+= log.(priors)'

    # apply a softmax transformation to convert Pr to a probability matrix
    Pr .-= maximum(Pr, dims=2)
    Pr  .= exp.(Pr)
    Pr ./= sum(Pr, dims=2)

    return MMI.UnivariateFinite(classes_seen, Pr)
end

####
#### SubspaceLDA
####

"""
    SubspaceLDA(; kwargs...)

$SubspaceLDA_DESCR

## Parameters

* `normalize=true`:   Option to normalize the between class variance for the number
                      of observations in each class, one of `true` or `false`.
* `out_dim`:        the dimension of the transformed space
                    to be used by `predict` and `transform` methods, automatically set if 0 is given (default).
* `dist=SqEuclidean`: the distance metric to use when performing classification
                      (to compare the distance between a new point and
                      centroids in the transformed space), an alternative
                      choice can be the `CosineDist`.

See also the [package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Howland & Park (2006), "Generalizing discriminant analysis using the generalized singular value decomposition", IEEE Trans. Patt. Anal. & Mach. Int., 26: 995-1006.
"""
@mlj_model mutable struct SubspaceLDA <: MMI.Probabilistic
    normalize::Bool=true
    out_dim::Int   = 0::(_ ≥ 0)
    dist::SemiMetric = SqEuclidean()
end

function MMI.fit(model::SubspaceLDA, ::Int, X, y)
    class_list = MMI.classes(y[1]) #class list containing entries in pool of y
    classes_seen  = filter(in(unique(y)), class_list) #class list containing unique entries in y
    nc   = length(classes_seen) #number of classes in pool of y
    nclasses = length(class_list)
    nc == nclasses || (integers_seen = MMI.int(classes_seen))

    # NOTE: copy/transpose
    Xm_t   = MMI.matrix(X, transpose=true) # now p x n matrix
    yplain = MMI.int(y) # vector of n ints in {1,..., nclasses}
    p, n   = size(Xm_t)
    #recode yplain to be in {1,..., nc}
    nc == nclasses ||  replace!(yplain, (integers_seen .=> 1:nc)...)

    #check to make sure we have more than one class in training sample
    #This is to prevent Sb from being a zero matrix
    nc >= 2 || throw(ArgumentError("The number of unique classes in traning sample"*
                                    " has to be greater than one"))
    #check to make sure we have more samples than classes
    #This is to prevent Sw from being the zero matrix
    n > nc || throw(ArgumentError("The number of training samples has to be greater"*
                                    " than the number of unique classes in `y`"))

    # check output dimension default is  min(p, nc-1)
    def_outdim = min(p, nc - 1)
    # if unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # check if the given one is sensible
    out_dim ≤ nc - 1 || throw(ArgumentError("`out_dim` must not be larger than `nc - 1`"*
                                            " where  `nc` is the number of unique classes in `y`"))

    core_res = MS.fit(MS.SubspaceLDA, Xm_t,
                        Int.(yplain),
                        nc;
                        normalize = model.normalize)

    λ = core_res.λ # λ is a (nc -1) x 1 vector containing the eigen values sorted in descending order.
    explained_variance_ratio = λ ./ sum(λ) #proportions of variance

    cache     = nothing
    report    = (explained_variance_ratio  = explained_variance_ratio,
                 classes       = classes_seen,
                 class_means   = MS.classmeans(core_res),
                 mean          = MS.mean(core_res),
                 class_weights = MS.classweights(core_res),
                 nc            = nc)
    fitresult = (core_res, out_dim, classes_seen)

    return fitresult, cache, report
end

function MMI.fitted_params(::SubspaceLDA, (core_res, _))
    return (class_means       = MS.classmeans(core_res),
            projection_matrix = MS.projection(core_res))
end

function MMI.predict(m::SubspaceLDA, (core_res, out_dim, classes_seen), Xnew)
     # projection of Xnew, XWt is nt x o  where o = number of out dims
    # nt = number ot test samples
    proj = core_res.projw * view(core_res.projLDA, :, 1:out_dim) #proj is the projection_matrix
    XWt = MMI.matrix(Xnew) * proj
    # centroids in the transformed space, nc x o
    centroids = permutedims(proj' * core_res.cmeans)

    # compute the distances in the transformed space between pairs of rows
    # the probability matrix is `nt x nc` and normalised accross rows
    P = pairwise(m.dist, XWt, centroids, dims=1)

    # apply a softmax transformation
    P .-= maximum(P, dims=2)
    P  .= exp.(-P)
    P ./= sum(P, dims=2)

    return MMI.UnivariateFinite(classes_seen, P)
end

function MMI.transform(m::T, (core_res, out_dim, _), X) where T<:Union{SubspaceLDA, BayesianSubspaceLDA}
    # projection of X, XWt is nt x o  where o = out dims
    proj = core_res.projw * view(core_res.projLDA, :, 1:out_dim) #proj is the projection_matrix
    XWt = MMI.matrix(X) * proj
    return MMI.table(XWt, prototype = X)
end

function MMI.transform(m::T, (core_res,), X) where T<:Union{LDA, BayesianLDA}
    # projection of X, XWt is nt x o  where o = out dims
    proj = core_res.projw * core_res.projLDA #proj is the projection_matrix
    XWt = MMI.matrix(X) * proj
    return MMI.table(XWt, prototype = X)
end

####
#### METADATA
####

metadata_pkg.(
    (RidgeRegressor, PCA, KernelPCA, ICA, LDA,
     BayesianLDA, SubspaceLDA,  BayesianSubspaceLDA),
    name       = "MultivariateStats",
    uuid       = "6f286f6a-111f-5878-ab1e-185364afe411",
    url        = "https://github.com/JuliaStats/MultivariateStats.jl",
    license    = "MIT",
    julia      = true,
    is_wrapper = false)

metadata_model(RidgeRegressor,
    input   = Table(Continuous),
    target  = AbstractVector{Continuous},
    weights = false,
    descr   = RIDGE_DESCR)

metadata_model(PCA,
    input   = Table(Continuous),
    target  = Table(Continuous),
    weights = false,
    descr   = PCA_DESCR)

metadata_model(KernelPCA,
    input   = Table(Continuous),
    target  = Table(Continuous),
    weights = false,
    descr   = KPCA_DESCR)

metadata_model(ICA,
    input   = Table(Continuous),
    target  = Table(Continuous),
    weights = false,
    descr   = ICA_DESCR)

metadata_model(LDA,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    output  = Table(Continuous),
    descr   = LDA_DESCR)

metadata_model(BayesianLDA,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    output  = Table(Continuous),
    descr   = BayesianLDA_DESCR)

metadata_model(SubspaceLDA,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    output  = Table(Continuous),
    descr   = SubspaceLDA_DESCR)

metadata_model(BayesianSubspaceLDA,
    input   = Table(Continuous),
    target  = AbstractVector{<:Finite},
    weights = false,
    output  = Table(Continuous),
    descr   = BayesianSubspaceLDA_DESCR)

end


text:
using Documenter, JuliaFormatter

makedocs(
    sitename = "JuliaFormatter",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [JuliaFormatter],
    pages = [
        "Introduction" => "index.md",
        "How It Works" => "how_it_works.md",
        "Code Style" => "style.md",
        "Skipping Formatting" => "skipping_formatting.md",
        "Syntax Transforms" => "transforms.md",
        "Custom Alignment" => "custom_alignment.md",
        "Custom Styles" => "custom_styles.md",
        "YAS Style" => "yas_style.md",
        "Blue Style" => "blue_style.md",
        "Configuration File" => "config.md",
        "API Reference" => "api.md",
    ],
)

deploydocs(repo = "github.com/domluna/JuliaFormatter.jl.git", push_preview = true)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Distributions, PDMats

@testset "density_product" begin
    mvec = [-0.3, 0.3, 0.1]
    cmat = [1.0 1.5 0.0; 1.5 4.0 0.0; 0.0 0.0 1.0]
    Σ = @inferred PDMat(cmat)
    mvnorm = @inferred  MvNormal(mvec, Σ)
    mvn_density = @inferred BAT.GenericDensity(params -> LogDVal(logpdf(mvnorm, params)))

    cmat = [3.76748 0.446731 0.625418; 0.446731 3.9317 0.237361; 0.625418 0.237361 3.43867]
    mvec = [1., 2, 1.]
    Σ = @inferred PDMat(cmat)
    mvt = MvTDist(3, mvec, Σ)
    mvt_density = @inferred BAT.GenericDensity(params -> LogDVal(logpdf(mvt, params)))

    params = VectorOfSimilarVectors([0.0 -0.3; 0.0 0.3; 0.0 1.0])

    pb = BAT.HyperRectBounds([-2.0, -1.0, -0.5], [2.0, 3.0, 1.0],
                                BAT.reflective_bounds)
    pb2 = BAT.HyperRectBounds([-1.5, -2.0, -0.5], [2.0, 2.5, 1.5],
                                BAT.reflective_bounds)
    dp = @inferred BAT.DensityProduct((mvt_density, mvn_density), pb)
    dp1 = @inferred BAT.DensityProduct((mvt_density,), pb)
    dp2 = @inferred BAT.DensityProduct((mvn_density,), pb2)

    res = Array{Float64}(undef, size(params, 1))

    @testset "DensityProduct" begin
        @test typeof(dp) <: BAT.DensityProduct{2,
            Tuple{BAT.GenericDensity{typeof(mvt_density.f)},
                  BAT.GenericDensity{typeof(mvn_density.f)}},BAT.HyperRectBounds{Float64}}

        @test parent(dp)[1] == mvt_density
        @test BAT.var_bounds(dp) == pb
        @test totalndof(dp) == 3
    end

    @testset "unsafe_pod" begin
        up = @inferred BAT._unsafe_prod(mvt_density, mvn_density, pb)
        @test parent(up)[2] == mvn_density
        @test BAT.var_bounds(up) == pb

        prd = @inferred dp1*dp2
        @test typeof(prd) <: BAT.DensityProduct{2,
            Tuple{BAT.GenericDensity{typeof(mvt_density.f)},
                  BAT.GenericDensity{typeof(mvn_density.f)}},BAT.HyperRectBounds{Float64}}
        @test parent(prd)[2] == mvn_density
        prd_pb = @inferred BAT.var_bounds(prd)
        cut_pb = @inferred pb ∩ pb2
        @test prd_pb.vol.lo ≈ cut_pb.vol.lo
        @test prd_pb.vol.hi ≈ cut_pb.vol.hi

        up = @inferred BAT._unsafe_prod(dp1, mvn_density, pb)
        @test parent(up)[1] == mvt_density
        @test parent(up)[2] == mvn_density
        @test BAT.var_bounds(up) == pb

        up = @inferred BAT._unsafe_prod(mvt_density, dp2, pb)
        @test parent(up)[1] == mvt_density
        @test parent(up)[2] == mvn_density
        @test BAT.var_bounds(up) == pb

        up = @inferred BAT._unsafe_prod(dp1, dp2, pb)
        @test parent(up)[1] == mvt_density
        @test parent(up)[2] == mvn_density
        @test BAT.var_bounds(up) == pb
    end

    @testset "BAT.eval_logval_unchecked" begin
        @test BAT.eval_logval_unchecked(dp, params[1]) ≈ -8.8547305
    end
end


text:
module InferenceTests

using Test
using InplaceArrays.Inference
import InplaceArrays.Inference: testargs

bar(x,y) = x*y
foo(x) = sqrt(x-1)
hola(x,y) = sqrt(x-1)*y

testargs(::typeof(foo),T::DataType) = (one(T),)
testargs(::typeof(hola),Tx,Ty) = (one(Tx),testvalue(Ty))

@test 0 == testvalue(Int)
@test (0,0.0) == testvalues(Int,Float64)

@test return_type(sqrt, Int) == Float64
@test return_type(foo, Int) == Float64
@test return_type(bar, Int, Int) == Int
@test return_type(bar, Float64, Int) == Float64
@test return_type(bar, Matrix{Float64}, Int) == Matrix{Float64}
@test return_type(hola, Int, Vector{Int}) == Vector{Float64}

@test return_type_broadcast(bar, Matrix{Float64}, Int) == Matrix{Float64}
@test return_type_broadcast(foo, Matrix{Float64}) == Matrix{Float64}
@test return_type_broadcast(foo, Matrix{Float64}) == Matrix{Float64}
@test return_type_broadcast(bar, Float64, Int) == Float64
@test return_type_broadcast(foo, Float64) == Float64
@test return_type_broadcast(bar, Int, Int) == Int
@test return_type_broadcast(bar, Int, Vector{Int}) == Vector{Int}
@test return_type_broadcast(hola, Int, Vector{Int}) == Vector{Float64}
@test return_type_broadcast(hola, Int, Int) == Float64

end # module


text:
"Example of a more complicated periodic contractor"
function sawtooth(X::IntervalBox)

    x, y = X

    x = x ∩ (-1..1)
    y = y ∩ (-2..2)

    y = y ∩ (2x)
    x = x ∩ (y/2)

    return IntervalBox(x, y)
end

function constant_contractor(X, y_val)
    x, y = X
    y = y ∩ Interval(y_val)
    return IntervalBox(x, y)
end

function add_one(X)  # y = x + 1
    x, y = X

    y = y ∩ (x + 1)
    x = x ∩ (y - 1)

    return IntervalBox(x, y)
end

function cube0(X::IntervalBox)  # contractor for y=x^3, x>=0

    x, y = X

    x = x ∩ (0..∞)

    y = y ∩ (x ^ 3)
    x = x ∩ Interval(y.lo ^ (1/3), y.hi^(1/3))   # not rigorous!

    return x × y
end

"""
Make the "odd" version of a box
odd(X::IntervalBox) = ( (x,y) = X; IntervalBox(-x, -y) )
"""

cube_neg = symmetrise(cube0, odd)
cube = cube0 ∪ cube_neg

ff(x) = x^2 - x^3

function ff(X::IntervalBox)
    x, y = X

    a = x^2
    b = x^3

    y = y ∩ (a - b)  # y = a - b

    a = a ∩ (y + b)
    b = b ∩ (a - y)

    x, a = square(IntervalBox(x, a))
    x, b = cube(IntervalBox(x, b))

    return IntervalBox(x, y)

end


text:
using Test
using Iris

@testset "Common" begin
    include("common.jl")
end

@testset "Spectral" begin
    include("spectral.jl")
end

@testset "Floquet" begin
    include("floquet.jl")
end

@testset "Scattering" begin
    include("scattering.jl")
end

@testset "SMatrices" begin
    include("smatrices.jl")
end

@testset "Lasing" begin
    include("lasing.jl")
end

@testset "SaturableCPA" begin
    include("saturablecpa.jl")
end

@testset "TimeDomain" begin
    include("timedomain.jl")
end


text:
using ArcadeLearningEnvironment
using ArgParse
using Cambrian
using CartesianGeneticProgramming
using Dates
using IICGP
using Random

# function extension
import Cambrian.mutate
import Cambrian.evaluate

s = ArgParseSettings()
@add_arg_table! s begin
    "--cfg"
    help = "configuration script"
    default = "cfg/test_mono.yaml"
    "--game"
    help = "game rom name"
    default = "assault"
    "--seed"
    help = "random seed for evolution"
    arg_type = Int
    default = 0
    "--ind"
    help = "individual for evaluation"
    arg_type = String
    default = ""
end
args = parse_args(ARGS, s)
const game = args["game"]
const seed = args["seed"]
Random.seed!(seed)

main_cfg, cont_cfg, reducer, bootstrap = IICGP.monocgp_config(args["cfg"], game)

const max_frames = main_cfg["max_frames"]
const stickiness = main_cfg["stickiness"]
const grayscale = main_cfg["grayscale"]
const downscale = main_cfg["downscale"]
const logid = cont_cfg.id
const state_ref = get_state_ref(game, seed)

function play_atari(
    reducer::Reducer,
    controller::CGPInd,
    lck::ReentrantLock;
    rom=game,
    seed=seed,
    rom_state_ref=state_ref,
    max_frames=max_frames,
    grayscale=grayscale,
    downscale=downscale,
    stickiness=stickiness
)
    # Random.seed!(seed)
    mt = MersenneTwister(seed)
    game = Game(rom, seed, lck=lck, state_ref=rom_state_ref)
    reward = 0.0
    frames = 0
    prev_action = Int32(0)
    while ~game_over(game.ale)
        if rand(mt) > stickiness || frames == 0
            s = get_state(game, grayscale, downscale)
            output = IICGP.process(reducer, controller, s)
            action = game.actions[argmax(output)]
        else
            action = prev_action
        end
        reward += act(game.ale, action)
        frames += 1
        prev_action = action
        if frames > max_frames
            break
        end
    end
    close!(game)
    [reward]
end

if length(args["ind"]) > 0
    ind = CGPInd(cfg, read(args["ind"], String))
    ftn = fitness(ind, inps, outs)
    println("Fitness: ", ftn)
else
    mutate(ind::CGPInd) = goldman_mutate(cont_cfg, ind)
    lck = ReentrantLock()
    fit(controller::CGPInd) = play_atari(reducer, controller, lck)
    evaluate(e::CGPEvolution) = IICGP.fitness_evaluate(e, e.fitness)
    e = CartesianGeneticProgramming.CGPEvolution(cont_cfg, fit)
    init_backup(logid, args["cfg"])
    run!(e)
end


text:
module TestDerivatives

end # module