text:
using AdventOfCodeSolutions
using Test
using IterTools
using OffsetArrays

function input(puzzle::Puzzle{2020, 10, n}) where n
    io = openInput(puzzle)
    split(read(io, String), "\n", keepempty=false)
end

function parseInput(input)
    parse.(Int, input)
end

adjacentElements(c) = Iterators.zip(c, Iterators.drop(c, 1))
adjacentDifferences(c) = map(((left, right),) -> right - left, adjacentElements(c))
function solve(::Puzzle{2020, 10, 1}, input)
    adapters = parseInput(input)
    adapters = sort(vcat(0, adapters, maximum(adapters) + 3))
    differences = map(length , groupby(identity, adjacentDifferences(adapters) |> sort))

    return prod(differences)
end

function solve(::Puzzle{2020, 10, 2}, input)
    adapters = vcat(0, sort(parseInput(input)))

    n = length(adapters)
    combinations = zeros(Int64, n)
    combinations[1] = 1

    for i = 1:n
        if i + 1 > n
            continue
        end
        if adapters[i + 1] - adapters[i] <= 3
            combinations[i + 1] += combinations[i]
        end
        if i + 2 > n
            continue
        end
        if adapters[i + 2] - adapters[i] <= 3
            combinations[i + 2] += combinations[i]
        end
        if i + 3 > n
            continue
        end
        if adapters[i + 3] - adapters[i] <= 3
            combinations[i + 3] += combinations[i]
        end
    end

    return combinations[end]
end

function solve(::Puzzle{2020, 10, 3}, input)
    adapters = vcat(0, sort(parseInput(input)))
    paths = OffsetVector(zeros(adapters[end] + 4), -1)
    paths[0] = 1

    n = length(adapters)
    for adapter = adapters
        p = paths[adapter]
        paths[adapter + 1] += p
        paths[adapter + 2] += p
        paths[adapter + 3] += p
    end

    return paths[adapters[end]]
end

@testset "Day 10" begin
    TEST_INPUT_1 = """
16
10
15
5
1
11
7
19
6
12
4
"""

    TEST_INPUT_1_2 = """
28
33
18
42
31
14
46
20
48
47
24
23
49
45
19
38
39
11
1
32
25
35
8
17
7
9
4
2
34
10
3
"""

    @test(solve(Puzzle(10, 1), split(TEST_INPUT_1, '\n', keepempty=false)) == 7 * 5)
    @test(solve(Puzzle(10, 1), split(TEST_INPUT_1_2, '\n', keepempty=false)) == 22 * 10)
    @test(solve(Puzzle(10, 1), input(Puzzle(10, 1))) == 2432)

    TEST_INPUT_2 = """
"""

    for p = 2:3
        @test(solve(Puzzle(10, p), ["1"]) == 1)
        @test(solve(Puzzle(10, p), ["1", "2"]) == 2)
        @test(solve(Puzzle(10, p), ["1", "2", "3"]) == 4)
        @test(solve(Puzzle(10, p), ["1", "3"]) == 2)
        @test(solve(Puzzle(10, p), ["1", "4"]) == 1)
        @test(solve(Puzzle(10, p), ["1", "2", "3", "4"]) == 7)
        @test(solve(Puzzle(10, p), split(TEST_INPUT_1, '\n', keepempty=false)) == 8)
        @test(solve(Puzzle(10, p), split(TEST_INPUT_1_2, '\n', keepempty=false)) == 19208)
        @test(solve(Puzzle(10, p), input(Puzzle(10, 2))) == 453551299002368)
    end
end


text:
using BuildExecutable
mkdir("bin")
build_executable("eve-julia", "src/main.jl", "bin", "native")

ext: jl

max_stars_repo_name: barcharcraz/eve-julia-gui

max_issues_repo_name: barcharcraz/eve-julia-gui

max_forks_repo_name: barcharcraz/eve-julia-gui


text:
export Derivative,Integral,Laplacian,Volterra

abstract type CalculusOperator{S,OT,T}<:Operator{T} end

## Note that all functions called in calculus_operator must be exported

macro calculus_operator(Op)
    ConcOp=Meta.parse("Concrete"*string(Op))
    WrappOp=Meta.parse(string(Op)*"Wrapper")
    DefaultOp=Meta.parse("Default"*string(Op))
    return esc(quote
        # The SSS, TTT are to work around #9312
        abstract type $Op{SSS,OT,TTT} <: CalculusOperator{SSS,OT,TTT} end

        struct $ConcOp{S<:Space,OT,T} <: $Op{S,OT,T}
            space::S        # the domain space
            order::OT
        end
        struct $WrappOp{BT<:Operator,S<:Space,OT,T} <: $Op{S,OT,T}
            op::BT
            order::OT
        end

        @wrapper $WrappOp

        ## Constructors
        $ConcOp(sp::Space,k) = $ConcOp{typeof(sp),typeof(k),prectype(sp)}(sp,k)

        $Op(sp::UnsetSpace,k) = $ConcOp(sp,k)
        $Op(sp::UnsetSpace,k::Real) = $ConcOp(sp,k)
        $Op(sp::UnsetSpace,k::Integer) = $ConcOp(sp,k)

        function $DefaultOp(sp::Space,k)
            csp=canonicalspace(sp)
            if conversion_type(csp,sp)==csp   # Conversion(sp,csp) is not banded, or sp==csp
               error("Implement $(string($Op))($(string(sp)),$k)")
            end
            $WrappOp(TimesOperator([$Op(csp,k),Conversion(sp,csp)]),k)
        end

        $DefaultOp(d,k) = $Op(Space(d),k)

        $DefaultOp(sp) = $Op(sp,1)
        $DefaultOp() = $Op(UnsetSpace())
        $DefaultOp(k::Number) = $Op(UnsetSpace(),k)
        $DefaultOp(k::Vector) = $Op(UnsetSpace(),k)

        $Op(x...) = $DefaultOp(x...)
        $ConcOp(S::Space) = $ConcOp(S,1)

        function Base.convert(::Type{Operator{T}},D::$ConcOp) where T
            if T==eltype(D)
                D
            else
                $ConcOp{typeof(D.space),typeof(D.order),T}(D.space,D.order)
            end
        end

        $WrappOp(op::Operator,order) =
            $WrappOp{typeof(op),typeof(domainspace(op)),typeof(order),eltype(op)}(op,order)
        $WrappOp(op::Operator) = $WrappOp(op,1)

        function Base.convert(::Type{Operator{T}},D::$WrappOp) where T
            if T==eltype(D)
                D
            else
                # work around typeinfernece bug
                op=convert(Operator{T},D.op)
                $WrappOp{typeof(op),typeof(domainspace(op)),typeof(D.order),T}(op,D.order)
            end
        end

        ## Routines
        ApproxFun.domain(D::$ConcOp) = domain(D.space)
        ApproxFun.domainspace(D::$ConcOp) = D.space

        Base.getindex(::$ConcOp{UnsetSpace,OT,T},k::Integer,j::Integer) where {OT,T} =
            error("Spaces cannot be inferred for operator")
        ApproxFun.rangespace(D::$ConcOp{UnsetSpace,T}) where {T} = UnsetSpace()

        #promoting domain space is allowed to change range space
        # for integration, we fall back on existing conversion for now
        ApproxFun.promotedomainspace(D::$Op,sp::UnsetSpace) = D

        function ApproxFun.promotedomainspace(D::$Op,sp::Space)
            if isambiguous(domain(sp))
                $Op(typeof(sp)(domain(D)),D.order)
            else
                $Op(sp,D.order)
            end
        end
    end)
#     for func in (:rangespace,:domainspace,:bandwidths)
#         # We assume the operator wrapped has the correct spaces
#         @eval $func(D::$WrappOp)=$func(D.op)
#     end
end

choosedomainspace(M::CalculusOperator{UnsetSpace},sp::Space) =
    iswrapper(M) ? choosedomainspace(M.op,sp) : sp  # we assume the space itself will work

@calculus_operator(Derivative)
@calculus_operator(Integral)
@calculus_operator(Volterra)

## Overrideable

## Convenience routines

integrate(d::Domain) = Integral(d,1)

# Default is to use ops
differentiate(f::Fun)=Derivative(space(f))*f
function integrate(f::Fun)
    d=domain(f)
    cd=canonicaldomain(d)
    if typeof(d)==typeof(cd)  || isa(d,PeriodicDomain)
        Integral(space(f))*f
    else
        # map to canonical domain
        setdomain(integrate(setdomain(f,cd)*fromcanonicalD(f)),d)
    end
end

function Base.sum(f::Fun)
    d=domain(f)
    cd=canonicaldomain(d)
    if typeof(cd)==typeof(d)  || isa(d,PeriodicDomain)
        g=integrate(f)
        last(g)-first(g)
    else
        # map first
        sum(setdomain(f,cd)*fromcanonicalD(f))
    end
end

function linesum(f::Fun)
    cd=canonicaldomain(f)
    d=domain(f)

    if isreal(d)
        a,b=leftendpoint(d),rightendpoint(d)
        sign(last(b)-first(a))*sum(f)
    elseif typeof(cd)==typeof(d)  || isa(d,PeriodicDomain)
        error("override linesum for $(f.space)")
    else
        # map first
        linesum(setdomain(f,cd)*abs(fromcanonicalD(f)))
    end
end

# Multivariate

@calculus_operator(Laplacian)

## Map to canonical
function DefaultDerivative(sp::Space,k::Integer)
    if typeof(canonicaldomain(sp)).name==typeof(domain(sp)).name
        # this is the normal default constructor
        csp=canonicalspace(sp)
        if conversion_type(csp,sp)==csp   # Conversion(sp,csp) is not banded, or sp==csp
           error("Implement Derivative($(string(sp)),$k)")
        end
        DerivativeWrapper(TimesOperator([Derivative(csp,k),Conversion(sp,csp)]),k)
    else
        D1=invfromcanonicalD(sp)*Derivative(setdomain(sp,canonicaldomain(sp)))
        D=DerivativeWrapper(SpaceOperator(D1,sp,setdomain(rangespace(D1),domain(sp))),1)
        if k==1
            D
        else
            DerivativeWrapper(TimesOperator(Derivative(rangespace(D),k-1),D),k)
        end
    end
end

function DefaultIntegral(sp::Space,k::Integer)
    if typeof(canonicaldomain(sp)).name==typeof(domain(sp)).name
        # this is the normal default constructor
        csp=canonicalspace(sp)
        if conversion_type(csp,sp)==csp   # Conversion(sp,csp) is not banded, or sp==csp
            # we require that Integral is overridden
            error("Implement Integral($(string(sp)),$k)")
        end
        IntegralWrapper(TimesOperator([Integral(csp,k),Conversion(sp,csp)]),k)
    elseif k > 1
        Q=Integral(sp,1)
        IntegralWrapper(TimesOperator(Integral(rangespace(Q),k-1),Q),k)
    else # k==1
        csp=setdomain(sp,canonicaldomain(sp))

        x=Fun(identity,domain(csp))
        M=Multiplication(fromcanonicalD(sp,x),csp)
        Q=Integral(rangespace(M))*M
        IntegralWrapper(SpaceOperator(Q,sp,setdomain(rangespace(Q),domain(sp))),1)
    end
end

for TYP in (:Derivative,:Integral,:Laplacian)
    @eval begin
        function *(D1::$TYP,D2::$TYP)
            @assert domain(D1) == domain(D2)

            $TYP(domainspace(D2),D1.order+D2.order)
        end
    end
end

"""
`Derivative(sp::Space,k::Int)` represents the `k`-th derivative on `sp`.
"""
Derivative(::Space,::Int)

"""
`Derivative(sp::Space,k::Vector{Int})` represents a partial derivative on a multivariate space.
For example,
```julia
Dx = Derivative(Chebyshev()^2,[1,0]) # ∂/∂x
Dy = Derivative(Chebyshev()^2,[0,1]) # ∂/∂y
```
"""
Derivative(::Space,::Vector{Int})

"""
`Derivative(sp::Space)` represents the first derivative `Derivative(sp,1)`.
"""
Derivative(::Space)

"""
`Derivative(k)` represents the `k`-th derivative, acting on an unset space.
Spaces will be inferred when applying or manipulating the operator.
"""
Derivative(k)

"""
`Derivative()` represents the first derivative on an unset space.
Spaces will be inferred when applying or manipulating the operator.
"""
Derivative()

"""
`Integral(sp::Space,k::Int)` represents a `k`-th integral on `sp`.
There is no guarantee on normalization.
"""
Integral(::Space,::Int)

