text:
using NCDatasets, ArchGDAL, GeoData, Test, Statistics, Dates, CFTime, Plots, GeoFormatTypes
using GeoData: name, window, mode, span, sampling, val, Ordered
include(joinpath(dirname(pathof(GeoData)), "../test/test_utils.jl"))

ncexamples = "https://www.unidata.ucar.edu/software/netcdf/examples/"
ncsingle = maybedownload(joinpath(ncexamples, "tos_O1_2001-2002.nc"))
ncmulti = maybedownload(joinpath(ncexamples, "test_echam_spectral.nc"))

stackkeys = (
    :abso4, :aclcac, :aclcov, :ahfcon, :ahfice, :ahfl, :ahfliac, :ahfllac,
    :ahflwac, :ahfres, :ahfs, :ahfsiac, :ahfslac, :ahfswac, :albedo, :albedo_nir,
    :albedo_nir_dif, :albedo_nir_dir, :albedo_vis, :albedo_vis_dif, :albedo_vis_dir,
    :alsobs, :alsoi, :alsol, :alsom, :alsow, :ameltdepth, :ameltfrac, :amlcorac,
    :ao3, :apmeb, :apmegl, :aprc, :aprl, :aprs, :aps, :az0i, :az0l, :az0w,
    :barefrac, :dew2, :drain, :evap, :evapiac, :evaplac, :evapwac, :fage, :friac,
    :geosp, :glac, :gld, :hyai, :hyam, :hybi, :hybm, :lsp, :q, :qres, :qvi, :relhum,
    :runoff, :sd, :seaice, :siced, :sicepdi, :sicepdw, :sicepres, :slm, :sn, :snacl,
    :snc, :sni, :snifrac, :snmel, :sofliac, :sofllac, :soflwac, :srad0, :srad0d,
    :srad0u, :sradl, :srads, :sradsu, :sraf0, :srafl, :srafs, :st, :svo, :t2max,
    :t2min, :temp2, :thvsig, :topmax, :tpot, :trad0, :tradl, :trads, :tradsu,
    :traf0, :trafl, :trafs, :trfliac, :trfllac, :trflwac, :tropo, :tsi, :tsicepdi,
    :tslm1, :tsurf, :tsw, :u10, :ustr, :ustri, :ustrl, :ustrw, :v10, :vdis, :vdisgw,
    :vstr, :vstri, :vstrl, :vstrw, :wimax, :wind10, :wl, :ws, :wsmx, :xi, :xivi,
    :xl, :xlvi
)

@testset "NCDarray" begin
    ncarray = NCDarray(ncsingle)

    @testset "open" begin
        @test all(Open(A -> A[Lat=1], ncarray) .=== ncarray[:, 1, :])
    end

    @testset "array properties" begin
        @test size(ncarray) == (180, 170, 24)
        @test ncarray isa NCDarray
        @test val(dims(ncarray, Ti())) == DateTime360Day(2001, 1, 16):Month(1):DateTime360Day(2002, 12, 16)
        @test val(dims(ncarray, Lat())) == -79.5:89.5
        @test val(dims(ncarray, Lon())) == 1.0:2:359.0
        @test bounds(ncarray) == (
            (0.0, 360.0),
            (-80.0, 90.0),
            (DateTime360Day(2001, 1, 16), DateTime360Day(2002, 12, 16)),
        )
    end

    @testset "dimensions" begin
        @test ndims(ncarray) == 3
        @test length.(dims(ncarray)) == (180, 170, 24)
        @test dims(ncarray) isa Tuple{<:Lon,<:Lat,<:Ti}
        @test refdims(ncarray) == ()
        # TODO detect the time span, and make it Regular
        @test mode(ncarray) ==
            (Mapped(Ordered(), Regular(2.0), Intervals(Center()), EPSG(4326), EPSG(4326)),
             Mapped(Ordered(), Regular(1.0), Intervals(Center()), EPSG(4326), EPSG(4326)),
             Sampled(Ordered(), Irregular(), Points()))
        @test bounds(ncarray) == ((0.0, 360.0), (-80.0, 90.0), (DateTime360Day(2001, 1, 16), DateTime360Day(2002, 12, 16)))
    end

    @testset "other fields" begin
        @test ismissing(missingval(ncarray))
        @test metadata(ncarray) isa NCDarrayMetadata
        @test name(ncarray) == :tos
    end

    @testset "indexing" begin
        @test ncarray[Ti(1)] isa GeoArray{<:Any,2}
        @test ncarray[Lat(1), Ti(1)] isa GeoArray{<:Any,1}
        @test ncarray[Lon(1), Ti(1)] isa GeoArray{<:Any,1}
        @test ncarray[Lon(1), Lat(1), Ti(1)] isa Missing
        @test ncarray[Lon(30), Lat(30), Ti(1)] isa Float32
        # Russia
        @test ncarray[Lon(50), Lat(100), Ti(1)] isa Missing
        # Alaska
        @test ncarray[Lat(Near(64.2008)), Lon(Near(149.4937)), Ti(1)] isa Missing
        @test ncarray[Ti(2), Lon(At(59.0)), Lat(At(-50.5))] == ncarray[30, 30, 2] === 278.47168f0
    end

    @testset "indexing with reverse lat" begin
        if !haskey(ENV, "CI") # CI downloads fail. But run locally
            ncrevlat = maybedownload("ftp://ftp.cdc.noaa.gov/Datasets/noaa.ersst.v5/sst.mon.ltm.1981-2010.nc")
            ncrevlatarray = NCDstack(ncrevlat; childkwargs=(missingval=-9.96921f36,))[:sst]
            @test order(dims(ncrevlatarray, Lat)) == Ordered(ReverseIndex(), ReverseArray(), ForwardRelation())
            @test ncrevlatarray[Lat(At(40)), Lon(At(100)), Ti(1)] == missingval(ncrevlatarray)
            @test ncrevlatarray[Lat(At(-40)), Lon(At(100)), Ti(1)] == ncrevlatarray[51, 65, 1] == 14.5916605f0
            @test val(span(ncrevlatarray, Ti)) == Month(1)
        end
    end

    @testset "selectors" begin
        a = ncarray[Lon(At(21.0)), Lat(Between(50, 52)), Ti(Near(DateTime360Day(2002, 12)))];
        @test bounds(a) == ((50.0, 52.0),)
        x = ncarray[Lon(Near(150)), Lat(Near(30)), Ti(1)]
        @test x isa Float32
        # TODO make sure we are getting the right cell.
        @test size(ncarray[Lat(Between(-80, 90)), Lon(Between(0, 360)),
            Ti(Between(DateTime360Day(2001, 1, 16), DateTime360Day(2003, 01, 16)))
        ]) == (180, 170, 24)
        nca = ncarray[Lat(Between(-80, -25)), Lon(Between(0, 180)),
                      Ti(Near(DateTime360Day(2002, 02, 20)))]
        @test size(nca) == (90, 55)
    end

    @testset "conversion to GeoArray" begin
        geoarray = ncarray[Lon(1:50), Lat(20:20), Ti(1)]
        @test size(geoarray) == (50, 1)
        @test eltype(geoarray) <: Union{Missing,Float32}
        @time geoarray isa GeoArray{Float32,1}
        @test dims(geoarray) isa Tuple{<:Lon,<:Lat}
        @test refdims(geoarray) isa Tuple{<:Ti}
        @test metadata(geoarray) == metadata(ncarray)
        @test ismissing(missingval(geoarray))
        @test name(geoarray) == :tos
    end

    @testset "save" begin
        @testset "to netcdf" begin
            # TODO save and load subset
            geoarray = GeoArray(ncarray)
            metadata(geoarray)
            @test size(geoarray) == size(ncarray)
            filename = tempname()
            write(filename, NCDarray, geoarray)
            saved = GeoArray(NCDarray(filename))
            @test size(saved) == size(geoarray)
            @test refdims(saved) == refdims(geoarray)
            @test missingval(saved) === missingval(geoarray)
            @test_broken metadata(saved) == metadata(geoarray)
            @test_broken all(metadata.(dims(saved)) .== metadata.(dims(geoarray)))
            @test GeoData.name(saved) == GeoData.name(geoarray)
            @test all(mode.(dims(saved)) .== mode.(dims(geoarray)))
            @test typeof(dims(saved)) <: typeof(dims(geoarray))
            @test val(dims(saved)[3]) == val(dims(geoarray)[3])
            @test all(val.(dims(saved)) .== val.(dims(geoarray)))
            @test all(data(saved) .=== data(geoarray))
            @test saved isa typeof(geoarray)
            # TODO test crs
        end
        @testset "to gdal" begin
            gdalfilename = tempname() * ".tif"
            nccleaned = replace_missing(ncarray[Ti(1)], -9999.0)
            write(gdalfilename, GDALarray, nccleaned)
            gdalarray = GDALarray(gdalfilename)
            # gdalarray WKT is missing one AUTHORITY
            # @test_broken crs(gdalarray) == convert(WellKnownText, EPSG(4326))
            # But the Proj representation is the same
            @test convert(ProjString, crs(gdalarray)) == convert(ProjString, EPSG(4326))
            @test bounds(gdalarray) == (bounds(nccleaned)..., (1, 1))
            # Tiff locus = Start, Netcdf locus = Center
            @test reverse(val(dims(gdalarray, Lat))) .+ 0.5 ≈ val(dims(nccleaned, Lat))
            @test val(dims(gdalarray, Lon)) .+ 1.0  ≈ val(dims(nccleaned, Lon))
            @test reverse(GeoArray(gdalarray); dims=Lat()) ≈ nccleaned
        end
        @testset "to grd" begin
            nccleaned = replace_missing(ncarray[Ti(1)], -9999.0)
            write("testgrd", GRDarray, nccleaned)
            grdarray = GRDarray("testgrd");
            @test crs(grdarray) == convert(ProjString, EPSG(4326))
            @test bounds(grdarray) == (bounds(nccleaned)..., (1, 1))
            @test val(dims(grdarray, Lat)) ≈ val(dims(nccleaned, Lat)) .- 0.5
            @test val(dims(grdarray, Lon)) ≈ val(dims(nccleaned, Lon)) .- 1.0
            @test GeoArray(grdarray) ≈ reverse(nccleaned; dims=Lat)
        end
    end

    @testset "show" begin
        sh = sprint(show, ncarray)
        # Test but don't lock this down too much
        @test occursin("NCDarray", sh)
        @test occursin("Latitude", sh)
        @test occursin("Longitude", sh)
        @test occursin("Time", sh)
    end

    @testset "plot" begin
        ncarray |> plot
        ncarray[Ti(1)] |> plot
        ncarray[Lat(100), Ti(1)] |> plot
    end

end

@testset "NCDstack" begin
    ncstack = NCDstack(ncmulti)

    @testset "load ncstack" begin
        @test ncstack isa NCDstack{String}
        @test ismissing(missingval(ncstack))
        @test metadata(ncstack) isa NCDstackMetadata
        @test refdims(ncstack) == ()
        # Loads child as a regular GeoArray
        @test_throws NCDatasets.NetCDFError ncstack[:not_a_key]
        @test ncstack[:albedo] isa GeoArray{<:Any,3}
        @test ncstack[:albedo, 2, 3, 1] isa Float32
        @test ncstack[:albedo, :, 3, 1] isa GeoArray{<:Any,1}
        @test dims(ncstack, :albedo) isa Tuple{<:Lon,<:Lat,<:Ti}
        @test keys(ncstack) isa NTuple{131,Symbol}
        @test keys(ncstack) == stackkeys
        @test first(keys(ncstack)) == :abso4
        @test metadata(ncstack) isa NCDstackMetadata
        @test metadata(ncstack)["institution"] == "Max-Planck-Institute for Meteorology"
        @test metadata(ncstack, :albedo) isa NCDarrayMetadata
        @test metadata(ncstack, :albedo)["long_name"] == "surface albedo"
        # Test some DimensionalData.jl tools work
        # Time dim should be reduced to length 1 by mean
        @test axes(mean(ncstack[:albedo, Lat(1:20)] , dims=Ti)) ==
              (Base.OneTo(192), Base.OneTo(20), Base.OneTo(1))
        geoarray = ncstack[:albedo][Ti(4:6), Lon(1), Lat(2)]
        @test geoarray == ncstack[:albedo, Ti(4:6), Lon(1), Lat(2)]
        @test size(geoarray) == (3,)
    end

    if VERSION > v"1.1-"
        @testset "copy" begin
            geoarray = ncstack[:albedo]
            copy!(geoarray, ncstack, :albedo);
            # First wrap with GeoArray() here or == loads from disk for each cell.
            # we need a general way of avoiding this in all disk-based sources
            @test geoarray == GeoArray(ncstack[:albedo])
        end
    end

    @testset "indexing" begin
        ncmultistack = NCDstack(ncsingle)
        @test dims(ncmultistack, :tos) isa Tuple{<:Lon,<:Lat,<:Ti}
        @test ncmultistack[:tos] isa GeoArray{<:Any,3}
        @test ncmultistack[:tos, Ti(1)] isa GeoArray{<:Any,2}
        @test ncmultistack[:tos, Lat(1), Ti(1)] isa GeoArray{<:Any,1}
        @test ncmultistack[:tos, 8, 30, 10] isa Float32
    end

    @testset "window" begin
        windowedstack = NCDstack(ncmulti; window=(Lat(1:5), Lon(1:5), Ti(1)))
        @test window(windowedstack) == (Lat(1:5), Lon(1:5), Ti(1))
        windowedarray = windowedstack[:albedo]
        @test size(windowedarray) == (5, 5)
        @test windowedarray[1:3, 2:2] == reshape([0.84936917f0, 0.8776228f0, 0.87498736f0], 3, 1)
        @test windowedarray[1:3, 2] == [0.84936917f0, 0.8776228f0, 0.87498736f0]
        @test windowedarray[1, 2] == 0.84936917f0
        windowedstack = NCDstack(ncmulti; window=(Lat(1:5), Lon(1:5), Ti(1:1)))
        windowedarray = windowedstack[:albedo]
        @test windowedarray[1:3, 2:2, 1:1] == reshape([0.84936917f0, 0.8776228f0, 0.87498736f0], 3, 1, 1)
        @test windowedarray[1:3, 2:2, 1] == reshape([0.84936917f0, 0.8776228f0, 0.87498736f0], 3, 1)
        @test windowedarray[1:3, 2, 1] == [0.84936917f0, 0.8776228f0, 0.87498736f0]
        @test windowedarray[1, 2, 1] == 0.84936917f0
        windowedstack = NCDstack(ncmulti; window=(Ti(1),))
        windowedarray = windowedstack[:albedo]
        @test windowedarray[1:3, 2:2] == reshape([0.84936917f0, 0.8776228f0, 0.87498736f0], 3, 1)
        @test windowedarray[1:3, 2] == [0.84936917f0, 0.8776228f0, 0.87498736f0]
        @test windowedarray[1, 2] ==  0.84936917f0
    end

    @testset "conversion to GeoStack" begin
        geostack = GeoStack(ncstack)
        @test Symbol.(Tuple(keys(geostack))) == keys(ncstack)
        smallstack = GeoStack(ncstack; keys=(:albedo, :evap, :runoff))
        @test keys(smallstack) == (:albedo, :evap, :runoff)
    end

    @testset "save" begin
        metadata(ncstack)
        geostack = GeoStack(ncstack);
        metadata(geostack)
        length(dims(geostack[:aclcac]))
        ndims(geostack[:aclcac])
        filename = tempname()
        write(filename, NCDstack, geostack);
        saved = GeoStack(NCDstack(filename))
        @test keys(saved) == keys(geostack)
        @test metadata(saved)["advection"] == "Lin & Rood"
        @test metadata(saved) == metadata(geostack) == metadata(ncstack)
        @test all(first(values(saved)) .== first(values(geostack)))
    end

