text:
push!(LOAD_PATH, ".")

module Population

include("utility.jl")

using Individual
using Fitness

IndFitType{IndType <: Individual.AbstractIndividual} = Tuple{IndType, NTuple{N, Float64} where N}

abstract type AbstractPopulation{IndType <: Individual.AbstractIndividual} <: AbstractArray{Tuple{IndType, NTuple{N, Float64} where N}, 1} end

mutable struct PopulationType{IndType} <: AbstractPopulation{IndType}
  ind_args::Tuple
  pop::Vector{IndType}
  fitness::Fitness.AbstractFitness{N} where N
  pop_fit::Vector{NTuple{N, Float64} where N}
  fit_refresh::Vector{Bool}
  function PopulationType{IndType}(x...) where {IndType}
    args = build(IndType, x...)
    new(args[1], args[2], args[3], args[4], args[5])
  end
end

function build(ind_type::Type, ind_args::Tuple, f::Fitness.AbstractFitness{N}) where N
  return ind_args, Vector{ind_type}(), f, Vector{NTuple{N, Float64}}(), Vector{Bool}()
end

function getPopSize(this::AbstractPopulation)
  return length(this.pop)
end

function insertIndividual!(this::AbstractPopulation{IndType}, ind::IndType, fit::NTuple{N, Float64} = (NaN,)) where {IndType, N}
  if N != Fitness.getSize(this.fitness)
    fit = tuple(fill(NaN, Fitness.getSize(this.fitness))...)
  end
  push!(this.pop, ind)
  push!(this.pop_fit, fit)
  push!(this.fit_refresh, fit[1] === NaN)
end

function getIndividual(this::AbstractPopulation, pos::Int64)
  return (this.pop[pos], this.pop_fit[pos])
end

function refreshIndividual!(this::AbstractPopulation, pos::Int64)
  if this.fit_refresh[pos]
    this.fit_refresh[pos] = false
    this.pop_fit[pos] = this.fitness(this.pop[pos])
  end
end

function evalFitness!(this::AbstractPopulation)
  for i = 1 : length(this.pop)
    refreshIndividual!(this, i)
  end
end

function getIndArgs(this::AbstractPopulation)
  return this.ind_args
end

function getFitFunction(this::AbstractPopulation)
  return this.fitness
end

function populateRandom!(this::AbstractPopulation{IndType}, num_ind::Int64) where {IndType}
  for i = 1 : num_ind
    new_ind = IndType(this.ind_args...)
    Individual.generateRandom!(new_ind)
    insertIndividual!(this, new_ind)
  end
end

function Base.show(io::IO, this::AbstractPopulation)
  println("Population ", getPopSize(this), "\n")
  for i = 1 : getPopSize(this)
    print("Individual ", i, ": [")
    for gene in this.pop[i]
      print(" ", gene)
    end
    println(" ] -> ", this.pop_fit[i])
  end
  print("\n")
end

function Base.size(this::AbstractPopulation)
  return (getPopSize(this), )
end

function Base.similar(this::AbstractPopulation, ::Type{IndFitType}, sz::Int64)
  return Vector{IndFitType}(undef, sz)
end

function Base.IndexStyle(::Type{<:AbstractPopulation})
  return IndexLinear()
end

function Base.getindex(this::AbstractPopulation, pos::Int64)
  return (this.pop[pos], this.pop_fit[pos])
end

function Base.setindex!(this::AbstractPopulation, value::IndFitType, pos::Int64)
  this.pop[pos], this.pop_fit[pos] = value
end

function Base.firstindex(this::AbstractPopulation)
  return 1
end

function Base.lastindex(this::AbstractPopulation)
  return getPopSize(this)
end

function Base.:<(x::IndFitType{IndType}, y::IndFitType{IndType}) where {IndType}
  fit_size = length(x[2])
  fit_size > 1 || return x[2][1] < y[2][1]
  x[2] != y[2] || return false
  for i in eachindex(x[2])
    x[2][i] <= y[2][i] || return false
  end
  return true
end

function Base.isless(x::IndFitType, y::IndFitType)
  return x < y
end

function Base.maximum(this::AbstractPopulation)
  res = this[1]
  for ind in @view this[2:end]
    res = max(ind, res)
  end
  return res
end

end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


text:
using FreeTypeAbstraction, Colors, ColorVectorSpace, GeometryBasics
using Test
using FreeTypeAbstraction: boundingbox, Vec, glyph_rects, get_extent, FTFont, kerning, glyph_ink_size
using FreeType

face = FreeTypeAbstraction.findfont("hack"; additional_fonts=@__DIR__)

@test repr(face) == "FTFont (family = Hack, style = Regular)"

bb = boundingbox("asdasd", face, 64)
@test round.(Int, minimum(bb)) == Vec(4, -1)
@test round.(Int, widths(bb)) == Vec2(221, 50)

FA = FreeTypeAbstraction

FA.set_pixelsize(face, 64) # should be the default
img, extent = renderface(face, 'C', 64)
@test size(img) == (30, 49)
@test typeof(img) == Array{UInt8,2}

a = renderstring!(zeros(UInt8, 20, 100), "helgo", face, 10, 10, 10)

@test any(a[3:12, :] .!= 0)
@test all(a[vcat(1:2, 13:20), :] .== 0)
@test any(a[:, 11:40] .!= 0)
@test all(a[:, vcat(1:10, 41:100)] .== 0)
a = renderstring!(zeros(UInt8, 20, 100), "helgo", face, 10, 15, 70)
@test any(a[8:17, :] .!= 0)
@test all(a[vcat(1:7, 18:20), :] .== 0)
@test any(a[:, 71:100] .!= 0)
@test all(a[:, 1:70] .== 0)

a = renderstring!(
    zeros(UInt8, 20, 100),
    "helgo",
    face,
    10,
    10,
    50,
    valign = :vtop,
)

@test all(a[1:10, :] .== 0)
@test any(a[11:20, :] .!= 0)
a = renderstring!(
    zeros(UInt8, 20, 100),
    "helgo",
    face,
    10,
    10,
    50,
    valign = :vcenter,
)
@test all(a[vcat(1:5, 16:end), :] .== 0)
@test any(a[6:15, :] .!= 0)
a = renderstring!(
    zeros(UInt8, 20, 100),
    "helgo",
    face,
    10,
    10,
    50,
    valign = :vbaseline,
)
@test all(a[vcat(1:2, 13:end), :] .== 0)
@test any(a[3:12, :] .!= 0)
a = renderstring!(
    zeros(UInt8, 20, 100),
    "helgo",
    face,
    10,
    10,
    50,
    valign = :vbottom,
)
@test any(a[1:10, :] .!= 0)
@test all(a[11:20, :] .== 0)
a = renderstring!(
    zeros(UInt8, 20, 100),
    "helgo",
    face,
    10,
    10,
    50,
    halign = :hleft,
)
@test all(a[:, 1:50] .== 0)
@test any(a[:, 51:100] .!= 0)
a = renderstring!(
    zeros(UInt8, 20, 100),
    "helgo",
    face,
    10,
    10,
    50,
    halign = :hcenter,
)
@test all(a[:, vcat(1:35, 66:end)] .== 0)
@test any(a[:, 36:65] .!= 0)
a = renderstring!(
    zeros(UInt8, 20, 100),
    "helgo",
    face,
    10,
    10,
    50,
    halign = :hright,
)
@test any(a[:, 1:50] .!= 0)
@test all(a[:, 51:100] .== 0)

a = renderstring!(
    zeros(UInt8, 20, 100),
    "helgo",
    face,
    10,
    10,
    50,
    fcolor = 0x80,
)
@test maximum(a) <= 0x80
a = renderstring!(
    zeros(UInt8, 20, 100),
    "helgo",
    face,
    10,
    10,
    50,
    fcolor = 0x80,
    bcolor = 0x40,
)
@test any(a .== 0x40)
a = renderstring!(
    fill(0x01, 20, 100),
    "helgo",
    face,
    10,
    10,
    50,
    bcolor = nothing,
)
@test !any(a .== 0x00)

a = renderstring!(zeros(Float32, 20, 100), "helgo", face, 10, 10, 50)
@test maximum(a) <= 1.0
a = renderstring!(zeros(Float64, 20, 100), "helgo", face, 10, 10, 50)
@test maximum(a) <= 1.0

renderstring!(zeros(Gray, 20, 100), "helgo", face, 10, 10, 50)
renderstring!(
    zeros(Gray{Float64}, 20, 100),
    "helgo",
    face,
    10,
    10,
    50,
    fcolor = Gray(0.5),
)

renderstring!(
    zeros(UInt8, 20, 100),
    "helgo",
    face,
    10,
    0,
    0,
    halign = :hcenter,
    valign = :vcenter,
)
renderstring!(zeros(UInt8, 20, 100), "helgo", face, 10, 25, 80)

# Find fonts
# these fonts should be available on all platforms:

# debug travis... does it even have fonts?
fontpaths = FreeTypeAbstraction.fontpaths()

isempty(fontpaths) && println("OS doesn't have any font folder")

if Sys.islinux()
    fonts = ["dejavu sans"]
    # apple on gh-actions doesn't seem to have any fonts...
elseif Sys.isapple()
    fonts = []
else # windows have some more fonts installed per default
    fonts = [
        "Times New Roman",
        "Arial",
        "Comic Sans MS",
        "Impact",
        "Tahoma",
        "Trebuchet MS",
        "Verdana",
        "Courier New",
    ]
end

@testset "finding fonts" begin

    for font in fonts
        @testset "finding $font" begin
            @test findfont(font) != nothing
        end
    end
    @testset "find in additional dir" begin
        @test findfont("Hack") == nothing
        @test findfont("Hack", additional_fonts = @__DIR__) != nothing
    end
end

@testset "loading lots of fonts" begin
    for i = 1:10
        for font in fonts
            @time findfont(font)
        end
    end
    @test "No Error" == "No Error"
end

ext: jl


text:
using GLVisualize, GeometryTypes, GLAbstraction, Colors

if !isdefined(:runtests)
	window = glscreen()
end
const static_example = true

# one could also create a GLNormalAttributeMesh manually:
sphere = Sphere{Float32}(Point3f0(0), 2f0)
# decompose decomposes a mesh or geometry into the primitive of the first argument
vertices = decompose(Point3f0, sphere, 50)
faces = decompose(GLTriangle, sphere, 50)
# create some colors for each vertex, must be Float32
colors = RGBA{Float32}[RGBA{Float32}(rand(), rand(), rand(), 1.) for i=1:length(vertices)]

sphere_mesh = GLNormalVertexcolorMesh(
    vertices=vertices, faces=faces,
    color=colors
)

_view(visualize(sphere_mesh), window)

if !isdefined(:runtests)
	renderloop(window)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""
    ObjsUserProfile(;
        always_active=nothing,
        api_app_id=nothing,
        avatar_hash=nothing,
        bot_id=nothing,
        display_name=nothing,
        display_name_normalized=nothing,
        is_app_user=nothing,
        is_custom_image=nothing,
        last_avatar_image_hash=nothing,
        memberships_count=nothing,
        phone=nothing,
        pronouns=nothing,
        real_name=nothing,
        real_name_normalized=nothing,
        skype=nothing,
        status_default_emoji=nothing,
        status_default_text=nothing,
        status_emoji=nothing,
        status_expiration=nothing,
        status_text=nothing,
        team=nothing,
        title=nothing,
        updated=nothing,
        user_id=nothing,
    )

    - always_active::Bool
    - api_app_id::DefsOptionalAppId
    - avatar_hash::String
    - bot_id::DefsBotId
    - display_name::String
    - display_name_normalized::String
    - is_app_user::Bool
    - is_custom_image::Bool
    - last_avatar_image_hash::String
    - memberships_count::Int32
    - phone::String
    - pronouns::String
    - real_name::String
    - real_name_normalized::String
    - skype::String
    - status_default_emoji::String
    - status_default_text::String
    - status_emoji::String
    - status_expiration::Int32
    - status_text::String
    - team::DefsWorkspaceId
    - title::String
    - updated::Int32
    - user_id::String