"""
`Integral(sp::Space)` represents the first integral `Integral(sp,1)`.
"""
Integral(::Space)

"""
Integral(k)` represents the `k`-th integral, acting on an unset space.
Spaces will be inferred when applying or manipulating the operator.
"""
Integral(k)

"""
`Intergral()` represents the first integral on an unset space.
Spaces will be inferred when applying or manipulating the operator.
"""
Integral()

"""
`Laplacian(sp::Space)` represents the laplacian on space `sp`.
"""
Laplacian(::Space)

"""
`Laplacian()` represents the laplacian on an unset space.
Spaces will be inferred when applying or manipulating the operator.
"""
Laplacian()


text:
"""
Fetch the "natural" block size of this band.

GDAL contains a concept of the natural block size of rasters so that
applications can organized data access efficiently for some file formats.
The natural block size is the block size that is most efficient for accessing
the format. For many formats this is simple a whole scanline in which case
`*pnXSize` is set to `GetXSize()`, and *pnYSize is set to 1.

However, for tiled images this will typically be the tile size.

Note that the X and Y block sizes don't have to divide the image size evenly,
meaning that right and bottom edge blocks may be incomplete. See `ReadBlock()`
for an example of code dealing with these issues.
"""
function getblocksize(rb::RasterBand)
    xy = Array{Cint}(2); x = pointer(xy); y = x + sizeof(Cint)
    GDAL.getblocksize(rb.ptr, x, y)
    xy
end

"Fetch the pixel data type for this band."
getdatatype(rb::RasterBand) = _JLTYPE[GDAL.getrasterdatatype(rb.ptr)]

"Fetch the width in pixels of this band."
width(rb::RasterBand) = GDAL.getrasterbandxsize(rb.ptr)

"Fetch the height in pixels of this band."
height(rb::RasterBand) = GDAL.getrasterbandysize(rb.ptr)

"Find out if we have update permission for this band."
getaccess(rb::RasterBand) = GDALAccess(GDAL.getrasteraccess(rb.ptr))

"""Fetch the band number (1+) within its dataset, or 0 if unknown.

This method may return a value of 0 to indicate overviews, or free-standing
`GDALRasterBand` objects without a relationship to a dataset.
"""
getnumber(rb::RasterBand) = GDAL.getbandnumber(rb.ptr)

"""
Fetch the handle to its dataset handle, or `NULL` if this cannot be determined.

Note that some `GDALRasterBands` are not considered to be a part of a dataset,
such as overviews or other "freestanding" bands.
"""
# GDAL wrapper checks null by default, but it is a valid result in this case
getdataset(rb::RasterBand) =
    Dataset(@gdal(GDALGetBandDataset::GDALDataset, rb.ptr::GDALRasterBand))

"""
Return a name for the units of this raster's values. For instance, it might be
"m" for an elevation model in meters, or "ft" for feet.
"""
getunittype(rb::RasterBand) = GDAL.getrasterunittype(rb.ptr)

"""
Set unit type of `rb` to `unittype`.

Values should be one of \"\" (the default indicating it is unknown), \"m\"
indicating meters, or \"ft\" indicating feet, though other nonstandard values
are allowed.
"""
function setunittype!(rb::RasterBand, unitstring::AbstractString)
    result = GDAL.setrasterunittype(rb.ptr, unitstring)
    @cplerr result "Failed to set unit type"
    rb
end

"""
Fetch the raster value offset.

This (in combination with `GetScale()`) is used to transform raw pixel values
into the units returned by `GetUnits()`. For e.g. this might be used to store
elevations in `GUInt16` bands with a precision of 0.1, starting from -100.

    Units value = (raw value * scale) + offset

For file formats that don't know this intrinsically, a value of 0 is returned.
"""
getoffset(rb::RasterBand) = GDAL.getrasteroffset(rb.ptr, C_NULL)

"Set scaling offset."
function setoffset!(rb::RasterBand, value::Real)
    result = GDAL.setrasteroffset(rb.ptr, value)
    @cplerr result "Failed to set scaling offset."
    rb
end

"""
Fetch the raster value scale.

This value (in combination with the `GetOffset()` value) is used to transform
raw pixel values into the units returned by `GetUnits()`. For example this
might be used to store elevations in GUInt16 bands with a precision of 0.1,
and starting from -100.

    Units value = (raw value * scale) + offset

For file formats that don't know this intrinsically a value of one is returned.
"""
getscale(rb::RasterBand) = GDAL.getrasterscale(rb.ptr, C_NULL)

"Set scaling ratio."
function setscale!(rb::RasterBand, ratio::Real)
    result = GDAL.setrasterscale(rb.ptr, ratio)
    @cplerr result "Failed to set scaling ratio"
end

"""
Fetch the no data value for this band.

If there is no out of data value, an out of range value will generally be
returned. The no data value for a band is generally a special marker value
used to mark pixels that are not valid data. Such pixels should generally
not be displayed, nor contribute to analysis operations.

### Returns
the nodata value for this band.
"""
# ### Parameters
# * `pbSuccess`   pointer to a boolean to use to indicate if a value is actually
# associated with this layer. May be `NULL` (default).
getnodatavalue(rb::RasterBand) = GDAL.getrasternodatavalue(rb.ptr, C_NULL)

"Set the no data value for this band."
function setnodatavalue!(rb::RasterBand, value::Real)
    result = GDAL.setrasternodatavalue(rb.ptr, value)
    @cplerr result "Could not set nodatavalue"
    value
end

function deletenodatavalue!(rb::RasterBand)
    result = GDAL.deleterasternodatavalue(rb.ptr)
    @cplerr result "Could not delete nodatavalue"
end

"Set the category names for this band."
function setcategorynames!(rb::RasterBand, names)
    result = @gdal(GDALSetRasterCategoryNames::GDAL.CPLErr,
        rb.ptr::GDALRasterBand,
        names::StringList
    )
    @cplerr result "Failed to set category names"
end

"Fetch the minimum value for this band."
getminimum(rb::RasterBand) = GDAL.getrasterminimum(rb.ptr, C_NULL)

"Fetch the maximum value for this band."
getmaximum(rb::RasterBand) = GDAL.getrastermaximum(rb.ptr, C_NULL)

"Fetch default Raster Attribute Table."
getdefaultRAT(rb::RasterBand) = GDAL.getdefaultrat(rb.ptr)

"""
Set default Raster Attribute Table.

Associates a default RAT with the band. If not implemented for the format a
CPLE_NotSupported error will be issued. If successful a copy of the RAT is made,
the original remains owned by the caller.
"""
function setdefaultRAT!(rb::RasterBand, rat::RasterAttrTable)
    result = GDAL.setdefaultrat(rb.ptr, rat.ptr)
    @cplerr result "Failed to set default raster attribute table"
    rat
end

"""
Copy all raster band raster data.

This function copies the complete raster contents of one band to another
similarly configured band. The source and destination bands must have the same
width and height. The bands do not have to have the same data type.

It implements efficient copying, in particular "chunking" the copy in
substantial blocks.

Currently the only `options` value supported is : "COMPRESSED=YES" to
force alignment on target dataset block sizes to achieve best compression.
More options may be supported in the future.

### Parameters
* `source`        the source band
* `dest`          the destination band
* `options`       transfer hints in "StringList" Name=Value format.
* `progressfunc`  progress reporting function.
* `progressdata`  callback data for progress function.
"""
function copywholeraster!(
        source::RasterBand,
        dest::RasterBand;
        options                 = StringList(C_NULL),
        progressdata            = C_NULL,
        progressfunc::Function  = GDAL.C.GDALDummyProgress
    )
    result = @gdal(GDALRasterBandCopyWholeRaster::GDAL.CPLErr,
        source.ptr::GDALRasterBand,
        dest.ptr::GDALRasterBand,
        options::StringList,
        @cplprogress(progressfunc)::GDALProgressFunc,
        progressdata::Ptr{Void}
    )
    @cplerr result "Failed to copy whole raster"
end

"Return the number of overview layers available, zero if none."
noverview(rb::RasterBand) = GDAL.getoverviewcount(rb.ptr)

"Fetch overview raster band object."
getoverview(rb::RasterBand, i::Integer) =
    RasterBand(GDAL.getoverview(rb.ptr, i))

"""
Fetch best overview satisfying `nsamples` number of samples.

Returns the most reduced overview of the given band that still satisfies the
desired number of samples `nsamples`. This function can be used with zero as the
number of desired samples to fetch the most reduced overview. The same band as
was passed in will be returned if it has not overviews, or if none of the
overviews have enough samples.
"""
function getsampleoverview(rb::RasterBand, nsamples::Integer)
    result = @gdal(GDALGetRasterSampleOverviewEx::GDALRasterBand,
        rb.ptr::GDALRasterBand,
        nsamples::GDAL.GUIntBig
    )
    GDAL.checknull(result)
    RasterBand(result)
end

"Color Interpretation value for band"
getcolorinterp(rb::RasterBand) =
    GDALColorInterp(GDAL.getrastercolorinterpretation(rb.ptr))

"Set color interpretation of a band."
function setcolorinterp!(rb::RasterBand, color::GDALColorInterp)
    result = @gdal(GDALSetRasterColorInterpretation::GDAL.CPLErr,
        rb.ptr::GDALRasterBand,
        color::GDAL.GDALColorInterp
    )
    @cplerr result "Failed to set color interpretation"
    color
end

"""
Fetch the color table associated with band.

If there is no associated color table, the return result is `NULL`. The
returned color table remains owned by the `GDALRasterBand`, and can't be
depended on for long, nor should it ever be modified by the caller.
"""
getcolortable(rb::RasterBand) = ColorTable(GDAL.getrastercolortable(rb.ptr))

"""
Set the raster color table.

The driver will make a copy of all desired data in the colortable. It remains
owned by the caller after the call.

### Parameters
* `colortable` color table to apply (where supported).
"""
function setcolortable!(rb::RasterBand, colortable::ColorTable)
    result = GDAL.setrastercolortable(rb.ptr, colortable.ptr)
    @cplwarn result "CPLError $(result): action is unsupported by the driver"
    colortable
end

function clearcolortable!(rb::RasterBand)
    result = GDAL.setrastercolortable(rb.ptr, GDALColorTable(C_NULL))
    @cplwarn result "CPLError $(result): action is unsupported by the driver"
end

"""
Generate downsampled overviews.

This function will generate one or more overview images from a base image using
the requested downsampling algorithm. Its primary use is for generating
overviews via BuildOverviews(), but it can also be used to generate downsampled
images in one file from another outside the overview architecture.

### Parameters
* `rb`              the source (base level) band.
* `overviewbands`   the list of downsampled bands to be generated.
* `resampling`      (optional) Resampling algorithm (eg. "AVERAGE").
* `progressfunc`    (optional) progress report function.
* `progressdata`    (optional) progress function callback data.

### Additional Remarks
The output bands need to exist in advance.

This function will honour properly `NODATA_VALUES` tuples (special dataset
metadata) so that only a given RGB triplet (in case of a RGB image) will be
considered as the nodata value and not each value of the triplet independantly
per band.
"""
function regenerateoverviews!(
        rb::RasterBand,
        overviewbands::Vector{RasterBand},
        resampling::AbstractString  = "NEAREST",
        progressfunc::Function      = GDAL.C.GDALDummyProgress,
        progressdata                = C_NULL
    )
    result = @gdal(GDALRegenerateOverviews::GDAL.CPLErr,
        rb.ptr::GDALRasterBand,
        length(overviewbands)::Cint,
        GDALRasterBand[band.ptr for band in overviewbands]::Ptr{GDALRasterBand},
        resampling::Cstring,
        @cplprogress(progressfunc)::GDALProgressFunc,
        progressdata::Ptr{Void}
    )
    @cplerr result "Failed to regenerate overviews"
    result
end

# "Advise driver of upcoming read requests."
# _rasteradviseread(hRB::GDALRasterBandH,
#                   nDSXOff::Integer,
#                   nDSYOff::Integer,
#                   nDSXSize::Integer,
#                   nDSYSize::Integer,
#                   nBXSize::Integer,
#                   nBYSize::Integer,
#                   eBDataType::GDALDataType,
#                   papszOptions::Ptr{Ptr{UInt8}}) =
#     GDALRasterAdviseRead(hRB, nDSXOff, nDSYOff, nDSXSize, nDSYSize, nBXSize,
#                          nBYSize, eBDataType, papszOptions)::CPLErr