end

@testset "NCD series" begin
    series = GeoSeries([ncmulti, ncmulti], (Ti,); childtype=NCDstack)
    geoarray = GeoArray(NCDarray(ncmulti, :albedo; name=:test))
    @test series[Ti(1)][:albedo] == geoarray
    @test typeof(series[Ti(1)][:albedo]) == typeof(geoarray)
    modified_series = modify(Array, series)
    @test typeof(modified_series) <: GeoSeries{<:GeoStack{<:NamedTuple{stackkeys,<:Tuple{<:GeoArray{Float32,3,<:Tuple,<:Tuple,<:Array{Float32,3}},Vararg}}}}
end

nothing


text:
using Krylov
using LinearOperators
using MatrixMarket
# using ProfileView

mtx = "data/1138bus.mtx";
# mtx = "data/bcsstk09.mtx";
# mtx = "data/bcsstk18.mtx";

A = MatrixMarket.mmread(mtx);
n = size(A, 1);
b = ones(n); b_norm = norm(b);

# Define a linear operator with preallocation.
Ap = zeros(n);
op = LinearOperator(n, n, true, true, p -> A_mul_B!(1.0,  A, p, 0.0, Ap))

# Solve Ax=b.
(x, stats) = cg(op, b);
# @profile x = cg(A, b);
@time (x, stats) = cg(op, b);
show(stats);
r = b - A * x;
@printf("Relative residual: %8.1e\n", norm(r)/b_norm);

# ProfileView.view()

ext: jl


text:
## legacy, not sure if needed now, JJS 12/5/19
# function TimeSeriesCheckDemand(loads::Array{T}) where {T<:WaterDemand}
#     t = length(loads[1].demand)
#     for l in loads
#         if t == length(l.demand)
#             continue
#         else
#             error("Inconsistent load scaling factor time series length")
#         end
#     end
#     return t
# end

"""
Create a System from an epanet ".inp" file.
"""
function parse_inp_file(inp_file::String)
    data = make_dict(inp_file) # data is a dictionary
    junctions, arcs, res, demands, tanks, patterns, curves, pipes, pumps, valves =
        dict_to_struct(data)

    # need to work on the constructor(s) for the Water System, JJS 200103
    return System()
end


text:
module MultidimensionalTools

using LinearAlgebra

# type wrappers
export AbstractIndex, AbstractIndices, AbstractIndexOrIndices
# helper functions
export n_adjacencies, extrema_indices, tryindex, get_directions
# main functions
export append_n_times, append_n_times_backwards, promote_to_nD, promote_to_3D,
    expand_as_required, adjacencies, expanded_adjacencies, n_adjacent_to,
    expanded_n_adjacent_to, global_adjacencies_indices, global_adjacencies,
    global_n_adjacent_to, reduce_dims

_IterableList{T} = Union{AbstractVector{T}, NTuple{N, T}} where N
AbstractIndex{N} = Union{NTuple{N, T}, CartesianIndex{N}} where {T <: Integer}
AbstractIndices{N}  = Union{AbstractVector{I}, NTuple{M, I}, CartesianIndices{N, NTuple{N, Base.OneTo{T}}}} where {I <: AbstractIndex{N}, T <: Integer, M}
AbstractIndexOrIndices{N} = Union{AbstractIndex{N}, AbstractIndices{N}}

"""
    𝟘(n::Integer) -> NTuple{n, Integer}
    \\bbzero<tab>(n::Integer) -> NTuple{n, Integer}

Simply returns the zero tuple of a given dimension.
"""
𝟘(n::T) where {T <: Integer} = ntuple(_ -> zero(T), n)

_tryindex(M::AbstractArray, inds::CartesianIndex...) where {T, N} = Tuple((checkbounds(Bool, M, i) ? M[i] : nothing) for i in inds)
_tryindex(M::AbstractArray, inds::AbstractArray{CartesianIndex}) where {T, N} = _tryindex(M, inds...)
_tryindex(M::AbstractArray, inds::NTuple...) = _tryindex(M, CartesianIndex.(inds)...)
_tryindex(M::AbstractArray, inds::AbstractArray{NTuple}) where {T, N} = _tryindex(M, inds...)

"""
    tryindex(M::AbstractArray{T, N}, inds::AbstractIndex{N}}...)
    tryindex(M::AbstractArray{T, N}, inds::AbstractIndices{N}})

For each index specified, gets the index or missing if the index is unavailable.
"""
tryindex(M::AbstractArray, inds::AbstractIndex...) = _tryindex(M, inds...)
tryindex(M::AbstractArray, inds::AbstractIndices) = _tryindex(M, inds...)

function _extrema_indices(I::CartesianIndex...)
    Tuple(extrema(reduce(vcat, permutedims(collect(i)) for i in I...), dims = 1))
end
function _extrema_indices(I::AbstractIndexOrIndices{N}...) where {N}
    return Tuple(extrema(reduce(vcat, permutedims(collect(i)) for i in I...), dims = 1))
end

"""
    extrema_indices(I::AbstractIndices{N}) -> NTuple{N, NTuple{N, Integer}}
    extrema_indices(i::AbstractIndex{N}...) -> NTuple{N, NTuple{N, Integer}}
    extrema_indices(A::AbstractArray{T, N}...) -> NTuple{N, NTuple{N, Integer}}

Given indices or arrays, returns an `NTuple`.  Each element in the `NTuple` represents a (min, max) tuple for each dimension.

```julia
julia> T = [(3, 1), (4, 5), (6, 2), (1, 1)]
4-element Array{Tuple{Int64,Int64},1}:
 (3, 1)
 (4, 5)
 (6, 2)
 (1, 1)

julia> extrema_indices(T)
1×2 Array{Tuple{Int64,Int64},2}:
(1, 6)  (1, 5)

julia> # that is, the (min, max) for dim 2 is (1, 5); etc.
```
"""
extrema_indices(I::AbstractIndices{N}) where {N} = _extrema_indices(I)
extrema_indices(i::AbstractIndex{N}...) where {N} = _extrema_indices(i)
extrema_indices(A::AbstractArray{T, N}) where {T, N} = size(A)
extrema_indices(A::AbstractArray{T, N}...) where {T, N} = _extrema_indices(size.(A))

@inline function allequal(A::_IterableList{T}) where T
    length(A) < 2 && return true

    @inbounds for i in 2:length(A)
        first(A) ≠ A[i] && return false
    end

    return true
end
@inline allequal(a::T...) where {T} = allequal(a)

"""
Check whether the given array is a cube (i.e., all its dimensions are the same size).  Despite the function name, this will work for dimension < 3
"""
ishypercube(::Number) = true # this is the case where there is zero dimensions: the point in space
ishypercube(A::AbstractVector{T}) where {T} = length(A) == 1 ? true : false # a vector of length one is a square, hence it is a lower-dimensional cube
ishypercube(A::AbstractArray{T, N}) where {T, N} = allequal(size(A))
function ishypercube_alt(A::AbstractArray{T, N}) where {T, N}
    ax1 = axes(A, 1)
    for d in 2:ndims(A)
        axes(A, d) == ax1 || return false
    end
    return true
end
ishypercube(A::AbstractArray{T, N}...) where {T, N} = Bool[ishypercube(a) for a in A]

"""
Check whether the 2-dimensional aspect of a given array is a square.
"""
is2Dsquare(::Union{Number, AbstractVector{T}}) where {T} = false # there is no 2-dimensional aspect of zero- or one-dimensional objects
is2Dsquare(A::AbstractArray{T, N}) where {T, N} = size(A, 1) == size(A, 2) ? true : false
is2Dsquare(A::AbstractArray{T, N}...) where {T, N} = Bool[is2Dsquare(a) for a in A]

include("generic.jl")
include("adjacent.jl")
# include("multilinalg.jl")

end # end module


text:
# This file was generated, do not modify it. # hide
using MLJ
using PyPlot
using PrettyPrinting
using StableRNGs
import DataFrames: DataFrame, describe
MLJ.color_off() # hide
ioff() # hide

X, y = @load_boston
sch = schema(X)
p = length(sch.names)
describe(y)  # From DataFrames

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# Run modified JSON checker tests

const JSON_DATA_DIR = joinpath(dirname(dirname(@__FILE__)), "data")

for i in 1:38
    file = "fail$(lpad(i, 2, 0)).json"
    filepath = joinpath(JSON_DATA_DIR, "jsonchecker", file)

    @test_throws ErrorException JSON.parsefile(filepath)
end

for i in 1:3
    # Test that the files parse successfully and match streaming parser
    tf = joinpath(JSON_DATA_DIR, "jsonchecker", "pass$(lpad(i, 2, 0)).json")
    @test JSON.parsefile(tf) == open(JSON.parse, tf)
end

# Run JSON roundtrip tests (check consistency of .json)

roundtrip(data) = JSON.json(JSON.Parser.parse(data))

for i in 1:27
    file = "roundtrip$(lpad(i, 2, 0)).json"
    filepath = joinpath(JSON_DATA_DIR, "roundtrip", file)

    rt = roundtrip(readstring(filepath))
    @test rt == roundtrip(rt)
end


text:
using DiffEqJump, DiffEqBase, SafeTestsets

@time begin
  @time @safetestset "Spatial utilities" begin include("utils_test.jl") end
  @time @safetestset "Spatial A + B <--> C" begin include("ABC.jl") end
  @time @safetestset "Pure diffusion" begin include("diffusion.jl") end
end


text:
using IncrementalInference
using Test

##

@testset "test _evalFactorTemporary" begin
## test utility to build a temporary graph

fct = EuclidDistance(Normal(10,1))
varTypes = (ContinuousScalar,ContinuousScalar);
varPts = ([0;],[9.5;])

##

dfg, _dfgfct = IIF._buildGraphByFactorAndTypes!(fct, varTypes, varPts)

@test length(intersect(ls(dfg), [:x1; :x2])) == 2
@test lsf(dfg) == [:x1x2f1;]

## test  the evaluation of factor without

B = IIF._evalFactorTemporary!(EuclidDistance(Normal(10,1)), varTypes, 2, [[10;]], varPts );

@test_broken B isa Vector{Vector{Float64}}
@test isapprox( B[1], [10.0;], atol=1e-6)

##
end

@testset "test residual slack prerequisite for numerical factor gradients, Euclidean(1)" begin
##

fct = EuclidDistance(Normal(10,1))
measurement = [[10;]]
varTypes = (ContinuousScalar,ContinuousScalar)
pts = ([0;],[9.5;])

##

slack_resid = calcFactorResidualTemporary(fct, varTypes, measurement[1], pts)

##

coord_1 = IIF._evalFactorTemporary!(fct, varTypes, 1, measurement, pts, _slack=slack_resid )
@test length(coord_1) == 1
@test isapprox( coord_1[1], [0.0], atol=1e-6)

coord_2 = IIF._evalFactorTemporary!(fct, varTypes, 2, measurement, pts, _slack=slack_resid )
@test length(coord_2) == 1
@test isapprox( coord_2[1], [9.5], atol=1e-6)

##

coord_1 = IIF._evalFactorTemporary!(fct, varTypes, 1, measurement, pts )
@test length(coord_1) == 1
@test isapprox( coord_1[1], [-0.5], atol=1e-6)

coord_2 = IIF._evalFactorTemporary!(fct, varTypes, 2, measurement, pts )
@test length(coord_2) == 1
@test isapprox( coord_2[1], [10.0], atol=1e-6)

##
end

@testset "test residual slack prerequisite for numerical factor gradients, Euclidean(2)" begin
##

fct = LinearRelative(MvNormal([10;0],[1 0; 0 1]))
measurement = [[10.0;0.0]]
varTypes = (ContinuousEuclid{2},ContinuousEuclid{2})
pts = ([0;0.0] ,[9.5;0])