"""
mutable struct ObjsUserProfile <: SwaggerModel
    always_active::Any # spec type: Union{ Nothing, Bool } # spec name: always_active
    api_app_id::Any # spec type: Union{ Nothing, DefsOptionalAppId } # spec name: api_app_id
    avatar_hash::Any # spec type: Union{ Nothing, String } # spec name: avatar_hash
    bot_id::Any # spec type: Union{ Nothing, DefsBotId } # spec name: bot_id
    display_name::Any # spec type: Union{ Nothing, String } # spec name: display_name
    display_name_normalized::Any # spec type: Union{ Nothing, String } # spec name: display_name_normalized
    is_app_user::Any # spec type: Union{ Nothing, Bool } # spec name: is_app_user
    is_custom_image::Any # spec type: Union{ Nothing, Bool } # spec name: is_custom_image
    last_avatar_image_hash::Any # spec type: Union{ Nothing, String } # spec name: last_avatar_image_hash
    memberships_count::Any # spec type: Union{ Nothing, Int32 } # spec name: memberships_count
    phone::Any # spec type: Union{ Nothing, String } # spec name: phone
    pronouns::Any # spec type: Union{ Nothing, String } # spec name: pronouns
    real_name::Any # spec type: Union{ Nothing, String } # spec name: real_name
    real_name_normalized::Any # spec type: Union{ Nothing, String } # spec name: real_name_normalized
    skype::Any # spec type: Union{ Nothing, String } # spec name: skype
    status_default_emoji::Any # spec type: Union{ Nothing, String } # spec name: status_default_emoji
    status_default_text::Any # spec type: Union{ Nothing, String } # spec name: status_default_text
    status_emoji::Any # spec type: Union{ Nothing, String } # spec name: status_emoji
    status_expiration::Any # spec type: Union{ Nothing, Int32 } # spec name: status_expiration
    status_text::Any # spec type: Union{ Nothing, String } # spec name: status_text
    team::Any # spec type: Union{ Nothing, DefsWorkspaceId } # spec name: team
    title::Any # spec type: Union{ Nothing, String } # spec name: title
    updated::Any # spec type: Union{ Nothing, Int32 } # spec name: updated
    user_id::Any # spec type: Union{ Nothing, String } # spec name: user_id

    function ObjsUserProfile(;always_active=nothing, api_app_id=nothing, avatar_hash=nothing, bot_id=nothing, display_name=nothing, display_name_normalized=nothing, is_app_user=nothing, is_custom_image=nothing, last_avatar_image_hash=nothing, memberships_count=nothing, phone=nothing, pronouns=nothing, real_name=nothing, real_name_normalized=nothing, skype=nothing, status_default_emoji=nothing, status_default_text=nothing, status_emoji=nothing, status_expiration=nothing, status_text=nothing, team=nothing, title=nothing, updated=nothing, user_id=nothing)
        o = new()
        validate_property(ObjsUserProfile, Symbol("always_active"), always_active)
        setfield!(o, Symbol("always_active"), always_active)
        validate_property(ObjsUserProfile, Symbol("api_app_id"), api_app_id)
        setfield!(o, Symbol("api_app_id"), api_app_id)
        validate_property(ObjsUserProfile, Symbol("avatar_hash"), avatar_hash)
        setfield!(o, Symbol("avatar_hash"), avatar_hash)
        validate_property(ObjsUserProfile, Symbol("bot_id"), bot_id)
        setfield!(o, Symbol("bot_id"), bot_id)
        validate_property(ObjsUserProfile, Symbol("display_name"), display_name)
        setfield!(o, Symbol("display_name"), display_name)
        validate_property(ObjsUserProfile, Symbol("display_name_normalized"), display_name_normalized)
        setfield!(o, Symbol("display_name_normalized"), display_name_normalized)
        validate_property(ObjsUserProfile, Symbol("is_app_user"), is_app_user)
        setfield!(o, Symbol("is_app_user"), is_app_user)
        validate_property(ObjsUserProfile, Symbol("is_custom_image"), is_custom_image)
        setfield!(o, Symbol("is_custom_image"), is_custom_image)
        validate_property(ObjsUserProfile, Symbol("last_avatar_image_hash"), last_avatar_image_hash)
        setfield!(o, Symbol("last_avatar_image_hash"), last_avatar_image_hash)
        validate_property(ObjsUserProfile, Symbol("memberships_count"), memberships_count)
        setfield!(o, Symbol("memberships_count"), memberships_count)
        validate_property(ObjsUserProfile, Symbol("phone"), phone)
        setfield!(o, Symbol("phone"), phone)
        validate_property(ObjsUserProfile, Symbol("pronouns"), pronouns)
        setfield!(o, Symbol("pronouns"), pronouns)
        validate_property(ObjsUserProfile, Symbol("real_name"), real_name)
        setfield!(o, Symbol("real_name"), real_name)
        validate_property(ObjsUserProfile, Symbol("real_name_normalized"), real_name_normalized)
        setfield!(o, Symbol("real_name_normalized"), real_name_normalized)
        validate_property(ObjsUserProfile, Symbol("skype"), skype)
        setfield!(o, Symbol("skype"), skype)
        validate_property(ObjsUserProfile, Symbol("status_default_emoji"), status_default_emoji)
        setfield!(o, Symbol("status_default_emoji"), status_default_emoji)
        validate_property(ObjsUserProfile, Symbol("status_default_text"), status_default_text)
        setfield!(o, Symbol("status_default_text"), status_default_text)
        validate_property(ObjsUserProfile, Symbol("status_emoji"), status_emoji)
        setfield!(o, Symbol("status_emoji"), status_emoji)
        validate_property(ObjsUserProfile, Symbol("status_expiration"), status_expiration)
        setfield!(o, Symbol("status_expiration"), status_expiration)
        validate_property(ObjsUserProfile, Symbol("status_text"), status_text)
        setfield!(o, Symbol("status_text"), status_text)
        validate_property(ObjsUserProfile, Symbol("team"), team)
        setfield!(o, Symbol("team"), team)
        validate_property(ObjsUserProfile, Symbol("title"), title)
        setfield!(o, Symbol("title"), title)
        validate_property(ObjsUserProfile, Symbol("updated"), updated)
        setfield!(o, Symbol("updated"), updated)
        validate_property(ObjsUserProfile, Symbol("user_id"), user_id)
        setfield!(o, Symbol("user_id"), user_id)
        o
    end
end # type ObjsUserProfile

const _property_map_ObjsUserProfile = Dict{Symbol,Symbol}(Symbol("always_active")=>Symbol("always_active"), Symbol("api_app_id")=>Symbol("api_app_id"), Symbol("avatar_hash")=>Symbol("avatar_hash"), Symbol("bot_id")=>Symbol("bot_id"), Symbol("display_name")=>Symbol("display_name"), Symbol("display_name_normalized")=>Symbol("display_name_normalized"), Symbol("is_app_user")=>Symbol("is_app_user"), Symbol("is_custom_image")=>Symbol("is_custom_image"), Symbol("last_avatar_image_hash")=>Symbol("last_avatar_image_hash"), Symbol("memberships_count")=>Symbol("memberships_count"), Symbol("phone")=>Symbol("phone"), Symbol("pronouns")=>Symbol("pronouns"), Symbol("real_name")=>Symbol("real_name"), Symbol("real_name_normalized")=>Symbol("real_name_normalized"), Symbol("skype")=>Symbol("skype"), Symbol("status_default_emoji")=>Symbol("status_default_emoji"), Symbol("status_default_text")=>Symbol("status_default_text"), Symbol("status_emoji")=>Symbol("status_emoji"), Symbol("status_expiration")=>Symbol("status_expiration"), Symbol("status_text")=>Symbol("status_text"), Symbol("team")=>Symbol("team"), Symbol("title")=>Symbol("title"), Symbol("updated")=>Symbol("updated"), Symbol("user_id")=>Symbol("user_id"))
const _property_types_ObjsUserProfile = Dict{Symbol,String}(Symbol("always_active")=>"Bool", Symbol("api_app_id")=>"DefsOptionalAppId", Symbol("avatar_hash")=>"String", Symbol("bot_id")=>"DefsBotId", Symbol("display_name")=>"String", Symbol("display_name_normalized")=>"String", Symbol("is_app_user")=>"Bool", Symbol("is_custom_image")=>"Bool", Symbol("last_avatar_image_hash")=>"String", Symbol("memberships_count")=>"Int32", Symbol("phone")=>"String", Symbol("pronouns")=>"String", Symbol("real_name")=>"String", Symbol("real_name_normalized")=>"String", Symbol("skype")=>"String", Symbol("status_default_emoji")=>"String", Symbol("status_default_text")=>"String", Symbol("status_emoji")=>"String", Symbol("status_expiration")=>"Int32", Symbol("status_text")=>"String", Symbol("team")=>"DefsWorkspaceId", Symbol("title")=>"String", Symbol("updated")=>"Int32", Symbol("user_id")=>"String")
Base.propertynames(::Type{ ObjsUserProfile }) = collect(keys(_property_map_ObjsUserProfile))
Swagger.property_type(::Type{ ObjsUserProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsUserProfile[name]))}
Swagger.field_name(::Type{ ObjsUserProfile }, property_name::Symbol) =  _property_map_ObjsUserProfile[property_name]

function check_required(o::ObjsUserProfile)
    (getproperty(o, Symbol("avatar_hash")) === nothing) && (return false)
    (getproperty(o, Symbol("display_name")) === nothing) && (return false)
    (getproperty(o, Symbol("display_name_normalized")) === nothing) && (return false)
    (getproperty(o, Symbol("phone")) === nothing) && (return false)
    (getproperty(o, Symbol("real_name")) === nothing) && (return false)
    (getproperty(o, Symbol("real_name_normalized")) === nothing) && (return false)
    (getproperty(o, Symbol("skype")) === nothing) && (return false)
    (getproperty(o, Symbol("status_emoji")) === nothing) && (return false)
    (getproperty(o, Symbol("status_text")) === nothing) && (return false)
    (getproperty(o, Symbol("title")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsUserProfile }, name::Symbol, val)
end


text:
using NearestNeighborDescent
using BenchmarkTools
using Distances
using MLDatasets
using PkgBenchmark

include("benchutils.jl")

const SUITE = BenchmarkGroup()

# graph construction time
SUITE["graph"] = BenchmarkGroup(["graph", "construction"])

# randomly generated dataset benchmarks
rand_vecs = [2. .* rand(100) .- 1 for _ in 1:5000]
rand_mats = 2. .* rand(100, 5000) .- 1
ham_vecs = [rand([0, 1], 100) for _ in 1:5000]
ham_mats = rand([0, 1], 100, 5000)

SUITE["graph"]["random"] = BenchmarkGroup()
SUITE["graph"]["random"][:vectors, :euclidean] = @benchmarkable DescentGraph($rand_vecs, 10, Euclidean())
SUITE["graph"]["random"][:matrices, :euclidean] = @benchmarkable DescentGraph($rand_mats, 10, Euclidean())
SUITE["graph"]["random"][:vectors, :cosine] = @benchmarkable DescentGraph($rand_vecs, 10, CosineDist())
SUITE["graph"]["random"][:matrices, :cosine] = @benchmarkable DescentGraph($rand_mats, 10, CosineDist())
SUITE["graph"]["random"][:vectors, :hamming] = @benchmarkable DescentGraph($ham_vecs, 10, Hamming())
SUITE["graph"]["random"][:matrices, :hamming] = @benchmarkable DescentGraph($ham_mats, 10, Hamming())

# real-world dataset benchmarks
SUITE["graph"]["real"] = BenchmarkGroup()
FMNIST_data = FashionMNIST.convert2features(FashionMNIST.traintensor(Float64))
FMNIST_queries = FashionMNIST.convert2features(FashionMNIST.testtensor(Float64))

MNIST_data = MNIST.convert2features(MNIST.traintensor(Float64))
MNIST_queries = MNIST.convert2features(MNIST.testtensor(Float64))

SUITE["graph"]["real"][:mnist] = @benchmarkable DescentGraph($MNIST_data, 10, Euclidean())
SUITE["graph"]["real"][:fmnist] = @benchmarkable DescentGraph($FMNIST_data, 5, Euclidean())

# queries per second
SUITE["query"] = BenchmarkGroup(["query", "search"])
rand_query_vecs = [2. .* rand(100) .- 1 for _ in 1:1000]
rand_query_mats = 2. .* rand(100, 1000) .- 1
ham_query_vecs = [rand([0, 1], 100) for _ in 1:1000]
ham_query_mats = rand([0, 1], 100, 1000)

SUITE["query"]["random"] = BenchmarkGroup()
SUITE["query"]["random"][:vectors, :euclidean] = @benchmarkable search(graph, $rand_query_vecs, 10, 10) setup=(graph=DescentGraph($rand_vecs, 10, Euclidean()))
SUITE["query"]["random"][:matrices, :euclidean] = @benchmarkable search(graph, $rand_query_mats, 10, 10) setup=(graph=DescentGraph($rand_mats, 10, Euclidean()))
SUITE["query"]["random"][:vectors, :cosine] = @benchmarkable search(graph, $rand_query_vecs, 10, 10) setup=(graph=DescentGraph($rand_vecs, 10, CosineDist()))
SUITE["query"]["random"][:matrices, :cosine] = @benchmarkable search(graph, $rand_query_mats, 10, 10) setup=(graph=DescentGraph($rand_mats, 10, CosineDist()))
SUITE["query"]["random"][:vectors, :hamming] = @benchmarkable search(graph, $ham_query_vecs, 10, 10) setup=(graph=DescentGraph($ham_vecs, 10, Hamming()))
SUITE["query"]["random"][:matrices, :hamming] = @benchmarkable search(graph, $ham_query_mats, 10, 10) setup=(graph=DescentGraph($ham_mats, 10, Hamming()))

# real-world dataset query benchmarks
SUITE["query"]["real"] = BenchmarkGroup()

SUITE["query"]["real"][:mnist] = @benchmarkable search(graph, $MNIST_queries, 10, 10) setup=(graph=DescentGraph($MNIST_data, 10, Euclidean()))
SUITE["query"]["real"][:fmnist] = @benchmarkable search(graph, $FMNIST_queries, 5, 10) setup=(graph=DescentGraph($FMNIST_data, 5, Euclidean()))


text:
@inline function Phi_se(Cell,s,z,Def,ϕ_tf,D,res0)
"""
   Solid-Electrolyte Potential Transfer Function

   Phi_se(Cell,s,z,Def)