"""
Fetch the list of category names for this raster.

The return list is a "StringList" in the sense of the CPL functions. That is a
NULL terminated array of strings. Raster values without associated names will
have an empty string in the returned list. The first entry in the list is for
raster values of zero, and so on.
"""
getcategorynames(band::RasterBand) =
    unsafe_loadstringlist(@gdal(GDALGetRasterCategoryNames::Ptr{Cstring},
        band.ptr::GDALRasterBand
    ))

"Set the category names for this band."
function setcategorynames!(band::RasterBand, names::Vector{String})
    result = @gdal(GDALSetRasterCategoryNames::GDAL.CPLErr,
        band.ptr::GDALRasterBand,
        names::StringList
    )
    @cplerr result "Failed to set category names for this band"
    result
end

"""
Flush raster data cache.

This call will recover memory used to cache data blocks for this raster band,
and ensure that new requests are referred to the underlying driver.
"""
function flushcache!(band::RasterBand)
    result = GDAL.flushrastercache(band.ptr)
    @cplerr result "Failed to flush raster data cache"
    result
end

"""
Fill this band with a constant value.

GDAL makes no guarantees about what values pixels in newly created files are
set to, so this method can be used to clear a band to a specified "default"
value. The fill value is passed in as a double but this will be converted to
the underlying type before writing to the file. An optional second argument
allows the imaginary component of a complex constant value to be specified.

### Parameters
* `realvalue`         Real component of fill value
* `imagvalue`    Imaginary component of fill value, defaults to zero
"""
function fillraster!(band::RasterBand, realvalue::Real, imagvalue::Real = 0)
    result = GDAL.fillraster(band.ptr, realvalue, imagvalue)
    @cplerr result "Failed to fill raster band"
end

"""
Return the mask band associated with the band.

The `GDALRasterBand` class includes a default implementation of `GetMaskBand()`
that returns one of four default implementations:

- If a corresponding .msk file exists it will be used for the mask band.
- If the dataset has a `NODATA_VALUES` metadata item, an instance of the new
GDALNoDataValuesMaskBand class will be returned. `GetMaskFlags()` will return
`GMF_NODATA | GMF_PER_DATASET`.
- If the band has a nodata value set, an instance of the new
`GDALNodataMaskRasterBand` class will be returned. `GetMaskFlags()` will
return `GMF_NODATA`.
- If there is no nodata value, but the dataset has an alpha band that seems to
apply to this band (specific rules yet to be determined) and that is of type
`GDT_Byte` then that alpha band will be returned, and the flags
`GMF_PER_DATASET` and `GMF_ALPHA` will be returned in the flags.
- If neither of the above apply, an instance of the new
`GDALAllValidRasterBand` class will be returned that has 255 values for all
pixels. The null flags will return `GMF_ALL_VALID`.

Note that the `GetMaskBand()` should always return a `GDALRasterBand` mask,
even if it is only an all 255 mask with the flags indicating `GMF_ALL_VALID`.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask

### Returns
a valid mask band.
"""
getmaskband(band::RasterBand) = RasterBand(GDAL.getmaskband(band.ptr))

"""
Return the status flags of the mask band associated with the band.

The GetMaskFlags() method returns an bitwise OR-ed set of status flags with the
following available definitions that may be extended in the future:

* `GMF_ALL_VALID` (`0x01`):    There are no invalid pixels, all mask values
will be 255. When used this will normally be the only flag set.
* `GMF_PER_DATASET` (`0x02`):  The mask band is shared between all bands on
the dataset.
- `GMF_ALPHA` (`0x04`):        The mask band is actually an alpha band and may
have values other than 0 and 255.
- `GMF_NODATA` (`0x08`):       Indicates the mask is actually being generated
from nodata values. (mutually exclusive of `GMF_ALPHA`)

The `GDALRasterBand` class includes a default implementation of `GetMaskBand()`
that returns one of four default implementations:

- If a corresponding .msk file exists it will be used for the mask band.
- If the dataset has a `NODATA_VALUES` metadata item, an instance of the new
`GDALNoDataValuesMaskBand` class will be returned. `GetMaskFlags()` will return
`GMF_NODATA | GMF_PER_DATASET`.
- If the band has a nodata value set, an instance of the new
`GDALNodataMaskRasterBand` class will be returned. `GetMaskFlags()` will return
`GMF_NODATA`.
- If there is no nodata value, but the dataset has an alpha band that seems to
apply to this band (specific rules yet to be determined) and that is of type
`GDT_Byte` then that alpha band will be returned, and the flags
`GMF_PER_DATASET` and `GMF_ALPHA` will be returned in the flags.
- If neither of the above apply, an instance of the new `GDALAllValidRasterBand`
class will be returned that has 255 values for all pixels. The null flags will
return `GMF_ALL_VALID`.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask

### Returns
a valid mask band.
"""
getmaskflags(band::RasterBand) = GDAL.getmaskflags(band.ptr)

"""
Adds a mask band to the current band.

The default implementation of the `CreateMaskBand()` method is implemented
based on similar rules to the `.ovr` handling implemented using the
`GDALDefaultOverviews` object. A `TIFF` file with the extension `.msk` will be
created with the same basename as the original file, and it will have as many
bands as the original image (or just one for `GMF_PER_DATASET`). The mask
images will be deflate compressed tiled images with the same block size as the
original image if possible.

If you got a mask band with a previous call to `GetMaskBand()`, it might be
invalidated by `CreateMaskBand()`. So you have to call `GetMaskBand()` again.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask
"""
function createmaskband!(band::RasterBand, nflags::Integer)
    result = GDAL.createmaskband(band.ptr, nflags)
    @cplerr result "Failed to create mask band"
end


text:
# Tests for Radial Multiplier DEA Models
@testset "RadialMultiplierDEAModel" begin

    ## Test Radial Multiplier DEA Models with FLS Book data
    X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17]
    Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12]

    # Input oriented CRS
    deamio = deam(X, Y, orient = :Input, rts = :CRS)

    @test typeof(deamio) == RadialMultiplierDEAModel

    @test nobs(deamio) == 11
    @test ninputs(deamio) == 2
    @test noutputs(deamio) == 1
    @test efficiency(deamio) ≈ [
        1.0000000000;
        0.6222896791;
        0.8198562444;
        1.0000000000;
        0.3103709311;
        0.5555555556;
        1.0000000000;
        0.7576690896;
        0.8201058201;
        0.4905660377;
        1.0000000000] atol = 1e-5

    @test multipliers(deamio, :X) ≈ [
        0.0901826  0.0422374
        0.0314397  0.0414137
        0.0354897  0.0166217
        0.0416228  0.0194942
        0.0274413  0.0361469
        0.0        0.166667
        0.0261969  0.0345077
        0.0178637  0.0235309
        0.0        0.0714286
        0.0133456  0.0175794
        0.2        2.61229e-17
    ] atol = 1e-5
    @test multipliers(deamio, :Y) ≈ [
        0.08333333333333336
        0.04444926279271466
        0.032794249775381853
        0.038461538461538464
        0.03879636638909917
        0.06172839506172836
        0.03703703703703704
        0.02525563631883701
        0.026455026455026464
        0.018867924528301886
        0.0833333333333333
    ] atol = 1e-5

    @test efficiency(deam(targets(deamio, :X), targets(deamio, :Y), orient = :Input, rts = :CRS)) ≈ ones(11) atol = 1e-5
    @test rts(deamio) ≈ zeros(11)

    # Otuput oriented CRS
    deamoo = deam(X, Y, orient = :Output, rts = :CRS)

    @test nobs(deamoo) == 11
    @test ninputs(deamoo) == 2
    @test noutputs(deamoo) == 1
    @test efficiency(deamoo) ≈ [
        1.0000000000;
        1.606968641;
        1.219726027;
        1.0000000000;
        3.221951220;
        1.800000000;
        1.0000000000;
        1.319837398;
        1.219354839;
        2.038461538;
        1.0000000000] atol = 1e-5

    @test multipliers(deamoo, :X) ≈ [
        0.0901826  0.0422374
        0.0505226  0.0665505
        0.0432877  0.020274
        0.0272045  0.0358349
        0.0884146  0.116463
        0.0        0.3
        0.0261969  0.0345077
        0.0235772  0.0310569
        0.0        0.0870968
        0.0272045  0.0358349
        0.2        0.0
    ] atol = 1e-5
    @test multipliers(deamoo, :Y) ≈ [
        0.08333333333333329
        0.07142857142857142
        0.04000000000000001
        0.03846153846153846
        0.12499999999999999
        0.11111111111111112
        0.037037037037037035
        0.03333333333333334
        0.03225806451612903
        0.03846153846153846
        0.08333333333333318
    ] atol = 1e-5

    @test efficiency(deam(targets(deamoo, :X), targets(deamoo, :Y), orient = :Output, rts = :CRS)) ≈ ones(11) atol = 1e-5
    @test rts(deamoo) ≈ zeros(11)

    # Input oriented VRS
    deamiovrs = deam(X, Y, orient = :Input, rts = :VRS)

    @test nobs(deamiovrs) == 11
    @test ninputs(deamiovrs) == 2
    @test noutputs(deamiovrs) == 1
    @test efficiency(deamiovrs) ≈ [
        1.0000000000;
        0.8699861687;
        1.0000000000;
        1.0000000000;
        0.7116402116;
        1.0000000000;
        1.0000000000;
        1.0000000000;
        1.0000000000;
        0.4931209269;
        1.0000000000]

    @test multipliers(deamiovrs, :X) ≈ [
        0.0901826  0.0422374
        0.0197095  0.0570539
        0.0612245  0.000784929
        0.0416228  0.0194942
        0.0185185  0.047619
        0.0        0.166667
        0.0261969  0.0345077
        0.0227671  0.0175131
        0.0        0.0714286
        0.013034   0.0181028
        0.2        2.61229e-17
    ] atol = 1e-5
    @test multipliers(deamiovrs, :Y) ≈ [
        0.08333333333333336
        0.014868603042876911
        0.05259026687598115
        0.038461538461538464
        0.0
        0.03703703703703698
        0.03703703703703704
        0.08756567425569173
        0.07142857142857138
        0.013758146270818254
        0.0833333333333333
    ] atol = 1e-5

    @test efficiency(deam(targets(deamiovrs, :X), targets(deamiovrs, :Y), orient = :Input, rts = :VRS)) ≈ ones(11) atol = 1e-5
    @test rts(deamiovrs) ≈ [0.0; -0.661826; 0.314757; 0.0; -0.711640; -0.666667; 0.0; 1.626970; 1.214286; -0.135409; 0.0] atol = 1e-5

    # Output oriented VRS
    deamoovrs = deam(X, Y, orient = :Output, rts = :VRS)

    @test nobs(deamoovrs) == 11
    @test ninputs(deamoovrs) == 2
    @test noutputs(deamoovrs) == 1
    @test efficiency(deamoovrs) ≈ [
        1.0000000000;
        1.507518797;
        1.0000000000;
        1.0000000000;
        3.203947368;
        1.000000000;
        1.0000000000;
        1.000000000;
        1.000000000;
        1.192307692;
        1.0000000000]

    @test multipliers(deamoovrs, :X) ≈ [
        0.0901826   0.0422374
        0.0676692   0.093985
        0.0465672   0.000597015
        0.0272045   0.0358349
        0.118421    0.164474
        0.0         0.5
        0.0261969   0.0345077
        0.00866667  0.00666667
        0.0         0.0322581
        0.0         0.0
        0.2         0.0
    ] atol = 1e-5
    @test multipliers(deamoovrs, :Y) ≈ [
        0.08333333333333329
        0.07142857142857144
        0.040000000000000015
        0.03846153846153846
        0.125
        0.11111111111111109
        0.037037037037037035
        0.03333333333333333
        0.03225806451612903
        0.038461538461538464
        0.08333333333333318
    ] atol = 1e-5

    @test efficiency(deam(targets(deamoovrs, :X), targets(deamoovrs, :Y), orient = :Output, rts = :VRS)) ≈ ones(11) atol = 1e-5
    @test rts(deamoovrs) ≈ [0.0; -0.703008; 0.239403; 0.0; -1.230263; -2.0; 0.0; 0.619333; 0.548387; 1.192308; 0.0] atol = 1e-5

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    deamio_ref_eff = zeros(size(X, 1))
    deamoo_ref_eff = zeros(size(X, 1))

    deamiovrs_ref_eff = zeros(size(X, 1))
    deamoovrs_ref_eff = zeros(size(X, 1))

    deamiovrs_ref_multX = zeros(size(X))
    deamiovrs_ref_multY = zeros(size(Y))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        deamio_ref_eff[i] = efficiency(deam(Xeval, Yeval, orient = :Input, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        deamoo_ref_eff[i] = efficiency(deam(Xeval, Yeval, orient = :Output, rts = :CRS, Xref = Xref, Yref = Yref))[1]

        deamiovrs_ref_eff[i] = efficiency(deam(Xeval, Yeval, orient = :Input, rts = :VRS, Xref = Xref, Yref = Yref))[1]
        deamoovrs_ref_eff[i] = efficiency(deam(Xeval, Yeval, orient = :Output, rts = :VRS, Xref = Xref, Yref = Yref))[1]

        deamiovrs_ref_multX[i,:] = multipliers(deam(Xeval, Yeval, orient = :Input, rts = :VRS, Xref = Xref, Yref = Yref), :X)
        deamiovrs_ref_multY[i,:] = multipliers(deam(Xeval, Yeval, orient = :Input, rts = :VRS, Xref = Xref, Yref = Yref), :Y)
    end

    @test deamio_ref_eff ≈ efficiency(deamio)
    @test deamoo_ref_eff ≈ efficiency(deamoo)

    @test deamiovrs_ref_eff ≈ efficiency(deamiovrs)
    @test deamoovrs_ref_eff ≈ efficiency(deamoovrs)

    @test deamiovrs_ref_multX ≈ multipliers(deamiovrs, :X)
    @test deamiovrs_ref_multY ≈ multipliers(deamiovrs, :Y)

    # Print
    show(IOBuffer(), deamio)

    # Test errors
    @test_throws DimensionMismatch deam([1; 2 ; 3], [4 ; 5]) #  Different number of observations
    @test_throws DimensionMismatch deam([1; 2], [4 ; 5], Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws DimensionMismatch deam([1 1; 2 2], [4 4; 5 5], Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws DimensionMismatch deam([1 1; 2 2], [4 4; 5 5], Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ArgumentError deam([1; 2; 3], [4; 5; 6], orient = :Error) # Invalid orientation
    @test_throws ArgumentError deam([1; 2; 3], [4; 5; 6], rts = :Error) # Invalid returns to scale

    @test_throws ArgumentError targets(deamio, :Error)    # Invalid target
    @test_throws ErrorException slacks(deamio, :X)        # No slacks in multiplier model
    @test_throws ArgumentError multipliers(deamio, :Error) # Invalid slacks

    # ------------------
    # Test Vector and Matrix inputs and outputs
    # ------------------
    # Tests against results in R

    # Inputs is Matrix, Outputs is Vector
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6	8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]

    @test efficiency(deam(X, Y, orient = :Input)) ≈ [1; 1; 1; 0.6; 0.4; 1; 0.6666666667; 0.625]

    # Inputs is Vector, Output is Matrix
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]

    @test efficiency(deam(X, Y, orient = :Output)) ≈ [1; 1; 1; 1.555555556; 2.333333333; 1; 1.272727273; 1.6]

    # Inputs is Vector, Output is Vector
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]

    @test efficiency(deam(X, Y, orient = :Input)) ≈ [0.4; 1; 0.8; 0.6; 0.4; 0.4; 0.5142857143; 0.2]

end


text:
using PyPlot
using PyCall
mpl = pyimport("tikzplotlib")
using JLD2
using NNFEM
include("../nnutil.jl")
nntype= "ae_scaled"
stress_scale = 100.0

tid = 3
@load "../Data/domain$tid.jld2" domain
@load "../Data/learn_domain$(tid)_te.jld2" domain_te

u1 = hcat(domain.history["state"]...)
u2 = hcat(domain_te.history["state"]...)
err = zeros(8, 101)
for i = 1:8
    u1_ = [u1[i,:] u1[i+8,:]]
    u2_ = [u2[i,:] u2[i+8,:]]
    err[i,:] = sqrt.(sum((u1_-u2_).^2,dims=2)[:])
end

T = 0.5
NT = 100
t = LinRange(0.0,T, NT+1)
for i = 1:8
    plot(t, err[i,:])
end
xlabel("t")
ylabel("\$u_{ref}-u_{est}\$")
mpl.save("truss2d_learn_loc_diff$tid.tex")

@load "Data/domain3.jld2" domain
X, Y = prepare_strain_stress_data1D(domain)
y = squeeze(nn(constant(X[:,1]), constant(X[:,2]), constant(X[:,3])))
sess = Session(); init(sess)
close("all")
ADCME.load(sess, "Data/learned_nn.mat")
out = run(sess, y)
plot(X[:,1], out,"+", label="NN")
plot(X[:,1], Y, ".", label="Exact")
legend()


text:
@testmodule ThresholdSamplingTest begin
  using MLLibrary
  function test_ThresholdTrainingIndices()

    resources=["abalone.data","3percent-miscategorization.csv","10percent-miscatergorization.csv","Seperable.csv"]
    for res in resources
      path = joinpath(Pkg.dir("MLLibrary"),"resources","datasets",res)
      count = MLLibrary.DataParser.DataStream.getDataSetCount(path)

      threshold=.1
      srand(123)
      trainingset = MLLibrary.DataSampling.ThresholdSampling.getThresholdTrainingSetIndicies(path,count,threshold)
      @test length(trainingset)==round(count*(1-threshold))
    end
  end

  function test_ThresholdTestIndicies()
    resources=["abalone.data","3percent-miscategorization.csv","10percent-miscatergorization.csv","Seperable.csv"]
    for res in resources
      path = joinpath(Pkg.dir("MLLibrary"),"resources","datasets",res)
      count = MLLibrary.DataParser.DataStream.getDataSetCount(path)

      threshold=.1
      srand(123)
      testset = MLLibrary.DataSampling.ThresholdSampling.getThresholdTestSetIndices(path,count,threshold)
      @test length(testset)==round(count*(threshold))
    end

  end

  function test_unionAndInterSection()
    resources=["abalone.data","3percent-miscategorization.csv","10percent-miscatergorization.csv","Seperable.csv"]
    for res in resources
      path = joinpath(Pkg.dir("MLLibrary"),"resources","datasets",res)
      count = MLLibrary.DataParser.DataStream.getDataSetCount(path)

      threshold=.1
      srand(123)
      trainingset = MLLibrary.DataSampling.ThresholdSampling.getThresholdTrainingSetIndicies(path,count,threshold)
      srand(123)
      testset = MLLibrary.DataSampling.ThresholdSampling.getThresholdTestSetIndices(path,count,threshold)

      @test length(intersect(trainingset,testset))==0
      @test length(union(trainingset,testset))==count
      @test in(0,trainingset)==false
      @test in(0,testset)==false

    end
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
"""
    update_runs!(array, runs, n)

Add a sequence which start with the first element `runs` and decreases
by `dec` to the first `n` elements of the vector `array`.
"""
function update_runs!(array :: AbstractVector{T}, runs, n) where T <: Number
    array[1:n] .+= take(countfrom(runs, -1), n)
end

function update_runs_periodic!(array :: AbstractVector{T},
                               left  :: Integer,
                               right :: Integer,
                               len   :: Integer) where T <: Number
    sum = left + right
    nupdate = min(len, sum)
    f(n) = min(n - 1, left, right, sum - (n - 1))
    array[1:nupdate] .+= (f(n) for n in 1:nupdate)
end

struct CountRuns{A, T}
    array :: A
    phase :: T
end

CountRuns(array :: AbstractArray{T}, phase :: T) where T =
    CountRuns{typeof(array), T}(array, phase)

function Base.iterate(counter :: CountRuns, state = nothing)
    idx = isnothing(state) ? 1 : state
    array = counter.array
    phase = counter.phase
    len = length(array)

    start = findfirst(x -> x == phase, view(array, idx:len))
    isnothing(start) && return nothing
    start = start + idx - 1

    stop = findfirst(x -> x != phase, view(array, start:len))
    stop = isnothing(stop) ? len + 1 : stop + start - 1

    return stop - start, stop
end

Base.IteratorSize(:: CountRuns) = Base.SizeUnknown()

"""
    l2(array, phase; [len][, directions,] periodic = false)

Calculate `L₂` (lineal path) correlation function for one-, two- or
three-dimensional multiphase system.

`L₂(x)` equals to probability that all elements of a line segment with
length `x` cut from the array belong to the same phase. This
implementation calculates `L₂(x)` for all `x`es in the range from `1`
to `len` which defaults to half of the minimal dimension of the array.

# Examples
```jldoctest
julia> l2([1,1,1,0,1,1], 1; len = 6)[:x]
6-element Array{Float64,1}:
 0.8333333333333334
 0.6
 0.25
 0.0
 0.0
 0.0
```

For a list of possible dimensions, see also: [`direction1Dp`](@ref),
[`direction2Dp`](@ref), [`direction3Dp`](@ref).
"""
function l2(array      :: AbstractArray,
            phase;
            len        :: Integer        = (array |> size |> minimum) ÷ 2,
            directions :: Vector{Symbol} = array |> default_directions,
            periodic   :: Bool           = false,
            plans      :: Any            = nothing)
    cd = CorrelationData(len, check_directions(directions, size(array), periodic))

    for direction in directions
        slicer = slice_generators(array, periodic, Val(direction))

        for slice in slicer
            slen = length(slice)
            firstrun = 0
            lastrun = 0

            # Update count of slices which satisfy "all elements
            # belong to phase" condition
            for run in CountRuns(slice, phase)
                firstrun = (firstrun == 0) ? run : firstrun
                lastrun = run
                update_runs!(cd.success[direction], run, min(run, len))
            end

            total_updates = min(slen, len)
            if periodic
                if (slice[begin] == slice[end] == phase)
                    update_runs_periodic!(cd.success[direction],
                                          firstrun, lastrun,
                                          total_updates)
                end

                # Update total number of slices
                cd.total[direction][1:total_updates] .+= slen
            else
                # Calculate total number of slices with lengths from 1
                # to len
                update_runs!(cd.total[direction], slen, total_updates)
            end
        end
    end

    return cd
end


text:
# Need https://github.com/JuliaLang/julia/pull/33970
# and  https://github.com/JuliaLang/julia/pull/34043
if VERSION < v"1.4.0-DEV.666"
    exit()
end

using CUDA
if capability(device()) < v"7.0"
    exit()
end

is_debug = ccall(:jl_is_debugbuild, Cint, ()) != 0
(is_debug && VERSION < v"1.5.0-DEV.437") && exit()

### START
using Test

using CUDA

a     = rand(Float16, (16, 16))
b     = rand(Float16, (16, 16))
c     = rand(Float32, (16, 16))

a_dev = CuArray(a)
b_dev = CuArray(b)
c_dev = CuArray(c)
d_dev = similar(c_dev)

function kernel(a_dev, b_dev, c_dev, d_dev)
    conf = WMMA.Config{16, 16, 16, Float32}

    a_frag = WMMA.load_a(pointer(a_dev), 16, WMMA.ColMajor, conf)
    b_frag = WMMA.load_b(pointer(b_dev), 16, WMMA.ColMajor, conf)
    c_frag = WMMA.load_c(pointer(c_dev), 16, WMMA.ColMajor, conf)

    c_frag = 0.5f0 .* c_frag

    d_frag = WMMA.mma(a_frag, b_frag, c_frag, conf)

    WMMA.store_d(pointer(d_dev), d_frag, 16, WMMA.ColMajor, conf)

    return
end

@cuda threads=32 kernel(a_dev, b_dev, c_dev, d_dev)
d = Array(d_dev)

@test all(isapprox.(a * b + 0.5 * c, d; rtol=0.01))
### END


text:
function setup_rhs(bulkconstrains::BulkPartition{Nsys}, boundary::Boundary,
                   bulkderivs::BulkPartition{Nsys},
                   systems::SystemPartition, integration::Integration) where {Nsys}
    # function to solve the nested system
    nested = Nested(systems, bulkconstrains, bulkderivs)

    function (ff_t::EvolVars, ff::EvolVars, evoleq::EvolutionEquations, t)
        bulkevols_t = getbulkevolvedpartition(ff_t)
        bulkevols   = getbulkevolvedpartition(ff)

        # solve nested system for the constrained variables
        nested(bulkevols, boundary, evoleq)

        # first u-domain
        sys           = systems[1]
        bulkevol_t    = bulkevols_t[1]
        bulkevol      = bulkevols[1]
        bulkconstrain = bulkconstrains[1]
        deriv         = bulkderivs[1]
        compute_bulkevolved_t_1st!(bulkevol_t, bulkconstrain, bulkevol, deriv,
                                   sys, evoleq)
        # remaining u-domains
        @inbounds for aa in 2:Nsys
            sys           = systems[aa]
            bulkevol_t    = bulkevols_t[aa]
            bulkevol      = bulkevols[aa]
            bulkconstrain = bulkconstrains[aa]
            deriv         = bulkderivs[aa]

            compute_bulkevolved_t!(bulkevol_t, bulkconstrain, bulkevol, deriv,
                                   sys, evoleq)
        end
        sync_bulkevolved!(bulkevols_t, bulkconstrains, systems, evoleq)

        nothing
    end
end


text:
"""
    inverse_flash_update!(storage, eos, c, V)

Update internal matrix of partial derivatives for a converged flash result.
"""
function inverse_flash_update!(storage, eos, c, V)
    x, y = storage.x, storage.y
    J_p, = update_primary_jacobian!(storage, eos, c, storage.forces, x, y, V)
    J_s, = update_secondary_jacobian!(storage, eos, c, x, y, V)
    if isa(J_p, MMatrix)
        tmp = SMatrix(J_p)\SMatrix(J_s)
        J_s .= tmp
    else
        F = lu!(J_p)
        ldiv!(F, J_s)
    end
    return J_s
end

"""
    set_partials(v, storage, eos, c, index)

Modify a value `v` of ForwardDiff.Dual type to get the correct derivatives.

If the mixture is two-phase, the partial derivatives of phase molar fractions and vapor fraction
with respect to the flash conditions (pressure, temperature and overall mole fractions) are not
trivial. This function is the main gateway for setting these values.

# Notes
Experimental interface, subject to change. You most likely want to use either
[`set_partials_phase_mole_fractions!`](@ref) or [`set_partials_vapor_fraction`](@ref)

In order for this routine to work, the storage must be initialized using `flash_storage`
with the following options enabled:
- `inc_jac = true`
- `diff_externals = true`
and `inverse_flash_update!` must be called after a successful flash. The partial derivatives
with respect to p, T, z is then contained in `storage.buf_inv` with a negative sign. This
function then performs the requisite chain rule operations for the input.
"""
function set_partials(v, storage, eos, c, index)
    M = storage.J_inv
    p, z, T = c.p, c.z, c.T
    buf = storage.buf_inv
    @. buf = 0
    is_dual(x) = typeof(x)<:ForwardDiff.Dual
    if is_dual(p)
        @inbounds @. buf -= M[index, 1].*p.partials
    end
    if is_dual(T)
        @inbounds @. buf -= M[index, 2].*T.partials
    end
    if eltype(z)<:ForwardDiff.Dual
        for (ind, zi) = enumerate(z)
            @inbounds ∂z = -M[index, ind+2]
            for i in eachindex(zi.partials)
                @inbounds buf[i] += ∂z*zi.partials[i]
            end
        end
    end
    ∂T = typeof(v)
    P = typeof(v.partials)

    return ∂T(v.value, P(Tuple(buf)))
end

"""
    set_partials_vapor_fraction(V, storage, eos, ∂c)

Set partial derivatives to a vapor mole fraction `ForwardDiff.Dual` instance with correct value, but
missing partial derivatives.
"""
set_partials_vapor_fraction(V, storage, eos, ∂c) = set_partials(V, storage, eos, ∂c, 2*number_of_components(eos) + 1)

"""
set_partials_phase_mole_fractions!(xy, storage, eos, ∂c, [phase_symbol])

Set partial derivatives to phase mole fraction vector with type `ForwardDiff.Dual` with correct values, but
missing partial derivatives.
"""
function set_partials_phase_mole_fractions!(xy, storage, eos, ∂c, phase = :liquid)
    n = length(xy)
    if phase == :liquid
        offset = 0
    else
        offset = n
    end
    for i = 1:n
        xy[i] = set_partials(xy[i], storage, eos, ∂c, offset + i)
    end
    return xy
end


text:
include("basic_setup.jl")
using Statistics

function estimator_variance(num_responses, integration_nodes, chain_length)
    algorithm = TIEstimate(1024, integration_nodes, chain_length)
    result = Trajectories.marginal_entropy(gen, algorithm=algorithm; num_responses=num_responses, duration=100.0)
    estimator_var = var(result["marginal_entropy"].Sample) / (num_responses - 1)
    time = sum(result["marginal_entropy"].TimeElapsed)
    (estimator_var, time)
end

integration_nodes = 4

chains = map(x -> round(Int, x), 2 .^ (8:0.25:11))

results = estimator_variance.(100, integration_nodes, chains)
times = map(x->x[2], results)
variance = map(x->x[1], results)

using DataFrames, GLM
data = DataFrame(N=chains, Time=times, Variance=variance)

ols = lm(@formula(Time ~ N), data)
print(ols)

fixtime = coef(ols)[1] / 100
tau_s = coef(ols)[2] / 100

nr(ns) = round(Int, 10 * 60 / (fixtime + ns * tau_s))

chains2 = map(x -> round(Int, x), 2 .^ (9:0.125:12))
nr.(chains2)

results2 = estimator_variance.(nr.(chains2), integration_nodes, chains2)

using Plots
plot(data.N, [data.Time, predict(ols)])

data2 = DataFrame(N=chains2, Time=map(x->x[2], results2), Variance=map(x->x[1], results2))

plot(data2.N, data2.Variance, xlim=(0, :auto), ylim=(0, :auto))

using DrWatson
using CSV

mkpath(projectdir("data", "estimator_variance"))
CSV.write(projectdir("data", "estimator_variance", "integration_nodes=4.csv"), data2)
CSV.write(projectdir("data", "estimator_variance", "timing_integration_nodes=4.csv"), data)


text:
function convert_fields(coalesced::Vector{<:LibClang.CXCursor}, indent::Int, func::Union{Function, Nothing} = nothing)
	comments = Dict{String, Comment}()
	common = nothing
	expr = nothing
	flds = map(coalesced) do fld
		fldName = jlname = convert_name(fld)
		if isnothing(func)
			fld.kind in (
				LibClang.CXCursor_FieldDecl,
			) || return nothing
		else
			fld.kind in (
				LibClang.CXCursor_VarDecl,
				LibClang.CXCursor_FunctionDecl,
			) || return nothing

			LibClang.clang_Cursor_getStorageClass(fld) in (
				LibClang.CX_SC_None,
				LibClang.CX_SC_Extern,
			) || return nothing

			jlname = func(fld)
			jlname === false && return nothing
			jlname = jlname isa AbstractString ? jlname : fldName
		end

		((pre, typ, post), fldComments) = fld.kind == LibClang.CXCursor_FunctionDecl ? convert_function(fld, indent) : convert_type(fld, LibClang.clang_getCursorType(fld), indent)

		isnothing(common) || Bool(LibClang.clang_equalTypes(common, typ)) || error("Unexpected nested type")
		if isnothing(common)
			common = typ
			cvt = convert_nested_type(typ, fld, indent)
			expr = cvt.expr
			merge_comments!(comments, cvt.comments)
		end
		merge_comments!(comments, fldComments)

		if !isnothing(func)
			merge_comment!(comments, convert_comment(fld, jlname))
			jlname = jlname == fldName ? jlname : "($(fldName) => $(jlname))"
		end

		return (fld, jlname, pre, post)
	end
	flds = filter(!isnothing, flds)

	if isnothing(expr)
		if coalesced[1].kind == LibClang.CXCursor_TypedefDecl
			convert_decl = convert_typedef
		elseif coalesced[1].kind == LibClang.CXCursor_EnumDecl
			convert_decl = convert_enum
		elseif coalesced[1].kind == LibClang.CXCursor_StructDecl
			convert_decl = convert_struct
		elseif coalesced[1].kind == LibClang.CXCursor_UnionDecl
			convert_decl = convert_union
		else
			return nothing
		end

		return convert_decl(coalesced[1], indent)
	end

	if length(flds) == 1
		(fld, jlname, pre, post) = flds[1]
		if startswith(expr, "𝐣𝐥.@") && (
			startswith(post, '[') ||                       # @cstruct {...}[N]
			(endswith(pre, '{') && startswith(post, ','))  # Cfunction{@cstruct {...}, ...}
		)
			expr = "($(expr))"
		end
		colons = fld.kind == LibClang.CXCursor_FunctionDecl ? "" : "::"
		expr = "$(jlname)$(colons)$(pre)$(expr)$(post)"
	else
		flds = map(flds) do (fld, jlname, pre, post)
			if !isempty(pre*post) && (
				startswith(post, '[') ||                          # _[N]
				(endswith(pre, "::") && isempty(post)) ||         # f()::_
				(endswith(pre, '(') && startswith(post, ')')) ||  # Cconst(_)
				(endswith(pre, '{') && startswith(post, '}')) ||  # Ptr{_}
				(endswith(pre, '{') && startswith(post, ','))     # Cfunction{_, ...}
			)
				pre = pre*"_"
			end
			colons = fld.kind == LibClang.CXCursor_FunctionDecl ? "" : "::"
			return isempty(pre*post) ? jlname : "$(jlname)$(colons)$(pre)$(post)"
		end
		expr = "($(join(flds, ", ")))::$(expr)"
	end

	if !isnothing(func)
		expr = "𝐣𝐥.@cextern $(expr)"
	end

	return Converted(
		expr,
		comments,
	)
end


text:
# constructors
voidptr_a = VEPtr{Cvoid}(Int(0xDEADBEEF))
@test reinterpret(Ptr{Cvoid}, voidptr_a) == Ptr{Cvoid}(Int(0xDEADBEEF))

# getters
@test eltype(voidptr_a) == Cvoid

# comparisons
voidptr_b = VEPtr{Cvoid}(Int(0xCAFEBABE))
@test voidptr_a != voidptr_b

@testset "conversions" begin

# between host and device pointers
@test_throws ArgumentError convert(Ptr{Cvoid}, voidptr_a)

# between VE pointers
intptr_a = VEPtr{Int}(Int(0xDEADBEEF))
@test convert(typeof(intptr_a), voidptr_a) == intptr_a

# convert back and forth from UInt
intptr_b = VEPtr{Int}(Int(0xDEADBEEF))
@test convert(UInt, intptr_b) == 0xDEADBEEF
@test convert(VEPtr{Int}, Int(0xDEADBEEF)) == intptr_b
@test Int(intptr_b) == Int(0xDEADBEEF)

# pointer arithmetic
intptr_c = VEPtr{Int}(Int(0xDEADBEEF))
intptr_d = 2 + intptr_c
@test isless(intptr_c, intptr_d)
@test intptr_d - intptr_c == 2
@test intptr_d - 2 == intptr_c

end

@testset "VE or CPU integration" begin

a = [1]
ccall(:clock, Nothing, (Ptr{Int},), a)
@test_throws Exception ccall(:clock, Nothing, (VEPtr{Int},), a)
ccall(:clock, Nothing, (PtrOrVEPtr{Int},), a)

end

@testset "reference values" begin
    # Ref

    @test typeof(Base.cconvert(Ref{Int}, 1)) == typeof(Ref(1))
    @test Base.unsafe_convert(Ref{Int}, Base.cconvert(Ref{Int}, 1)) isa Ptr{Int}

    ptr = reinterpret(Ptr{Int}, C_NULL)
    @test Base.cconvert(Ref{Int}, ptr) == ptr
    @test Base.unsafe_convert(Ref{Int}, Base.cconvert(Ref{Int}, ptr)) == ptr

    arr = [1]
    @test Base.cconvert(Ref{Int}, arr) isa Base.RefArray{Int, typeof(arr)}
    @test Base.unsafe_convert(Ref{Int}, Base.cconvert(Ref{Int}, arr)) == pointer(arr)

    # VERef

    @test typeof(Base.cconvert(VERef{Int}, 1)) == typeof(VERef(1))
    @test Base.unsafe_convert(VERef{Int}, Base.cconvert(VERef{Int}, 1)) isa VERef{Int}

    veptr = reinterpret(VEPtr{Int}, C_NULL)
    @test Base.cconvert(VERef{Int}, veptr) == veptr
    @test Base.unsafe_convert(VERef{Int}, Base.cconvert(VERef{Int}, veptr)) == Base.bitcast(VERef{Int}, veptr)

    # TODO
    #vearr = oneAPI.oneArray([1])
    #@test Base.cconvert(VERef{Int}, vearr) isa oneL0.VERefArray{Int, typeof(vearr)}
    #@test Base.unsafe_convert(VERef{Int}, Base.cconvert(VERef{Int}, vearr)) == Base.bitcast(VERef{Int}, pointer(vearr))

    # RefOrVERef

    @test typeof(Base.cconvert(RefOrVERef{Int}, 1)) == typeof(Ref(1))
    @test Base.unsafe_convert(RefOrVERef{Int}, Base.cconvert(RefOrVERef{Int}, 1)) isa RefOrVERef{Int}

    @test Base.cconvert(RefOrVERef{Int}, ptr) == ptr
    @test Base.unsafe_convert(RefOrVERef{Int}, Base.cconvert(RefOrVERef{Int}, ptr)) == Base.bitcast(RefOrVERef{Int}, ptr)

    @test Base.cconvert(RefOrVERef{Int}, veptr) == veptr
    @test Base.unsafe_convert(RefOrVERef{Int}, Base.cconvert(RefOrVERef{Int}, veptr)) == Base.bitcast(RefOrVERef{Int}, veptr)

    @test Base.cconvert(RefOrVERef{Int}, arr) isa Base.RefArray{Int, typeof(arr)}
    @test Base.unsafe_convert(RefOrVERef{Int}, Base.cconvert(RefOrVERef{Int}, arr)) == Base.bitcast(RefOrVERef{Int}, pointer(arr))

    # TODO
    #@test Base.cconvert(RefOrVERef{Int}, vearr) isa oneL0.VERefArray{Int, typeof(vearr)}
    #@test Base.unsafe_convert(RefOrVERef{Int}, Base.cconvert(RefOrVERef{Int}, vearr)) == Base.bitcast(RefOrVERef{Int}, pointer(vearr))
end


text:
using Random
Random.seed!(7)

dt = 0.01
integration = :rk4

###################
## Parallel Park ##
###################

opts = SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-4
opts.resolve_feasible = true
opts.outer_loop_update_type = :default
opts.R_infeasible = 10

# Set up model, objective, and solver
model, obj = TrajectoryOptimization.Dynamics.dubinscar_parallelpark
solver = Solver(model, obj, integration=integration, dt=dt, opts=opts)
U0 = ones(solver.model.m,solver.N-1)
X0 = line_trajectory(solver)

results, stats = TrajectoryOptimization.solve(solver,U0)
results_inf, stats_inf = TrajectoryOptimization.solve(solver,X0,U0)

@test norm(results.X[end]-obj.xf) < 1e-4
@test TrajectoryOptimization.max_violation(results) < 1e-4
@test norm(results_inf.X[end]-obj.xf) < 1e-4
@test TrajectoryOptimization.max_violation(results_inf) < 1e-4

# # Parallel Park (boxed)
# x_min = obj.x_min
# x_max = obj.x_max
# plt = plot(title="Parallel Park")#,aspect_ratio=:equal)
# plot!(x_min[1]*ones(1000),collect(range(x_min[2],stop=x_max[2],length=1000)),color=:red,width=2,label="")
# plot!(x_max[1]*ones(1000),collect(range(x_min[2],stop=x_max[2],length=1000)),color=:red,width=2,label="")
# plot!(collect(range(x_min[1],stop=x_max[1],length=1000)),x_min[2]*ones(1000),color=:red,width=2,label="")
# plot!(collect(range(x_min[1],stop=x_max[1],length=1000)),x_max[2]*ones(1000),color=:red,width=2,label="")
# plot_trajectory!(to_array(results.X),width=2,color=:green,label="Constrained",legend=:bottomright)
# plot_trajectory!(to_array(results_inf.X),width=2,color=:yellow,label="Constrained (infeasible)",legend=:bottomright)

########################
## Obstacle Avoidance ##
########################

model, obj_con_obstacles, circles = TrajectoryOptimization.Dynamics.dubinscar_obstacles
model, obj_con_obstacles_control, = TrajectoryOptimization.Dynamics.dubinscar_obstacles_control_limits
solver_con_obstacles = Solver(model, obj_con_obstacles_control, integration=integration, dt=dt, opts=opts)

# -Initial state and control trajectories
U0 = ones(solver_con_obstacles.model.m,solver_con_obstacles.N-1)
X0 = line_trajectory(solver_con_obstacles)
results_inf, stats_inf = TrajectoryOptimization.solve(solver_con_obstacles,X0,U0)

# plt = plot(title="Obstacle Avoidance")
# plot_obstacles(circles)
# plot_trajectory!(to_array(results_con_obstacles.X),width=2,color=:orange,label="Infeasible")
# plot(to_array(results.U)')

@test norm(results_inf.X[end]-obj_con_obstacles.xf) < 1e-4
@test TrajectoryOptimization.max_violation(results_inf) < 1e-4

############
## Escape ##
############

N_escape = 101
model, obj_escape, circles_escape = TrajectoryOptimization.Dynamics.dubinscar_escape
solver_escape = Solver(model, obj_escape, integration=integration, N=N_escape, opts=opts)
X_guess = [2.5 2.5 0.;4. 5. .785;5. 6.25 0.;7.5 6.25 -.261;9 5. -1.57;7.5 2.5 0.]
X0 = TrajectoryOptimization.interp_rows(solver_escape.N,solver_escape.obj.tf,Array(X_guess'))
U0 = ones(solver_escape.model.m,solver_escape.N-1)

solver_escape.opts.R_infeasible = 1e-1
solver_escape.opts.resolve_feasible = false
solver_escape.opts.cost_tolerance = 1e-6
solver_escape.opts.cost_tolerance_intermediate = 1e-3
solver_escape.opts.constraint_tolerance = 1e-4#1e-5
solver_escape.opts.constraint_tolerance_intermediate = 0.01
solver_escape.opts.penalty_scaling = 100.0
solver_escape.opts.penalty_initial = 100.0
solver_escape.opts.outer_loop_update_type = :default
solver_escape.opts.iterations_outerloop = 20
solver_escape.opts.use_penalty_burnin = false
solver_escape.opts.verbose = false
solver_escape.opts.live_plotting = false
results_escape, stats_escape = solve(solver_escape,X0,U0)

# plt = plot(title="Escape",aspect_ratio=:equal)
# plot_obstacles(circles_escape)
# plot_trajectory!(to_array(results_escape.X),width=2,color=:purple,label="Infeasible",aspect_ratio=:equal,xlim=[-1,11],ylim=[-1,11])
# plot!(X0[1,:],X0[2,:],label="Infeasible Initialization",width=1,color=:purple,linestyle=:dash)
# display(plt)
# plot(to_array(results_escape.U[1:solver_escape.N-1])',label="")

@test norm(results_escape.X[end]-solver_escape.obj.xf) < 1e-4
@test TrajectoryOptimization.max_violation(results_escape) < 1e-4


text:
# load packages
using Plots
#pyplot()
gr()

# load modules
# select a model later
include("HexaLatticeModel02a.jl")
include("SquareLatticeModel02a.jl")
# MCMC
include("IsingMCMC02a.jl")

###################################################################################################
# view functions
function drawEnergy(sim, trial, T, d=d, N=N)
    xs = [t for t in 1:1000:trial]
    ys = [E(sim[t],N) for t in xs]
    p = Plots.plot(xs, ys, label="Energy history with T=$(T)")
end

function drawSim(sim, f, d=d, N=N)
    snap = sim[f]

    p=Plots.plot([],[],legend=false, size=(500,500))

    p0xs = [];p0ys = []
    p1xs = [];p1ys = []
    for i in 1:N
        # bits
        (xs,ys) = genPoly(i, d)
        if (snap[i] == 1) # up-spin
            p1xs = vcat(p1xs, NaN, xs)
            p1ys = vcat(p1ys, NaN, ys)
        else # down-spin
            p0xs = vcat(p0xs, NaN, xs)
            p0ys = vcat(p0ys, NaN, ys)
        end
    end
    poly0 = Shape(p0xs, p0ys) # creating polygon
    p = Plots.plot!(poly0, fillcolor = plot_color(:blue)  ,legend=false)
    poly1 = Shape(p1xs, p1ys)
    p = Plots.plot!(poly1, fillcolor = plot_color(:yellow),legend=false)
end

###################################################################################################
# prepare
# setting for energy function
function getJ(i, j, d=d)
    J0 = 1
    (j in getNeighbors(i,d)) ? J0 : 0
end

function getB(i, d=d)
    B0 = 0
    B0
end

###################################################################################################
# do simulation

# setting for model
d = SquareLatticeModel.ParamDict(64,64,0)
#d = SquareLatticeModel.ParamDict(32,32,0)
#d = SquareLatticeModel.ParamDict(8,8,0)
N = d.nrows * d.ncols
d.N = N

# main function
# mdlno: 1(square), 2(hexa)
# T: temperature e.g. T=2.26
function main(mdlno, T)
    # setting for simulation
    trial = 100000 # MCMC trial

    # gen functors as global
    global getNeighbors, genPoly
    if mdlno == 1
        (getNeighbors, genPoly)= SquareLatticeModel.genFunctor(d)
    elseif mdlno == 2
        (getNeighbors, genPoly)= HexaLatticeModel.genFunctor(d)
    end

    # simulation
    global E, MCMC
    (E, MCMC) = IsingMCMC.genFunctor(getJ, getB)
    print("start MCMC...")
    @time sim = MCMC(T, N, trial);

    # save images
    # energy history
    print("start drawEnergy...")
    @time p = drawEnergy(sim, trial, T)
    Plots.savefig(p, "Ising_energy_M$(mdlno)_T$(T).png")

    # cell state
    print("start drawSim...")
    @time p = drawSim(sim, trial, d)
    Plots.savefig(p, "Ising_sim_M$(mdlno)_T$(T).png")
    print("done.")
end
###################################################################################################


text:
function getEndpoints(distr::UnivariateDistribution, epsilon::Real)
    (left,right) = map(x -> quantile(distr,x), (0,1))
    leftEnd = left!=-Inf ? left : quantile(distr, epsilon)
    rightEnd = right!=-Inf ? right : quantile(distr, 1-epsilon)
    (leftEnd, rightEnd)
end

function expectation(distr::ContinuousUnivariateDistribution, g::Function, epsilon::Real)
    f = x->pdf(distr,x)
    (leftEnd, rightEnd) = getEndpoints(distr, epsilon)
    integrate(x -> f(x)*g(x), leftEnd, rightEnd)
end

## Assuming that discrete distributions only take integer values.
function expectation(distr::DiscreteUnivariateDistribution, g::Function, epsilon::Real)
    f = x->pmf(distr,x)
    (leftEnd, rightEnd) = getEndpoints(distr, epsilon)
    sum(map(x -> f(x)*g(x), leftEnd:rightEnd))
end

function expectation(distr::UnivariateDistribution, g::Function)
    expectation(distr, g, 1e-10)
end

## Leave undefined until we've implemented a numerical integration procedure
# function entropy(distr::UnivariateDistribution)
#     pf = typeof(distr)<:ContinuousDistribution ? pdf : pmf
#     f = x -> pf(distr, x)
#     expectation(distr, x -> -log(f(x)))
# end

function kldivergence(P::UnivariateDistribution, Q::UnivariateDistribution)
    expectation(P, x -> log(pdf(P,x)/pdf(Q,x)))
end


text:
#Common functions and data structures
include("common.jl")
include("orthogonalize.jl")

# Hessenberg problem solver
include("hessenberg.jl")

#Stationary solvers
include("stationary.jl")

#Conjugate gradients
include("cg.jl")

#BiCGStab(l)
include("bicgstabl.jl")

#MINRES
include("minres.jl")

#GMRES
include("gmres.jl")

#IDRS
include("idrs.jl")

# QMR
include("qmr.jl")

#Chebyshev
include("chebyshev.jl")

#Simple Eigensolvers
include("simple_eigensolvers.jl")
include("lobpcg.jl")

#Golub-Kahan-Lanczos singular values computation
include("svdl.jl")

include("lsqr.jl")
include("lsmr.jl")

#History data structure
include("history.jl")

# Test deprecations
include("deprecations.jl")


text:
"""
	EmptyNetwork <: AbstractNetwork

A network with no connections.

# Fields
- N              :: Integer -> Number of nodes
"""
struct EmptyNetwork <: AbstractNetwork
	N :: Integer
	_props :: NetworkProperties
end

"""
	EmptyNetwork(N::Integer)

Create an empty network with `N` nodes
"""
function EmptyNetwork(N::Integer)
	(N <= 0) && throw(ArgumentError("Number of nodes must be a positive integer!"))
	EmptyNetwork(N, NetworkProperties(N, false))
end

show(network::EmptyNetwork) = println("Empty Network\n- N = $(network.N)")


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct BayesianProbitRegressor_Gaussian <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function BayesianProbitRegressor_Gaussian(; kwargs...)
        obj = new(meta(BayesianProbitRegressor_Gaussian), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct BayesianProbitRegressor_Gaussian
const __meta_BayesianProbitRegressor_Gaussian = Ref{ProtoMeta}()
function meta(::Type{BayesianProbitRegressor_Gaussian})
    ProtoBuf.metalock() do
        if !isassigned(__meta_BayesianProbitRegressor_Gaussian)
            __meta_BayesianProbitRegressor_Gaussian[] = target = ProtoMeta(BayesianProbitRegressor_Gaussian)
            allflds = Pair{Symbol,Union{Type,String}}[:mean => Float64, :precision => Float64]
            meta(target, BayesianProbitRegressor_Gaussian, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_BayesianProbitRegressor_Gaussian[]
    end
end
function Base.getproperty(obj::BayesianProbitRegressor_Gaussian, name::Symbol)
    if name === :mean
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :precision
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

mutable struct BayesianProbitRegressor_FeatureValueWeight <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function BayesianProbitRegressor_FeatureValueWeight(; kwargs...)
        obj = new(meta(BayesianProbitRegressor_FeatureValueWeight), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct BayesianProbitRegressor_FeatureValueWeight
const __meta_BayesianProbitRegressor_FeatureValueWeight = Ref{ProtoMeta}()
function meta(::Type{BayesianProbitRegressor_FeatureValueWeight})
    ProtoBuf.metalock() do
        if !isassigned(__meta_BayesianProbitRegressor_FeatureValueWeight)
            __meta_BayesianProbitRegressor_FeatureValueWeight[] = target = ProtoMeta(BayesianProbitRegressor_FeatureValueWeight)
            allflds = Pair{Symbol,Union{Type,String}}[:featureValue => UInt32, :featureWeight => BayesianProbitRegressor_Gaussian]
            meta(target, BayesianProbitRegressor_FeatureValueWeight, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_BayesianProbitRegressor_FeatureValueWeight[]
    end
end
function Base.getproperty(obj::BayesianProbitRegressor_FeatureValueWeight, name::Symbol)
    if name === :featureValue
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :featureWeight
        return (obj.__protobuf_jl_internal_values[name])::BayesianProbitRegressor_Gaussian
    else
        getfield(obj, name)
    end
end

mutable struct BayesianProbitRegressor_FeatureWeight <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function BayesianProbitRegressor_FeatureWeight(; kwargs...)
        obj = new(meta(BayesianProbitRegressor_FeatureWeight), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct BayesianProbitRegressor_FeatureWeight
const __meta_BayesianProbitRegressor_FeatureWeight = Ref{ProtoMeta}()
function meta(::Type{BayesianProbitRegressor_FeatureWeight})
    ProtoBuf.metalock() do
        if !isassigned(__meta_BayesianProbitRegressor_FeatureWeight)
            __meta_BayesianProbitRegressor_FeatureWeight[] = target = ProtoMeta(BayesianProbitRegressor_FeatureWeight)
            allflds = Pair{Symbol,Union{Type,String}}[:featureId => UInt32, :weights => Base.Vector{BayesianProbitRegressor_FeatureValueWeight}]
            meta(target, BayesianProbitRegressor_FeatureWeight, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_BayesianProbitRegressor_FeatureWeight[]
    end
end
function Base.getproperty(obj::BayesianProbitRegressor_FeatureWeight, name::Symbol)
    if name === :featureId
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :weights
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{BayesianProbitRegressor_FeatureValueWeight}
    else
        getfield(obj, name)
    end
end

mutable struct BayesianProbitRegressor <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function BayesianProbitRegressor(; kwargs...)
        obj = new(meta(BayesianProbitRegressor), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct BayesianProbitRegressor
const __meta_BayesianProbitRegressor = Ref{ProtoMeta}()
function meta(::Type{BayesianProbitRegressor})
    ProtoBuf.metalock() do
        if !isassigned(__meta_BayesianProbitRegressor)
            __meta_BayesianProbitRegressor[] = target = ProtoMeta(BayesianProbitRegressor)
            fnum = Int[1,2,3,10,11,12,13,20,21,22,23]
            allflds = Pair{Symbol,Union{Type,String}}[:numberOfFeatures => UInt32, :bias => BayesianProbitRegressor_Gaussian, :features => Base.Vector{BayesianProbitRegressor_FeatureWeight}, :regressionInputFeatureName => AbstractString, :optimismInputFeatureName => AbstractString, :samplingScaleInputFeatureName => AbstractString, :samplingTruncationInputFeatureName => AbstractString, :meanOutputFeatureName => AbstractString, :varianceOutputFeatureName => AbstractString, :pessimisticProbabilityOutputFeatureName => AbstractString, :sampledProbabilityOutputFeatureName => AbstractString]
            meta(target, BayesianProbitRegressor, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_BayesianProbitRegressor[]
    end
end
function Base.getproperty(obj::BayesianProbitRegressor, name::Symbol)
    if name === :numberOfFeatures
        return (obj.__protobuf_jl_internal_values[name])::UInt32
    elseif name === :bias
        return (obj.__protobuf_jl_internal_values[name])::BayesianProbitRegressor_Gaussian
    elseif name === :features
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{BayesianProbitRegressor_FeatureWeight}
    elseif name === :regressionInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :optimismInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :samplingScaleInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :samplingTruncationInputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :meanOutputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :varianceOutputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :pessimisticProbabilityOutputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :sampledProbabilityOutputFeatureName
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

export BayesianProbitRegressor_Gaussian, BayesianProbitRegressor_FeatureValueWeight, BayesianProbitRegressor_FeatureWeight, BayesianProbitRegressor


text:
# Load Julia packages (libraries) needed for clip

using StatisticalRethinking
using CSV, DataFrames
using StanSample, MonteCarloMeasurements

ProjDir = @__DIR__

# ### snippet 5.1

df = DataFrame!(CSV.File(rel_path("..", "data", "WaffleDivorce.csv"), delim=';'));
scale!(df, [:Marriage, :MedianAgeMarriage, :Divorce])

# Define the Stan language model

m_5_2 = "
data {
  int N;
  vector[N] divorce_s;
  vector[N] marriage_s;
}
parameters {
  real a;
  real bM;
  real<lower=0> sigma;
}
model {
  vector[N] mu = a + bM * marriage_s;
  a ~ normal( 0 , 0.2 );
  bM ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  divorce_s ~ normal( mu , sigma );
}
";

# Define the SampleModel
m5_2s = SampleModel("m5.2", m_5_2);

# Input data

m5_2_data = Dict(
  "N" => size(df, 1),
  "divorce_s" => df[:, :Divorce_s],
  "marriage_s" => df[:, :Marriage_s]
);

# Sample using cmdstan

rc = stan_sample(m5_2s, data=m5_2_data);

if success(rc)

  # Describe the draws

  dfa2 = read_samples(m5_2s; output_format=:dataframe)

  # Result rethinking

  rethinking = "
          mean   sd  5.5% 94.5%
    a     0.00 0.11 -0.17  0.17
    bM    0.35 0.13  0.15  0.55
    sigma 0.91 0.09  0.77  1.05
  "

  Particles(dfa2)

end


text:
module AbstractPatterns
export spec_gen, RedyFlavoured, TypeObject
export and, or, literal, and, wildcard, decons,
       guard, effect
export PatternCompilationError, Target, PatternImpl, PComp
export APP, NoncachablePre, NoPre
export ChainDict, for_chaindict, child, for_chaindict_dup
export BasicPatterns
export P_bind, P_tuple, P_type_of, P_vector, P_capture, P_vector3, P_slow_view, P_fast_view
export P_svec, P_svec3
export SimpleCachablePre, see_captured_vars, see_captured_vars!
export CFGSpec, CFGJump, CFGLabel, CFGItem, init_cfg

mutable struct CFGSpec
    exp :: Expr
end

struct CFGItem
    kind :: Symbol
    name :: Symbol
end

CFGJump(x::Symbol) = CFGItem(Symbol("@goto"), x)
CFGLabel(x::Symbol) = CFGItem(Symbol("@label"), x)

init_cfg(ex::Expr) = init_cfg(CFGSpec(ex))
function init_cfg(cfg::CFGSpec)
    exp = copy(cfg.exp)
    cfg_info = Dict{Symbol, Symbol}()
    init_cfg!(exp, cfg_info)
    exp
end

const _const_lineno = LineNumberNode(32, "<codegen>")
function init_cfg!(ex::Expr, cf_info::Dict{Symbol, Symbol})
    args = ex.args
    for i in eachindex(args)
        @inbounds arg = args[i]
        if arg isa CFGItem
            label = get!(cf_info, arg.name) do
                gensym(arg.name)
            end
            @inbounds args[i] = Expr(:macrocall, arg.kind, _const_lineno, label)
        elseif arg isa Expr
            init_cfg!(arg, cf_info)
        elseif arg isa CFGSpec
            @inbounds args[i] = init_cfg(arg)
        end
    end
end

include("DataStructure.jl")
include("Target.jl")
include("PatternSignature.jl")
# include("Print.jl")            # maybe use one day
# include("structures/Print.jl") # maybe use one day
include("structures/TypeTagExtraction.jl")
include("ADT.jl")
include("CaseMerge.jl")
include("UserSignature.jl")
include("Retagless.jl")
include("impl/RedyFlavoured.jl")
include("impl/BasicPatterns.jl")
using .BasicPatterns

const _points_of_view = Dict{Function, Int}(tag_extract => 1, untagless => 2)
function spec_gen(branches :: Vector{Pair{Function, Tuple{LineNumberNode, Int}}})
    cores = Branch[]
    for (tf, ln_and_cont) in branches
        impls = (tag_extract(_points_of_view), untagless(_points_of_view))
        type, pat = tf(impls)
        push!(cores, PatternInfo(pat::TagfulPattern, type::TypeObject) => ln_and_cont)
    end
    split_cores = Branch[]
    case_split!(split_cores, cores)
    case_merge(split_cores)
end

end # module


text:
"""
    histogram(data; kwargs...)

Description
============

Draws a horizontal histogram of the given `data`. The positional
parameter `data` can either be a `StatsBase.Histogram`, or some
`AbstractArray`. In the later case the `Histogram` will be
fitted automatically.

Note internally that `histogram` is a simply wrapper for
[`barplot`](@ref), which means that it supports the same keyword
arguments.

Usage
======

    histogram(x; nbins, closed = :left, kwargs...)

    histogram(hist; xscale = identity, title = "", xlabel = "", ylabel = "", labels = true, border = :barplot, margin = 3, padding = 1, color = :green, width = 40, symb = "▇")

Arguments
==========

- **`x`** : Array of numbers for which the histogram should be computed.

- **`nbins`** : The approximate number of bins that should be used.

- **`closed`** : If `:left` (the default), the bin intervals are
  left-closed ``[a,b)``; if `:right`, intervals are right-closed
  ``(a,b]``.

- **`hist`** : A fitted `StatsBase.Histogram` that should be plotted.

- **`xscale`** : Function to transform the bar length before plotting.
  This effectively scales the x-axis without influencing the captions
  of the individual bars. e.g. use `xscale = log10` for logscale.

$DOC_PLOT_PARAMS

- **`symb`** : Specifies the character that should be used to
  render the individual bars.

Returns
========

A plot object of type `Plot{BarplotGraphics}`

Author(s)
==========

- Iain Dunning (Github: https://github.com/IainNZ)
- Christof Stocker (Github: https://github.com/Evizero)
- Kenta Sato (Github: https://github.com/bicycle1885)

Examples
=========

```julia-repl
julia> histogram(randn(1000) * 0.1, closed = :right, nbins = 15)
                  ┌                                        ┐
   (-0.3 , -0.25] ┤▇ 4
   (-0.25, -0.2 ] ┤▇▇ 12
   (-0.2 , -0.15] ┤▇▇▇▇▇▇▇▇▇ 48
   (-0.15, -0.1 ] ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 114
   (-0.1 , -0.05] ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 143
   (-0.05, -0.0 ] ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 177
   ( 0.0 ,  0.05] ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 192
   ( 0.05,  0.1 ] ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 152
   ( 0.1 ,  0.15] ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 83
   ( 0.15,  0.2 ] ┤▇▇▇▇▇▇▇▇ 45
   ( 0.2 ,  0.25] ┤▇▇▇▇ 23
   ( 0.25,  0.3 ] ┤▇ 6
   ( 0.3 ,  0.35] ┤ 1
                  └                                        ┘
                                  Frequency
```

See also
=========

[`Plot`](@ref), [`barplot`](@ref), [`BarplotGraphics`](@ref)
"""
function histogram(
        hist::Histogram;
        symb = "▇",
        xscale = identity,
        xlabel = transform_name(xscale, "Frequency"),
        kw...)
    edges, counts = hist.edges[1], hist.weights
    labels = Vector{String}(undef, length(counts))
    binwidth = typeof(edges.step) == Float64 ? edges.step : edges.step.hi
    # compute label padding based on all labels.
    # this is done to make all decimal points align.
    pad_left, pad_right = 0, 0
    for I in eachindex(edges)
        val1 = float_round_log10(edges[I], binwidth)
        val2 = float_round_log10(val1+binwidth, binwidth)
        a1 = Base.alignment(IOBuffer(), val1)
        a2 = Base.alignment(IOBuffer(), val2)
        pad_left = max(pad_left, a1[1], a2[1])
        pad_right = max(pad_right, a1[2], a2[2])
    end
    # compute the labels using the computed padding
    l_str = hist.closed == :right ? "(" : "["
    r_str = hist.closed == :right ? "]" : ")"
    for i in 1:length(counts)
        val1 = float_round_log10(edges[i], binwidth)
        val2 = float_round_log10(val1+binwidth, binwidth)
        a1 = Base.alignment(IOBuffer(), val1)
        a2 = Base.alignment(IOBuffer(), val2)
        labels[i] =
            "\e[90m" * l_str * "\e[0m" *
            repeat(" ", pad_left - a1[1]) *
            string(val1) *
            repeat(" ", pad_right - a1[2]) *
            "\e[90m, \e[0m" *
            repeat(" ", pad_left - a2[1]) *
            string(val2) *
            repeat(" ", pad_right - a2[2]) *
            "\e[90m" * r_str * "\e[0m"
    end
    barplot(labels, counts; symb = symb, xlabel = xlabel, xscale = xscale, kw...)
end

function histogram(x; bins = nothing, closed = :left, kw...)
    if bins != nothing
        Base.depwarn("The keyword parameter `bins` is deprecated, use `nbins` instead", :histogram)
        hist = fit(Histogram, x; nbins = bins, closed = closed)
    else
        hargs = filter(p -> p.first == :nbins, kw)
        hist = fit(Histogram, x; closed = closed, hargs...)
    end
    pargs = filter(p -> p.first != :nbins, kw)
    histogram(hist; pargs...)
end

@deprecate histogram(x::AbstractArray, nbins::Int; kw...) histogram(x; nbins = nbins, kw...)


text:
# Practical Machine learning
# Artificial Neural Network
# Chapter 11

function randInitializeWeights(L_in, L_out)
   W = RANDINITIALIZEWEIGHTS(L_in, L_out)
  return W
end


text:
# Script for formatting the raw data

using DataFrames, StatFiles

# ------------- loading data -----------------------------------------------------------------------------------------------------------------------

df = DataFrame(load("raw/wiot_full.dta")) # entire 2013 version of WIOD (1995-2011) in long format (Stata .dta file)
# can be obtained from http://www.wiod.org/database/wiots13 under "WIOT tables in STATA format"
df2 = copy(df) # make a copy in case code does not work so we dont need to load entire dataset every time!
#df = df[end-3_000_000:end, :] # work with smaller df for now

rename!(df, :row_item => :row_sector, :col_item => :col_sector)
transform!(df, [:row_country, :col_country] .=> ByRow(string) .=> [:row_country, :col_country], renamecols=false) # convert columns to strings

transform!(df, [:year, :row_sector, :col_sector] .=> ByRow(Int64) .=> [:year, :row_sector, :col_sector], renamecols=false) # convert columns to Float64
sort!(df, names(df[:, Not([:col_country, :value])])) # now sorted alphabetically and numerically except for row_country where the additional variables
# TOT, GO, VA are at the bottom of each year

# -------------- data wrangling -----------------------------------------------------------------------------------------------------------------------

years = 1995:2011
ctrys = sort(unique(df.col_country)) # all country iso3 codes
goods_sectors = 1:16 # numbers associated with goods sectors
service_sectors = 17:35 # numbers associated with service sectors

n_ctrys = length(ctrys) # 40 + ROW = 41 countries in WIOD Release 2013
n_sectors = length([goods_sectors; service_sectors]) # 35 different sectors in WIOD Release 2013

# The function "make matrix" accepts is designed to use the input data from WIOD in stata format (above), it returns the specified year's
# 1. ID ... intermediate demand, N*S×N*S
# 2. FD ... final demand, N*S×N
# 3. GO ... total final output, N*S×1
# 4. VA ... value added, N*S×1
# 5. IV ... net inventory changes (notice negative values, i.e. reduction in inventories), N*S×N
# 6. other ... total other expenses, N*S×1
# Notes: use "sort!" exessively in order to make sure the matrices are sorted accordingly to country-sector pairs
# Notes: ID/FD_corrected which are corrected at the country level (should correct at country-industry level instead)

function obtain_matrices(df::DataFrame, year::Int64)

    # subset df at start to make difference in subsets below clearer
    df1 = subset(df, :year => ByRow(x -> x == year), :row_country => ByRow(x -> x in ctrys), :col_country => ByRow(x -> x in ctrys))

    # ---------------- Intermediate demand matrix, N*S×N*S
    ID_long = subset(df1, :row_sector => ByRow(x -> x in 1:35), :col_sector => ByRow(x -> x in 1:35)) # subsetting for correct sectors

    sort!(ID_long, [:row_country, :row_sector, :col_country, :col_sector]) # sort according to country and sector
    ID_long[:, :id_col] .= ID_long.col_country .* "___" .* string.(ID_long.col_sector) # identifier column for wide format

    ID_wide = unstack(ID_long, [:row_country, :row_sector], :id_col, :value) # transform into wide format
    ID = Matrix(convert.(Float64, ID_wide[:, Not([:row_country, :row_sector])])) # N*S×N*S, matrix object to perform matrix algebra

    # ---------------- Final demand matrix, N*S×N
    FD_long = subset(df1, :row_sector => ByRow(x -> x in 1:35), :col_sector => ByRow(x -> x in 37:41)) # subsetting for correct sectors

    gdf = groupby(FD_long, [:row_country, :row_sector, :col_country]) # sum the four sources of final demand into one
    FD_long = combine(gdf, :value => sum => :total)

    sort!(FD_long, [:row_country, :col_country, :row_sector]) # sort according to country and sector
    FD_wide = unstack(FD_long, [:row_country, :row_sector], :col_country, :total) # transform into wide format
    FD = Matrix(convert.(Float64, FD_wide[:, Not([:row_country, :row_sector])])) # N*S×N, matrix object to perform matrix algebra

    # ---------------- Inventory demand matrix, N*S×N
    IV_long = subset(df1, :row_sector => ByRow(x -> x in 1:35), :col_sector => ByRow(x -> x == 42)) # subsetting for correct sectors

    sort!(IV_long, [:row_country, :row_sector, :col_country]) # sort according to country and sector
    gdf = groupby(IV_long, [:row_country, :row_sector, :col_country])
    IV_long = combine(gdf, :value => sum => :total)
    IV_wide = unstack(IV_long, [:row_country, :row_sector], :col_country, :total) # transform into wide format

    IV = Matrix(convert.(Float64, IV_wide[:, Not([:row_country, :row_sector])])) # N*S×N, matrix object to perform matrix algebra

    # ---------------- Gross output, N*S×1
    GO_long = subset(df, :year => ByRow(x -> x == year), :row_country => ByRow(x -> x == "GO"), :col_country => ByRow(x -> x in ctrys),
    :col_sector => ByRow(x -> x in 1:35)) # subsetting for GO (gross output), countries and sectors (notice using df not df1 anymore)

    sort!(GO_long, [:col_country, :col_sector])
    GO = GO_long.value

    # ---------------- Value added matrix, N*S×1
    VA_long = subset(df, :year => ByRow(x -> x == year), :row_country => ByRow(x -> x == "VA"), :col_country => ByRow(x -> x in ctrys),
    :col_sector => ByRow(x -> x in 1:35)) # subsetting for VA, countries and sectors (notice using df not df1 anymore)

    sort!(VA_long, [:col_country, :col_sector])
    VA = VA_long.value

    # ---------------- Other records matrix, N*S×1
    other = ["CIF", "ITM", "PUA", "PUF", "TXP"] # subsetting for [CIF, ITM, PUA, PUF, TXP], countries and sectors
    other_long = subset(df, :year => ByRow(x -> x == year), :row_country => ByRow(x -> x in other), :col_country => ByRow(x -> x in ctrys),
    :col_sector => ByRow(x -> x in 1:35)) # (notice using df not df1 anymore)

    gdf = groupby(other_long, [:col_country, :col_sector]) # sum the four sources of final demand into one
    other_long = combine(gdf, :value => sum => :total)

    sort!(other_long, [:col_country, :col_sector])
    other = other_long.total

    return ID, FD, GO, VA, IV, other

end

#ID, FD, GO, VA, IV, other = obtain_matrices(df, 2011)


text:
path = string(cd(pwd, ".."), "/images/Large_rep1")
dcms = dcmdir_parse(path)
dcm_array = load_dcm_array(dcms)
header = dcms[1].meta
PixelSpacing = get_pixel_size(header)

@testset ExtendedTestSet "mask_heart" begin
    global masked_array
    global center_insert
    global masked_array
    masked_array, center_insert, mask = mask_heart(
        header, dcm_array, size(dcm_array, 3) ÷ 2
    )
    @testset ExtendedTestSet "mask_heart" begin
        global rows
        global cols
        rows, cols = Int(header[(0x0028, 0x0010)]), Int(header[(0x0028, 0x0011)])
        Y, X = collect(1:512), collect(1:512)'
        dist_from_center = @. sqrt((X - center_insert[2])^2 + (Y - center_insert[1])^2)
        answer = dist_from_center .<= 97.3360655737705
        @test mask == answer
    end
end

@testset ExtendedTestSet "get_calcium_slices" begin
    global slice_dict
    global large_index
    slice_dict, large_index = get_calcium_slices(masked_array, header)
    @testset ExtendedTestSet "get_calcium_slices" begin
        answer_dict = Dict(25 => 2, 26 => 2, 24 => 1)
        @test slice_dict == answer_dict
    end

    @testset ExtendedTestSet "get_calcium_slices" begin
        answer_index = [10, 11, 12, 13, 14, 15, 16]
        @test large_index == answer_index
    end
end

@testset ExtendedTestSet "get_calcium_center_slices" begin
    global slice_dict2
    global flipped
    global flipped_index
    slice_dict2, flipped, flipped_index = get_calcium_center_slices(
        masked_array, slice_dict, large_index
    )
    @testset ExtendedTestSet "get_calcium_center_slices" begin
        answer_dict = Dict(25 => 2, 26 => 2, 24 => 1)
        @test slice_dict2 == answer_dict
    end
end

@testset ExtendedTestSet "mask_rod" begin
    global calcium_image
    global slice_CCI
    global quality_slice
    global cal_rod_slice
    calcium_image, slice_CCI, quality_slice, cal_rod_slice = mask_rod(masked_array, header)
    @testset ExtendedTestSet "mask_rod" begin
        @test slice_CCI == 25
    end

    @testset ExtendedTestSet "mask_rod" begin
        @test quality_slice == 32
    end

    @testset ExtendedTestSet "mask_rod" begin
        @test cal_rod_slice == 15
    end
end

@testset ExtendedTestSet "calc_output" begin
    global output
    output = calc_output(masked_array, header, slice_CCI)
    @testset ExtendedTestSet "calc_output" begin
        @test output[1] == 7
    end
end

@testset ExtendedTestSet "center_points" begin
    global center
    global center1
    global center2
    global center3
    center, center1, center2, center3 = center_points(
        dcm_array, output, header, center_insert, slice_CCI
    )
    @testset ExtendedTestSet "center_points" begin
        @test center == [217, 257]
    end
end

@testset ExtendedTestSet "calc_centers" begin
    dict = calc_centers(dcm_array, output, header, center_insert, slice_CCI)
    @testset ExtendedTestSet "calc_centers" begin
        @test dict[:Large_MD] == [151, 292]
    end
end

@testset ExtendedTestSet "mask_inserts" begin
    @testset ExtendedTestSet "mask_inserts" begin
        mask_L_HD, mask_M_HD, mask_S_HD, mask_L_MD, mask_M_MD, mask_S_MD, mask_L_LD, mask_M_LD, mask_S_LD = mask_inserts(
            dcm_array, masked_array, header, slice_CCI, center_insert
        )
        dict = calc_centers(dcm_array, output, header, center_insert, slice_CCI)
        mask_L_HD_answer = create_circular_mask(
            cols, rows, dict[:Large_HD], (round(5 / PixelSpacing[1], RoundUp) / 2) + 1
        )
        @test mask_L_HD == transpose(mask_L_HD_answer)
    end
end