## test the building of factor graph to be correct

_fg,_ = IIF._buildGraphByFactorAndTypes!(fct, varTypes, pts);

@test length(getVal(_fg[:x1])) == 1
@test length(getVal(_fg[:x1])[1]) == 2

@test length(getVal(_fg[:x2])) == 1
@test length(getVal(_fg[:x2])[1]) == 2

##

_fg = initfg()

slack_resid = calcFactorResidualTemporary(fct, varTypes, measurement[1], pts, tfg=_fg)

@test length(getVal(_fg[:x1])) == 1
@test length(getVal(_fg[:x1])[1]) == 2

@test length(getVal(_fg[:x2])) == 1
@test length(getVal(_fg[:x2])[1]) == 2

## Manually provide a common temp graph and force no factor and same variables via keywords

tfg,_ = IIF._buildGraphByFactorAndTypes!(fct, varTypes, pts);
coord_1 = IIF._evalFactorTemporary!(fct, varTypes, 1, measurement, pts, _slack=slack_resid, tfg=tfg, newFactor=false, currNumber=0 )

##

@test length(coord_1) == 1
@test length(coord_1[1]) == 2

@test isapprox( coord_1[1], [0;0.0], atol=1e-6)

coord_2 = IIF._evalFactorTemporary!(fct, varTypes, 2, measurement, pts, _slack=slack_resid )
@test length(coord_2) == 1
@test length(coord_2[1]) == 2

@test isapprox( coord_2[1], [9.5; 0], atol=1e-6)

## Repeat the same test but allow _evalFactorTemporary to self construct internal temporary graph

coord_1 = IIF._evalFactorTemporary!(fct, varTypes, 1, measurement, pts )

##
@test length(coord_1) == 1
@test length(coord_1[1]) == 2

@test isapprox( coord_1[1], [-0.5;0], atol=1e-6)

coord_2 = IIF._evalFactorTemporary!(fct, varTypes, 2, measurement, pts )
@test length(coord_2) == 1
@test length(coord_2[1]) == 2

@test isapprox( coord_2[1], [10.0;0], atol=1e-6)

##
end

#


text:
import PolynomialApproximations: Remez
using SpecialFunctions: erfcx
using PyPlot

T = Float64

function g(y::T) where T <: AbstractFloat
    if -one(T) ≤ y ≤ -one(T) + eps(T)
        return one(T)
    elseif one(T) - eps(T) ≤ y ≤ one(T)
        return zero(T)
    else
        x = (1+y) / (1-y)
        return erfcx(x)
    end
end

function erfcx_approx(x)
    if x < one(T)
        y = (x-1) / (x+1)
    else
        rx = 1 / x
        y = (1-rx) / (1+rx)
    end
    return g_approx(y)
end

n = 15
iterations = 6
clamp = :both
ddp, pt, zmax, zmin = Remez.minimax(T, g, n, iterations, clamp)

g_approx(y) = Remez.Newton_poly(ddp, pt, y)

figure(1)
y = range(-one(T), one(T), length=201)

plot(y, g.(y))
grid(true)
xlabel(L"$y$")
ylabel(L"$g(y)$")

figure(2)
plot(y, g.(y)-g_approx.(y))
grid(true)
xlabel(L"$y$")
ylabel(L"error in $g$")

figure(3)
x = Float64[ 10 * (n/400)^2 for n = 0:400 ]
plot(x, erfcx_approx.(x) - erfcx.(x))
grid(true)
xlabel(L"$x$")
ylabel("error in erfcx")


text:
struct WEsampler{TM, TS, TR, TA} <: AbstractWEsampler
    mutation!::TM
    selection!::TS
    rebin!::TR
    analysis!::TA
end

struct DistributedWEsampler{TM, TS, TR, TA} <: AbstractWEsampler
    mutation::TM
    selection!::TS
    rebin!::TR
    analysis!::TA
end

function WEsampler(mutation!::TM, selection!::TS, rebin!::TR; analysis!::TA = trivial_analysis!) where {TM<:Function, TS<:Function, TR<:Function, TA<:Function}
    return WEsampler(mutation!, selection!, rebin!, analysis!)
end

function DistributedWEsampler(mutation::TM, selection!::TS, rebin!::TR; analysis!::TA = trivial_analysis!) where {TM<:Function, TS<:Function, TR<:Function, TA<:Function}
    return DistributedWEsampler(mutation, selection!, rebin!, analysis!)
end


text:
# Now that you have the structure of your transmission decoded, you can
# calculate the value of the expression it represents.
#
# Literal values (type ID 4) represent a single number as described above. The
# remaining type IDs are more interesting:
#
# - Packets with type ID 0 are sum packets - their value is the sum of the
#   values of their sub-packets. If they only have a single sub-packet, their
#   value is the value of the sub-packet.
# - Packets with type ID 1 are product packets - their value is the result of
#   multiplying together the values of their sub-packets. If they only have a
#   single sub-packet, their value is the value of the sub-packet.
# - Packets with type ID 2 are minimum packets - their value is the minimum of
#   the values of their sub-packets.
# - Packets with type ID 3 are maximum packets - their value is the maximum of
#   the values of their sub-packets.
# - Packets with type ID 5 are greater than packets - their value is 1 if the
#   value of the first sub-packet is greater than the value of the second
#   sub-packet; otherwise, their value is 0. These packets always have exactly
#   two sub-packets.
# - Packets with type ID 6 are less than packets - their value is 1 if the
#   value of the first sub-packet is less than the value of the second
#   sub-packet; otherwise, their value is 0. These packets always have exactly
#   two sub-packets.
# - Packets with type ID 7 are equal to packets - their value is 1 if the value
#   of the first sub-packet is equal to the value of the second sub-packet;
#   otherwise, their value is 0. These packets always have exactly two
#   sub-packets.
#
# Using these rules, you can now work out the value of the outermost packet in
# your BITS transmission.
#
# For example:
#
# C200B40A82 finds the sum of 1 and 2, resulting in the value 3.
# 04005AC33890 finds the product of 6 and 9, resulting in the value 54.
# 880086C3E88112 finds the minimum of 7, 8, and 9, resulting in the value 7.
# CE00C43D881120 finds the maximum of 7, 8, and 9, resulting in the value 9.
# D8005AC2A8F0 produces 1, because 5 is less than 15.
# F600BC2D8F produces 0, because 5 is not greater than 15.
# 9C005AC2F8F0 produces 0, because 5 is not equal to 15.
# 9C0141080250320F1802104A08 produces 1, because 1 + 3 = 2 * 2.
#
# What do you get if you evaluate the expression represented by your
# hexadecimal-encoded BITS transmission?

# read input
line = readline("input.txt")
bits = parse(BigInt, line; base=16)
pos = ndigits(bits; base=2)

# utility function to read n bits
function read_bits(n)
  global pos -= n
  pos < 0 && throw("pos is less than zero")
  mask = (1 << n) - 1
  (bits >> pos) & mask
end

# utility function - returns true if we're done reading bits
function iseob()
  # ensure 1 << pos will fit in an int (assuming 64-bit) - most likely, if
  # there are at least 63 bits remaining, we probably aren't done anyway
  pos > 63 && return false
  mask = (1 << pos) - 1
  iszero(bits & mask)
end

# do eet
function get_packet_value()
  start_pos = pos
  read_bits(3) # version - unneeded
  id = read_bits(3)
  if id == 4
    # literal number
    f = read_bits(1)
    v = BigInt(read_bits(4))
    while f > 0
      f = read_bits(1)
      v = (v << 4) | read_bits(4)
    end
    return (v, start_pos - pos)
  end

  # process subpackets
  stack = BigInt[]
  lentype = read_bits(1)
  if lentype == 0
    # len is total bits of all subpackets
    len = read_bits(15)
    while len > 0
      (v, nread) = get_packet_value()
      push!(stack, v)
      len -= nread
    end
  else
    # len is total number of subpackets
    len = read_bits(11)
    for _ in 1:len
      (v, _) = get_packet_value()
      push!(stack, v)
    end
  end

  # compute result
  result = 0
  if id == 0
    result = sum(stack)
  elseif id == 1
    result = prod(stack)
  elseif id == 2
    result = minimum(stack)
  elseif id == 3
    result = maximum(stack)
  elseif id == 5
    result = BigInt(stack[1] > stack[2])
  elseif id == 6
    result = BigInt(stack[1] < stack[2])
  elseif id == 7
    result = BigInt(stack[1] == stack[2])
  end
  return (result, start_pos - pos)
end

(result, _) = get_packet_value()
println("Result: ", result)


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() # HIDE
let
    x = LinRange(0.01,30π,3000)
    y = cos.(x)
    lines(y, x, color = :orangered, figure = (resolution = (700,450),),
        axis = (yscale = log10, xlabel = "x", ylabel = "y",xgridstyle=:dash,
            ygridstyle=:dash, yminorticksvisible = true,
            yminorticks = IntervalsBetween(9)))
    xlims!(-1,1)
    current_figure()
    save(joinpath(@__DIR__, "output", "LineLogy.png"), current_figure(), px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
# Error for aliases used in recipes
function warn_on_recipe_aliases!(plotattributes, recipe_type, args...)
    for k in keys(plotattributes)
        if !is_default_attribute(k)
            dk = get(_keyAliases, k, k)
            if k !== dk
                @warn "Attribute alias `$k` detected in the $recipe_type recipe defined for the signature $(signature_string(Val{recipe_type}, args...)). To ensure expected behavior it is recommended to use the default attribute `$dk`."
            end
            plotattributes[dk] = pop_kw!(plotattributes, k)
        end
    end
end
function warn_on_recipe_aliases!(v::AbstractVector, recipe_type, args...)
    foreach(x -> warn_on_recipe_aliases!(x, recipe_type, args...), v)
end
function warn_on_recipe_aliases!(rd::RecipeData, recipe_type, args...)
    warn_on_recipe_aliases!(rd.plotattributes, recipe_type, args...)
end

function signature_string(::Type{Val{:user}}, args...)
    return string("(::", join(string.(typeof.(args)), ", ::"), ")")
end
signature_string(::Type{Val{:type}}, T) = "(::Type{$T}, ::$T)"
signature_string(::Type{Val{:plot}}, st) = "(::Type{Val{:$st}}, ::AbstractPlot)"
signature_string(::Type{Val{:series}}, st) = "(::Type{Val{:$st}}, x, y, z)"

# ------------------------------------------------------------------
# preprocessing

function series_idx(kw_list::AVec{KW}, kw::AKW)
    Int(kw[:series_plotindex]) - Int(kw_list[1][:series_plotindex]) + 1
end

function _expand_seriestype_array(plotattributes::AKW, args)
    sts = get(plotattributes, :seriestype, :path)
    if typeof(sts) <: AbstractArray
        reset_kw!(plotattributes, :seriestype)
        rd = Vector{RecipeData}(undef, size(sts, 1))
        for r in axes(sts, 1)
            dc = copy(plotattributes)
            dc[:seriestype] = sts[r:r,:]
            rd[r] = RecipeData(dc, args)
        end
        rd
    else
        RecipeData[RecipeData(copy(plotattributes), args)]
    end
end

function _preprocess_args(plotattributes::AKW, args, still_to_process::Vector{RecipeData})
    # the grouping mechanism is a recipe on a GroupBy object
    # we simply add the GroupBy object to the front of the args list to allow
    # the recipe to be applied
    if haskey(plotattributes, :group)
        args = (extractGroupArgs(plotattributes[:group], args...), args...)
    end

    # if we were passed a vector/matrix of seriestypes and there's more than one row,
    # we want to duplicate the inputs, once for each seriestype row.
    if !isempty(args)
        append!(still_to_process, _expand_seriestype_array(plotattributes, args))
    end

    # remove subplot and axis args from plotattributes... they will be passed through in the kw_list
    if !isempty(args)
        for (k,v) in plotattributes
            if k in _all_subplot_args || k in _all_axis_args
                reset_kw!(plotattributes, k)
            end
        end
    end

    args
end

# ------------------------------------------------------------------
# user recipes

function _process_userrecipes(plt::Plot, plotattributes::AKW, args)
    still_to_process = RecipeData[]
    args = _preprocess_args(plotattributes, args, still_to_process)

    # for plotting recipes, swap out the args and update the parameter dictionary
    # we are keeping a stack of series that still need to be processed.
    # each pass through the loop, we pop one off and apply the recipe.
    # the recipe will return a list a Series objects... the ones that are
    # finished (no more args) get added to the kw_list, the ones that are not
    # are placed on top of the stack and are then processed further.
    kw_list = KW[]
    while !isempty(still_to_process)
        # grab the first in line to be processed and either add it to the kw_list or
        # pass it through apply_recipe to generate a list of RecipeData objects (data + attributes)
        # for further processing.
        next_series = popfirst!(still_to_process)
        # recipedata should be of type RecipeData.  if it's not then the inputs must not have been fully processed by recipes
        if !(typeof(next_series) <: RecipeData)
            error("Inputs couldn't be processed... expected RecipeData but got: $next_series")
        end
        if isempty(next_series.args)
            _process_userrecipe(plt, kw_list, next_series)
        else
            rd_list = RecipesBase.apply_recipe(
                next_series.plotattributes,
                next_series.args...
            )
            warn_on_recipe_aliases!(rd_list, :user, next_series.args...)
            prepend!(still_to_process,rd_list)
        end
    end

    # don't allow something else to handle it
    plotattributes[:smooth] = false
    kw_list
end

function _process_userrecipe(plt::Plot, kw_list::Vector{KW}, recipedata::RecipeData)
    # when the arg tuple is empty, that means there's nothing left to recursively
    # process... finish up and add to the kw_list
    kw = recipedata.plotattributes
    preprocessArgs!(kw)
    _preprocess_userrecipe(kw)
    warnOnUnsupported_scales(plt.backend, kw)

    # add the plot index
    plt.n += 1
    kw[:series_plotindex] = plt.n

    push!(kw_list, kw)
    _add_errorbar_kw(kw_list, kw)
    _add_smooth_kw(kw_list, kw)
    return
end

function _preprocess_userrecipe(kw::AKW)
    _add_markershape(kw)

    # if there was a grouping, filter the data here
    _filter_input_data!(kw)

    # map marker_z if it's a Function
    if isa(get(kw, :marker_z, nothing), Function)
        # TODO: should this take y and/or z as arguments?
        kw[:marker_z] = isa(kw[:z], Nothing) ? map(kw[:marker_z], kw[:x], kw[:y]) : map(kw[:marker_z], kw[:x], kw[:y], kw[:z])
    end

    # map line_z if it's a Function
    if isa(get(kw, :line_z, nothing), Function)
        kw[:line_z] = isa(kw[:z], Nothing) ? map(kw[:line_z], kw[:x], kw[:y]) : map(kw[:line_z], kw[:x], kw[:y], kw[:z])
    end

    # convert a ribbon into a fillrange
    if get(kw, :ribbon, nothing) !== nothing
        make_fillrange_from_ribbon(kw)
    end
    return
end

function _add_errorbar_kw(kw_list::Vector{KW}, kw::AKW)
    # handle error bars by creating new recipedata data... these will have
    # the same recipedata index as the recipedata they are copied from
    for esym in (:xerror, :yerror)
        if get(kw, esym, nothing) !== nothing
            # we make a copy of the KW and apply an errorbar recipe
            errkw = copy(kw)
            errkw[:seriestype] = esym
            errkw[:label] = ""
            errkw[:primary] = false
            push!(kw_list, errkw)
        end
    end
end

function _add_smooth_kw(kw_list::Vector{KW}, kw::AKW)
    # handle smoothing by adding a new series
    if get(kw, :smooth, false)
        x, y = kw[:x], kw[:y]
        β, α = convert(Matrix{Float64}, [x ones(length(x))]) \ convert(Vector{Float64}, y)
        sx = [ignorenan_minimum(x), ignorenan_maximum(x)]
        sy = β .* sx .+ α
        push!(kw_list, merge(copy(kw), KW(
            :seriestype => :path,
            :x => sx,
            :y => sy,
            :fillrange => nothing,
            :label => "",
            :primary => false,
        )))
    end
end

# ------------------------------------------------------------------
# plot recipes

# Grab the first in line to be processed and pass it through apply_recipe
# to generate a list of RecipeData objects (data + attributes).
# If we applied a "plot recipe" without error, then add the returned datalist's KWs,
# otherwise we just add the original KW.
function _process_plotrecipe(plt::Plot, kw::AKW, kw_list::Vector{KW}, still_to_process::Vector{KW})
    if !isa(get(kw, :seriestype, nothing), Symbol)
        # seriestype was never set, or it's not a Symbol, so it can't be a plot recipe
        push!(kw_list, kw)
        return
    end
    try
        st = kw[:seriestype]
        st = kw[:seriestype] = get(_typeAliases, st, st)
        datalist = RecipesBase.apply_recipe(kw, Val{st}, plt)
        warn_on_recipe_aliases!(datalist, :plot, st)
        for data in datalist
            preprocessArgs!(data.plotattributes)
            if data.plotattributes[:seriestype] == st
                error("Plot recipe $st returned the same seriestype: $(data.plotattributes)")
            end
            push!(still_to_process, data.plotattributes)
        end
    catch err
        if isa(err, MethodError)
            push!(kw_list, kw)
        else
            rethrow()
        end
    end
    return
end

# ------------------------------------------------------------------
# setup plot and subplot

function _plot_setup(plt::Plot, plotattributes::AKW, kw_list::Vector{KW})
    # merge in anything meant for the Plot
    for kw in kw_list, (k,v) in kw
        haskey(_plot_defaults, k) && (plotattributes[k] = pop!(kw, k))
    end

    # TODO: init subplots here
    _update_plot_args(plt, plotattributes)
    if !plt.init
        plt.o = Base.invokelatest(_create_backend_figure, plt)

        # create the layout and subplots from the inputs
        plt.layout, plt.subplots, plt.spmap = build_layout(plt.attr)
        for (idx,sp) in enumerate(plt.subplots)
            sp.plt = plt
            sp.attr[:subplot_index] = idx
        end

        plt.init = true
    end

    # handle inset subplots
    insets = plt[:inset_subplots]
    if insets !== nothing
        if !(typeof(insets) <: AVec)
            insets = [insets]
        end
        for inset in insets
            parent, bb = is_2tuple(inset) ? inset : (nothing, inset)
            P = typeof(parent)
            if P <: Integer
                parent = plt.subplots[parent]
            elseif P == Symbol
                parent = plt.spmap[parent]
            else
                parent = plt.layout
            end
            sp = Subplot(backend(), parent=parent)
            sp.plt = plt
            push!(plt.subplots, sp)
            push!(plt.inset_subplots, sp)
            sp.attr[:relative_bbox] = bb
            sp.attr[:subplot_index] = length(plt.subplots)
        end
    end
    plt[:inset_subplots] = nothing
end

function _subplot_setup(plt::Plot, plotattributes::AKW, kw_list::Vector{KW})
    # we'll keep a map of subplot to an attribute override dict.
    # Subplot/Axis attributes set by a user/series recipe apply only to the
    # Subplot object which they belong to.
    # TODO: allow matrices to still apply to all subplots
    sp_attrs = Dict{Subplot,Any}()
    for kw in kw_list
        # get the Subplot object to which the series belongs.
        sps = get(kw, :subplot, :auto)
        sp = get_subplot(plt, _cycle(sps == :auto ? plt.subplots : plt.subplots[sps], series_idx(kw_list,kw)))
        kw[:subplot] = sp

        # extract subplot/axis attributes from kw and add to sp_attr
        attr = KW()
        for (k,v) in collect(kw)
            if is_subplot_attr(k) || is_axis_attr(k)
                attr[k] = pop!(kw, k)
            end
            if is_axis_attr_noletter(k)
                v = pop!(kw, k)
                for letter in (:x,:y,:z)
                    attr[Symbol(letter,k)] = v
                end
            end
            for k in (:scale,), letter in (:x,:y,:z)
                # Series recipes may need access to this information
                lk = Symbol(letter,k)
                if haskey(attr, lk)
                    kw[lk] = attr[lk]
                end
            end
        end
        sp_attrs[sp] = attr
    end

    # override subplot/axis args.  `sp_attrs` take precendence
    for (idx,sp) in enumerate(plt.subplots)
        attr = if !haskey(plotattributes, :subplot) || plotattributes[:subplot] == idx
            merge(plotattributes, get(sp_attrs, sp, KW()))
        else
            get(sp_attrs, sp, KW())
        end
        _update_subplot_args(plt, sp, attr, idx, false)
    end

    # do we need to link any axes together?
    link_axes!(plt.layout, plt[:link])
end

# getting ready to add the series... last update to subplot from anything
# that might have been added during series recipes
function _prepare_subplot(plt::Plot{T}, plotattributes::AKW) where T
    st::Symbol = plotattributes[:seriestype]
    sp::Subplot{T} = plotattributes[:subplot]
    sp_idx = get_subplot_index(plt, sp)
    _update_subplot_args(plt, sp, plotattributes, sp_idx, true)

    st = _override_seriestype_check(plotattributes, st)

    # change to a 3d projection for this subplot?
    if is3d(st)
        sp.attr[:projection] = "3d"
    end

    # initialize now that we know the first series type
    if !haskey(sp.attr, :init)
        _initialize_subplot(plt, sp)
        sp.attr[:init] = true
    end
    sp
end

# ------------------------------------------------------------------
# series types

function _override_seriestype_check(plotattributes::AKW, st::Symbol)
    # do we want to override the series type?
    if !is3d(st) && !(st in (:contour,:contour3d))
        z = plotattributes[:z]
        if !isa(z, Nothing) && (size(plotattributes[:x]) == size(plotattributes[:y]) == size(z))
            st = (st == :scatter ? :scatter3d : :path3d)
            plotattributes[:seriestype] = st
        end
    end
    st
end

function _prepare_annotations(sp::Subplot, plotattributes::AKW)
    # strip out series annotations (those which are based on series x/y coords)
    # and add them to the subplot attr
    sp_anns = annotations(sp[:annotations])
    # series_anns = annotations(pop!(plotattributes, :series_annotations, []))
    # if isa(series_anns, SeriesAnnotations)
    #     series_anns.x = plotattributes[:x]
    #     series_anns.y = plotattributes[:y]
    # elseif length(series_anns) > 0
    #     x, y = plotattributes[:x], plotattributes[:y]
    #     nx, ny, na = map(length, (x,y,series_anns))
    #     n = max(nx, ny, na)
    #     series_anns = [(x[mod1(i,nx)], y[mod1(i,ny)], text(series_anns[mod1(i,na)])) for i=1:n]
    # end
    # sp.attr[:annotations] = vcat(sp_anns, series_anns)
end

function _expand_subplot_extrema(sp::Subplot, plotattributes::AKW, st::Symbol)
    # adjust extrema and discrete info
    if st == :image
        xmin, xmax = ignorenan_extrema(plotattributes[:x]); ymin, ymax = ignorenan_extrema(plotattributes[:y])
        expand_extrema!(sp[:xaxis], (xmin, xmax))
        expand_extrema!(sp[:yaxis], (ymin, ymax))
    elseif !(st in (:pie, :histogram, :bins2d, :histogram2d))
        expand_extrema!(sp, plotattributes)
    end
    # expand for zerolines (axes through origin)
    if sp[:framestyle] in (:origin, :zerolines)
        expand_extrema!(sp[:xaxis], 0.0)
        expand_extrema!(sp[:yaxis], 0.0)
    end
end

function _add_the_series(plt, sp, plotattributes)
    warnOnUnsupported_args(plt.backend, plotattributes)
    warnOnUnsupported(plt.backend, plotattributes)
    series = Series(plotattributes)
    push!(plt.series_list, series)
    push!(sp.series_list, series)
    _series_added(plt, series)
end

# -------------------------------------------------------------------------------

# this method recursively applies series recipes when the seriestype is not supported
# natively by the backend
function _process_seriesrecipe(plt::Plot, plotattributes::AKW)
    #println("process $(typeof(plotattributes))")
    # replace seriestype aliases
    st = Symbol(plotattributes[:seriestype])
    st = plotattributes[:seriestype] = get(_typeAliases, st, st)

    # shapes shouldn't have fillrange set
    if plotattributes[:seriestype] == :shape
        plotattributes[:fillrange] = nothing
    end

    # if it's natively supported, finalize processing and pass along to the backend, otherwise recurse
    if is_seriestype_supported(st)
        sp = _prepare_subplot(plt, plotattributes)
        _prepare_annotations(sp, plotattributes)
        _expand_subplot_extrema(sp, plotattributes, st)
        _update_series_attributes!(plotattributes, plt, sp)
        _add_the_series(plt, sp, plotattributes)

    else
        # get a sub list of series for this seriestype
        x, y, z = plotattributes[:x], plotattributes[:y], plotattributes[:z]
        datalist = RecipesBase.apply_recipe(plotattributes, Val{st}, x, y, z)
        warn_on_recipe_aliases!(datalist, :series, st)

        # assuming there was no error, recursively apply the series recipes
        for data in datalist
            if isa(data, RecipeData)
                preprocessArgs!(data.plotattributes)
                if data.plotattributes[:seriestype] == st
                    error("The seriestype didn't change in series recipe $st.  This will cause a StackOverflow.")
                end
                _process_seriesrecipe(plt, data.plotattributes)
            else
                @warn("Unhandled recipe: $(data)")
                break
            end
        end
    end
    nothing
end


text:
# %% imports
using Cambrian
using EvolutionaryStrategies
using JLD
using PyCall
using Statistics
gym = pyimport("gym")
include("../../Modules/Networks/net_arch.jl")
include("../../Modules/Poisson/PoissonStateSpike.jl")
include("../../Actors/cartpole_fitness.jl")

function create_net(elites::Array{AbstractESIndividual,1}, indiv::Int64)
    params = YAML.load(open("Julia/Actors/cartpole_cfglif.yml"))
    arch = params["arch"]

    genes = elites[indiv].genes

    # Input linear combination
    genes_alpha = map(x->borne(x,-1.0, 1.0),genes[1:arch[1]*4])
    matalpha = reshape(genes_alpha, (arch[1],4))

    # Weights of the network
    genes_w = genes[arch[1]*4+1:end]

    # Initialize network
    n = Network(arch, params)
    exc, inh = ei_idx(arch)
    set_weight!(n,genes_w, exc, inh)

    return n, matalpha
end

function init_teacher_student(filepath::String, indiv::Int64)
    elites = load(filepath, "elites")
    n_teach, matalpha = create_net(elites, indiv)
    n_stud = deepcopy(n_teach)
    n_stud.s = Weights(n_stud.arch, n_stud.param)

    return n_teach, n_stud, matalpha
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
struct CompositeMap{T, As<:LinearMapTuple} <: LinearMap{T}
    maps::As # stored in order of application to vector
    function CompositeMap{T, As}(maps::As) where {T, As}
        N = length(maps)
        for n in 2:N
            check_dim_mul(maps[n], maps[n-1])
        end
        for TA in Base.Generator(eltype, maps)
            # like lazy map; could use Base.Iterators.map in Julia >= 1.6
            promote_type(T, TA) == T ||
                error("eltype $TA cannot be promoted to $T in CompositeMap constructor")
        end
        new{T, As}(maps)
    end
end
CompositeMap{T}(maps::As) where {T, As<:LinearMapTuple} = CompositeMap{T, As}(maps)

# basic methods
Base.size(A::CompositeMap) = (size(A.maps[end], 1), size(A.maps[1], 2))
Base.isreal(A::CompositeMap) = all(isreal, A.maps) # sufficient but not necessary

# the following rules are sufficient but not necessary
for (f, _f, g) in ((:issymmetric, :_issymmetric, :transpose),
                    (:ishermitian, :_ishermitian, :adjoint),
                    (:isposdef, :_isposdef, :adjoint))
    @eval begin
        LinearAlgebra.$f(A::CompositeMap) = $_f(A.maps)
        $_f(maps::Tuple{}) = true
        $_f(maps::Tuple{<:LinearMap}) = $f(maps[1])
        $_f(maps::Tuple{Vararg{LinearMap}}) =
            maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps)))
        # since the introduction of ScaledMap, the following cases cannot occur
        # function $_f(maps::Tuple{Vararg{LinearMap}}) # length(maps) >= 2
            # if maps[1] isa UniformScalingMap{<:RealOrComplex}
            #     return $f(maps[1]) && $_f(Base.tail(maps))
            # elseif maps[end] isa UniformScalingMap{<:RealOrComplex}
            #     return $f(maps[end]) && $_f(Base.front(maps))
            # else
                # return maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps)))
            # end
        # end
    end
end

# scalar multiplication and division (non-commutative case)
function Base.:(*)(α::Number, A::LinearMap)
    T = promote_type(typeof(α), eltype(A))
    return CompositeMap{T}(tuple(A, UniformScalingMap(α, size(A, 1))))
end
function Base.:(*)(α::Number, A::CompositeMap)
    T = promote_type(typeof(α), eltype(A))
    Alast = last(A.maps)
    if Alast isa UniformScalingMap
        return CompositeMap{T}(tuple(Base.front(A.maps)..., α * Alast))
    else
        return CompositeMap{T}(tuple(A.maps..., UniformScalingMap(α, size(A, 1))))
    end
end
# needed for disambiguation
function Base.:(*)(α::RealOrComplex, A::CompositeMap)
    T = Base.promote_op(*, typeof(α), eltype(A))
    return ScaledMap{T}(α, A)
end
function Base.:(*)(A::LinearMap, α::Number)
    T = promote_type(typeof(α), eltype(A))
    return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A))
end
function Base.:(*)(A::CompositeMap, α::Number)
    T = promote_type(typeof(α), eltype(A))
    Afirst = first(A.maps)
    if Afirst isa UniformScalingMap
        return CompositeMap{T}(tuple(Afirst * α, Base.tail(A.maps)...))
    else
        return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A.maps...))
    end
end
# needed for disambiguation
function Base.:(*)(A::CompositeMap, α::RealOrComplex)
    T = Base.promote_op(*, typeof(α), eltype(A))
    return ScaledMap{T}(α, A)
end

Base.:(\)(α::Number, A::LinearMap) = inv(α) * A
Base.:(/)(A::LinearMap, α::Number) = A * inv(α)

# composition of linear maps
"""
    *(A::LinearMap, B::LinearMap)::CompositeMap

Construct a (lazy) representation of the product of the two operators.
Products of `LinearMap`/`CompositeMap` objects and `LinearMap`/`CompositeMap`
objects are reduced to a single `CompositeMap`. In products of `LinearMap`s and
`AbstractMatrix`/`UniformScaling` objects, the latter get promoted to `LinearMap`s
automatically.

# Examples
```jldoctest; setup=(using LinearAlgebra, LinearMaps)
julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;

julia> LinearMap(ones(Int, 3, 3)) * CS * I * rand(3, 3);
```
"""
function Base.:(*)(A₁::LinearMap, A₂::LinearMap)
    T = promote_type(eltype(A₁), eltype(A₂))
    return CompositeMap{T}(tuple(A₂, A₁))
end
function Base.:(*)(A₁::LinearMap, A₂::CompositeMap)
    T = promote_type(eltype(A₁), eltype(A₂))
    return CompositeMap{T}(tuple(A₂.maps..., A₁))
end
function Base.:(*)(A₁::CompositeMap, A₂::LinearMap)
    T = promote_type(eltype(A₁), eltype(A₂))
    return CompositeMap{T}(tuple(A₂, A₁.maps...))
end
function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap)
    T = promote_type(eltype(A₁), eltype(A₂))
    return CompositeMap{T}(tuple(A₂.maps..., A₁.maps...))
end
# needed for disambiguation
Base.:(*)(A₁::ScaledMap, A₂::CompositeMap) = A₁.λ * (A₁.lmap * A₂)
Base.:(*)(A₁::CompositeMap, A₂::ScaledMap) = (A₁ * A₂.lmap) * A₂.λ

# special transposition behavior
LinearAlgebra.transpose(A::CompositeMap{T}) where {T} =
    CompositeMap{T}(map(transpose, reverse(A.maps)))
LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} =
    CompositeMap{T}(map(adjoint, reverse(A.maps)))

# comparison of CompositeMap objects
Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps)

# multiplication with vectors
_unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector) =
    _compositemul!(y, A, x)
function _compositemul!(y::AbstractVecOrMat,
                        A::CompositeMap{<:Any,<:Tuple{LinearMap}},
                        x::AbstractVector,
                        source = nothing,
                        dest = nothing)
    return _unsafe_mul!(y, A.maps[1], x)
end
function _compositemul!(y::AbstractVecOrMat,
                        A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x::AbstractVector,
                        source = similar(y, size(A.maps[1], 1)),
                        dest = nothing)
    _unsafe_mul!(source, A.maps[1], x)
    _unsafe_mul!(y, A.maps[2], source)
    return y
end
function _compositemul!(y::AbstractVecOrMat,
                        A::CompositeMap,
                        x::AbstractVector,
                        source = similar(y, size(A.maps[1], 1)),
                        dest = similar(y, size(A.maps[2], 1)))
    N = length(A.maps)
    _unsafe_mul!(source, A.maps[1], x)
    for n in 2:N-1
        try
            resize!(dest, size(A.maps[n], 1))
        catch err
            if err == ErrorException("cannot resize array with shared data")
                dest = similar(y, size(A.maps[n], 1))
            else
                rethrow(err)
            end
        end
        _unsafe_mul!(dest, A.maps[n], source)
        dest, source = source, dest # alternate dest and source
    end
    _unsafe_mul!(y, A.maps[N], source)
    return y
end


text:
#=
Created on Sunday 19 April 2020
Last update: Sunday 16 April 2020

@author: Michiel Stock
michielfmstock@gmail.com

Utilities for the TSP problem.
=#

module TSP

import STMO:dist
using STMO
using Plots: plot!, scatter, scatter!, plot

export TravelingSalesmanProblem, cities, computecost, coordinates
export plot_cities, plot_cities!, plot_tour, plot_tour!
export swap!, deltaswapcost, deltaflipcost, flip!
export totoro_tsp, got_coords

struct TravelingSalesmanProblem{Tc}
    coordinates::Matrix{Tc}
    distance::Matrix{Float64}
end

TravelingSalesmanProblem(coordinates) = TravelingSalesmanProblem(coordinates,
                                dist(coordinates, coordinates))

"""Returns the number of cities."""
Base.length(tsp::TravelingSalesmanProblem) = size(tsp.coordinates, 1)

"""Returns the coordinates of the cities."""
coordinates(tsp::TravelingSalesmanProblem) = tsp.coordinates

"""Returns the cities"""
cities(tsp::TravelingSalesmanProblem) = collect(1:length(tsp))

"""Returns the distance matrix of the cities."""
dist(tsp::TravelingSalesmanProblem) = tsp.distance

"""Returns the distance (or cost) of going from city `ci` to city `cj`"""
dist(tsp::TravelingSalesmanProblem, ci, cj) = dist(tsp)[ci,cj]

Base.isvalid(tsp::TravelingSalesmanProblem, tour) = length(tour) == length(tsp) &&
    Set(tour) == Set(cities(tsp))

"""
    computecost(tsp::TravelingSalesmanProblem, tour)

Computes the cost of travessing a tour.
"""
function computecost(tsp::TravelingSalesmanProblem, tour)
    !isvalid(tsp, tour) && throw(AssertionError("invalid tour provided"))
    c = 0.0
    for (i, j) in zip(tour[1:end-1], tour[2:end])
        c += dist(tsp, i, j)
    end
    # complete tour
    c += dist(tsp, tour[end], tour[1])
    return c
end

split_coord(X) = X[:,1], X[:,2]

plot_cities(tsp::TravelingSalesmanProblem; kwargs...) = scatter(split_coord(coordinates(tsp))...,
                color=myblue, markersize=1, label="", aspect_ratio=:equal; kwargs...)

plot_cities!(tsp::TravelingSalesmanProblem; kwargs...) = scatter!(split_coord(coordinates(tsp))...,
                color=myblue, markersize=1, label=""; kwargs...)

coords_tour(tsp, tour) = [coordinates(tsp)[tour,:];coordinates(tsp)[[tour[1]],:]]

plot_tour(tsp::TravelingSalesmanProblem, tour; kwargs...) = plot(
                split_coord(coords_tour(tsp, tour))...,
                color=myred, label="", aspect_ratio=:equal; kwargs...)

plot_tour!(tsp::TravelingSalesmanProblem, tour; kwargs...) = plot!(
                split_coord(coords_tour(tsp, tour))...,
                color=myred, label=""; kwargs...)

"""
    deltaswapcost(tsp, tour, i, j)

Compute the change in tour cost if the cities at positions `i` and `j` are
swapped.
"""
function deltaswapcost(tsp, tour, i, j)
    n = length(tsp)
    i == j && return 0.0
    # put in order
    i, j = i < j ? (i, j) : (j, i)
    # choose indices respecting cyclic  invariance
    i₋₁ = i == 1 ? n : i - 1
    i₊₁ = i == n ? 1 : i + 1
    j₋₁ = j == 1 ? n : j - 1
    j₊₁ = j == n ? 1 : j + 1
    ci₋₁, ci, ci₊₁, cj₋₁, cj, cj₊₁ = tour[[i₋₁, i, i₊₁, j₋₁, j, j₊₁]]
    if j - i == 1  # i and j are neighbors
        Δc = ((dist(tsp, ci₋₁, cj) + dist(tsp, ci, cj₊₁)) -
                    (dist(tsp, ci₋₁, ci) + dist(tsp, cj, cj₊₁)))
    elseif (i==1 && j==n)
        Δc = ((dist(tsp, tour[end-1], tour[1]) + dist(tsp, tour[end], tour[2])) -
                    (dist(tsp, tour[end-1], tour[end]) + dist(tsp, tour[1], tour[2])))
    else
        Δc = ((dist(tsp, ci₋₁, cj) + dist(tsp, cj, ci₊₁) +
                    dist(tsp, cj₋₁, ci) + dist(tsp, ci, cj₊₁)) -
                    (dist(tsp, ci₋₁, ci) + dist(tsp, ci, ci₊₁) +
                    dist(tsp, cj₋₁, cj) + dist(tsp, cj, cj₊₁)))
    end
    return Δc
end

"""
    swap!(tour, i, j)

Swaps the cities at positions `i` and `j` in `tour`.
"""
function swap!(tour, i, j)
    tour[i], tour[j] = tour[j], tour[i]
    return tour
end

"""
    deltaflipcost(tsp, tour, i, j)

Compute the change in tour cost if the subtour between the `i`th and `j`th city
are flipped.
"""
function deltaflipcost(tsp, tour, i, j)
    i == j && return 0.0
    n = length(tsp)
    # put in order
    i, j = i < j ? (i, j) : (j, i)
    (i, j) == (1, n) && return 0.0
    # choose indices respecting cyclic  invariance
    i₋₁ = i == 1 ? n : i - 1
    i₊₁ = i == n ? 1 : i + 1
    j₋₁ = j == 1 ? n : j - 1
    j₊₁ = j == n ? 1 : j + 1
    ci₋₁, ci, ci₊₁, cj₋₁, cj, cj₊₁ = tour[[i₋₁, i, i₊₁, j₋₁, j, j₊₁]]
    Δc = - (dist(tsp, ci₋₁, ci) + dist(tsp, cj, cj₊₁))
    Δc += dist(tsp, ci₋₁, cj) + dist(tsp, ci, cj₊₁)
    return Δc
end

function flip!(tour, i, j)
    i, j = i < j ? (i, j) : (j, i)
    reverse!(@view tour[i:j])
    return tour
end

# example network
include("targ.jl")
include("totoro.jl")

got_tsp() = TravelingSalesmanProblem(got_coords)
totoro_tsp() = TravelingSalesmanProblem(totoro_coords)

end  # module TSP


text:
@testset verbose=true "Test Java Code Generation" begin
    using JavaCall: JavaCodeGeneration

    @testset "Test Simple Object" begin
        obj_module = Module(:Object)
        Core.eval(obj_module, :(using JavaCall))
        Core.eval(
            obj_module,
            JavaCodeGeneration.loadclass(obj_module, Symbol("java.lang.Object"))
        )
        @test isdefined(obj_module, :JObject)
        @test isdefined(obj_module, :JObjectJuliaImpl)
        @test isdefined(obj_module, :JString)
        @test isdefined(obj_module, :JStringJuliaImpl)
        @test isdefined(obj_module, :j_equals)
        @test isdefined(obj_module, :j_to_string)
    end

    @testset "Test Static Methods" begin
        @testset "Test Arrays" begin
            arrays_mod = Module(:Arrays)
            Core.eval(arrays_mod, :(using JavaCall))
            Core.eval(arrays_mod, JavaCodeGeneration.loadclass(
                arrays_mod,
                Symbol("java.util.Arrays")
            ))

            a = [1, 2, 3, 4, 5]
            # Java Indexes at 0
            @test 1 == arrays_mod.j_binary_search(arrays_mod.JArrays, a, 2)
        end

        @testset "Test Dates" begin
            dates_mod = Module(:Dates)
            Core.eval(dates_mod, :(using JavaCall))
            Core.eval(dates_mod, JavaCodeGeneration.loadclass(
                dates_mod,
                Symbol("java.time.LocalDate")
            ))
            a = dates_mod.j_of(dates_mod.JLocalDate, Int32(2000), Int32(1), Int32(1))
            b = dates_mod.j_plus_days(a, 1)
            c = dates_mod.j_plus_months(b, 1)
            @test dates_mod.j_get_year(a) == 2000
            @test dates_mod.j_get_year(b) == 2000
            @test dates_mod.j_get_year(c) == 2000
            @test dates_mod.j_get_month_value(a) == 1
            @test dates_mod.j_get_month_value(b) == 1
            @test dates_mod.j_get_month_value(c) == 2
            @test dates_mod.j_get_day_of_month(a) == 1
            @test dates_mod.j_get_day_of_month(b) == 2
            @test dates_mod.j_get_day_of_month(c) == 2
        end
    end

    @testset "Test Constructors" begin
        cons_mod = Module(:Constructors)
        Core.eval(cons_mod, :(using JavaCall))
        Core.eval(
            cons_mod,
            JavaCodeGeneration.loadclass(cons_mod, Symbol("java.lang.String"))
        )

        hellochars = ['H', 'e', 'l', 'l', 'o']
        hello1 = cons_mod.JString(hellochars)
        hello2 = cons_mod.JString(hello1)

        @test cons_mod.j_length(hello1) == 5
        @test cons_mod.j_length(hello2) == 5

        for (i, c) in enumerate(hellochars)
            @test cons_mod.j_char_at(hello1, Int32(i-1)) == c
            @test cons_mod.j_char_at(hello2, Int32(i-1)) == c
        end

        helloworldchars = ['H', 'e', 'l', 'l', 'o', ' ', 'W', 'o', 'r', 'l', 'd']

        finalworld = cons_mod.JString(helloworldchars, Int32(5), Int32(6))

        helloworld = cons_mod.j_concat(hello1, finalworld)
        for (i, c) in enumerate(helloworldchars)
            @test cons_mod.j_char_at(helloworld, Int32(i-1)) == c
        end

        @test cons_mod.j_equals(hello1, hello2)
        @test_false cons_mod.j_equals(hello1, helloworld)
        @test hello1 == hello2
        @test hello1 != helloworld

        # Test string special constructor
        helloworld2 = cons_mod.JString("Hello World")
        @test cons_mod.j_equals(helloworld, helloworld2)
        @test helloworld == helloworld2

        # Test string method to convert JString
        @test cons_mod.string(helloworld2) == "Hello World"
    end

    @testset "Test Load Superclass" begin
        super_mod = Module(:Constructors)
        Core.eval(super_mod, :(using JavaCall))
        Core.eval(
            super_mod,
            JavaCodeGeneration.loadclass(super_mod, Symbol("java.lang.Integer"))
        )
        @test isdefined(super_mod, :JInteger)
        @test isdefined(super_mod, :JIntegerJuliaImpl)
        @test isdefined(super_mod, :JNumber)
        @test isdefined(super_mod, :JNumberJuliaImpl)

        a = super_mod.JInteger(Int32(1))
        b = super_mod.JInteger(Int32(1))
        @test 1 == super_mod.j_long_value(a)
        @test a == b
    end

    @testset "Test Load Return Value" begin
        net_mod = Module(:Net)
        Core.eval(net_mod, :(using JavaCall))
        Core.eval(net_mod, JavaCodeGeneration.loadclass(
            net_mod,
            Symbol("javax.net.SocketFactory")
        ))
        socketfactory = net_mod.j_get_default(net_mod.JSocketFactory)
        socket = net_mod.j_create_socket(socketfactory)
        @test isdefined(net_mod, :j_bind) # Socket.bind
    end
end


text:
function print_results(m, frac, nodeLabelSize, edgeLabelSize, T)
    x, y = m[:x], m[:y]

    println(termination_status(m))
    println("$T : ", JuMP.value(x[T]))
    println("Obj : ", round(objective_value(m), digits=3), "\n")

    returnedProducts = []
    returnedRecipes = []

    println("Products : ")
    for p in subtypes(Product)
        production = getvalue(x[p])
        if production > 0
            println(p, " : ", round(production, digits=3))
            push!(returnedProducts, (p, production))
        end
    end

    println("\nRecipes : ")
    for r in allRecipes
        used = getvalue(y[r]) * frac
        if (used > 1e-9)
            println(r)
            if !isempty(r.in)
                print("\tinput : ")
                for p in r.in
                    print(p[1], " : ", round(p[2] * used, digits=4), "/min ", )
                end
                println()
            end
            for p in r.out
                println("\t", p[1], " : ", round(p[2] * used, digits=4), "/min")
            end
            println("\t", round(used, digits=3), " ", r.building)
            push!(returnedRecipes, (r, used))
        end
    end

    buildGraph(returnedRecipes, nodeLabelSize, edgeLabelSize, T, objective_value(m))
end

intOrRound(x, d = 2) = isinteger(round(x, digits=3)) ? string(round(Int, x)) : string(round(x, digits = d))

function buildGraph(returnedRecipes, nodeLabelSize, edgeLabelSize, T, objValue)
    isempty(returnedRecipes) && return
    g = SimpleDiGraph()
    add_vertices!(g, length(returnedRecipes))
    indices = Dict(r => i for (i, r) in enumerate(first.(returnedRecipes)))

    # sinks = filter(rec -> any(x -> x[1] == T, rec[1].out), returnedRecipes)

    edgeLabels = Dict{Tuple{Int, Int}, Float64}()
    edgeProducts = Dict{Tuple{Int, Int}, String}()
    removedEdges = Tuple{Int, Int}[]

    for (r1, used1) in returnedRecipes, (r2, used2) in returnedRecipes
        for (p1, qty1) in r1.out, (p2, qty2) in r2.in
            if p1 == p2
                add_edge!(g, indices[r1], indices[r2])
                push!(edgeLabels, (indices[r1], indices[r2]) => min(qty1 * used1, qty2 * used2))
                push!(edgeProducts, (indices[r1], indices[r2]) => string(p1))
                if LightGraphs.is_cyclic(g)
                    rem_edge!(g, indices[r1], indices[r2])
                    println("removed acyclic edge ", r1, " => ", r2, " : ", intOrRound(min(qty1 * used1, qty2 * used2)))
                    push!(removedEdges, (indices[r1], indices[r2]))
                end
            end
        end
    end

    labels = [r.name  * "(" * intOrRound(used) * ")" for (r, used) in returnedRecipes]
    for (i, (r, used)) in enumerate(returnedRecipes)
        x = findfirst(y -> y[1] == T, r.out)
        if x !== nothing
            labels[i] = labels[i] * " => " * intOrRound(r.out[x][2] * used)
        end
    end
    xs, ys, edge2path = LayeredLayouts.solve_positions(LayeredLayouts.Zarate(time_limit=LayeredLayouts.Dates.Second(0)), g)
    p = plot(showaxis=false, ticks=false, xlims = (0.5, maximum(xs) + 2.), dpi = 300)

    for (edge, path) in edge2path
        x, y = path
        ann = text(edgeProducts[(edge.src, edge.dst)] * "\n" * intOrRound(edgeLabels[(edge.src, edge.dst)]), pointsize=edgeLabelSize, valign=y[end] + 1 < y[end-1] ? :bottom : :top)
        curves!(p, x, y, arrow=Plots.arrow(:closed), label=nothing)
        annotate!(p, [((x[end-1] + 0.4 * (x[end] - x[end-1])), (y[end-1] + 0.4 * (y[end] - y[end-1])), ann)])
    end

    for e in removedEdges
        x = [xs[e[1]], (xs[e[1]] + xs[e[2]]) / 2, xs[e[2]]]
        y = [ys[e[1]], (ys[e[1]] + ys[e[2]]) / 2 + 0.15, ys[e[2]]]
        ann = text(intOrRound(edgeLabels[e[1], e[2]]), pointsize=edgeLabelSize, valign=:top)
        curves!(p, x, y, arrow=Plots.arrow(:closed), label=nothing)
        annotate!(p, [((x[end-1] + x[end]) / 2, (y[end-1] + y[end]) / 2, ann)])
    end

    scatter!(p, xs, ys, text=text.(labels, pointsize=nodeLabelSize, valign=:bottom, halign=:left), label=nothing,)

    p |> display
end


text:
#ДАНО: Робот - в юго-западном углу поля, на котором расставлено некоторое количество маркеров
#РЕЗУЛЬТАТ: Функция вернула значение средней температуры всех замаркированных клеток

function stemp(r::Robot, side::HorizonSide) #возврат значения суммы температур всех клеток в направлении
    stemp=0
    while isborder(r,side)==false
        stemp+=check_temp(r)
        move!(r,side)
    end
    stemp+=check_temp(r)
    return stemp
end

s=0
function glavnaia(r::Robot)
    global s=0
    overal_temp=0
    side=Ost
    while isborder(r,Nord)==false
        overal_temp+=stemp(r,side)
        move!(r,Nord)
        side=inverse(side)
    end
    overal_temp+=stemp(r,side)
    average_temp=overal_temp/s
    return average_temp
end

function check_temp(r::Robot)
    global s
    if ismarker(r)==true
        s+=1
        return temperature(r)
    else
        return 0
    end
end

inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4))

max_stars_repo_name: Aminulik/mirea-progs

max_issues_repo_name: Aminulik/mirea-progs

max_forks_repo_name: Aminulik/mirea-progs


text:
using StatsBase

include("discrete.jl")

const LEFT = 1
const DOWN = 2
const RIGHT = 3
const UP = 4

MAPS = Dict(
    :fourxfour => [
        "SFFF",
        "FHFH",
        "FFFH",
        "HFFG"
    ],

    :eightxeight => [
        "SFFFFFFF",
        "FFFFFFFF",
        "FFFHFFFF",
        "FFFFFHFF",
        "FFFHFFFF",
        "FHHFFFHF",
        "FHFFHFHF",
        "FFFHFFFG"
    ]
)

"""
    generate_random_map(size::Integer, prob::AbstractFloat)
Generates a random valid map of size `size` and `prob` determines the probability
of a tile being frozen.
"""
function generate_random_map(shape=8, p=8f-1)
    valid = false

    # BFS to check that it's a valid path...
    function is_valid(arr, r=1, c=1)
        if arr[r, c] == 'G'
            return true
        end

        tmp = arr[r, c]
        arr[r, c] = '#'

        if r+1 ≤ shape && !(arr[r+1, c] ∈ "#H")
            if is_valid(arr, r+1, c)
                arr[r, c] = tmp
                return true
            end
        end

        if c+1 ≤ shape && !(arr[r, c+1] ∈ "#H")
            if is_valid(arr, r, c+1)
                arr[r, c] = tmp
                return true
            end
        end

        if r-1 ≥ 1 && !(arr[r-1, c] ∈ "#H")
            if is_valid(arr, r-1, c)
                arr[r, c] = tmp
                return true
            end
        end

        if c-1 ≥ 1 && !(arr[r, c-1] ∈ "#H")
            if is_valid(arr, r, c-1)
                arr[r, c] = tmp
                return true
            end
        end
        arr[r, c] = tmp
        return false
    end

    res = Array{Char, 2}(undef, shape, shape)
    while !valid
        p = min(1, p)
        res = sample(['F', 'H'], ProbabilityWeights(Float32[p, 1-p]), (shape, shape), replace=true)
        res[1, 1] = 'S'
        res[end, end] = 'G'
        valid = is_valid(res)
    end
    return [join(res[x, :]) for x in 1:size(res, 1)]
end

mutable struct FrozenLakeEnv <: DiscreteEnv
    desc::Array{Char, 2}
    nrow::Int
    ncol::Int
    reward_range
    discenv_obj::DiscreteEnvObj
end

function FrozenLakeEnv(desc_::AbstractArray{<:AbstractString, 1}, is_slippery::Bool)
    desc = Char.(vcat([Int.(collect(row))' for row in desc_]...))

    nrow, ncol = size(desc)

    reward_range = [0, 1]

    nA = 4
    nS = nrow * ncol

    isd = desc .== 'S'
    isd ./= sum(isd)
    isd = reshape(isd, length(isd))

    P = Dict(s => Dict(a => [] for a=1:nA) for s=1:nS)

    @inline to_s(row, col) = (col-1)*nrow + row

    function inc(row, col, a)
        if a == LEFT
            col = max(col-1, 1)
        elseif a == RIGHT
            col = min(col+1, ncol)
        elseif a == DOWN
            row = min(row+1, nrow)
        elseif a == UP
            row = max(row-1, 1)
        end
        return (row, col)
    end

    for row=1:nrow
        for col=1:ncol
            s = to_s(row, col)
            for a=1:4
                li = P[s][a]
                letter = desc[row, col]
                if letter ∈ "GH"
                    push!(li, (1.0, s, 0, true))
                else
                    if is_slippery
                        for b ∈[((a-2)%4)+1, a, ((a%4)+1)]
                            newrow, newcol = inc(row, col, b)
                            newstate = to_s(newrow, newcol)
                            newletter = desc[newrow, newcol]
                            done = newletter == 'G' || newletter == 'H'
                            rew = Float32(newletter == 'G')
                            push!(li, (1f0/3f0, newstate, rew, done))
                        end
                    else
                        newrow, newcol = inc(row, col, a)
                        newstate = to_s(newrow, newcol)
                        newletter = desc[newrow, newcol]
                        done = newletter == 'G' || newletter == 'H'
                        rew = Float32(newletter == 'G')
                        push!(li, (1f0, newstate, rew, done))
                    end
                end
            end
        end
    end
    discenv_obj = DiscreteEnvObj(nS, nA, P, isd)
    FrozenLakeEnv(desc, nrow, ncol, reward_range, discenv_obj)
end

function FrozenLakeEnv(map_name::Symbol, is_slippery::Bool=true)
    desc = MAPS[map_name]
    FrozenLakeEnv(desc, is_slippery)
end

function FrozenLakeEnv(is_slippery::Bool=true)
    desc = generate_random_map()
    FrozenLakeEnv(desc, is_slippery)
end


text:
"""
This section relies heavily on AbstractTrees.jl
The goal is to build a tree structure starting from
an arbitrary type to all other defined types.
Once this tree is created, it is relatively easy to
go back and find where I want to get on the tree
(basically, the end type) and see the tables necessary
to get there.
"""

mutable struct FirmIdentifierNode

    data::Type{<:FirmIdentifier}
    parent::FirmIdentifierNode
    children::Vector{FirmIdentifierNode}

    # Root constructor
    FirmIdentifierNode(data) = new(data)
    # child node constructor
    FirmIdentifierNode(data, parent::FirmIdentifierNode) = new(data, parent)
end

function addchild(data, parent::FirmIdentifierNode)
    if !isdefined(parent, :children)
        parent.children = FirmIdentifierNode[]
    end
    push!(parent.children, FirmIdentifierNode(data, parent))
end

function AbstractTrees.children(node::FirmIdentifierNode)
    return node
    #isdefined(node, :children) && node.children
    # return node.children
end

AbstractTrees.printnode(io::IO, node::FirmIdentifierNode) = print(io, node.data)

function Base.iterate(node::FirmIdentifierNode, state::Int=1)
    isdefined(node, :children) && length(node.children) >= state && return (node.children[state], state+1)
    return nothing
end

Base.IteratorSize(::Type{FirmIdentifierNode}) = Base.SizeUnknown()
Base.eltype(::Type{FirmIdentifierNode}) = FirmIdentifierNode

Base.IteratorEltype(::Type{<:TreeIterator{FirmIdentifierNode}}) = Base.HasEltype()
AbstractTrees.parentlinks(::Type{FirmIdentifierNode}) = AbstractTrees.StoredParents()
AbstractTrees.siblinglinks(::Type{FirmIdentifierNode}) = AbstractTrees.StoredSiblings()

Base.parent(root::FirmIdentifierNode, node::FirmIdentifierNode) = isdefined(node, :parent) ? node.parent : nothing
Base.parent(node::FirmIdentifierNode) = isdefined(node, :parent) ? node.parent : nothing

Base.pairs(node::FirmIdentifierNode) = enumerate(node)

##

"""
These are the functions necessary to construct the tree
"""

function get_root(node::FirmIdentifierNode)
    if parent(node) !== nothing
        return get_root(parent(node))
    end
    return node
end

function existing_types(
    node::FirmIdentifierNode;
    out=Type{<:FirmIdentifier}[]
)

    push!(out, node.data)
    for x in node
        push!(out, x.data)
        existing_types(x; out)
    end

    out
end

function base_firm_identifiers_types(
    x::Type{<:FirmIdentifier};
    out=Type{<:FirmIdentifier}[]
)
    for y in subtypes(x)
        if isabstracttype(y)
            out=base_firm_identifiers_types(y; out)
        else
            push!(out, y)
        end
    end

    out |> unique
end

function build_tree_base(x::FirmIdentifierNode)
    for y in base_firm_identifiers_types(FirmIdentifier)

        if isabstracttype(y) || y == x.data || y ∈ existing_types(get_root(x))
            continue
        end
        if !isempty(methods(LinkTable, (Type{x.data}, Type{y})))
            addchild(y, x)
        end
    end
    if isdefined(x, :children)
        for y in x
            build_tree_base(y)
        end
    end
end

function build_tree(x::Type{<:FirmIdentifier})
    root = FirmIdentifierNode(x)
    build_tree_base(root)
    return root
end


text:
using Test
using BenchmarkTools

using LinearAlgebra
using Random
using SparseArrays
using StaticArrays
using FiniteDiff
using Rotations

using Dojo

@testset "Integrator"                    verbose=true begin include("integrator.jl") end
@testset "Minimal Coordinates"           verbose=true begin include("minimal.jl") end
@testset "Modified Rodrigues Parameters" verbose=true begin include("mrp.jl") end
@testset "Jacobian (Solution Matrix)"    verbose=true begin include("jacobian.jl") end
@testset "Data"                          verbose=true begin include("data.jl") end
@testset "Momentum"                      verbose=true begin include("momentum.jl") end
@testset "Energy"                        verbose=true begin include("energy.jl") end
@testset "Behavior"                      verbose=true begin include("behaviors.jl") end
@testset "Joint Limits"                  verbose=true begin include("joint_limits.jl") end
@testset "Impulse Map"                   verbose=true begin include("impulse_map.jl") end
@testset "Bodies"                        verbose=true begin include("bodies.jl") end
@testset "Mechanism (Miscellaneous)"     verbose=true begin include("mechanism.jl") end
@testset "Simulate"                      verbose=true begin include("simulate.jl") end
@testset "Visuals"                       verbose=true begin include("visuals.jl") end
@testset "Utilities"                     verbose=true begin include("utilities.jl") end
@testset "Mechanisms"                    verbose=true begin include("mechanisms.jl") end
@testset "Environments"                  verbose=true begin include("environments.jl") end


text:
"Module for tracing a ray."
module trace_simple

export trace_pixel

using ..model
using ..constants
using ..geometry

# Simply use midpoint formula and adaptive steps to measure tau

const TAU_THRESH = 8.0

"Adaptive stepper using Midpoint Method. https://en.wikipedia.org/wiki/Midpoint_method "
function integrate_tau(zstart::Real, v0::Real, pars::AbstractParameters, interp, gpre; h_tau::Real=0.1, nsteps::Int=1000, tau_start::Real=0.0, max_ds::Real=-5*AU)
    # Initialize an array of taus
    zps = Array{Float64}(nsteps)
    alpha_mps = Array{Float64}(nsteps) # alpha evaluated at the midpoints
    taus = Array{Float64}(nsteps)

    tot_intensity = 0.0

    zps[1] = zstart
    taus[1] = tau_start

    # Write in coordinate conversions for querying alpha and sfunc
    rcyl, z, vlos = get_coords(zstart, gpre)

    # Calculate Delta v at this position
    Deltav = v0 * 1.e5 - vlos # get_vlos(xprime, yprime, zprime, pars)

    # Look up RT quantities from nearest-neighbor interp.
    DV2, sfunc, Upsilon = interp(rcyl, z)

    # Evaluate alpha at current zstart position
    alpha = Upsilon * exp(-Deltav^2/DV2)
    alpha_mps[1] = alpha
    zp2 = 0.0

    println("Max ds: ", max_ds)

    for i=2:nsteps
        # Based upon current alpha value, calculate how much of a dz we would need to get dz * alpha = h_tau
        # Because these are negative numbers, the smaller step is actually the maximum (closer to 0)
        dzp = maximum([-h_tau / alpha_mps[i-1], max_ds])

        # Update current z position
        zps[i] = dzp + zps[i-1]

        # Midpoint step position
        zp2 = zps[i-1] + 0.5 * dzp

        # Write in coordinate conversions for querying alpha and sfunc
        rcyl, z, vlos = get_coords(zp2, gpre)

        # Calculate Delta v at this position
        Deltav = v0 * 1.e5 - vlos # get_vlos(xprime, yprime, zprime, pars)

        # Look up RT quantities from nearest-neighbor interp.
        # sfunc is evaluated at the midpoint
        DV2, sfunc, Upsilon = interp(rcyl, z)

        # Calculate alpha at new midpoint location
        alpha_mps[i] = Upsilon * exp(-Deltav^2/DV2)
        taus[i] = taus[i-1] - dzp * alpha_mps[i]

        # dtau = - alpha * dzp
        # tau_(n+1) = tau_n + dzp * alpha(zp_n + h/2)
        # since we have a simple ODE, explicit and implicit techniques are the same.

        # Calculate S at the new midpoint location
        # tot_intensity += -exp(-0.5(taus[i] + taus[i-1])) * sfunc * alpha_mps[i] * dzp # why is this the most accurate?
        # tot_intensity += -exp(-taus[i]) * sfunc * alpha_mps[i] * dzp
        # the reason it wasn't working before is because  sfunc is actually already evaluated at the midpoint.
        tot_intensity += (exp(-taus[i-1]) - exp(-taus[i])) * sfunc # this is also very accurate

        # Update trailing sfuncs
        # sfunc_ = sfunc

        # Query alpha at current z position
        if taus[i] > TAU_THRESH
            # total intensity
            return (zps[1:i], taus[1:i], tot_intensity)
        end

    end

    return (zps, taus, tot_intensity)
end

"v0 is the central velocity of the channel, relative to the disk systemic velocity."
function trace_pixel(xprime, yprime, v0, mol::Molecule, pars::AbstractParameters, interp)

    rmax = 700.0 * AU

    # Calculate DeltaVmax for this ray
    rcyl_min = abs(xprime)
    DeltaVmax = sqrt(2 * kB * temperature(rcyl_min, pars)/mol.mol_weight + (pars.ksi * 1e5)^2) * 1e-5 # km/s

    # Calculate the bounding positions for emission
    # If it's two, just trace it.
    # If it's four, break it up into two integrals.

    bounds = get_bounding_zps(xprime, yprime, pars, v0, DeltaVmax, rmax)
    if length(bounds) == 4
        z1start, z1end, z2start, z2end = bounds
    elseif length(bounds) == 2
        z1start, z1end = bounds
    end

    gpre = precalc_geo(xprime, yprime, pars)

    ans = integrate_tau(z1start, v0, pars, interp, gpre)

    return ans

    # if (tau_final < TAU_THRESH) & (length(bounds) == 4)
    #     #We didn't reach the threshold in the first bounding region, so continue in the second.
    #     tspan = (0, z2start - z2end)
    #
    #     # Start the ode to integrate the first path from the end of the second.
    #     u0 = Float64[tau_final; int_final]
    #     args = (z2start, v0, pars, interp, gpre)
    #     prob = ODEProblem(f, u0, tspan, args)
    #     sol = solve(prob, alg, callback=cb, reltol=rt, abstol=at, dense=false, save_everystep=false, dtmax=dtmax, dt=dt)
    #     tau_final, int_final = sol.u[end]
    #     return int_final
    #
    # else
    #     return int_final
    # end
end

end # module


text:
module BeliefUpdaters

using POMDPs
import POMDPs: Updater, update, initialize_belief, pdf, mode, updater, iterator
import Base: rand, rand!, mean, ==

using POMDPToolbox.ordered_states
using StatsBase

export
    VoidUpdater
include("void.jl")

export
    DiscreteBelief,
    DiscreteUpdater,
    uniform_belief,
    product                     # Remove because deprecated
include("discrete.jl")

export
    PreviousObservationUpdater,
    FastPreviousObservationUpdater,
    PrimedPreviousObservationUpdater

include("previous_observation.jl")

end


text:
array1 = [(2n+1)^2 for n in 1:5]
array2 = [sqrt(i) for i in array1]
println(typeof(1:5), "  ", typeof(array1), "  ", typeof(array2))
1:5, array1, array2

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# Observations
# observation(pomdp, a, sp)
function POMDPs.observation(pomdp::FireWorld, a::Array{Int64,1}, sp::FireState)
    neighbors = FireObs[]
    probabilities = Array{Float64,1}(undef,0)

    burning = sp.burning
    fuels = sp.fuels

    P_xy = fire_spread(pomdp.wind, sp)
    s_new = update_burn_probs(sp, P_xy)
    burn_probs_new = s_new.burn_probs

    # idea is
    # update burning in observations
    # only change the cells that were not acted upon
    # because if action was applied, fp and fn rates are 0 (eyes on fire)
    # otherwise, some delay happens and observation may not be the same

    # update the burn_probs map for actions taken to be 1.0 (so it won't be sensitive to threshold we pick)
    burn_probs_for_obs = deepcopy(burn_probs_new)
    for a_i in a
        burn_probs_for_obs[a_i] = 1.0
    end

    thresholds = collect(12:19)/20
    individual_prob = 1/length(thresholds)
    burns = []
    for threshold_k in thresholds
        burning_k = burn_probs_for_obs .> threshold_k
        # need to check if already in
        if burning_k in burns
            b_idx = findall(x->x==burning_k, burns)
            prob = probabilities[b_idx[1]]
            probabilities[b_idx[1]] = prob + individual_prob
        else
            push!(burns, burning_k)
            push!(probabilities, individual_prob)
            push!(neighbors, FireObs(burning_k))
        end
    end
    return SparseCat(neighbors, probabilities)
end


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.2
#   kernelspec:
#     display_name: Julia 1.5.0
#     language: julia
#     name: julia-1.5
# ---

# # fixed_arrays.jl

using StaticArrays

function unit(::Type{T}, i::Integer) where T<:StaticVector
    T(ntuple(Val(length(T))) do j
        ifelse(i==j, 1,0)
    end)
end

macro fixed_vector(name, parent)
    esc(quote
        """
            S: Size
            T: Type
        """
        struct $(name){S, T} <: $(parent){S, T}
                data::NTuples{S, T}
                function $(name){S, T}(x::NTuple{S, T}) where {S, T}
                   new{S, T}(x)
                end

                function $(name){S, T}(x::NTuple{S, Any}) where {S, T}
                   new{S, T}(StaticArrays.convert_ntuple(T, x))
                end
            end
        end
    )
end

@macroexpand @fixed_vector goma kyu


text:
function Multiply()
    return Multiply(())
end

function value(obj::Multiply, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "value", jdouble, (jdouble, jdouble), arg0, arg1)
end


text:
"""
    jpeg_encode(filename::AbstractString, img; kwargs...) -> Int
    jpeg_encode(io::IO, img; kwargs...) -> Int
    jpeg_encode(img; kwargs...) -> Vector{UInt8}

Encode 2D image `img` as JPEG byte sequences and write to given I/O stream or file. The
return value is number of bytes. If output is not specified, the encoded result is stored
in memory as return value.

# Parameters

- `transpose::Bool`: whether we need to permute the image's width and height dimension
  before encoding. The default value is `false`.
- `quality::Int`: Constructs JPEG quantization tables appropriate for the indicated quality
  setting. The quality value is expressed on the 0..100 scale recommended by IJG. The
  default value is `92`. Pass `quality=nothing` to let libjpeg-turbo dynamicly guess a
  value.

!!! info "Custom compression parameters"
    JPEG has a large number of compression parameters that determine how the image is
    encoded. Most applications don't need or want to know about all these parameters. For
    more detailed information and explaination, please refer to the "Compression parameter
    selection" in [1]. Unsupported custom parameters might cause Julia segmentation fault.

# Examples

```jldoctest
julia> using JpegTurbo, TestImages

julia> img = testimage("cameraman");

julia> jpeg_encode("out.jpg", img) # write to file
51396

julia> buf = jpeg_encode(img); length(buf) # directly write to memory
51396
```

# References

- [1] [libjpeg API Documentation (libjpeg.txt)](https://raw.githubusercontent.com/libjpeg-turbo/libjpeg-turbo/main/libjpeg.txt)
"""
function jpeg_encode(img::AbstractMatrix{T}; transpose=false, kwargs...) where T<:Union{Real, Colorant}
    # quantilized into 8bit sequences first
    CT = T <: Colorant ? n0f8(eltype(img)) : Gray{N0f8}
    AT = Array{CT, ndims(img)}
    clamp01nan!(img)
    # jpegturbo is a C library and assumes row-major memory order, thus `collect` the data into
    # contiguous memeory layout already makes a transpose.
    img = transpose ? convert(AT, img) : convert(AT, PermutedDimsArray(img, (2, 1)))

    return _encode(img; kwargs...)
end

function jpeg_encode(filename::AbstractString, img; kwargs...)
    open(filename, "w") do io
        jpeg_encode(io, img; kwargs...)
    end
end
# TODO(johnnychen94): further improve the performance via asynchronously IO and buffer reuse.
jpeg_encode(io::IO, img; kwargs...) = write(io, jpeg_encode(img; kwargs...))

function _encode(
    img::Matrix{<:Colorant};
    colorspace::Union{Nothing,Type} = nothing,
    # ImageMagick: "the default is to use the estimated quality of your input image if it can be determined, otherwise 92."
    quality::Union{Nothing,Int} = 92,
    arith_code::Union{Nothing,Bool} = nothing,
    optimize_coding::Union{Nothing,Bool} = nothing,
    smoothing_factor::Union{Nothing,Int} = nothing,
    write_JFIF_header::Union{Nothing,Bool} = nothing,
    JFIF_version::Union{Nothing,VersionNumber} = nothing,
    density_unit::Union{Nothing,Int} = nothing,
    X_density::Union{Nothing,Int} = nothing,
    Y_density::Union{Nothing,Int} = nothing,
    write_Adobe_marker::Union{Nothing,Bool} = nothing
)
    cinfo = LibJpeg.jpeg_compress_struct()
    cinfo_ref = Ref(cinfo)
    jerr = Ref{LibJpeg.jpeg_error_mgr}()
    cinfo.err = LibJpeg.jpeg_std_error(jerr)
    LibJpeg.jpeg_create_compress(cinfo_ref)

    # set input image information
    cinfo.image_width = size(img, 1)
    cinfo.image_height = size(img, 2)
    cinfo.input_components = jpeg_components(img)
    cinfo.in_color_space = jpeg_color_space(img)

    # set compression keywords
    # it's recommended to call `jpeg_set_defaults` first before setting custom parameters
    # as it's more likely to provide a working parameters and is more likely to be working
    # correctly in the future.
    LibJpeg.jpeg_set_defaults(cinfo_ref)
    isnothing(colorspace) || LibJpeg.jpeg_set_colorspace(cinfo_ref, jpeg_color_space(colorspace))
    isnothing(quality) || LibJpeg.jpeg_set_quality(cinfo_ref, quality, true)
    isnothing(arith_code) || (cinfo.arith_code = arith_code)
    isnothing(optimize_coding) || (cinfo.optimize_coding = optimize_coding)
    isnothing(smoothing_factor) || (cinfo.smoothing_factor = smoothing_factor)
    isnothing(write_JFIF_header) || (cinfo.write_JFIF_header = write_JFIF_header)
    if !isnothing(JFIF_version)
        cinfo.JFIF_major_version = UInt8(JFIF_version.major)
        cinfo.JFIF_minor_version = UInt8(JFIF_version.minor)
    end
    isnothing(density_unit) || (cinfo.density_unit = density_unit)
    isnothing(X_density) || (cinfo.X_density = X_density)
    isnothing(Y_density) || (cinfo.Y_density = Y_density)
    isnothing(write_Adobe_marker) || (cinfo.write_Adobe_marker = write_Adobe_marker)

    # set destination
    # TODO(johnnychen94): allow pre-allocated buffer
    bufsize = Ref{Culong}(0)
    buf_ptr = Ref{Ptr{UInt8}}(C_NULL)
    LibJpeg.jpeg_mem_dest(cinfo_ref, buf_ptr, bufsize)

    # compression stage
    LibJpeg.jpeg_start_compress(cinfo_ref, true)
    row_stride = size(img, 1) * jpeg_components(img)
    row_pointer = Ref{Ptr{UInt8}}(0)
    while (cinfo.next_scanline < cinfo.image_height)
        row_pointer[] = pointer(img) + cinfo.next_scanline * row_stride
        LibJpeg.jpeg_write_scanlines(cinfo_ref, row_pointer, 1);
    end
    LibJpeg.jpeg_finish_compress(cinfo_ref)
    LibJpeg.jpeg_destroy_compress(cinfo_ref)

    return unsafe_wrap(Array, buf_ptr[], bufsize[]; own=true)
end


text:
∑ = sum
using StatsFuns:logistic
sigmoid = logistic
using CSV:File
using DataFrames
using Printf:@sprintf

# DataTypes
using LinearAlgebra:Transpose
Numeric = Union{Int64,Float64}
NumericV = Union{Array{Int64,1},Array{Float64,1}}
NumericM = Union{Array{Int64,2},Array{Float64,2},Transpose{Float64,Array{Float64,2}},Transpose{Int64,Array{Int64,2}}}
arrayNumericV = Union{Array{Array{Int64,1},1},Array{Array{Float64,1},1}}
arrrayNumericArray = Union{Array{Array{Array{Float64,2},1},1},Array{Array{Array{Int64,2},1},1}}
TupleNumericM = Tuple{NumericM, NumericM};

function maskBiasUnits(array::NumericV)::NumericV
    mask = ones(size(array));
    mask[1,1] = 0;
    maskedArray = array.*mask;
end

function maskBiasUnits(array::NumericM)::NumericM
    mask = ones(size(array));
    mask[:,1] .= 0;
    maskedArray = array.*mask;
end

function powerElements(array::Union{NumericV,NumericM}, power::Int)::NumericM
    # function defined to be broadcasted to elements of a tuple
    array = array.^power
end

function getLabels(y::NumericV)::NumericV
    labels = collect(BitSet(y));
    sort!(labels)
end

function sortArrayByCol(array::Array{Any,2}, col::Int)
    array = DataFrame(array);
    sort!(array, col)
    array = Matrix(array)
end

function getLayersQ(all_θ)::Int
    θ_matricesQ = length(collect(all_θ));
    layersQ = θ_matricesQ+1;
end

function nnClasificationFeedForward(all_θ, X, y, labels, layersQ::Int=getLayersQ(all_θ))
    labelsQ = length(labels);
    hx_and_indices = reshape([],0,labelsQ+1);
    #zDict = Dict(i => Dict(j => Array{Float64,2}  for j in 1:layersQ) for i in 1:labelsQ);
    all_z = [];
    for lb in labels
        push!(all_z, [])
        ly = 1;
        lb_idx = findall(==(lb), y);
        a_lb_ly = X[lb_idx,:];
        ly += 1;
        while ly ≤ layersQ
            # add bias units to a_lb_ly
            a_lb_ly = [ones(size(a_lb_ly)[1],1) a_lb_ly];
            a_previous_ly = a_lb_ly;
            θ_previous_ly = all_θ[ly-1];
            z_lb_ly = a_previous_ly*θ_previous_ly';
            push!(all_z[lb],z_lb_ly);
            a_lb_ly = sigmoid.(z_lb_ly);
            ly += 1;
        end
        ly = 1;
        hx_lb = a_lb_ly;
        hx_and_indices_lb = hcat(hx_lb, lb_idx);
        hx_and_indices = vcat(hx_and_indices, hx_and_indices_lb);
    end
    hx_and_indices = sortArrayByCol(hx_and_indices, labelsQ+1);
    hx = hx_and_indices[:,1:labelsQ];
    hx = convert(Array{Float64,2}, hx)
    return hx, all_z
end

function nnCostClasification(all_θ::TupleNumericM, hx, y::NumericV, labels::NumericV, λ::Numeric)::Numeric
    J = 0;
    for lb in labels
        lb_idx = findall(==(lb), y)
        hx_lb = hx[lb_idx,:];
        y_lb = vec(zeros(size(labels))); y_lb[lb] = 1;
        J += ∑(-((log.(hx_lb))*y_lb)-(log.(1 .- hx_lb)*(1 .- y_lb)));
    end
    casesQ = size(hx)[1];
    J = (1/casesQ)*J;
    # add regularization cost
    J += (λ/(2*casesQ))* ∑(∑.(powerElements.(maskBiasUnits.(all_θ),2)));
end

function calcSigmoidDerivate(z::Numeric)::Numeric
    sigmoid(z)*(1-sigmoid(z));
end

function calcSigmoidDerivate(z::NumericM)::NumericM
    sigmoid.(z).*(1 .- sigmoid.(z));
end

function calcSigmoidDerivate(z::Array{Any,1})
    calcSigmoidDerivate.(z);
end

function calcNiceInitϵ(layerInputsQ::Int, layerOutputsQ::Int)::Numeric
    initϵ = √6/√(layerInputsQ+layerOutputsQ)
end

function genRandomInitialθ(layerInputsQ::Int, layerOutputsQ::Int, initϵ::Numeric=calcNiceInitϵ(layerInputsQ,layerOutputsQ))::NumericM
    initϵInt = floor(initϵ*100);
    randInitialθ = rand(-initϵInt:initϵInt,layerOutputsQ,layerInputsQ+1);
    randInitialθ = randInitialθ./100;
end

function nnClasificationBackPropagation(hx::NumericM, all_θ::TupleNumericM, all_SigmoidDerivates::arrrayNumericArray, layersQ::Int, y::NumericV, labels::NumericV)
    all_δ = [];
    labelsQ = length(labels);
    casesQ = size(hx)[1];
    ly = layersQ;
    lb_idx_v = collect(1:labelsQ);
    outputLayer = map(idx_lb -> (vec(hx[:,idx_lb])), lb_idx_v); #Matrix of subarrays where each of them correspond to each row of hx'
    δ_ly =  map(idx_lb -> outputLayer[idx_lb] .- y, lb_idx_v);
    push!(all_δ, δ_ly)
    ly -= 1;
    while ly ≥ 2
        δ_previous_ly = δ_ly;
        θ_ly = all_θ[ly][:,2:size(all_θ[ly])[2]];
        δ_ly = map(idx_lb -> (δ_previous_ly[idx_lb][findall(==(idx_lb), y)]*(θ_ly[idx_lb,:])').*all_SigmoidDerivates[idx_lb][ly-1], lb_idx_v);
        # note, there is one matrix less than the total number of layers in the all_SigmoidDerivates[idx_lb] array, hence the [ly-1] slicing
        push!(all_δ, δ_ly)
        ly -= 1;
    end
    reverse!(all_δ)
end

#parameters
# csv files -> DataFrames -> arrays
X =  Matrix{Float64}(File("./X.csv", header=false)
    |> DataFrame);
y = vec(Matrix(File("./y.csv", header=false) |> DataFrame));
θ1 = Matrix{Float64}(File("./theta1.csv", header=false)
    |> DataFrame);
θ2 = Matrix{Float64}(File("./theta2.csv", header=false) |> DataFrame);
λ = 1;

#implementation
all_θ = (θ1, θ2);
labels = getLabels(y);
layersQ = getLayersQ(all_θ);
hx, all_z = nnClasificationFeedForward(all_θ, X, y, labels, layersQ);
J = nnCostClasification(all_θ, hx, y, labels, λ);
all_SigmoidDerivates= calcSigmoidDerivate.(all_z);
all_δ = nnClasificationBackPropagation(hx, all_θ, all_SigmoidDerivates, layersQ, y, labels);

## consider next

# collect(Iterators.flatten(all_δ[1]))
# check that 125000/25 = 5000
    # 25 is the output nodes from layer 1

# collect(Iterators.flatten(all_δ[2]))
# check that 50000/10 = 5000
    # 10 is the output nodes from layer 2

# what is what?
    # how to slice / reshape?
        # check original output from all_δ[1] and all_δ[2]

max_stars_repo_name: luis-perez-one/courseStandfordML

max_issues_repo_name: luis-perez-one/courseStandfordML

max_forks_repo_name: luis-perez-one/courseStandfordML


text:
facts("Account type constructor works") do

    context("Good things work") do
        @fact 1 => 1
    end

    context("Bad things don't work") do
        @fact_throws Account("foo")
    end
end