"""

if Def == "Pos"
   Electrode = Cell.Pos #Electrode Length
else
   Electrode = Cell.Neg #Electrode Length
end

κ_eff = Cell.Const.κ*Electrode.ϵ_e^Electrode.κ_brug #Effective Electrolyte Conductivity
σ_eff = Electrode.σ*Electrode.ϵ_s^Electrode.σ_brug #Effective Electrode Conductivity

#Defining SOC
θ = Cell.Const.SOC * (Electrode.θ_100-Electrode.θ_0) + Electrode.θ_0

#Prepare for j0
ce0 = Cell.Const.ce0
cs_max = Electrode.cs_max
cs0 = cs_max * θ
α = Electrode.α

#Current Flux Density
if Cell.Const.CellTyp == "Doyle_94"
   κ = Electrode.k_norm/Electrode.cs_max/ce0^(1-α)
   j0 = κ*(ce0*(cs_max-cs0))^(1-α)*cs0^α
else
   j0 = Electrode.k_norm*(Cell.Const.ce0*cs0*(Electrode.cs_max-cs0))^(1-Electrode.α)
end

#Resistance
Rtot = R*Cell.Const.T/(j0*F^2) + Electrode.RFilm
#Rtot = R*Cell.Const.T/(j0*Cell.Const.CC_A*F) + Electrode.RFilm

∂Uocp_elc = Cell.Const.∂Uocp(Def,θ)/cs_max #Open Circuit Potential Partial
res0 .= @. -3*∂Uocp_elc/(Electrode.as*F*Electrode.L*Cell.Const.CC_A*Electrode.Rs) # residual for pole removal

ν = @. Electrode.L*sqrt((Electrode.as/σ_eff+Electrode.as/κ_eff)/(Rtot+∂Uocp_elc*(Electrode.Rs/(F*Electrode.Ds))*(tanh(Electrode.β)/(tanh(Electrode.β)-Electrode.β)))) #Condensing Variable - eq. 4.13
ν_∞ = @. Electrode.L*sqrt(Electrode.as*((1/κ_eff)+(1/σ_eff))/(Rtot))

ϕ_tf .= @. Electrode.L/(Cell.Const.CC_A*ν*sinh(ν))*((1/κ_eff)*cosh(ν*z)+(1/σ_eff)*cosh(ν*(z-1)))-res0/s #Transfer Function - eq. 4.14

zero_tf = @. (6*(5*Electrode.Ds*F*Rtot-∂Uocp_elc*Electrode.Rs)*σ_eff)/(30*Cell.Const.CC_A*Electrode.as*Electrode.Ds*F*σ_eff*Electrode.L) + (5*Electrode.as*Electrode.Ds*F*Electrode.L^2*(σ_eff*(-1+3*z^2)+κ_eff*(2-6*z+3*z^2)))/(30*Cell.Const.CC_A*Electrode.as*Electrode.Ds*F*σ_eff*κ_eff*Electrode.L)
D .= @. Electrode.L/(Cell.Const.CC_A*ν_∞*sinh(ν_∞))*((1/κ_eff)*cosh(ν_∞*z)+(1/σ_eff)*cosh(ν_∞*(z-1))) # Contribution to D as G->∞
ϕ_tf[:,findall(s.==0)] .= zero_tf[:,findall(s.==0)]
res0 .= zeros(length(z))

if Def == "Pos" #Double check this implementation
   ϕ_tf .= -ϕ_tf
   D .= -D
end

end

max_stars_repo_name: BradyPlanden/LiBRA

max_forks_repo_name: BradyPlanden/LiBRA


text:
import Base: @propagate_inbounds, @inbounds # compiler inlines a function, while retaining the caller's inbounds context
import BenchmarkTools: @btime, @time, @benchmark

#euclid fans have a right to throw tomato at me
@propagate_inbounds sumInterval(a, b) = abs(a + b) +1> # Range (min … max):  0.001 ns … 0.100 ns

@propagate_inbounds function sumInterval(a, b) #sumInterval subtracts 1 (complies with julia logic)
    cond = abs(a + b) - 1
    @inbounds if cond >= 1
        return cond
    end
end

"""isEven1 experimential Performance: always below isEven2
edit: there's nothing called as `always`
"""
@propagate_inbounds function isEven1(a = 1, b = 10) #3 if-statements
    try
        @inbounds if a >= 0 && b >= 0
            distance = sumInterval(a, b) #1. distance
            _isPositive = isPositive(distance) #2. positive distance
            @inbounds if _isPositive #positive?
                isEven = distance % 2 == 0 ? true : false
                @inbounds if isEven           # even?
                    return true
                elseif !isEven # || isPositive
                    return false
                else
                    println(UnexpMsg)
                end
            else
                throw(error("Positive number error"))
            end
            #end
        end
    catch notPositiveInputError
        #Exception
        @error positiveMsg exception = (notPositiveInputError, catch_backtrace())
        println(positiveMsg) #positive arguments error
    end
end
isEven1(-1, 10)
isEven1(1, 10)

#=
@inbounds if lenV < lenA # first assumption
    @inbounds a[i:lenV] = v[i:lenV]
    #elseif lenV > lenA #v can't be larger than original vector a (throws an outofBoundsError )
else
    throw(error("Positive number error"))
    println(UnexpMsg)
end
#return a
#end =#
 #computeRange!(a, α, β) # # TODO: middle(arr, a, b) #or # TODO: middle(a, b, arr)
@propagate_inbounds function replaceVector2(v = [1, 2], a = [2, 3, 4, 5]; i = 1)

"""divideConquer """

#  mit 6.006  #unused #orphaned #to-be-deleted
#choosing a,b is choosing arr , this is a mean, not a goal
# divideConquer(a, range[1], length(range[1])) #check
@propagate_inbounds function divideConquer(a = [1, 2, 3, 4]; α = arr[1], β = length(a) - 1)  ## arr is an initial value, while the goal would be , 1  length(other size of bound to walk on  ) a code that never walks..
    count = 1
    mid = middle(α, β) # first Domino #this presumes middle returns 1 middle # immature
    #check valid Interval
    cond = sumInterval(mid, count) > 0 ? true : false #= mid - count >= 0=#

    q = []

    if cond
        count += 1
        if a[mid] < a[mid-1] # only look at left half 1 . . . n/2 − − − 1 to look for peak
            pushfirst!(q, goleft(arr, a, mid))

        elseif a[mid] < a[mid+1] # only look at right half n/2 + 1 . . . n to look for peak
            pushfirst!(q, goright(a))

        else # a[mid] =?= a[mid+1]
            pushfirst!(q, mid) # the peak!
        end
    end
end
#no errors

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

cond = isEven(1, 3)
res = -1

#---test run tryout

#1,3 total range
# 1+3/2 = 2
#yeilds 1:2 , 2:3
#check dist

abs(2 - 3) == 1 # bounds are apart by 1 #sweet!

# final comparison
#doCompare(2,3,arr)
#compare current ranges with their value , do comparison when required
a[1]

#=Plain comparison - flip if lower index has a higher value, otherwise return
given bound values of indicies, find their corresponding value that they weigh,
then compare them together, as well

```input:
st: start bound
ed: ending bound
a: corresponding vector array to traverse
```

```output:
```

=#
#inputs, outputs

"""
```input:
```

```output:
```
"""
compareVector()

#--- q Zone /
a = buildAboveSoBelow(below, above, arr)
#= [1, 2]
 [2, 3]
 [3, 4, 5]=#
l = length(q)
count = 1 #starting from 1
#=
0 -> 3 (3)
1-> 4 (3)
this correctly counts 1 -> 3 (sumRanges 3+1 = 4 but max is 3 )
=#

l = length(q[2]) # 2
array = []; # [firstvalue]
count = 1
temp = deepcopy(q[2])

typeof(q[2][1]) #Vector #int64

res = q[2] # pop!(q[2])
typeof(q[2])

#---repeat
for i = 1:l #enumerate(q) # in 0:l#4 Always!    #enumerate(q) also works fine as Base.Iterator.Enumerate
    #count would be 3 = length(q) A lways
    #count<=l ? count += 1 : break; # count ==3 max
    #array[i] =i
    push!(arr, i)
    item = Base.popfirst!(q[2]) #pop!(collection) -> item # Remove the last item in collection and return it.
    #i<=l  ? count+= 1 : break;
    #count+=1 # counts 4 ( size = length(q) + 1 ) #either need to count -1

end # will count 4

#array arr
count -= 1 # learnd heuristic # length adjustment [of the Interval ] # q. why not measure length(v) directly (after its adjustments)?
#count = 0->3 0 index or 1 ->4 1-index

#q = buildAboveSoBelow(a, below, above, b) # instead: make adjustments on the array arr directtly

#--- end of q
#--- end q Zone

#=
"""a valid middle - Classic #old #depreciated ,why it uses q """
function middle(α, β) #
    condition = isEven(α, β) #  not | b-a
    #q = []
    if condition == true
        #return true #a.s. #eucledian Distance divided by 2 returing a whole integer
        check = Int(sumInterval(α, β) // 2) # | b + a | // 2 isa Integer #sumInterval -to-> ϟ#5
        #middleExtraction(condition, check) # Here we didn't get anything ! <------------- # check not defined here
        #return condition, check
        #push!(q, check)
    elseif !condition # == false
        #return false #a.s.# check = sumInterval(a,b)//2*1.0
        #GET Ceil & Floor

        #call middle
        #sumInterval(α,β)/2
        #does other things as well : # warning  #old-thinking
        check = sumInterval(α , β)/2 # floating-point division sumInterval(a, b) / 2 * 1.0 # freely allowing floats, to be ceiled & floored
        if isEven2(check)
            above = Int(ceil(check)) #nearest index above
            below = Int(floor(check))

      #  push!(q, below)#old thinking detected # to be changed
       # push!(q, above)

    else # faulty Input or Unexpected Error Occured
        #    return check  # nothing
        return #q #condition, check, above, below
    end
    return #q #condition, check, above, below
end
resultvector = middle(1,2)
#check length
length(resultvector) > 1 #there's atleast 1 middle

isEven2(length(resultvector )) #divisible /2
=#

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
export AbstractKet, Ket

"""
    AbstractKet{S <: AbstractTag}

Any type with `LHS` or `RHS` tag is a ket object. It should
support `dot` operation that raise an error when `LHS` and `RHS`
is in wrong order
"""
abstract type AbstractKet{S <: AbstractTag} end

import Base: show
import Compat.LinearAlgebra: dot

dot(rhs::AbstractKet{RHS}, lhs::AbstractKet{LHS}) =
    throw(TypeError(:dot, "wrong order of kets", LHS, RHS))

# A Naive Implementation
mutable struct Ket{T, S <: AbstractTag} <: AbstractKet{S}
    data::T
end

Ket(::Type{S}, data::T) where {T, S} = Ket{T, S}(data)

show(io::IO, ket::Ket{T, S}) where {T, S} = print(io, "Ket{$T, $S}")
dot(lhs::Ket{LT, LHS}, rhs::Ket{RT, RHS}) where {LT, RT} =
    dot(lhs.data, rhs.data)


text:
#!/usr/bin/env julia

function get_distance(coords::Array{Float64, 2}, c1::Int, c2::Int)
    # Returns the distance between two points
    return norm(coords[c1, :] - coords[c2, :])
end

function validate_fish_constraints(coords::Array{Float64, 2}, mds::MDSProblem)
    # checks that the FISH constraints are sort of close...
    # TODO: not yet implemented
    c = size(coords, 1)/size(mds.dist, 1)
end

function autoscale(coords::Array{Float64, 2}, mds::MDSProblem)
    # finds the scale that makes the FISH inferred distances correct
    d1 = mds.c1_c2_distance
    return d1/get_distance(coords, mds.control_locus_1, mds.control_locus_2)
end

function radius_scale(coords::Array{Float64, 2}, radius::Real)
    # finds the scale such that the radius of the chromosome becomes the
    # given radius.
    coords = center_data(coords)
    max_radius = maximum([norm(coords[i, :]) for i in 1:size(coords, 1)])
    return radius/max_radius
end

function run_mds(filename::AbstractString; scale::Number=1, constraint::AbstractString="",
            interp::Bool=true, exponent::Number=1/3, output_name::AbstractString="",
            auto_scale::Bool=false, shortest_paths::Bool=false,
            starting_points_file::AbstractString="")
    # loading data
    data = load_data(filename, scale=1, exponent=exponent)
    if shortest_paths
        data = contact_freq_to_distance_matrix(filename, scale=1,
                exponent=exponent, remove_inf=false)
    end
    mds = MDSProblem(data)
    # auto-creating output file name
    if output_name == ""
        prefix = split(filename, ".")[1]
        output_name = string(prefix, "_mds_coords.txt")
    end
    # initializing MDS problem
    mds = MDSProblem(data)
    zero_indices = get_removed_indices(mds)
    mds = remove_infs(mds)
    println("Running MDS on ", filename, "\n")
    println("Removed indices: ", zero_indices)
    println()
    ipopt_problem = make_ipopt_problem(mds, radius_constraint=false,
            fish_constraint=false)
    # dealing with constraints
    if constraint != ""
        mds = new_mds(data, 1.0, constraint)
        mds = remove_infs(mds)

        print_constraints(mds)
        println()
        ipopt_problem = make_ipopt_problem(mds, radius_constraint=false,
                fish_constraint=true)
    end
    # initial point
    if length(starting_points_file) > 0
        start_points = float64(readdlm(starting_points_file)[2:end,:])
        starting_x = reshape(start_points', length(start_points))
        if length(starting_x) != length(ipopt_problem.x)
            l1 = length(starting_x)
            l2 = length(ipopt_problem.x)
            if l1 < l2
                println("Error: not enough starting coords. Filling with random starting points.")
                # TODO: error?
                starting_x = cat(1, starting_x, ipopt_problem.x[l1+1:end])
            else
                println("Error: too many starting coords. Truncating.")
                println(l1)
                println(l2)
                starting_x = starting_x[1:l2]
            end
        else
            ipopt_problem.x = starting_x
        end
    end
    # solve problem
    solveProblem(ipopt_problem)
    coords = reshape(ipopt_problem.x, 3, round(Int, length(ipopt_problem.x)/3))'
    # provided scale
    if !auto_scale && scale != 1
        scale = radius_scale(coords, scale);
    end
    # autoscale
    if auto_scale && scale == 1 && constraint != ""
        scale = autoscale(coords, mds)
    end
    println("\nInferred scaling factor: ", scale, "\n")
    coords = center_data(coords)
    coords = scale * coords
    # interpolation
    if interp
        coords = interp_3d(coords, 10)
    end
    # file output
    println("\nWriting to file ", output_name)
    output_txt(coords, output_name)
    return coords
end

function parse_cl()
    # uses ArgParse to parse the command line
    s = ArgParseSettings()

    @add_arg_table s begin
        "filename"
            help = "CSV file containing HiC counts"
            required = true
        "--fish", "-f"
            help = "FISH constraints file"
            default = ""
        "--radius", "-r"
            help = "Radius of chromosome (determined via imaging)"
            arg_type = Float64
            default = 1.0
        "--output", "-o"
            help = "Output file name"
            default = ""
        "--init", "-i"
            help = "Starting coords file for optimization"
            default = ""
        "--interp"
            help = "Flag: use interpolation"
            action = :store_true
        "--auto-scale"
            help = "Flag: automatically infer the scaling factor (requires FISH constraints)"
            action = :store_true
        "--shortest-paths"
            help = "Flag: use shortest paths reconstruction"
            action = :store_true
    end

    return parse_args(s)
end

function mds_main()
    args = parse_cl()
    run_mds(args["filename"],
            constraint = args["fish"],
            scale = args["radius"],
            output_name = args["output"],
            auto_scale = args["auto-scale"],
            interp = args["interp"],
            shortest_paths = args["shortest-paths"],
            starting_points_file = args["init"])
end


text:
# Run tests by typing `Pkg.test("SPN")`.

using SPN
using Base.Test

"""
Creates the example SPN from the Poon paper.
"""
function create_toy_spn()
	s = SumNode()

	p1 = ProdNode()
	p2 = ProdNode()
	p3 = ProdNode()

	s1 = SumNode()
	s2 = SumNode()
	s3 = SumNode()
	s4 = SumNode()

	i1 = IndicatorNode(1,1)
	i2 = IndicatorNode(1,0)
	i3 = IndicatorNode(2,1)
	i4 = IndicatorNode(2,0)

	connect!(s,p1,weight=0.5)
	connect!(s,p2,weight=0.2)
	connect!(s,p3,weight=0.3)

	connect!(p1,s1)
	connect!(p1,s3)
	connect!(p2,s1)
	connect!(p2,s4)
	connect!(p3,s2)
	connect!(p3,s4)

	connect!(s1,i1,weight=0.6)
	connect!(s1,i2,weight=0.4)
	connect!(s2,i1,weight=0.9)
	connect!(s2,i2,weight=0.1)

	connect!(s3,i3,weight=0.3)
	connect!(s3,i4,weight=0.7)
	connect!(s4,i3,weight=0.2)
	connect!(s4,i4,weight=0.8)

	return s, p1, p2, p3, s1, s2, s3, s4, i1, i2, i3, i4
end

include("testnodes.jl")
include("testnetwork.jl")
include("testinference.jl")
include("testutils.jl")
include("learn/testpoon.jl")


text:
include("SentinelMissings.jl")
import .SentinelMissings
import YAXArrays.DAT: DATConfig
import YAXArrays.YAXTools: PickAxisArray

struct CubeIterator{R,ART,ARTBC,LAX,ILAX,S}
    dc::DATConfig
    r::R
    inars::ART
    inarsBC::ARTBC
    loopaxes::LAX
end
Base.IteratorSize(::Type{<:CubeIterator})=Base.HasLength()
Base.IteratorEltype(::Type{<:CubeIterator})=Base.HasEltype()
Base.eltype(i::Type{<:CubeIterator{A,B,C,D,E,F}}) where {A,B,C,D,E,F} = F

tuplelen(::Type{<:NTuple{N,<:Any}}) where N=N

lift64(::Type{Float32})=Float64
lift64(::Type{Int32})=Int64
lift64(T)=T

defaultval(t::Type{<:AbstractFloat})=convert(t,NaN)
defaultval(t::Type{<:Signed})=typemin(t)+1
defaultval(t::Type{<:Unsigned})=typemax(t)-1

function CubeIterator(dc,r;varnames::Tuple=ntuple(i->Symbol("x$i"),length(dc.incubes)),include_loopvars=())
    loopaxes = ntuple(i->dc.LoopAxes[i],length(dc.LoopAxes))
    inars, _ = getCubeCache(dc)
    length(varnames) == length(dc.incubes) || error("Supplied $(length(varnames)) varnames and $(length(dc.incubes)) cubes.")
    inarsbc = map(dc.incubes, inars) do ic, cache
      allax = getindsall(geticolon(ic), ic.loopinds, i->true)
      if ic.colonperm === nothing
        pa = PickAxisArray(cache,allax)
      else
        pa = PickAxisArray(cache,allax,ic.colonperm)
      end
    end
    et = map(inarsbc) do ibc
      ti = eltype(ibc)
      if ti <: AbstractArray
        ti
      else
        ti = Base.nonmissingtype(ti)
        SentinelMissings.SentinelMissing{ti,defaultval(ti)}
      end
    end
    if include_loopvars == true
      include_loopvars = map(axname,(loopaxes...,))
    end
    if isa(include_loopvars,Tuple) && !isempty(include_loopvars)
      ilax = map(i->findAxis(i,collect(loopaxes)),include_loopvars)
      any(isequal(nothing),ilax) && error("Axis not found in cubes")
      et=(et...,map(i->eltype(loopaxes[i]),ilax)...)
    else
      ilax=()
    end

    elt = NamedTuple{(map(Symbol,varnames)...,map(Symbol,include_loopvars)...),Tuple{et...}}
    CubeIterator{typeof(r),typeof(inars),typeof(inarsbc),typeof(loopaxes),ilax,elt}(dc,r,inars,inarsbc,loopaxes)
end
iternames(::CubeIterator{<:Any,<:Any,<:Any,<:Any,<:Any,E}) where E = E
tuplenames(t::Type{<:NamedTuple{N}}) where N = string.(N)
function Base.show(io::IO,ci::CubeIterator{<:Any,<:Any,<:Any,<:Any,<:Any,E}) where E
  print(io,"Datacube iterator with ", length(ci), " elements with fields: ",tuplenames(E))
end
Base.length(ci::CubeIterator)=prod(length.(ci.loopaxes))
function Base.iterate(ci::CubeIterator)
    rnow,blockstate = iterate(ci.r)
    updatears(ci.dc.incubes,rnow,:read,ci.inars)
    innerinds = CartesianIndices(length.(rnow))
    indnow, innerstate = iterate(innerinds)
    offs = map(i->first(i)-1,rnow)
    getrow(ci,ci.inarsBC,indnow,offs),
    (rnow=rnow,
      blockstate=blockstate,
      innerinds = innerinds,
      innerstate=innerstate
    )
end
function Base.iterate(ci::CubeIterator,s)
    t1 = iterate(s.innerinds,s.innerstate)
    N = tuplelen(eltype(ci.r))
    if t1 === nothing
        t2 = iterate(ci.r,s.blockstate)
        if t2 === nothing
            return nothing
        else
            rnow = t2[1]
            blockstate = t2[2]
            updatears(ci.dc.incubes,rnow,:read,ci.inars)
            innerinds = CartesianIndices(length.(rnow))
            indnow,innerstate = iterate(innerinds)

        end
    else
        rnow, blockstate = s.rnow, s.blockstate
        innerinds = s.innerinds
        indnow, innerstate = iterate(innerinds,s.innerstate)
    end
    offs = map(i->first(i)-1,rnow)
    getrow(ci,ci.inarsBC,indnow,offs),(rnow=rnow::NTuple{N,UnitRange{Int64}},
      blockstate=blockstate::Int64,
      innerinds=innerinds::CartesianIndices{N,NTuple{N,Base.OneTo{Int64}}},
      innerstate=innerstate::CartesianIndex{N})
end
abstract type CubeRow
end
abstract type CubeRowAx<:CubeRow
end

function getrow(ci::CubeIterator{<:Any,<:Any,<:Any,<:Any,ILAX,S},inarsBC,indnow,offs) where {ILAX,S}
   axvalsall = map((ax,i,o)->ax.values[i+o],ci.loopaxes,indnow.I,offs)
   axvals = map(i->axvalsall[i],ILAX)
   cvals  = map(i->i[indnow],inarsBC)
   S((cvals...,axvals...))
end
function getrow(ci::CubeIterator{<:Any,<:Any,<:Any,<:Any,(),S},inarsBC,indnow,offs) where S
   cvals  = map(i->i[indnow],inarsBC)
   S(cvals)
end

function Base.show(io::IO,s::CubeRow)
  print(io,"Cube Row: ")
  for n in propertynames(s)
    print(io,string(n), "=",getproperty(s,n)," ")
  end
end
function Base.show(io::IO,s::CubeRowAx)
  print(io,"Cube Row: ")
  for n in propertynames(s)
    print(io,string(n), "=",getproperty(s,n)," ")
  end
end
function Base.show(io::IO,s::Type{<:CubeRow})
  foreach(fieldnames(s)) do fn
    print(io,fn,"::",fieldtype(s,fn),", ")
  end
end
function Base.iterate(s::CubeRow,state=1)
  allnames = propertynames(s)
  if state<=length(allnames)
    (getproperty(s,allnames[state]),state+1)
  else
    nothing
  end
end

export CubeTable
"""
    CubeTable()

Function to turn a DataCube object into an iterable table. Takes a list of as arguments,
specified as a `name=cube` expression. For example
`CubeTable(data=cube1,country=cube2)` would generate a Table with the entries `data` and `country`,
where `data` contains the values of `cube1` and `country` the values of `cube2`. The cubes
are matched and broadcasted along their axes like in `mapCube`.

In addition, one can specify
`include_axes=(ax1,ax2...)` when one wants to include the values of certain axes in the table. For example
the command `(CubeTable(tair=cube1 axes=("lon","lat","time"))` would produce an iterator over a data structure
with entries `tair`, `lon`, `lat` and `time`.
"""
function CubeTable(;include_axes=(),expandaxes=(),cubes...)
  c = (map((k,v)->v,keys(cubes),values(cubes))...,)
  all(i->isa(i,Union{YAXArray,AbstractArray}),c) || throw(ArgumentError("All inputs must be DataCubes"))
  varnames = map(string,keys(cubes))
  expandaxes = isa(expandaxes,Tuple) ? expandaxes : (expandaxes,)
  inaxnames = Set{String}()
  indims = if isempty(expandaxes)
    map(i->InDims(),c)
  else
    map(c) do i
      axn = filter(collect(expandaxes)) do ax
        findAxis(ax,i) !== nothing
      end
      foreach(j->push!(inaxnames,axname(getAxis(j,i))), axn)
      InDims(axn...)
    end
  end
  axnames = map(i->axname.(caxes(i)),c)
  allvars = union(axnames...)
  allnums = collect(1:length(allvars))

  configiter = mapCube(identity,c,debug=true,indims=indims,outdims=(),ispar=false);
  # if inax !== nothing
  #   linax = length(inax)
  #   foreach(configiter.incubes) do ic1
  #     if !isempty(ic1.axesSmall)
  #       empty!(ic1.axesSmall)
  #       map!(i->i+1,ic1.loopinds,ic1.loopinds)
  #       pushfirst!(ic1.loopinds,1)
  #     else
  #       map!(i->i+1,ic1.loopinds,ic1.loopinds)
  #     end
  #   end
  # end
  r = collect(distributeLoopRanges((configiter.loopcachesize...,),(map(length,configiter.LoopAxes)...,),getchunkoffsets(configiter)))
  ci = CubeIterator(configiter,r,include_loopvars=include_axes,varnames=varnames)
end

import Tables
Tables.istable(::Type{<:CubeIterator}) = true
Tables.rowaccess(::Type{<:CubeIterator}) = true
Tables.rows(x::CubeIterator) = x
Tables.schema(x::CubeIterator) = Tables.schema(typeof(x))
Tables.schema(x::Type{<:CubeIterator}) = Tables.Schema(fieldnames(eltype(x)),map(s->fieldtype(eltype(x),s),fieldnames(eltype(x))))

Tables.istable(::Type{<:YAXArray}) = true
Tables.rowaccess(::Type{<:YAXArray}) = true
Tables.rows(x::YAXArray) = CubeTable(value = x, include_axes=true)
Tables.schema(x::YAXArray) = Tables.schema(Tables.rows(x))


text:
using Base.Test
using FluxMO.Validation

# http://swl.htwsaar.de/lehre/ss17/ml/slides/2017-vl-ml-ch4-1-clustering.pdf
function naive_intra_cluster_weights(C::AbstractArray)
    W_in = 0.0
    N_in = 0
    for (i,c) in enumerate(C)
        W_in += weights(c,c)
        N_in += length(c) * (length(c)-1)
    end
    0.5 * W_in, convert(Int64, 0.5 * N_in)
end

# http://swl.htwsaar.de/lehre/ss17/ml/slides/2017-vl-ml-ch4-1-clustering.pdf
function naive_n_in(C::AbstractArray)
    N_in = 0
    for (i,c) in enumerate(C)
        N_in += length(c) * (length(c)-1)
    end
    convert(Int64, 0.5 * N_in)
end

# http://swl.htwsaar.de/lehre/ss17/ml/slides/2017-vl-ml-ch4-1-clustering.pdf
function naive_inter_cluster_weights(C::AbstractArray)
    W_out = 0.0
    N_out = 0
    for (i,S) in enumerate(C)
        for (j,R) in enumerate(C)
            if i != j
                W_out += weights(S,R)
                N_out += length(S) * length(R)
            end
        end
    end
    0.5 * W_out, convert(Int64, 0.5 * N_out)
end

# http://swl.htwsaar.de/lehre/ss17/ml/slides/2017-vl-ml-ch4-1-clustering.pdf
function naive_n_out(C::AbstractArray)
    N_out = 0
    for (i, S) in enumerate(C)
        for (j, R) in enumerate(C)
            if i != j
                N_out += length(S) * length(R)
            end
        end
    end
    convert(Int64, 0.5 * N_out)
end

@testset "naive vs optimized for loops" begin

    # low BetaCV
    C = [
        [[1.0, 1.0], [1.1, 1.1]], # 1
        [[2.0, 2.0], [2.1, 2.1]], # 2
        [[3.0, 3.0], [3.1, 3.1]], # 3
        [[4.0, 4.0], [4.1, 4.1]], # 4
        [[5.0, 5.0], [5.1, 5.1]], # 5
        [[6.0, 6.0], [6.1, 6.1]], # 6
    ]

    # intra
    expected = naive_intra_cluster_weights(C)
    result   = intra_cluster_weights(C)
    @test expected[1] ≈ result[1]
    @test expected[2] ≈ result[2]
    @test naive_n_in(C) == 6
    @test naive_n_in(C) == result[2]

    # inter
    expected = naive_inter_cluster_weights(C)
    result   = inter_cluster_weights(C)
    @test expected[1] ≈ result[1]
    @test expected[1] ≈ 197.98989873223329
    @test expected[2] == result[2]
    @test expected[2] == 60
    @test naive_n_out(C) == 60

    @test betacv(C) ≈ 0.0428571428571428
    for _ in 1:10
        @time betacv(C)
    end

    @test betacv_fused(C) ≈ 0.0428571428571428
    for _ in 1:10
        @time betacv_fused(C)
    end

    # high BetaCV
    C = [
        [[1.0, 1.0], [6.1, 6.1]], # 1
        [[2.0, 2.0], [5.1, 5.1]], # 2
        [[3.0, 3.0], [4.1, 4.1]], # 3
        [[4.0, 4.0], [3.1, 3.1]], # 4
        [[5.0, 5.0], [2.1, 2.1]], # 5
        [[6.0, 6.0], [1.1, 1.1]], # 6
    ]

    @test betacv(C) ≈ 1.4681892332789561
    #               > 0.0428571428571428

end

@testset "tracked betacv with for loops" begin
    using Flux.Tracker

    C = [
        [param([1.0, 1.0]), param([1.1, 1.1])], # 1
        [param([2.0, 2.0]), param([2.1, 2.1])], # 2
        [param([3.0, 3.0]), param([3.1, 3.1])], # 3
        [param([4.0, 4.0]), param([4.1, 4.1])], # 4
        [param([5.0, 5.0]), param([5.1, 5.1])], # 5
        [param([6.0, 6.0]), param([6.1, 6.1])], # 6
    ]

    # intra
    expected = naive_intra_cluster_weights(C)
    result   = intra_cluster_weights(C)
    @test expected[1] ≈ result[1]
    @test expected[1].tracker.data ≈ 0.8485281374238982
    @test expected[2] == result[2]
    @test expected[2] == 6

    # inter
    expected = naive_inter_cluster_weights(C)
    result   = inter_cluster_weights(C)
    @test expected[1] ≈ result[1]
    @test expected[1] ≈ 197.98989873223329
    @test expected[2] == result[2]
    @test expected[2] == 60
    @test naive_n_out(C) == 60

    @test betacv(C).tracker.data ≈ 0.0428571428571428
    for _ in 1:10
        @time betacv(C)
    end
    @test betacv_fused(C).tracker.data ≈ 0.0428571428571428
    for _ in 1:10
        @time betacv_fused(C)
    end

    C[1][1].grad ≈ [0.0, 0.0]

    back!(betacv(C))

    @test C[1][1].grad[1] ≈ -0.034183673469387735
    @test C[1][1].grad[2] ≈ -0.034183673469387735
    @test C[1][2].grad[1] ≈ 0.03724489795918369
    @test C[1][2].grad[2] ≈ 0.03724489795918369

end

@testset "tracked betacv with distance-matrices" begin
    using Flux.Tracker

    C = [
        [1.0 1.1; 1.0 1.1], # 1
        [2.0 2.1; 2.0 2.1], # 2
        [3.0 3.1; 3.0 3.1], # 3
        [4.0 4.1; 4.0 4.1], # 4
        [5.0 5.1; 5.0 5.1], # 5
        [6.0 6.1; 6.0 6.1], # 6
    ]
    @test typeof(C[1]) == Array{Float64,2}

    C_tracked = [
        [param(1.0) param(1.1); param(1.0) param(1.1)], # 1
        [param(2.0) param(2.1); param(2.0) param(2.1)], # 2
        [param(3.0) param(3.1); param(3.0) param(3.1)], # 3
        [param(4.0) param(4.1); param(4.0) param(4.1)], # 4
        [param(5.0) param(5.1); param(5.0) param(5.1)], # 5
        [param(6.0) param(6.1); param(6.0) param(6.1)], # 6
    ]
    @test typeof(C_tracked[1]) == Array{Flux.Tracker.TrackedReal{Float64},2}

    # intra
    expected = intra_cluster_weights_pairwise(C)
    result   = intra_cluster_weights_pairwise(C_tracked)
    @test expected[1] ≈ result[1]
    @test expected[1] ≈ 0.8485281374238982
    @test expected[2] == result[2]
    @test expected[2] == 6

    # inter
    expected = inter_cluster_weights_pairwise(C)
    result   = inter_cluster_weights_pairwise(C_tracked)
    @test expected[1] ≈ result[1]
    @test expected[1] ≈ 197.98989873223329
    @test expected[2] == result[2]
    @test expected[2] == 60

    @time expected = betacv_pairwise(C)
    @time result = betacv_pairwise(C_tracked)
    @test expected            ≈ 0.0428571428571428
    @test result.tracker.data ≈ expected

    before = [t.tracker.grad for t in C_tracked[1]]
    @test isapprox(before, [0.0 0.0; 0.0 0.0], atol=1e-20)

    @time back!(betacv_pairwise(C_tracked))
    after = [t.tracker.grad for t in C_tracked[1]]
    @test isapprox(after,
        [-0.034183673469387735 0.03724489795918369;
         -0.034183673469387735 0.03724489795918369],
        atol=1e-20)

    # NOTE: it is not possible to compute the gradient
    # 2×2 Array{Float64,2}:
    # NaN  NaN
    # NaN  NaN

end

max_stars_repo_name: swiesend/FluxMO

max_issues_repo_name: swiesend/FluxMO

max_forks_repo_name: swiesend/FluxMO


text:
module BridgeSampling
using Turing, StatsBase, Printf
using LinearAlgebra: dot
import Base: show

include("iterative_algorithm.jl")
include("proposal.jl")
include("transform.jl")
include("model_interface.jl")
include("formatting.jl")
include("error_estimate.jl")

export bridgesampling, value, error_estimate, bayes_factor, posterior_probabilities, show

struct LogMarginalLikelihood{TV}
    value::Float64
    niter::Int
    p_posterior::TV
    g_posterior::TV
    p_proposal::TV
    g_proposal::TV
end

"""
    bridgesampling(samples, log_posterior, lb, ub; tol, maxiter, names)

    Estimates the log marginal likelihood with bridge sampling, given a matrix of posterior samples and an unnormalized log posterior function.

    # Arguments
    * `samples` : An ndims x nsamples matrix. The eltype can be either a scalar or an array of any dimension (e.g. multivariate parameters)
    * `log_posterior` : A function that takes one argument (one sample) and outputs the unnormalized log posterior probability (log joint)
    * `lb` and `ub` : respectively the lower and upper boundary of each dimension in the samples
    * `tol` : tolerance for the iterative scheme. Default = `1e-10`
    * `maxiter` : maximum number of iterations of the iterative scheme. Default = `1000`
    * `names` : A vector of names for each dimension in the samples. If specified, the input argument of `log_posterior` will be a `NamedTuple`. Default = `nothing`

    # Examples
    ```julia
        # Bivariate standard normal distribution
        samples = rand(MvNormal(ones(2)), 10_000)

        log_posterior(x) = -0.5 * x' * x

        Bridge = bridgesampling(samples, log_posterior, [-Inf, -Inf], [Inf, Inf])
    ```

    # Sources
        Quentin F. Gronau, Alexandra Sarafoglou, Dora Matzke, Alexander Ly, Udo Boehm, Maarten Marsman, David S. Leslie, Jonathan J. Forster, Eric-Jan Wagenmakers, Helen Steingroever,
        A tutorial on bridge sampling,
        Journal of Mathematical Psychology,
        Volume 81,
        2017,
        Pages 80-97,
        ISSN 0022-2496,
        https://doi.org/10.1016/j.jmp.2017.09.005.
"""
function bridgesampling(samples::AbstractMatrix, log_posterior::Function, lb, ub; tol=1e-10, maxiter=1_000, names=nothing)
    lb = informissing.(lb)
    ub = informissing.(ub)
    nd, ns = size(samples)
    nd < ns || @warn "The number of samples ($ns) is lower than the number of dimensions ($nd). Make sure the samples matrix is ndim x nsamples."

    ## Split samples
    smp1 = samples[:,1:2:end]
    smp2 = samples[:,2:2:end]
    n₁ = size(smp1, 2)
    n₂ = size(smp2, 2)

    ## Transform samples
    trans_smp1 = transform(smp1, lb, ub)
    trans_smp2 = transform(smp2, lb, ub)
    wide_smp1, name_map = check_samples(trans_smp1, names) # widen the transformed samples to match with the MvNormal proposal
    wide_smp2, _ = check_samples(trans_smp2, names)

    ## Fit the proposal
    prop_dist = proposal(wide_smp2)

    ## Sample from the proposal
    prop_samples = rand(prop_dist, n₂)
    comp_prop_smp = format_samples(prop_samples, name_map) # Reformat for posterior evaluation
    comp_prop_smp = invtransform(comp_prop_smp, lb, ub) # Inverse transform to evaluate the posterior

    ## Get the (log) pdf for each sample for both distributions
    p_post, g_post, l₁ = get_pdf_samples(smp1, wide_smp1, log_posterior, prop_dist, lb, ub, names)
    p_prop, g_prop, l₂ = get_pdf_samples(comp_prop_smp, prop_samples, log_posterior, prop_dist, lb, ub, names)

    ## Iterative algorithm
    logml, i = iterative_algorithm(l₁, l₂, n₁, n₂; tol=tol, maxiter=maxiter)
    return LogMarginalLikelihood(logml, i, p_post, g_post, p_prop, g_prop)
end

function bridgesampling(chn::MCMCChains.Chains, mdl::Tm, trans=missing; kwargs...) where Tm <: AbstractMCMC.AbstractModel

    samples, param_names = extract_samples(chn, mdl)
    log_posterior, dist = log_posterior_func(mdl, param_names)
    if isa(trans, NamedTuple)
        trans_vec = [trans[pn] for pn in param_names]
    else
        trans_vec = trues(length(param_names))
    end
    #all(isa.(dist, Distribution{Univariate})) || throw("Multivariate distributions are not supported")
    return bridgesampling(samples, log_posterior, dist, trans_vec; names=nothing, kwargs...)
end

function bridgesampling(chn::MCMCChains.Chains, args...; kwargs...)
    samples = permutedims(Array(MCMCChains.get_sections(chn, :parameters)), [2,1])
    return bridgesampling(samples, args...; kwargs...)
end

informissing(x) = x
informissing(x::T) where T<:Real = isinf(x) ? missing : x

value(LM::LogMarginalLikelihood; log::Bool=true) = log ? LM.value : exp(LM.value)

## Bayes factor
bayes_factor(LM1::LogMarginalLikelihood, LM2::LogMarginalLikelihood) =
    exp(value(LM1) - value(LM2))

## Posterior model probabilities
function posterior_probabilities(LM::TV, prior_prob=fill(1/length(LM), length(LM))) where {T <: LogMarginalLikelihood, TV <: AbstractVector{T}}
    isprobvec(prior_prob) || throw("Prior probabilities must be a proability vector")
    ml = value.(LM, log=false)
    post_prob = ml .* prior_prob
    post_prob ./= sum(post_prob)
    return post_prob
end

## Show
function Base.show(io::IO, ::MIME"text/plain", LML::LogMarginalLikelihood)
    @printf(io, "Log-marginal likelihood ≈ %.3f (%d iterations)", LML.value, LML.niter)
    err = error_estimate(LML)
    println(io)
    @printf(io, "Error estimate : %.3f %%", err.percent)
end

end


text:
# struct to model a cube

@tm1def mutable struct Cube
  Name::Union{String,Nothing}
  Rules::Union{String,Nothing}
  LastSchemaUpdate::Union{Dates.DateTime,Nothing}
  LastDataUpdate::Union{Dates.DateTime,Nothing}
  Attributes::Union{Dict,Nothing}
  Dimensions::Union{Vector{Dimension},Nothing}
  Views::Union{Vector,Nothing}
  ViewAttributes::Union{Vector,Nothing}
  PrivateViews::Union{Vector,Nothing}
  Annotations::Union{Vector,Nothing}
  LocalizedAttributes::Union{Vector,Nothing}
end

Cube(name::AbstractString) = Cube(Dict("Name" => name))

namefield(cube::Cube) = cube.Name

# functions for endpoints

@api_default function get_all_cubes(api::TM1API; options...)
  params = Dict("\$expand" => "Dimensions(\$select=Name)")
  result, page_data = tm1_get_paged_json(api, "Cubes"; params = params, options...)
  map(Cube, get(result, "value", [])), page_data
end

@api_default function get_all_model_cubes(api::TM1API; options...)
  params = Dict("\$expand" => "Dimensions(\$select=Name)")
  result, page_data = tm1_get_paged_json(api, "ModelCubes()"; params = params, options...)
  map(Cube, get(result, "value", [])), page_data
end

@api_default function get_all_control_cubes(api::TM1API; options...)
  params = Dict("\$expand" => "Dimensions(\$select=Name)")
  result, page_data = tm1_get_paged_json(api, "ControlCubes()"; params = params, options...)
  map(Cube, get(result, "value", [])), page_data
end

@api_default function get_cube(api::TM1API, cube_name::AbstractString; options...)
  params = Dict("\$expand" => "Dimensions(\$select=Name)")
  tm1_get_json(api, "Cubes('" * cube_name * " ')"; params = params, options...)
end

@api_default function delete_cube(api::TM1API, cube_name::AbstractString; options...)
  tm1_delete(api, "Cubes('" * cube_name * "')"; options...)
end

# is this the best way to achieve this?
@api_default function delete_cube(api::TM1API, cube::Cube; options...)
  delete_cube(api, name(cube); options...)
end

@api_default function create_cube(api::TM1API, cube::Cube; options...)
  # not working, not entirely sure this is how it should work
  # getting a weird odata error that doesn't help much
  tm1_put(api, "Cubes"; params = cube, options...)
end


text:
function horner(cfs::AbstractVector,z::T) where T
    ret = zero(T)
    for k = length(cfs):-1:1
        ret = ret .* z + cfs[k]
    end
    ret
end

## Evaluate at a point using Horner's method (above)
function evaluate(cfs::AbstractVector,S::MalmquistTakenaka,x)
    T = eltype(x)
    z = (-one(T)im)*(S.λ .- x)./(conj(S.λ) .- x)
    ret = (-one(T)im)*horner(cfs[1:2:end],z) ./ (conj(S.λ) .- x)
    ret = ret .+ horner(cfs[2:2:end],inv.(z)) ./ (S.λ .- x)
    ret *= sqrt(abs(imag(S.λ))/π)
    ret
end

function evaluate(cfs::AbstractVector,S::Weideman,x)
    T = eltype(x)
    z = (-one(T)im)*(S.λ .- x)./(conj(S.λ) .- x)
    ret = horner(cfs[1:2:end],z)
    ret = ret .+ horner([zero(T);cfs[2:2:end]],inv.(z))
    ret
end

## Reorders coefficients to be indexed by ℤ
function hatify(cfs::AbstractVector)
    newcfs = copy(cfs)
    n = length(cfs)
    newcfs[div(n,2)+mod(n,2)+1:n] = cfs[1:2:n]
    newcfs[div(n,2)+mod(n,2):-1:1] = cfs[2:2:n]
    -div(n,2)-mod(n,2):div(n,2)-1, newcfs
end

@recipe function f(F::Fun{Sp,T,Vector{T}}) where {Sp<:Union{MalmquistTakenaka,Weideman},T}
    S = space(F)
    λ = S.λ
    a = real(λ) - 30*imag(λ)
    b = real(λ) + 30*imag(λ)
    cfs = coefficients(F)
    if length(cfs) < 200
        pad!(cfs,200)
    end
    n = length(cfs)
    pts = points(S,n)
    pts = [pts[div(n,2)+2:n]; pts[1:div(n,2)+mod(n,2)]]
    vals = values(F)
    vals = [vals[div(n,2)+2:n]; vals[1:div(n,2)+mod(n,2)]]

    inds = findall(x -> a ≤ x ≤ b,pts)
    @series begin
        pts[inds],real(vals[inds])
    end
    if maximum(abs.(imag(vals[inds]))) ≥ 1e-6
        @series begin
            pts[inds],imag(vals[inds])
        end
    end
end


text:
using Distributions
using ChainRulesTestUtils
using Test

@testset "uniform.jl" begin
    # affine transformations
    test_affine_transformations(Uniform, rand(), 4 + rand())

    @testset "ChainRules" begin
        # run test suite for values in the support
        dist = Uniform(- 1 - rand(), 1 + rand())
        tangent = ChainRulesTestUtils.Tangent{Uniform{Float64}}(; a=randn(), b=randn())
        for x in (rand(), -rand())
            test_frule(logpdf, dist ⊢ tangent, x)
            test_rrule(logpdf, dist ⊢ tangent, x)
        end

        # check manually that otherwise derivatives are zero (FiniteDifferences returns NaN)
        for x in (-2, 2)
            # frule
            @test @inferred(
                ChainRulesTestUtils.frule(
                    (ChainRulesTestUtils.NoTangent(), tangent, randn()),
                    logpdf,
                    dist,
                    x,
                ),
            ) == (-Inf, 0.0)

            # rrule
            Ω, pullback = @inferred(ChainRulesTestUtils.rrule(logpdf, dist, x))
            @test Ω == -Inf
            @test @inferred(pullback(randn())) == (
                ChainRulesTestUtils.NoTangent(),
                ChainRulesTestUtils.Tangent{Uniform{Float64}}(; a=0.0, b=0.0),
                ChainRulesTestUtils.ZeroTangent(),
            )
        end
    end
end


text:
using Test, CornerTransferMethods, TensorOperations, TensorNetworkTensors, LinearAlgebra
@testset "CornerTransferMethods" begin
    @testset "tfisingmpo" begin
        β = randn()
        σz = [1 0; 0 -1]
        σx = [0 1; 1  0]
        id = [1 0; 0  1]
        @test ozofβ(β)[:,:] == exp(β * σz)
        oxx = oxxofβ(β, ComplexF64)
        @tensor res[1,2,3,4] := oxx[1,-1,3,-2] * oxx[2,-2,4,-1]
        ref = exp(2β*kron(σx,σx))
        resmat = fuselegs(res, ((1,2),(3,4)))[1]
        @test toarray(resmat) ≈ ref

        let β = 1e-5
            refe = exp(β * (kron(σz, id) + kron(id, σz) + 2kron(σx,σx)))
            ref = exp(β/2 * (kron(σz, id) + kron(id, σz))) *
                exp(2β * kron(σx,σx)) * #factor two probably edge-case
                exp(β/2 * (kron(σz, id) + kron(id, σz)))
            op = tfisingpropagator(β)
            @tensor op2[1,2,3,4] := op[1,-1,3,-2] *
                                    op[2,-2,4,-1];
            res = toarray(fuselegs(op2,((1,2),(3,4)))[1])
            @test ref ≈ refe && ref ≈ res
        end

        let β=1e-3
            refe = exp(β * (kron(σz,id,id) + kron(id,σz,id) + kron(id,id,σz) +
                            kron(σx,σx,id) + kron(σx,id,σx) + kron(id,σx,σx)))
            ref = exp(β/2 * (kron(σz,id,id) + kron(id,σz,id) + kron(id,id,σz))) *
                  exp(β *   (kron(σx,σx,id) + kron(id,σx,σx) + kron(σx,id,σx))) *
                  exp(β/2 * (kron(σz,id,id) + kron(id,σz,id) + kron(id,id,σz)))
            op = tfisingpropagator(β)
            @tensor op2[1,2,3,4,5,6] := op[1,-2,4,-1] *
                                        op[2,-3,5,-2] *
                                        op[3,-1,6,-3];
            res = toarray(fuselegs(op2,((1,2,3),(4,5,6)))[1])
            @test ref ≈ refe && ref ≈ res
        end
        op = tfisingpropagator(1e-5)
        @tensor op2[1,2,3,4,5,6] := op[1,-2,4,-1] * op[2,-3,5,-2] * op[3,-1,6,-3];
        prop = toarray(fuselegs(op2,((1,2,3),(4,5,6)))[1])
        ham = -(kron(σz,id,id) + kron(id,σz,id) + kron(id,id,σz) +
                kron(σx,σx,id) + kron(σx,id,σx) + kron(id,σx,σx))
        λp, vecsp = eigen(prop)
        λh, vecsh = eigen(ham)
        @test vecsp[:,argmax(λp)] ≈ vecsh[:,argmin(λh)]

        @test kron(-σx,σx) + kron(id,-σz) + kron(-σz,id) ≈ toarray(fuselegs(tfisinghamiltonian(true),((1,2),(3,4)))[1])

        sth0 = tfisingctm(1e-3, 8, 0, maxit=10^4, period = 10^2, verbose=false)
        @test mag(sth0) ≈ 1
        @test energy(sth0) ≈ -1/2

        sth∞ = tfisingctm(1e-3, 8, 1e5, maxit=10^4, period = 10^2, verbose=false)
        @test isapprox(mag(sth∞), 0; atol=1e-20)
        @test energy(sth∞) ≈ -1
    end

    @testset "algorithms equality" begin
        A = CornerTransferMethods.isingtensor(CornerTransferMethods.βc)
        χ = 4
        rotsymiter = rotsymctmiterable(A, χ)
        tciter = transconjctmiterable(A, χ)
        iter = ctmiterable(A, χ)
        rotsymstate = ctm_kernel(rotsymiter, maxit=10^6, verbose=false)
        tcstate = ctm_kernel(tciter, maxit=10^6, verbose=false)
        state = ctm_kernel(iter, maxit=10^6, verbose=false)
        eoes = cornereoe.([state.Cs[1], state.Cs[2], state.Cs[3], state.Cs[4],
                           rotsymstate.C, tcstate.C])
        @test all((x ≈ y for x in eoes, y in eoes))

        cls = clength.([state.Ts[1], state.Ts[2], state.Ts[3], state.Ts[4],
                        rotsymstate.T, tcstate.Ts[1], tcstate.Ts[2]])
        @test all((x ≈ y for x in cls, y in cls))
    end
end


text:
module ControlSystems

export  LTISystem,
        StateSpace,
        TransferFunction,
        ss,
        tf,
        tfg,
        zpk,
        ss2tf,
        # Linear Algebra
        balance,
        care,
        dare,
        dlyap,
        lqr,
        dlqr,
        kalman,
        dkalman,
        covar,
        norm,
        gram,
        ctrb,
        obsv,
        place,
        # Model Simplification
        reduce_sys,
        sminreal,
        minreal,
        balreal,
        # Stability Analysis
        isstable,
        pole,
        tzero,
        dcgain,
        zpkdata,
        damp,
        dampreport,
        markovparam,
        margin,
        delaymargin,
        gangoffour,
        # Connections
        append,
        series,
        parallel,
        feedback,
        feedback2dof,
        # Discrete
        c2d,
        # Time Response
        step,
        impulse,
        lsim,
        # Frequency Response
        freqresp,
        evalfr,
        bode,
        nyquist,
        sigma,
        # utilities
        numpoly,
        denpoly

import Plots
import Base: +, -, *, /, (./), (==), (.+), (.-), (.*), (!=), isapprox, call, convert, promote_op

include("types/lti.jl")
include("types/transferfunction.jl")
include("types/statespace.jl")
include("types/tf2ss.jl")

include("connections.jl")
include("discrete.jl")
include("matrix_comps.jl")
include("simplification.jl")
include("synthesis.jl")
include("analysis.jl")
include("timeresp.jl")
include("freqresp.jl")
include("utilities.jl")
include("plotting.jl")
include("pid_design.jl")

# The path has to be evaluated upon initial import
const __CONTROLSYSTEMS_SOURCE_DIR__ = dirname(Base.source_path())

end


text:
@testset "376.wiggle-subsequence.jl" begin
    @test wiggle_max_length([1, 2, 3, 4, 5, 6, 7, 8, 9]) == 2
    @test wiggle_max_length([1, 17, 5, 10, 13, 15, 10, 5, 16, 8]) == 7
    @test wiggle_max_length([1, 7, 4, 9, 2, 5]) == 6
    @test wiggle_max_length([0, 0]) == 1
    @test wiggle_max_length([0]) == 1
end


text:
"""
Abstract type used to represent genes on a `Record`. Subtypes must implement:
 * parent(gene)
 * index(gene)
 * locus(gene)
 * addgene!(chr, feature, locus)

This mainly exists to allow `Record`s and `Gene`s to refer to eachother,
and may be removed in the future.
"""
abstract type AbstractGene end

"""
Struct for storing information on genomic locations. `strand` can be '+', '-',
or '.' when the strand is irrelevant. `order` is used to store discontiguous
sequences, indicated in the GenBank file with the order() and join() operators.
"""
struct Locus
    position::UnitRange{Int}
    strand::Char
    complete_left::Bool
    complete_right::Bool
    order::Vector{UnitRange{Int}}
    join::Bool
end

"""
    Locus()
    Locus(position::UnitRange{Int})
    Locus(position::UnitRange{Int}, strand::Char)

"""
Locus() = Locus(1:1, '.', true, true, UnitRange{Int}[], false)
Locus(position::UnitRange{Int}) = Locus(position, '.', true, true, UnitRange{Int}[], false)
Locus(position::UnitRange{Int}, strand::Char) = Locus(position, strand, true, true, UnitRange{Int}[], false)
Locus(position::UnitRange{Int}, strand::Char, complete_left, complete_right) = Locus(position, strand, complete_left, complete_right, UnitRange{Int}[], false)

Base.convert(::Type{Locus}, x::UnitRange{Int}) = Locus(x)
function Base.convert(::Type{Locus}, x::StepRange{Int, Int})
    if x.step == -1
        return Locus(x.stop:x.start, '-')
    elseif x.step == 1
        return Locus(x.start:x.stop)
    end
    throw(DomainError(x, "`x` must have a step of 1 or -1"))
end

"""
Struct for storing annotations for a single chromosome, plasmid, contig, etc.
Contains five fields: `name`, `sequence`, `header`, `genes`, and `genedata`.
Annotations are stored as a `DataFrame` in `genedata`, but can be accessed
more easily through `genes` using the API provided in this module.
"""
mutable struct Record{G <: AbstractGene}
    name::String
    sequence::LongDNA{4}
    header::String
    genes::Vector{G}
    genedata::DataFrame
    circular::Bool
    Record{G}(name, sequence, header, genes::Vector{G}, genedata, circular) where G = new(name, sequence, header, genes, genedata, circular)
end

Record(args...) = Record{Gene}(args...)

struct Gene <: AbstractGene
    parent::Record{Gene}
    index::UInt
    locus::Locus
    feature::Symbol
end

Record{Gene}() = Record{Gene}("", dna"", "", Gene[], DataFrame(), false)

iscircular(c::Record{T}) where {T<:AbstractGene} = c.circular

"""
    addgene!(chr::Record, feature, locus; kw...)

Add gene to `chr`. `locus` can be a `Locus`, a UnitRange, or a StepRange (for
decreasing ranges, which will be annotated on the complementary strand).

# Example
```julia
addgene!(chr, "CDS", 1:756;
    locus_tag = "gene0001",
    product = "Chromosomal replication initiator protein dnaA")
```
"""
function addgene!(chr::Record{Gene}, feature, locus; kw...)
    locus = convert(Locus, locus)
    push!(chr.genedata, fill(missing, size(chr.genedata, 2)))
    index = UInt32(length(chr.genes) + 1)
    gene = Gene(chr, index, locus, feature)
    for (k, v) in kw
        Base.setproperty!(gene, k, v)
    end
    push!(chr.genes, gene)
    return gene
end

"""
    delete!(gene::AbstractGene)

Delete `gene` from `parent(gene)`. Warning: does not work when broadcasted! Use
`delete!(::AbstractVector{Gene}) instead.`
"""
function Base.delete!(gene::AbstractGene)
    i = index(gene)
    #deleteat!(parent(gene).genedata, i)
    deleteat!(parent(gene).genes, i)
    nothing
end

"""
    delete!(genes::AbstractArray{Gene, 1})

Delete all genes in `genes` from `parent(genes[1])`.

# Example
```julia
delete!(@genes(chr, length(gene) <= 60))
```
"""
function Base.delete!(genes::AbstractArray{Gene, 1})
    indices = index(gene)
    DataFrames.delete!(parent(genes).genedata, Int.(indices))
    lastindices = length(parent(genes).genes) - length(indices) + 1 : length(parent(genes).genes)
    deleteat!(parent(genes).genes, lastindices)
    nothing
end

function Base.propertynames(gene::G) where {G <: AbstractGene}
    propertynames(parent(gene).genedata)
end

"""
    vectorise(A::AbstractArray{Union{Missing, T}, 1}) where T

Convert an array of type Vector{Union{Missing, T}} to Vector{Union{Missing, Vector{T}}}.
"""
function vectorise(A::AbstractArray{Union{Missing, T}, 1}) where T
    B = Union{Missing, Vector{T}}[]
    for a in A
        if ismissing(a)
            push!(B, missing)
        else
            push!(B, [a])
        end
    end
    return B
end

"""
    pushproperty!(gene::AbstractGene, qualifier::Symbol, value::T)

Add a property to `gene`, similarly to `Base.setproperty!(::gene)`, but if the
property is not missing in `gene`, it will be transformed to store a vector
instead of overwriting existing data.

```julia
julia> eltype(chr.genedata[!, :EC_number])
Union{Missing,String}

julia> chr.genes[1].EC_number = "EC:1.2.3.4"
"EC:1.2.3.4"

julia> pushproperty!(chr.genes[1], :EC_number, "EC:4.3.2.1"); chr.genes[1].EC_number
2-element Array{String,1}:
 "EC:1.2.3.4"
 "EC:4.3.2.1"

julia> eltype(chr.genedata[!, :EC_number])
Union{Missing, Array{String,1}}
```
"""
function pushproperty!(gene::AbstractGene, qualifier::Symbol, value::T; forceany = true) where T
    if hasproperty(parent(gene).genedata, qualifier)
        C = eltype(parent(gene).genedata[!, qualifier])
        if T <: C
            if ismissing(parent(gene).genedata[index(gene), qualifier])
                parent(gene).genedata[index(gene), qualifier] = value
            else
                parent(gene).genedata[!, qualifier] = vectorise(parent(gene).genedata[!, qualifier])
                push!(parent(gene).genedata[index(gene), qualifier], value)
            end
        elseif Vector{T} <: C
            if ismissing(parent(gene).genedata[index(gene), qualifier])
                parent(gene).genedata[index(gene), qualifier] = [value]
            else
                push!(parent(gene).genedata[index(gene), qualifier], value)
            end
        elseif forceany && !(C <: AbstractVector)
            if ismissing(parent(gene).genedata[index(gene), qualifier])
                parent(gene).genedata[!, qualifier] = convert(Vector{Any}, parent(gene).genedata[!, qualifier])
                parent(gene).genedata[index(gene), qualifier] = value
            else
                parent(gene).genedata[!, qualifier] = vectorise(convert(Vector{Any}, parent(gene).genedata[! ,qualifier]))
                push!(parent(gene).genedata[index(gene), qualifier], value)
            end
        elseif forceany && C <: AbstractVector
            if ismissing(parent(gene).genedata[index(gene), qualifier])
                parent(gene).genedata[!, qualifier] = convert(Vector{Any}, parent(gene).genedata[!, qualifier])
                parent(gene).genedata[index(gene), qualifier] = [value]
            else
                @error "This shouldn't happen"
            end
        else
            @error "Tried to add a '$T' to '$qualifier::$(typeof(parent(gene).genedata[!, qualifier]))'"
        end
    else
        s = size(parent(gene).genedata, 1)
        parent(gene).genedata[!, qualifier] = Vector{Union{Missing, T}}(missing, s)
        parent(gene).genedata[index(gene), qualifier] = value
    end
    return value
end

"""
    get(g::AbstractGene, key, default)

Retrieve `key` from `g`. If `key` is missing, return `default`.
"""
function Base.get(g::AbstractGene, key, default)
    p = getproperty(g, key)
    if ismissing(p)
        return default
    else
        return p
    end
end

"""
    sequence(gene::AbstractGene; translate = false)

Return genomic sequence for `gene`. If `translate` is `true`, the sequence will be translated to a `LongAA`, excluding the stop, otherwise it will be returned as a `LongDNA{4}` (including the stop codon).
```
"""
function sequence(gene::AbstractGene; translate = false)
    if locus(gene).strand == '-'
        s = reverse_complement(parent(gene).sequence[locus(gene).position])
    else
        s = parent(gene).sequence[locus(gene).position]
    end
    translate ? BioSequences.translate(s)[1:end-1] : s
end

Base.length(gene::AbstractGene) = length(locus(gene))
Base.length(locus::Locus) = length(locus.position)

function sequence(record::Record)
    record.sequence
end

"""
    iscomplement(gene)

Return `true` if `locus(gene).compliment == '-'`, otherwise return `false`.
"""
iscomplement(gene::AbstractGene) = locus(gene).strand == '-'

"""
    iscomplete(gene)

Return `true` if `gene` is a complete gene, i.e. not a pseudo gene or partial.
"""
iscomplete(gene::AbstractGene) = !any(get(gene, :pseudo, false)) && !any(get(gene, :ribosomal_slippage, false)) && locus(gene).complete_right && locus(gene).complete_left

function appendstring(field, v::Bool)
    return "\n" * join(fill(' ', 21)) * "/$field"
end
function appendstring(field, v::Union{Number, Symbol})
    return "\n" * join(fill(' ', 21)) * "/$field=" * string(v)
end
function appendstring(field, v)
    v = multiline(v, field)
    v = replace(v, "\n" => "\n" * join(fill(' ', 21)))
    return "\n" * join(fill(' ', 21)) * "/$field=\"" * v * "\""
end

function Base.show(io::IO, gene::AbstractGene)
    buf = IOBuffer()
    print(buf, "     " * rpad(string(feature(gene)), 16, ' '))
    print(buf, locus(gene))
    for field in names(parent(gene).genedata)
        field in [:feature, :locus] && continue
        v = parent(gene).genedata[index(gene), field]
        if !ismissing(v)
            if v isa AbstractVector
                for i in eachindex(v)
                    print(buf, appendstring(field, v[i]))
                end
            else
                print(buf, appendstring(field, v))
            end
        end
    end
    println(io, String(take!(buf)))
end

Base.show(io::IO, mime::MIME"text/plain", genes::Array{G, 1}) where {G <: AbstractGene} = Base.show(io, genes)
function Base.show(io::IO, genes::Array{G, 1}) where {G <: AbstractGene}
    for gene in genes
        show(io, gene)
    end
end

function Base.show(io::IO, chr::Record)
    s = "Chromosome '$(chr.name)' ($(length(chr.sequence)) bp) with $(length(chr.genes)) annotations\n"
    print(io, s)
end

function Base.show(io::IO, genes::Vector{AbstractGene})
    print(join(genes, "\n"))
end

function Base.show(io::IO, locus::Locus)
    s = ""
    locus.strand == '-'     && (s *= "complement(")
    !locus.complete_left    && (s *= ">")
    if length(locus.order) > 0
        if locus.join
            s *= "join(" * join([join((r.start, r.stop), "..") for r in locus.order], ",") * ")"
        else
            s *= "order(" * join([join((r.start, r.stop), "..") for r in locus.order], ",") * ")"
        end
    else
        s *= join((locus.position.start, locus.position.stop), "..")
    end
    !locus.complete_right   && (s *= "<")
    locus.strand == '-'     && (s *= ")")
    print(io, s)
end

Base.copy(chr::Record) = deepcopy(chr)
Base.copy(gene::AbstractGene) = deepcopy(gene)

function Base.isless(l1::Locus, l2::Locus)
    if (l1.position.start < l2.position.start) || (l1.position.start == l2.position.start && l1.position.stop > l2.position.stop)
        return true
    end
    return false
end
Base.isless(g1::AbstractGene, g2::AbstractGene) = ((locus(g1) == locus(g2)) && (feature(g1) == "gene" && feature(g2) != "gene")) || (locus(g1) < locus(g2))

function Base.:(==)(x::Locus, y::Locus)
    (x.position == y.position) && (x.strand == y.strand) && (x.complete_left == y.complete_left) && (x.complete_right == y.complete_right) && (x.order == y.order)
end

Base.in(loc::Locus, r::UnitRange{Int}) = loc.position.start in r && loc.position.stop in r
Base.in(loc1::Locus, loc2::Locus) = loc1 in loc2.position
Base.intersect(loc1::Locus, loc2::Locus) = intersect(loc1.position, loc2.position)
Base.iterate(loc::Locus) = iterate(loc.position)

index(g::Gene) = getfield(g, :index)
locus(g::Gene) = getfield(g, :locus)
feature(g::Gene) = getfield(g, :feature)
Base.parent(g::Gene) = getfield(g, :parent)
function Base.parent(gs::AbstractVector{G}) where {G <: AbstractGene}
    @assert !isempty(gs) "Trying to get parent of empty Array{Gene}"
    getfield(gs[1], :parent)
end


text:
"""
    subspace_construction(model, cost, data, opt; T = 10, c = 1, M = 3, print_freq = 1)

Izmailov, P., Maddox, W. J., Kirichenko, P., Garipov, T., Vetrov, D., & Wilson, A. G. (2020, August). Subspace inference for Bayesian deep learning.
In Uncertainty in Artificial Intelligence (pp. 1169-1179). PMLR.

To construct subspace from pretrained weights.

# Input Arguments
- `model` 	 : Machine learning model. Eg: Chain(Dense(10,2)). Model should be created with Chain in Flux
- `cost`  	 : Cost function. Eg: L(x, y) = Flux.Losses.mse(m(x), y)
- `data` 	 : Inputs and outputs. Eg:	X = rand(10,100); Y = rand(2,100); data = DataLoader(X,Y);
- `opt`		 : Optimzer. Eg: opt = ADAM(0.1)

# Keyword Arguments
- `T` 		  : Number of steps for subspace calculation. Eg: T= 1
- `c` 		  : Moment update frequency. Eg: c = 1
- `M` 		  : Maximum number of columns in deviation matrix. Eg: M= 2
- `print_freq`: Loss printing frequency

# Outputs
- `W_swa`    : Mean weights
- `P` 		 : Projection Matrix
"""
function subspace_construction(model, cost, data, opt; T = 10, c = 1, M = 3, print_freq = 1
)
	training_loss = 0.0

	ps = Flux.params(model)
	W_swa = zeros(length(extract_params(ps)))
	all_len = length(W_swa)
	A = Array{eltype(W_swa)}(undef,0) #initialize deviation matrix

	# #initaize weights with mean
	all_weights = []
	for i in 1:T
		for d in data
			gs = gradient(ps) do
				training_loss = cost(model, d...)
				return training_loss
			end
			Flux.update!(opt, ps, gs)
			if mod(i,c) == 0
				W = extract_params(ps)
				n = i/c
				W_swa = (n.*W_swa + W)./(n+1)
				# if(length(A) >= M*all_len)
				# 	A = A[1:(end - all_len)]
				# end
				W_dev =  W - W_swa
				append!(A, W_dev)
			end
		end
		#print loss based o print frequency
		if (mod(i,print_freq) == 0 )|| (i == T)
			println("Traing loss: ", training_loss," Epoch: ", i)
		end
	end

	A = reshape(A, all_len, :)
	#take pseudo SVD and take only first M columns
	U,s,V = psvd(A)
	#generate projection matrix
	P = U[:,1:M]*LinearAlgebra.Diagonal(s[1:M])
	return W_swa, P
end

"""
	auto_encoder_subspace(model, cost, data, opt, encoder, decoder; T = 10, c = 1, M = 3, print_freq = 1
)

To construct subspace from pretrained weights using autoencoders.

# Input Arguments
- `model` 	 : Machine learning model. Eg: Chain(Dense(10,2)). Model should be created with Chain in Flux
- `cost`  	 : Cost function. Eg: L(x, y) = Flux.Losses.mse(m(x), y)
- `data` 	 : Inputs and outputs. Eg:	X = rand(10,100); Y = rand(2,100); data = DataLoader(X,Y);
- `opt`		 : Optimzer. Eg: opt = ADAM(0.1)
- `encoder`	 : Encoder to generate subspace from NN or Neural ODE parameters
- `decoder`	 : Decoder to generate NN or Neural ODE parameters from subspace

# Keyword Arguments
- `T` 		  : Number of steps for subspace calculation. Eg: T= 1
- `c` 		  : Moment update frequency. Eg: c = 1
- `M` 		  : Maximum number of columns in deviation matrix. Eg: M= 2
- `print_freq`: Loss printing frequency

# Outputs
- `W_swa`    : Mean weights
- `decoder`  : Trained decoder to generate NN or Neural ODE parameters from subspace
"""
function auto_encoder_subspace(model, cost, data, opt, encoder, decoder; T = 10, c = 1, M = 3, print_freq = 1
)
	training_loss = 0.0

	ps = Flux.params(model)
	W_swa = zeros(length(extract_params(ps)))
	all_len = length(W_swa)
	A = Array{eltype(W_swa)}(undef,0) #initialize deviation matrix

	# #initaize weights with mean
	all_weights = []
	for i in 1:T
		for d in data
			gs = gradient(ps) do
				training_loss = cost(model, d...)
				return training_loss
			end
			Flux.update!(opt, ps, gs)
			if mod(i,c) == 0
				W = extract_params(ps)
				n = i/c
				W_swa = (n.*W_swa + W)./(n+1)
				W_dev =  W - W_swa
				append!(A, W_dev)
			end
		end
		#print loss based o print frequency
		if (mod(i,print_freq) == 0 )|| (i == T)
			println("Traing loss: ", training_loss," Epoch: ", i)
		end
	end

	re_weight = reshape(A, all_len, :)
	# model
	sae = Chain(encoder,decoder)
	autops = Flux.params(sae);
	# loss
	autoloss(X) = Flux.mse(sae(X), X)
	+ sum(sqnorm, autops)

	# optimization
	opt = ADAM()
	# parameters

	cb() = @show autoloss(re_weight)
	#data
	wdata = DataLoader(re_weight, batchsize = 5, shuffle = true)
	Flux.train!(autoloss, autops, wdata, opt; cb=cb)
	return W_swa, decoder

end

"""
    diffusion_subspace(model, cost, data, opt; T = 10, c = 1, M = 3, print_freq = 1)

To construct subspace from pretrained weights using diffusion maps.

# Input Arguments
- `model` 	 : Machine learning model. Eg: Chain(Dense(10,2)). Model should be created with Chain in Flux
- `cost`  	 : Cost function. Eg: L(x, y) = Flux.Losses.mse(m(x), y)
- `data` 	 : Inputs and outputs. Eg:	X = rand(10,100); Y = rand(2,100); data = DataLoader(X,Y);
- `opt`		 : Optimzer. Eg: opt = ADAM(0.1)

# Keyword Arguments
- `T` 		  : Number of steps for subspace calculation. Eg: T= 1
- `c` 		  : Moment update frequency. Eg: c = 1
- `M` 		  : Maximum number of columns in deviation matrix. Eg: M= 2
- `print_freq`: Loss printing frequency

# Outputs
- `W_swa`    : Mean weights
- `P` 		 : Projection Matrix
"""
function diffusion_subspace(model, cost, data, opt; T = 10, c = 1, M = 3, print_freq = 1
)
	training_loss = 0.0

	ps = Flux.params(model)
	W_swa = zeros(length(extract_params(ps)))
	all_len = length(W_swa)
	A = Array{eltype(W_swa)}(undef,0) #initialize deviation matrix

	# #initaize weights with mean
	all_weights = []
	for i in 1:T
		for d in data
			gs = gradient(ps) do
				training_loss = cost(model, d...)
				return training_loss
			end
			Flux.update!(opt, ps, gs)
			if mod(i,c) == 0
				W = extract_params(ps)
				n = i/c
				W_swa = (n.*W_swa + W)./(n+1)
				# if(length(A) >= M*all_len)
				# 	A = A[1:(end - all_len)]
				# end
				W_dev =  W - W_swa
				append!(A, W_dev)
			end
		end
		#print loss based o print frequency
		if (mod(i,print_freq) == 0 )|| (i == T)
			println("Traing loss: ", training_loss," Epoch: ", i)
		end
	end

	A = reshape(A, all_len, :)
	#diffusion map
	diff_M = fit(DiffMap, A', maxoutdim=M) # construct diffusion map model
	P = transform(diff_M)
	return W_swa, P'
end


text:
a = 1:1:20
b = 21:1:40

iter_prod = Iterators.product(a,b)

Iterators.take(iter_prod,100)

size(iter_prod)

for i=1:10
    subset = Iterators.take(iter_prod,100)
    for x in subset
        println(x)
    end
    iter_prod = Iterators.drop(iter_prod,100)
end

length(iter_prod)

n_batch = 10
batch_size = length(iter_prod)÷n_batch

for i=1:n_batch
    subset = Iterators.take(iter_prod,batch_size)
    for x in subset
        println(x)
    end
    iter_prod = Iterators.drop(iter_prod,batch_size)
end

iter_prod = Iterators.product(a,b)
n_batch = 10
batch_size = length(iter_prod)÷n_batch
fid = h5open("../DHC/scratch_AKS/data/test.h5","w")
dset = create_dataset(fid, "A", datatype(Float64), dataspace(2,400), chunk=(2,40))
close(fid)
for i=1:n_batch
    subset = Iterators.take(iter_prod,batch_size)
    temp = zeros(2,batch_size)
    for (ind,x) in enumerate(subset)
        temp[:,ind] .= x[:]
    end
    fid = h5open("../DHC/scratch_AKS/data/test.h5","r+")
    data = open_dataset(fid, "A")
    data[:,(1:batch_size).+(i-1)*batch_size] = temp
    close(fid)
    iter_prod = Iterators.drop(iter_prod,batch_size)
end

A = h5read("../DHC/scratch_AKS/data/test.h5", "A")

iter_prod = Iterators.product(a,b)
n_batch = 10
batch_size = length(iter_prod)÷n_batch
fid = h5open("../DHC/scratch_AKS/data/testc.h5","w")
dset = create_dataset(fid, "A", datatype(Float64), dataspace(2,400), chunk=(2,40),compress=3)
close(fid)
for i=1:n_batch
    subset = Iterators.take(iter_prod,batch_size)
    temp = zeros(2,batch_size)
    for (ind,x) in enumerate(subset)
        temp[:,ind] .= x[:]
    end
    fid = h5open("../DHC/scratch_AKS/data/testc.h5","r+")
    data = open_dataset(fid, "A")
    data[:,(1:batch_size).+(i-1)*batch_size] = temp
    close(fid)
    iter_prod = Iterators.drop(iter_prod,batch_size)
end

A = h5read("../DHC/scratch_AKS/data/testc.h5", "A")

iter_prod = Iterators.product(a,b)
n_batch = 10
batch_size = length(iter_prod)÷n_batch
fid = h5open("../DHC/scratch_AKS/data/testc2.h5","w")
dset = create_dataset(fid, "A", datatype(Float64), dataspace(2,400), chunk=(2,40),compress=6)
close(fid)
for i=1:n_batch
    subset = Iterators.take(iter_prod,batch_size)
    temp = zeros(2,batch_size)
    for (ind,x) in enumerate(subset)
        temp[:,ind] .= x[:]
    end
    fid = h5open("../DHC/scratch_AKS/data/testc2.h5","r+")
    data = open_dataset(fid, "A")
    data[:,(1:batch_size).+(i-1)*batch_size] = temp
    close(fid)
    iter_prod = Iterators.drop(iter_prod,batch_size)
end

A = h5read("../DHC/scratch_AKS/data/testc2.h5", "A")

iter_prod = Iterators.product(a,b)
n_batch = 10
batch_size = length(iter_prod)÷n_batch
fid = h5open("../DHC/scratch_AKS/data/testc2.h5","w")
dset = create_dataset(fid, "A", datatype(Float64), dataspace(2,400), chunk=(2,40),suhffle=(),deflate=3)
close(fid)
for i=1:n_batch
    subset = Iterators.take(iter_prod,batch_size)
    temp = zeros(2,batch_size)
    for (ind,x) in enumerate(subset)
        temp[:,ind] .= x[:]
    end
    fid = h5open("../DHC/scratch_AKS/data/testc2.h5","r+")
    data = open_dataset(fid, "A")
    data[:,(1:batch_size).+(i-1)*batch_size] = temp
    close(fid)
    iter_prod = Iterators.drop(iter_prod,batch_size)
end

A = h5read("../DHC/scratch_AKS/data/testc2.h5", "A")

## Ok, so the summary is to just batch it and save with some showprogress on the outer on and deflate
# by = 3 (need to better understand what that number is) I think it is related to window width of
# zlib... tbd

iter_prod = Iterators.product(a,b)
n_batch = 10
batch_size = length(iter_prod)÷n_batch
fid = h5open("../DHC/scratch_AKS/data/testc2.h5","w")
dset = create_dataset(fid, "A", datatype(Float64), dataspace(2,400), chunk=(2,40),shuffle=(),deflate=3)
close(fid)
for i=1:n_batch
    subset = Iterators.take(iter_prod,batch_size)
    temp = zeros(2,batch_size)
    for (ind,x) in enumerate(subset)
        temp[:,ind] .= x[:]
    end
    fid = h5open("../DHC/scratch_AKS/data/testc2.h5","r+")
    data = open_dataset(fid, "A")
    data[:,(1:batch_size).+(i-1)*batch_size] = temp
    close(fid)
    iter_prod = Iterators.drop(iter_prod,batch_size)
end

A = h5read("../DHC/scratch_AKS/data/testc2.h5", "A")

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC