text:
# ----------------------------------------------------------------------------------------------------
# Xavier Gandibleux - Metaheuristiques - Algorithmes genetiques - Probleme du robot dans un labyrinthe

using Printf

# ----------------------------------------------------------------------------------------------------
# Contruit un labyrinthe de 10x10 pour un affichage style 1970
function construireLabyrinthe()

    carte = Matrix{Char}(undef,21,21)
    for i=1:21
        for j=1:21
            carte[i,j] = ' '
        end
    end

    for i=1:2:21
        for j=1:2:21
            carte[i,j] = '.'
        end
    end

    for j=2:2:20
        carte[1,j] = '-'
        carte[21,j] = '-'
    end

    for i=2:2:20
        carte[i,1] = '|'
        carte[i,21] = '|'
    end

    carte[2,1] = '>'
    carte[20,21] = ' '
    return carte
end

# ----------------------------------------------------------------------------------------------------
# Affichage alphanumerique du labyrinthe
function afficheLabyrinthe(carte)

    @printf("\n")
    for i=1:21
        @printf("    ")
        for j=1:21
            @printf("%c",carte[i,j])
        end
        @printf("\n")
    end
    @printf("\n")
end

# ----------------------------------------------------------------------------------------------------
# Affichage alphanumerique du labyrinthe
function afficheIndividu(indIn, carteIn)

    ind = indIn[1]
    fitness = indIn[2]
    carte = copy(carteIn)

    i=1; iplot = 2
    j=1; jplot = 2
    carte[iplot,jplot] = repr(0)[1]
    for longueur=2:fitness
            if     (ind[i,j] == 1) # Nord
                i-=1
            elseif (ind[i,j] == 2) # Est
                j+=1
            elseif (ind[i,j] == 3) # Sud
                i+=1
            elseif (ind[i,j] == 4) # Ouest
                j-=1
            end
            carte[2*i,2*j] = repr((longueur-1)%10)[1]
    end

    @printf("\n")
    for i=1:21
        @printf("    ")
        for j=1:21
            @printf("%c",carte[i,j])
        end
        @printf("\n")
    end
    @printf("\n")
end

# ----------------------------------------------------------------------------------------------------
# Ajoute quelques obstacles au labyrinthe
function contrainteSoftLabyrinthe(carte)

    for i=1:10
        j=rand(1:10)
        pattern=rand([1,2,4,8])
        if     (pattern == 1)
            carte[2*i,2*j-1] = '|'
        elseif (pattern == 2)
            carte[2*i+1,2*j] = '-'
        elseif (pattern == 4)
            carte[2*i,2*j+1] = '|'
        elseif (pattern == 8)
            carte[2*i-1,2*j] = '-'
        end
    end
    carte[2,1] = '>'
    carte[20,21] = ' '
end

# ----------------------------------------------------------------------------------------------------
# Creation d'une population d'individus
function creerPopulation(n,popSize,carte)

    NbRealisable = 0
    realisable = false
    maxFitness = 0
    population = 0
    ind=Array{Int}(undef, 10,10)

    pop=Vector(undef, popSize)

    for individu =1:popSize
    #while realisable ==false
        ind=rand(1:4,n,n)
        ind[1,1]=rand([1,2,3]) # interdit de sortir
        ind[10,10]=2 # force a sortir
        #individu +=1

        visite=Matrix{Bool}(undef,10,11)
        visite=fill(false,10,11)
        i=1; j=1; fitness=0; avance = true
        visite[i,j]=true
        while avance && j!=11
            if ind[i,j] ==1 && carte[2*i-1,2*j] ==' ' && visite[i-1,j]==false
                i-=1 ; fitness +=1; visite[i,j]=true #; println("Nord")
            elseif ind[i,j] ==2 && carte[2*i,2*j+1] ==' ' && visite[i,j+1]==false
                j+=1 ; fitness +=1; visite[i,j]=true #; println("Est")
            elseif ind[i,j] ==3 && carte[2*i+1,2*j] ==' ' && visite[i+1,j]==false
                i+=1 ; fitness +=1; visite[i,j]=true #; println("Sud")
            elseif ind[i,j] ==4 && carte[2*i,2*j-1] ==' ' && visite[i,j-1]==false
                j-=1 ; fitness +=1; visite[i,j]=true #; println("Ouest")
            else avance = false #; print(individu," "); println("fitness = ", fitness, "   i = ",i, " j = ",j);
                maxFitness = max(fitness, maxFitness)
            end
        end # while
        if j==11 println("REALISABLE")
            print(individu," "); println("fitness = ", fitness, "   i = ",i, " j = ",j);
            NbRealisable +=1
            realisable = true
        end
        pop[individu] = (ind , fitness , realisable)
    end # for individu
    println("Nbre Realisable = ", NbRealisable, " maxFitness = ", maxFitness)
    return pop
end

# ----------------------------------------------------------------------------------------------------
# Evaluation d'un individu
function evaluerIndividu(n,carte,ind)

    visite = fill(false,10,11)
    i=1; j=1; fitness = 0; avance = true; realisable = false
    visite[i,j]=true
    while (avance) && (!realisable)
            if     (ind[i,j] == 1) && (carte[2*i-1,2*j] ==' ') && (visite[i-1,j] == false)  # Nord
                i-=1 ; fitness +=1; visite[i,j]=true
            elseif (ind[i,j] == 2) && (carte[2*i,2*j+1] ==' ') && (visite[i,j+1] == false)  # Est
                j+=1 ; fitness +=1; visite[i,j]=true
            elseif (ind[i,j] == 3) && (carte[2*i+1,2*j] ==' ') && (visite[i+1,j] == false)  # Sud
                i+=1 ; fitness +=1; visite[i,j]=true
            elseif (ind[i,j] == 4) && (carte[2*i,2*j-1] ==' ') && (visite[i,j-1] == false) # Ouest
                j-=1 ; fitness +=1; visite[i,j]=true
            else
                avance = false
            end
        realisable = j==11
    end
    return fitness , realisable
end

# ----------------------------------------------------------------------------------------------------
# Selection d'un parent dans la population
function selectionParent(pop)

    # A ECRIRE

    ind1, fitness1, realisable1 = pop[rand(1:length(pop))]
    p1 = copy(ind1)
    return p1
 end

# ----------------------------------------------------------------------------------------------------
# crossover entre deux individus
function crossover(p1,p2)

    # A ECRIRE

    return p1,p2
end

# ----------------------------------------------------------------------------------------------------
# mutation d'un individu
function mutation( carte , individu )

    # A ECRIRE

    ind = copy(individu)
    return ind
end

# ----------------------------------------------------------------------------------------------------
# Selectionne un individu survivant entre deux individus
function survivantEnfant( carte, e , em )

    # A ECRIRE

    return e
end

# ----------------------------------------------------------------------------------------------------
# Recupere la nouvelle generation comme population de base
function changeGeneration(newGen, popSize)

    pop = Vector(undef, popSize)
    NbRealisable = 0
    maxFitness = 0
    minFitness = 100

    for i=1:popSize
        pop[i] = pop!(newGen)
        ind, fitness, realisable = pop[i]
        if realisable
            NbRealisable +=1
            minFitness = min(fitness,minFitness)
        end
        maxFitness = max(fitness,maxFitness)
    end
    println("Nbre Realisable = ", NbRealisable, " minFitnessRealisable = ", minFitness, " maxFitness = ", maxFitness)
    return pop
end

# ----------------------------------------------------------------------------------------------------
# Identifie deux individus elites membres de la population
function IdentifieMeilleur( pop , popSize )

    iElite1=0
    existe = false

    # Premier meilleurs
    minFitness=100
    for i=1:popSize
        ind, fitness, realisable = pop[i]
        if (realisable) && (fitness < minFitness)
            iElite1 = i
            minFitness = fitness
            existe = true
        end
    end
    if existe
        @printf(" Longueur de meilleure solution trouvee : %3d\n",minFitness)
    else
        println(" Pas de solution realisable trouvee !")
    end
    return iElite1
end

# ----------------------------------------------------------------------------------------------------
# Rock and Roll
function main()

    # partie "creation du labyrinthe"
    n = 10
    carte = Matrix{Char}( undef, 21, 21 )
    carte = construireLabyrinthe()
    afficheLabyrinthe( carte )
    contrainteSoftLabyrinthe( carte )
    afficheLabyrinthe( carte )

    # partie "algorithme genetique"
    popSize = 100 # multiple de 2
    pop = creerPopulation( n , popSize , carte )

    for generation=1:20
        newGen = []
        @printf("[%5d]  ",generation)

        # reproduction entre individus
        for reproduction=1:Int(popSize / 2)
            p1 = selectionParent( pop )
            p2 = selectionParent( pop )
            e1, e2 = crossover( p1, p2 )
            e1 = survivantEnfant( carte, e1 , mutation( carte, e1 ) )
            e2 = survivantEnfant( carte, e2 , mutation( carte, e2 ) )
            fitness, realisable = evaluerIndividu( n, carte, e1)
            push!( newGen, (e1, fitness, realisable) )
            fitness, realisable = evaluerIndividu( n, carte, e2)
            push!( newGen, (e2, fitness, realisable) )
        end # reproduction
        pop = changeGeneration( newGen, popSize )
     end # generation
     iBest = IdentifieMeilleur( pop , popSize )
     if iBest !=0
         afficheIndividu(pop[iBest],carte)
     end
end

max_stars_repo_name: xgandibleux/metaMAZEstu

max_issues_repo_name: xgandibleux/metaMAZEstu

max_forks_repo_name: xgandibleux/metaMAZEstu


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct HVDCLine <: DCBranch
        name::String
        available::Bool
        active_power_flow::Float64
        arc::Arc
        active_power_limits_from::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        active_power_limits_to::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        reactive_power_limits_from::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        reactive_power_limits_to::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        loss::NamedTuple{(:l0, :l1), Tuple{Float64, Float64}}
        services::Vector{Service}
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::InfrastructureSystemsInternal
    end

a High voltage DC line.

# Arguments
- `name::String`
- `available::Bool`
- `active_power_flow::Float64`
- `arc::Arc`
- `active_power_limits_from::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `active_power_limits_to::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `reactive_power_limits_from::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `reactive_power_limits_to::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `loss::NamedTuple{(:l0, :l1), Tuple{Float64, Float64}}`
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct HVDCLine <: DCBranch
    name::String
    available::Bool
    active_power_flow::Float64
    arc::Arc
    active_power_limits_from::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    active_power_limits_to::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    reactive_power_limits_from::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    reactive_power_limits_to::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    loss::NamedTuple{(:l0, :l1), Tuple{Float64, Float64}}
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function HVDCLine(name, available, active_power_flow, arc, active_power_limits_from, active_power_limits_to, reactive_power_limits_from, reactive_power_limits_to, loss, services=Device[], ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), )
    HVDCLine(name, available, active_power_flow, arc, active_power_limits_from, active_power_limits_to, reactive_power_limits_from, reactive_power_limits_to, loss, services, ext, time_series_container, InfrastructureSystemsInternal(), )
end

function HVDCLine(; name, available, active_power_flow, arc, active_power_limits_from, active_power_limits_to, reactive_power_limits_from, reactive_power_limits_to, loss, services=Device[], ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), internal=InfrastructureSystemsInternal(), )
    HVDCLine(name, available, active_power_flow, arc, active_power_limits_from, active_power_limits_to, reactive_power_limits_from, reactive_power_limits_to, loss, services, ext, time_series_container, internal, )
end

# Constructor for demo purposes; non-functional.
function HVDCLine(::Nothing)
    HVDCLine(;
        name="init",
        available=false,
        active_power_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        active_power_limits_from=(min=0.0, max=0.0),
        active_power_limits_to=(min=0.0, max=0.0),
        reactive_power_limits_from=(min=0.0, max=0.0),
        reactive_power_limits_to=(min=0.0, max=0.0),
        loss=(l0=0.0, l1=0.0),
        services=Device[],
        ext=Dict{String, Any}(),
        time_series_container=InfrastructureSystems.TimeSeriesContainer(),
    )
end

InfrastructureSystems.get_name(value::HVDCLine) = value.name
"""Get [`HVDCLine`](@ref) `available`."""
get_available(value::HVDCLine) = value.available
"""Get [`HVDCLine`](@ref) `active_power_flow`."""
get_active_power_flow(value::HVDCLine) = get_value(value, value.active_power_flow)
"""Get [`HVDCLine`](@ref) `arc`."""
get_arc(value::HVDCLine) = value.arc
"""Get [`HVDCLine`](@ref) `active_power_limits_from`."""
get_active_power_limits_from(value::HVDCLine) = get_value(value, value.active_power_limits_from)
"""Get [`HVDCLine`](@ref) `active_power_limits_to`."""
get_active_power_limits_to(value::HVDCLine) = get_value(value, value.active_power_limits_to)
"""Get [`HVDCLine`](@ref) `reactive_power_limits_from`."""
get_reactive_power_limits_from(value::HVDCLine) = get_value(value, value.reactive_power_limits_from)
"""Get [`HVDCLine`](@ref) `reactive_power_limits_to`."""
get_reactive_power_limits_to(value::HVDCLine) = get_value(value, value.reactive_power_limits_to)
"""Get [`HVDCLine`](@ref) `loss`."""
get_loss(value::HVDCLine) = value.loss
"""Get [`HVDCLine`](@ref) `services`."""
get_services(value::HVDCLine) = value.services
"""Get [`HVDCLine`](@ref) `ext`."""
get_ext(value::HVDCLine) = value.ext

InfrastructureSystems.get_time_series_container(value::HVDCLine) = value.time_series_container
"""Get [`HVDCLine`](@ref) `internal`."""
get_internal(value::HVDCLine) = value.internal

InfrastructureSystems.set_name!(value::HVDCLine, val) = value.name = val
"""Set [`HVDCLine`](@ref) `available`."""
set_available!(value::HVDCLine, val) = value.available = val
"""Set [`HVDCLine`](@ref) `active_power_flow`."""
set_active_power_flow!(value::HVDCLine, val) = value.active_power_flow = val
"""Set [`HVDCLine`](@ref) `arc`."""
set_arc!(value::HVDCLine, val) = value.arc = val
"""Set [`HVDCLine`](@ref) `active_power_limits_from`."""
set_active_power_limits_from!(value::HVDCLine, val) = value.active_power_limits_from = val
"""Set [`HVDCLine`](@ref) `active_power_limits_to`."""
set_active_power_limits_to!(value::HVDCLine, val) = value.active_power_limits_to = val
"""Set [`HVDCLine`](@ref) `reactive_power_limits_from`."""
set_reactive_power_limits_from!(value::HVDCLine, val) = value.reactive_power_limits_from = val
"""Set [`HVDCLine`](@ref) `reactive_power_limits_to`."""
set_reactive_power_limits_to!(value::HVDCLine, val) = value.reactive_power_limits_to = val
"""Set [`HVDCLine`](@ref) `loss`."""
set_loss!(value::HVDCLine, val) = value.loss = val
"""Set [`HVDCLine`](@ref) `services`."""
set_services!(value::HVDCLine, val) = value.services = val
"""Set [`HVDCLine`](@ref) `ext`."""
set_ext!(value::HVDCLine, val) = value.ext = val

InfrastructureSystems.set_time_series_container!(value::HVDCLine, val) = value.time_series_container = val

ext: jl


text:
using Pkg
pkg"activate .."

using Documenter
using DemoPackage_ak

##makedocs(
#    sitename = "DemoPackage_ak",
 #   Authors = "Ankit DS",
 #   format = Documenter.HTML(),
 #   modules = [DemoPackage_ak]
#)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "https://github.com/ankit48365/DemoPackage_ak.jl"
)=#

DEPLOYDOCS = (get(ENV, "DemoPackage_ak", nothing) == "true")

makedocs(
    sitename = "Demo Package Test by Ankit",
    format = Documenter.HTML(
        prettyurls = DEPLOYDOCS,
    ),
    authors = "Ankit DS",
    pages = ["index.md"]
)

# these were all available --pages = ["index.md", "semaphores.md", "sharedmemory.md", "reference.md"]

if DEPLOYDOCS
    deploydocs(
        repo = "https://github.com/ankit48365/DemoPackage_ak.jl.git",
    )
end


text:
module MakieCore

using Observables
using Observables: to_value
using Base: RefValue

include("types.jl")
include("attributes.jl")
include("recipes.jl")
include("basic_plots.jl")
include("conversion.jl")

end


text:
#=
Aodwt.jl
2019-02-23 Jeff Fessler, University of Michigan
=#

export Aodwt

using LinearMapsAA: LinearMapAA, LinearMapAM, LinearMapAO
using Wavelets: dwt!, idwt!, wavelet, WT

"""
    A, levels, mfun = Aodwt(dims ; level::Int=3, wt=wavelet(WT.haar))

Create orthogonal discrete wavelet transform (ODWT) `LinearMapAA`

in
- `dims::Dims` tuple of dimensions

option
- `level::Int` # of levels; default 3
- `wt` wavelet transform type (see `Wavelets` package); default Haar
- `operator::Bool=true` default to `LinearMapAO`

out
- `A` a `LinearMapAX` object
- `scales` array of size `dims` showing the scale of each coefficient
which is useful when imposing scale-dependent regularization
- `mfun` convenience function for A*X when X is a Matrix or Array (not vector)

2019-02-23 Jeff Fessler, University of Michigan
"""
function Aodwt(
	dims::Dims ;
	T::DataType = Float32,
	level::Int = 3,
	wt = wavelet(WT.haar),
	operator::Bool = true, # !
)

	function mfunA(lev)
		# todo: avoiding convert involves `dwt!` limitation; see:
		# https://github.com/JuliaDSP/Wavelets.jl/issues/53
		# https://github.com/JuliaDSP/Wavelets.jl/pull/54
	#	forw!(y,x) = dwt!(y, convert(AbstractArray{T},x), wt, level)
	#	back!(x,y) = idwt!(x, convert(AbstractArray{T},y), wt, level)
		forw!(y,x) = dwt!(y, x, wt, lev)
		back!(x,y) = idwt!(x, y, wt, lev)

		if operator
			mfun = (A, x) -> A * x
			return mfun, LinearMapAA(forw!, back!,
				(prod(dims), prod(dims)) ;
				prop = (wt=wt, level=lev), T=T,
				operator = true, idim=dims, odim=dims,
			)
		else
			mfun = (A, x) -> reshape(A * vec(x), dims)
			return mfun, LinearMapAA(
				(y,x) -> vec(forw!(reshape(y, dims), reshape(x, dims))),
				(x,y) -> vec(back!(reshape(x, dims), reshape(y, dims))),
			#	x -> vec(dwt(reshape(x, dims), wt, level)),
			#	y -> vec(idwt(reshape(y, dims), wt, level)),
				(prod(dims), prod(dims)) ;
				prop = (wt=wt, level=lev), T=T,
			)
		end
	end

	mfun, A = mfunA(level)

	scales = zeros(dims)
	for il=1:level
		_,Al = mfunA(il)
		tmp = mfun(Al, ones(dims)) .== 0
		scales += il * (tmp .& (scales .== 0))
	end

	return A, scales, mfun
end


text:
using OMEinsum
# tensor for classical 2-d model
"""
    tensorfromclassical(h::Matrix)

given a classical 2-body hamiltonian `h`, return the corresponding tensor
for use in e.g. `trg` for a two-dimensional square-lattice.

# Example
```julia
julia> model_tensor(Ising(),β) ≈ tensorfromclassical([β -β; -β β])

true
```
"""
function tensorfromclassical(ham::Matrix)
    wboltzmann = exp.(ham)
    q = sqrt(wboltzmann)
    ein"ij,ik,il,im -> jklm"(q,q,q,q)
end

"""
    model_tensor(::Ising,β)

return the isingtensor at inverse temperature `β` for a two-dimensional
square lattice tensor-network.
"""
function model_tensor(model::Ising)
    β = model.β
    a = reshape(Float64[1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1], 2,2,2,2)
    cβ, sβ = sqrt(cosh(β)), sqrt(sinh(β))
    q = 1/sqrt(2) * [cβ+sβ cβ-sβ; cβ-sβ cβ+sβ]
    ein"abcd,ai,bj,ck,dl -> ijkl"(a,q,q,q,q)
end

"""
    mag_tensor(::Ising,β)

return the operator for the magnetisation at inverse temperature `β`
at a site in the two-dimensional ising model on a square lattice in tensor-network form.
"""
function mag_tensor(model::Ising)
    β = model.β
    a = reshape(Float64[1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 -1] , 2,2,2,2)
    cβ, sβ = sqrt(cosh(β)), sqrt(sinh(β))
    q = 1/sqrt(2) * [cβ+sβ cβ-sβ; cβ-sβ cβ+sβ]
    ein"abcd,ai,bj,ck,dl -> ijkl"(a,q,q,q,q)
end

"""
    energy_tensor(::Ising,β)

return the operator for the energy at inverse temperature `β`
at a site in the two-dimensional ising model on a square lattice in tensor-network form.
"""
function energy_tensor(model::Ising)
    β = model.β
    ham = [-1 1;1 -1]
    wboltzmann = exp.(-β .* ham)
    wenergy = ham .* wboltzmann
    qb = sqrt(wboltzmann)
    ein"as,si,bi,ci,di -> abcd"(qb^(-1),wenergy,qb,qb,qb)
end


text:
using FdeSolver
using Plots

## inputs
tSpan = [0, 50] # time span

h = 0.1 # time step

N = 10 # number of species

β = ones(N) # order of derivatives

X0 = 2*rand(N) # initial abundances

## System definition

# parametrisation
par = [2,
      2*rand(N),
      rand(N),
      4*rand(N,N),
       N]

function F(t, n, β, x, par)
    l = par[1] # Hill coefficient
    b = par[2] # growth rates
    k = par[3] # death rates
    K = par[4] # inhibition matrix
    N = par[5] # number of species

# ODE
    Fun = zeros(N)
    for i in 1:N
    # inhibition functions
    f = prod(K[i,1:end .!= i] .^ l ./
             (K[i,1:end .!= i] .^ l .+ x[n, 1:end .!= i] .^l))
    # System of equations
    Fun[i] = x[n, i] .* (b[i] .* f .- k[i] .* x[n, i])
    end

    return Fun

end

# numerical solution
t, Xapp = FDEsolver(F, tSpan, X0, β, nothing, par, h = h, nc = 3, tol = 10^(-8))

# plot
plot(t, Xapp, linewidth = 5,
     title = "Dynamic of microbial intercation model",
     xaxis = "Time (t)")
     yaxis!("Log Abundance", :log10, minorgrid = true)


text:
using LogicCircuits.IO: constant

#####################
# Compilers to ProbabilisticCircuits data structures starting from already parsed line objects
#####################

# reuse some internal infrastructure of LogicCircuits' IO module
using LogicCircuits.IO: CircuitFormatLines, CircuitFormatLine, VtreeFormatLines, CircuitHeaderLine, UnweightedLiteralLine, WeightedLiteralLine, DecisionLine, LCElement, BiasLine, WeightedNamedConstantLine, PSDDElement, CircuitCommentLine, ID,
compile_smooth_struct_logical_m, compile_smooth_logical_m

"""
Compile lines into a probabilistic circuit.
"""
function compile_prob(lines::CircuitFormatLines)::ProbΔ
    # first compile a logical circuit
    logical_circuit, id2lognode = compile_smooth_logical_m(lines)
    decorate_prob(lines, logical_circuit, id2lognode)
end

"""
Compile lines into a logistic circuit.
"""
function compile_logistic(lines::CircuitFormatLines, classes::Int)::LogisticΔ
    # first compile a logical circuit
    logical_circuit, id2lognode = compile_smooth_logical_m(lines)
    decorate_logistic(lines, logical_circuit, classes, id2lognode)
end

"""
Compile circuit and vtree lines into a structured probabilistic circuit (one whose logical circuit origin is structured).
"""
function compile_struct_prob(circuit_lines::CircuitFormatLines, vtree_lines::VtreeFormatLines)
    logical_circuit, vtree, id2vtree, id2lognode = compile_smooth_struct_logical_m(circuit_lines, vtree_lines)
    prob_circuit = decorate_prob(circuit_lines, logical_circuit, id2lognode)
    return prob_circuit, vtree
end

function decorate_prob(lines::CircuitFormatLines, logical_circuit::LogicalΔ, id2lognode::Dict{ID,<:LogicalΔNode})::ProbΔ
    # set up cache mapping logical circuit nodes to their probabilistic decorator
    lognode2probnode = ProbCache()
    # build a corresponding probabilistic circuit
    prob_circuit = ProbΔ(logical_circuit,lognode2probnode)
    # map from line node ids to probabilistic circuit nodes
    id2probnode(id) = lognode2probnode[id2lognode[id]]

    # go through lines again and update the probabilistic circuit node parameters

    function compile(ln::CircuitFormatLine)
        error("Compilation of line $ln into probabilistic circuit is not supported")
    end
    function compile(::Union{CircuitHeaderLine,CircuitCommentLine,UnweightedLiteralLine})
        # do nothing
    end
    function compile(ln::WeightedNamedConstantLine)
        @assert constant(ln) == true
        node = id2probnode(ln.node_id)::Prob⋁
        node.log_thetas .= [ln.weight, log1p(-exp(ln.weight)) ]
    end
    function compile(ln::DecisionLine{<:PSDDElement})
        node = id2probnode(ln.node_id)::Prob⋁
        node.log_thetas .= [x.weight for x in ln.elements]
    end
    for ln in lines
        compile(ln)
    end

    prob_circuit
end

function decorate_logistic(lines::CircuitFormatLines, logical_circuit::LogicalΔ,
                            classes::Int, id2lognode::Dict{ID,<:LogicalΔNode})::LogisticΔ

    # set up cache mapping logical circuit nodes to their logistic decorator
    log2logistic = LogisticCache()
    # build a corresponding probabilistic circuit
    logistic_circuit = LogisticΔ(logical_circuit, classes, log2logistic)
    # map from line node ids to probabilistic circuit nodes
    id2logisticnode(id) = log2logistic[id2lognode[id]]

    # go through lines again and update the probabilistic circuit node parameters

    function compile(ln::CircuitFormatLine)
        error("Compilation of line $ln into logistic circuit is not supported")
    end
    function compile(::Union{CircuitHeaderLine,CircuitCommentLine,UnweightedLiteralLine})
        # do nothing
    end

    function compile(ln::CircuitHeaderLine)
        # do nothing
    end

    function compile(ln::WeightedLiteralLine)
        node = id2logisticnode(ln.node_id)::Logistic⋁
        node.thetas[1, :] .= ln.weights
    end

    function compile(ln::DecisionLine{<:LCElement})
        node = id2logisticnode(ln.node_id)::Logistic⋁
        for (ind, elem) in enumerate(ln.elements)
            node.thetas[ind, :] .= elem.weights
        end
    end

    function compile(ln::BiasLine)
        node = id2logisticnode(ln.node_id)::Logistic⋁
        # @assert length(node.thetas) == 1
        node.thetas[1,:] .= ln.weights
    end

    for ln in lines
        compile(ln)
    end

    logistic_circuit
end


text:
using Test, Statistics
using EfficientGlobalOptimization
using KrigingModel, GaussianProcesses

include("test_testfunction.jl")
include("test_ego.jl")
include("test_cons.jl")
include("test_pareto.jl")


text:
using ExactOptimalTransport

using Distances
using PythonOT: PythonOT
using Tulip
using MathOptInterface
using Distributions
using HCubature

using LinearAlgebra
using Random
using SparseArrays

const MOI = MathOptInterface
const POT = PythonOT

Random.seed!(100)

@testset "exact.jl" begin
    @testset "Earth-Movers Distance" begin
        M = 200
        N = 250
        μ = normalize!(rand(M), 1)
        ν = normalize!(rand(N), 1)

        @testset "example" begin
            # create random cost matrix
            C = pairwise(SqEuclidean(), rand(1, M), rand(1, N); dims=2)

            # compute optimal transport map and cost with POT
            pot_P = POT.emd(μ, ν, C)
            pot_cost = POT.emd2(μ, ν, C)

            # compute optimal transport map and cost with Tulip
            lp = Tulip.Optimizer()
            P = emd(μ, ν, C, lp)
            @test size(C) == size(P)
            @test MOI.get(lp, MOI.TerminationStatus()) == MOI.OPTIMAL
            @test maximum(abs, P .- pot_P) < 1e-2

            lp = Tulip.Optimizer()
            cost = emd2(μ, ν, C, lp)
            @test dot(C, P) ≈ cost atol = 1e-5
            @test MOI.get(lp, MOI.TerminationStatus()) == MOI.OPTIMAL
            @test cost ≈ pot_cost atol = 1e-5
        end

        @testset "pre-computed plan" begin
            # create random cost matrix
            C = pairwise(SqEuclidean(), rand(1, M), rand(1, N); dims=2)

            # compute optimal transport map
            P = emd(μ, ν, C, Tulip.Optimizer())

            # do not use μ and ν to ensure that provided map is used
            cost = emd2(similar(μ), similar(ν), C, Tulip.Optimizer(); plan=P)
            @test cost ≈ emd2(μ, ν, C, Tulip.Optimizer())
        end

        # https://github.com/JuliaOptimalTransport/OptimalTransport.jl/issues/71
        @testset "cost matrix with integers" begin
            C = pairwise(SqEuclidean(), rand(1:10, 1, M), rand(1:10, 1, N); dims=2)
            emd2(μ, ν, C, Tulip.Optimizer())
        end
    end

    @testset "1D Optimal Transport for Convex Cost" begin
        @testset "continuous distributions" begin
            # two normal distributions (has analytical solution)
            μ = Normal(randn(), rand())
            ν = Normal(randn(), rand())

            # compute OT plan
            γ = ot_plan(sqeuclidean, μ, ν)
            x = randn()
            @test γ(x) ≈ quantile(ν, cdf(μ, x))

            # compute OT cost
            c = ot_cost(sqeuclidean, μ, ν)
            @test c ≈ (mean(μ) - mean(ν))^2 + (std(μ) - std(ν))^2

            # do not use ν to ensure that the provided plan is used
            @test ot_cost(sqeuclidean, μ, Normal(randn(), rand()); plan=γ) ≈ c
        end

        @testset "semidiscrete case" begin
            μ = Normal(randn(), rand())
            νprobs = normalize!(rand(30), 1)
            ν = Categorical(νprobs)

            # compute OT plan
            γ = ot_plan(euclidean, μ, ν)
            x = randn()
            @test γ(x) ≈ quantile(ν, cdf(μ, x))

            # compute OT cost, without and with provided plan
            # do not use ν in the second case to ensure that the provided plan is used
            c = ot_cost(euclidean, μ, ν)
            @test ot_cost(euclidean, μ, Categorical(reverse(νprobs)); plan=γ) ≈ c

            # check that OT cost is consistent with OT cost of a discretization
            m = 500
            xs = rand(μ, m)
            μdiscrete = fill(1 / m, m)
            C = pairwise(Euclidean(), xs', (1:length(νprobs))'; dims=2)
            c2 = emd2(μdiscrete, νprobs, C, Tulip.Optimizer())
            @test c2 ≈ c rtol = 1e-1
        end

        @testset "discrete case" begin
            # random source and target marginal
            m = 30
            μprobs = normalize!(rand(m), 1)
            μsupport = randn(m)
            μ = DiscreteNonParametric(μsupport, μprobs)

            n = 50
            νprobs = normalize!(rand(n), 1)
            νsupport = randn(n)
            ν = DiscreteNonParametric(νsupport, νprobs)

            # compute OT plan
            γ = @inferred(ot_plan(euclidean, μ, ν))
            @test γ isa SparseMatrixCSC
            @test size(γ) == (m, n)
            @test vec(sum(γ; dims=2)) ≈ μ.p
            @test vec(sum(γ; dims=1)) ≈ ν.p

            # consistency checks
            I, J, W = findnz(γ)
            @test all(w > zero(w) for w in W)
            @test sum(W) ≈ 1
            @test sort(unique(I)) == 1:m
            @test sort(unique(J)) == 1:n
            @test sort(I .+ J) == 2:(m + n)

            # compute OT cost
            c = @inferred(ot_cost(euclidean, μ, ν))

            # compare with computation with explicit cost matrix
            # DiscreteNonParametric sorts the support automatically, here we have to sort
            # manually
            C = pairwise(Euclidean(), μsupport', νsupport'; dims=2)
            c2 = emd2(μprobs, νprobs, C, Tulip.Optimizer())
            @test c2 ≈ c rtol = 1e-5

            # compare with POT
            # disabled currently since https://github.com/PythonOT/POT/issues/169 causes bounds
            # error
            # @test γ ≈ POT.emd_1d(μ.support, ν.support; a=μ.p, b=μ.p, metric="euclidean")
            # @test c ≈ POT.emd2_1d(μ.support, ν.support; a=μ.p, b=μ.p, metric="euclidean")

            # do not use the probabilities of μ and ν to ensure that the provided plan is
            # used
            μ2 = DiscreteNonParametric(μsupport, reverse(μprobs))
            ν2 = DiscreteNonParametric(νsupport, reverse(νprobs))
            c2 = @inferred(ot_cost(euclidean, μ2, ν2; plan=γ))
            @test c2 ≈ c
            c2 = @inferred(ot_cost(euclidean, μ2, ν2; plan=Matrix(γ)))
            @test c2 ≈ c
        end
    end

    @testset "Multivariate Gaussians" begin
        @testset "translation with constant covariance" begin
            m = randn(100)
            τ = rand(100)
            Σ = Matrix(Hermitian(rand(100, 100) + 100I))
            μ = MvNormal(m, Σ)
            ν = MvNormal(m .+ τ, Σ)
            @test ot_cost(SqEuclidean(), μ, ν) ≈ norm(τ)^2

            x = rand(100, 10)
            T = ot_plan(SqEuclidean(), μ, ν)
            @test pdf(ν, mapslices(T, x; dims=1)) ≈ pdf(μ, x)
        end

        @testset "comparison to grid approximation" begin
            μ = MvNormal([0, 0], [1 0; 0 2])
            ν = MvNormal([10, 10], [2 0; 0 1])
            # Constructing circular grid approximation
            # Angular grid step
            θ = collect(0:0.2:(2π))
            θx = cos.(θ)
            θy = sin.(θ)
            # Radius grid step
            δ = collect(0:0.2:1)
            μsupp = [0.0 0.0]
            νsupp = [10.0 10.0]
            for i in δ[2:end]
                a = [θx .* i θy .* i * 2]
                b = [θx .* i * 2 θy .* i] .+ [10 10]
                μsupp = vcat(μsupp, a)
                νsupp = vcat(νsupp, b)
            end

            # Create discretized distribution
            μprobs = normalize!(pdf(μ, μsupp'), 1)
            νprobs = normalize!(pdf(ν, νsupp'), 1)
            C = pairwise(SqEuclidean(), μsupp', νsupp'; dims=2)
            @test emd2(μprobs, νprobs, C, Tulip.Optimizer()) ≈ ot_cost(SqEuclidean(), μ, ν) rtol =
                1e-3

            # Use hcubature integration to perform ``\\int c(x,T(x)) d\\mu``
            T = ot_plan(SqEuclidean(), μ, ν)
            c_hcubature, _ = hcubature([-10, -10], [10, 10]) do x
                return sqeuclidean(x, T(x)) * pdf(μ, x)
            end
            @test ot_cost(SqEuclidean(), μ, ν) ≈ c_hcubature rtol = 1e-3
        end
    end
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    Base.precompile(Tuple{Core.kwftype(typeof(FixedEffects.normalize!)),NamedTuple{(:tol, :maxiter),Tuple{Float64,Int64}},typeof(FixedEffects.normalize!),Array{Array{Float64,1},1},Array{FixedEffects.FixedEffect,1}})
    Base.precompile(Tuple{Core.kwftype(typeof(FixedEffects.normalize!)),NamedTuple{(:tol, :maxiter),Tuple{Float64,Int64}},typeof(FixedEffects.normalize!),Array{Array{Float64,1},1},Array{FixedEffect{Array{UInt32,1},UnitWeights{Float64}},1}})
    Base.precompile(Tuple{typeof(FixedEffects.full),Array{Array{Float64,1},1},Array{FixedEffects.FixedEffect,1}})

    Base.precompile(Tuple{Type{FixedEffects.AbstractFixedEffectSolver{Float64}},Array{FixedEffects.FixedEffect,1},Weights{Float64,Float64,Array{Float64,1}},Type{Val{:cpu}}})
    Base.precompile(Tuple{Type{FixedEffects.AbstractFixedEffectSolver{Float64}},Array{FixedEffects.FixedEffect,1},UnitWeights{Float64},Type{Val{:cpu}}})
    Base.precompile(Tuple{Type{FixedEffects.AbstractFixedEffectSolver{Float32}},Array{FixedEffects.FixedEffect,1},StatsBase.Weights{Float64,Float64,Array{Float64,1}},Type{Val{:cpu}},Int64})
    Base.precompile(Tuple{typeof(FixedEffects.update_weights!),FixedEffects.FixedEffectSolverCPU{Float64},Weights{Int64,Int64,Array{Int64,1}}})

    Base.precompile(Tuple{Core.kwftype(typeof(FixedEffects.solve_coefficients!)),NamedTuple{(:tol, :maxiter),Tuple{Float64,Int64}},typeof(FixedEffects.solve_coefficients!),Array{Float64,1},FixedEffects.FixedEffectSolverCPU{Float64}})

    Base.precompile(Tuple{typeof(solve_residuals!),Array{Float64,1},Array{FixedEffect{Array{UInt32,1},UnitWeights{Float64}},1}})
    Base.precompile(Tuple{typeof(solve_residuals!),Array{Float64,2},Array{FixedEffect{Array{UInt32,1},UnitWeights{Float64}},1}})
    Base.precompile(Tuple{Core.kwftype(typeof(FixedEffects.solve_residuals!)),NamedTuple{(:maxiter, :tol),Tuple{Int64,Float64}},typeof(solve_residuals!),Array{Float64,1},FixedEffects.FixedEffectSolverCPU{Float64}})

    Base.precompile(Tuple{Core.kwftype(typeof(FixedEffects.lsmr!)),NamedTuple{(:atol, :btol, :maxiter),Tuple{Float64,Float64,Int64}},typeof(FixedEffects.lsmr!),FixedEffects.FixedEffectCoefficients{Array{Float64,1}},FixedEffects.FixedEffectLinearMapCPU{Float64},Array{Float64,1},FixedEffects.FixedEffectCoefficients{Array{Float64,1}},FixedEffects.FixedEffectCoefficients{Array{Float64,1}},FixedEffects.FixedEffectCoefficients{Array{Float64,1}}})
    Base.precompile(Tuple{Core.kwftype(typeof(FixedEffects.lsmr!)),NamedTuple{(:atol, :btol, :maxiter),Tuple{Float64,Float64,Int64}},typeof(FixedEffects.lsmr!),FixedEffects.FixedEffectCoefficients{Array{Float32,1}},FixedEffects.FixedEffectLinearMapCPU{Float32},Array{Float32,1},FixedEffects.FixedEffectCoefficients{Array{Float32,1}},FixedEffects.FixedEffectCoefficients{Array{Float32,1}},FixedEffects.FixedEffectCoefficients{Array{Float32,1}}})

    Base.precompile(Tuple{typeof(group),Array{Union{Missing, String},1}})
    Base.precompile(Tuple{typeof(group),Array{Union{Missing, String},1},Vararg{Any,N} where N})
end


text:
""" plotvaf(VAF, [subclonefreqs=[]]; <keyword arguments>)
"""
# plotvaf

@userplot PlotVAF

@recipe function f(pv::PlotVAF; cumulative = false, fstep = 0.01, sampled=true)
    VAFresult = pv.args[1]
    VAF = sampled ? VAFresult.sampledVAF : VAFresult.trueVAF
    df = gethist(VAF, fstep = fstep)
    VAF = (df[!,:VAF] .*2 .- fstep) ./ 2 #set x values to middle of bins
    freq = cumulative ? df[:,:cumfreq] : df[:,:freq]

    if get(plotattributes, :yscale, nothing) == :log10
        VAF, freq = VAF[freq .> 0], freq[freq .>0]
    end
    yguide = cumulative ? "Cumulative number of mutations" : "Number of mutations"

    @series begin
        seriestype --> :bar
        linecolor --> :white
        linewidth --> 0.2
        fillcolor --> :darkslategrey
        markerstrokecolor --> :white
        VAF, freq
    end

    if length(VAFresult.subclonefreq) > 0
        subclonefreq = VAFresult.subclonefreq ./VAFresult.input.ploidy * VAFresult.cellularity

        @series begin
            seriestype --> :vline
            legend --> false
            fillcolor --> :darkred
            linewidth --> 2
            legend --> false
            grid --> false
            subclonefreq
        end
    end
    yguide --> yguide
    xguide --> "VAF"
    legend --> false
    grid --> false
    ()
end

""" plotinversevaf(VAF; <keyword arguments>)
"""
plotvaf

@userplot PlotInverseVAF

@recipe function f(pv::PlotInverseVAF; fmin = 0.12, fmax = 0.24, fit_fmax=fmax, sampled=true,
                    fstep = 0.001, fitcoef = nothing, cumulative = true, dataseries = :line)

    VAFresult = pv.args[1]
    VAF = sampled ? VAFresult.sampledVAF : VAFresult.trueVAF
    df = gethist(VAF, fmin = fmin, fmax = fmax, fstep = fstep)
    VAF = df[!,:VAF]
    freq = cumulative ? df[!,:cumfreq] : df[!,:freq]
    yguide = cumulative ? "Cumulative numberof mutations" : "Number of mutations"

    @series begin
        yguide --> yguide
        xguide --> "Inverse VAF"
        seriestype --> dataseries
        1 ./ VAF, freq
    end

    if fitcoef !== nothing
        if length(fitcoef) == 1
            m = fitcoef[1]
            c = - m/fit_fmax
        else
            c, m = fitcoef
        end
        x, y = 1 ./ VAF, m ./ VAF .+ c
        x, y = x[y .>= 0], y[y .>= 0]
        @series begin
            yguide --> yguide
            xguide --> "Inverse VAF"
            seriestype --> :line
            linecolor --> :red
            linestyle --> :dash
            x, y
        end
    end

    yguide --> yguide
    xguide --> "Inverse VAF"
    fvals = [df[end, :VAF], 2/(1/df[1, :VAF]+ 1/df[end, :VAF]), df[1, :VAF]]
    xticks --> (1 ./ fvals, map(x -> "1/" * x, string.(round.(fvals, digits = 2))))
    legend --> false
    grid --> false
    ()
end

@recipe function f(output::ModuleTracker)
    @series begin
        yguide --> "Population size"
        xguide --> "Time"
        seriestype --> :line
        output.tvec, output.Nvec
    end
end

@recipe function f(multisim::MultiSimulation; plottype=:popsize, tstep=nothing)
    if plottype == :modulesize
        yguide --> "Module size"
        xguide --> "Time"
        for moduletracker in multisim
            @series begin
                seriestype --> :line
                legend --> false
                moduletracker.tvec, moduletracker.Nvec
            end
        end
    elseif plottype == :popsize
        yguide --> "Number of modules"
        xguide --> "Time"
        @series begin
            seriestype --> :line
            legend --> false
            newmoduletimes(multisim), 1:length(multisim)
        end
    elseif plottype == :cellpopsize
        tstep = isnothing(tstep) ? 1 / multisim.input.bdrate : tstep
        yguide --> "Number of cells"
        xguide --> "Time"
        @series begin
            seriestype --> :line
            legend --> false
            cellpopulationsize(multisim, tstep)
        end

    end
end

@userplot PairwisePlot

@recipe function f(pp::PairwisePlot; sampleids=nothing)
    if length(pp.args) == 1
        z = pp.args[1]
        n = isnothing(sampleids) ? size(z)[1] : length(sampleids)
        x, y = 1:n, 1:n
    elseif length(pp.args) == 3
        x, y, z = pp.args
    end
    if sampleids !== nothing
        z = z[sampleids, sampleids]
    end
    yguide --> "Module"
    xguide --> "Module"
    title --> "Pairwise fixed differences"
    @series begin
        grid --> false
        seriestype := :heatmap
        seriescolor --> cgrad(:bone, rev=true)
        xticks --> x
        yticks --> y
        aspectratio --> 1
        size --> (500,510)
        x, y, z
    end
end

@userplot ModuleSharedPlot

@recipe function f(msp::ModuleSharedPlot)
    @series begin
        grid --> false
        seriestype := :bar
        xguide --> "Number of modules"
        yguide --> "Number fixed mutations"
        legend --> false
        msp.args[1]
    end
end

@userplot PairwiseDistributionPlot

@recipe function f(pdp::PairwiseDistributionPlot; sampleids=nothing, samplesize=nothing,
    showclonal=false, xticks=:auto)

    if typeof(pdp.args[1]) <: Dict
        pfddata = reduce(vcat, [fill(key, val) for (key,val) in pdp.args[1]])
    else
        pfdmatrix = pdp.args[1]
        n = size(pfdmatrix)[1]
        if sampleids === nothing && samplesize !== nothing
            sampleids = sample(1:n, samplesize, replace=false)
        end
        if sampleids !== nothing
            pfdmatrix = pfdmatrix[sampleids, sampleids]
            n = length(sampleids)
        end
        pfddata = convert_pfdmatrix_to_vector(pfdmatrix)
    end
    @series begin
        grid --> false
        seriestype := :hist
        xguide --> "Number of pairwise fixed differences"
        yguide --> "Frequency"
        color --> :dimgrey
        linecolor --> :dimgrey
        legend --> false
        titlefontsize -->10
        titlelocation --> :left
        pfddata
    end
    if showclonal && length(pdp.args) > 1
        @series begin
            legend --> false
            grid --> false
            xguide --> "Number of pairwise fixed differences"
            yguide --> "Frequency"
            titlefontsize -->10
            titlelocation --> :left
            seriestype --> :vline
            fillcolor --> :red
            linewidth --> 2
            [mean(pdp.args[2])]
        end
    end
end

function convert_pfdmatrix_to_vector(pfdmatrix)
    n = size(pfdmatrix)[1]
    vals = Int64[]
    for i in 1:n
        for j in i+1:n
            push!(vals, pfdmatrix[j,i])
        end
    end
    return vals
end


text:
struct SomeGroup <: Group end

struct SomeGroupElement <: GroupElement
    elts::Vector{Int} # SomeGroupElement is not isbits anymore
end

@testset "SomeGroup: No interface implemented" begin

    @testset "InterfaceNotImplemented exception" begin
        @test GroupsCore.InterfaceNotImplemented(
            :Castle,
            "Aaaaarghhhhh....",
        ) isa Exception
        ex = GroupsCore.InterfaceNotImplemented(:Castle, "Aaaaarghhhhh....")
        @test ex isa GroupsCore.InterfaceNotImplemented

        @test sprint(showerror, ex) ==
              "Missing method from Castle interface: `Aaaaarghhhhh....`"
    end

    INI = GroupsCore.InterfaceNotImplemented

    @testset "Group Interface" begin

        G = SomeGroup()

        # Iteration
        @test_throws INI eltype(G)
        @test_throws INI iterate(G)
        @test_throws INI iterate(G, 1)

        # Assumption 1: Groups are finite unless claimed otherwise
        @test Base.IteratorSize(G) == Base.HasLength()
        @test Base.isfinite(G)

        Base.IteratorSize(::Type{SomeGroup}) = Base.HasShape{1}()
        @test Base.isfinite(G)

        Base.IteratorSize(::Type{SomeGroup}) = Base.SizeUnknown()
        @test !Base.isfinite(G)

        Base.IteratorSize(::Type{SomeGroup}) = Base.IsInfinite()
        @test !Base.isfinite(G)

        # return to the default:
        Base.IteratorSize(::Type{SomeGroup}) = Base.HasLength()

        # Assumption 2: Groups have generators:
        @test hasgens(G)

        # Group Interface
        @test_throws INI one(G)
        @test_throws INI order(G)
        @test_throws INI gens(G)

        Base.eltype(::Type{SomeGroup}) = SomeGroupElement
        @test_throws INI rand(G, 2)

        @test_throws INI gens(G, 1)
        @test_throws INI ngens(G)

        @test_throws INI GroupsCore.pseudo_rand(G, 2, 2)
    end

    @testset "GroupElem Interface" begin

        g = SomeGroupElement(Int[])

        @test_throws INI parent(g)
        @test_throws INI GroupsCore.parent_type(g)
        @test_throws INI istrulyequal(g, g)

        @test_throws INI hasorder(g)
        @test_throws INI deepcopy(g)

        @test_throws INI inv(g)
        @test_throws INI g * g
    end

end


text:
# # MISMIP+ example
#
# This example demonstrates the Marine Ice Sheet Model Intercomparison (MISMIP) + ice0 experiment (doi: 10.5194/tc-14-2283-2020).
# This intercomparison exercise considers a plan-view ice sheet, in which the grounding line can stabilize on a section of bed which has a locally positive slope in the flow direction.
# Such configurations are theoretically impossible in one horizontal dimension (see doi: 10.1029/2006JF000664), demonstrating the importance of variations in the second dimension for buttressing ice sheets.
# This example demonstrates how to apply boundary conditions in `WAVI.jl`, and control the number of iterations in the velocity solve, and zip the output into a friendly format.

# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# using Pkg
# Pkg.add("https://github.com/RJArthern/WAVI.jl"), Pkg.add(Plots)
using WAVI, Plots, NetCDF

# ## Basal Topography
#
# The MISMIP+ domain is 640km in the x-direction and 80km in the y-direction, centred around $y = 0$.
# The basal topography is given by $z_b = \max [B_x(x) + B_y(y), -720] where $B_x(x)$ is a sixth order, even polynomial and $B_y(y)$ introduces two bumps in the domain:
function mismip_plus_bed(x,y)
    xbar = 300000.0
    b0 = -150.0; b2 = -728.8; b4 = 343.91; b6 = -50.75;
    wc = 24000.0; fc = 4000.0; dc = 500.0;
    bx(x)=b0+b2*(x/xbar)^2+b4*(x/xbar)^4+b6*(x/xbar)^6;
    by(y)= dc*( (1+exp(-2(y-wc)/fc))^(-1) + (1+exp(2(y+wc)/fc))^(-1) );
    b = max(bx(x) + by(y), -720.0);
    return b;
end

# Let's take a look at this bed. First we define the grid sizes and build some arrays, so we can plot. We'll use a high resolution to get a nice plot:
dx = 1.e3;
dy = 1.e3;
nx = round(Int, 640*1e3/dx);
ny = round(Int, 80*1e3/dx);
xx=[i*dx for i=1:nx, j=1:ny];
yy=[j*dy for i=1:nx, j=1:ny] .- 42000;
x = xx[:,1];
y = yy[1,:];

# Now we can plot
plt =  Plots.heatmap(x/1e3, y/1e3, mismip_plus_bed.(xx,yy)',
                    xlabel = "x (km)",
                    ylabel = "y (km)",
                    colorbar_title = "bed depth (m)")
plot!(size = (800,400))

# ## Boundary Conditions
#
# The MISMIP+ experiment specifies no slip (zero velocity in both directions) boundary conditions at $x = 0$, and free-slip boundary conditions (zero velocity in the direction normal to the walls) at the lateral boundaries at $y = 0$km and $y = 84$km.
# First, let's redefine the grid size to be lower resolution (to make the later simulations quicker)
dx = 4.e3;
dy = 4.e3;
nx = round(Int, 640*1e3/dx);
ny = round(Int, 80*1e3/dx);

# Velocity boundary conditions are controlled by specifying zeros in appropriate entries in arrays, which are then passed to the grid:
u_iszero = falses(nx+1,ny); #build x-direction velocity boundary condition matrix with no zero boundary conditions anywhere
u_iszero[1,:].=true;        #set the x-direction velocity to zero at x = 0.
v_iszero=falses(nx,ny+1);   #build x-direction velocity boundary condition matrix with no zero boundary conditions anywhere
v_iszero[:,1].=true;        #set the y-direction velocity to zero at y = 0 (free slip)
v_iszero[:,end].=true;       #set the y-direction velocity to zero at y = 84km (free slip)
v_iszero[1,:].=true;         #set the y-direction velocity to zero at x = 0km (no slip in combination with u_iszero)

# Now we build the grid as usual, passing the arrays we just constructed via optional arguments.
grid = Grid(nx = nx,
            ny = ny,
            dx = dx,
            dy = dy,
            u_iszero = u_iszero,
            v_iszero = v_iszero);

# ## Solver Parameters
# We're interested in the steady state reached at late times, rather than the solution along the way. We don't need to get the velocity right along the way, just have it correct eventually.
# We therefore set the number of iterations in the velocity solve to be small: at each timestep, the solver just does a small number of iterations, and the velocity is only approximate. But, since we do a lot of iterations getting to steady state, the velocity gets to the right thing eventually.
# This number, and other parameters relating to how the equations are solved, are set via a `SolverParams` object:
solver_params = SolverParams(maxiter_picard = 1);
# Explicitly, we set the number of iterations (formally, Picard iterations) to be as small as possible, i.e. one iteration.

# ## Make the model
# Now we have our grid, bed, and solver parameters, we just need to set the appropriate initial conditions and physical parameters for MISMIP+, and then we can build our model.
# In MISMIP+, the initial condition specifies 100m thick ice everywhere
initial_conditions = InitialConditions(initial_thickness = 100 .* ones(nx,ny));

# And the accumulation rate is set to 0.3m/a (all other defaults in WAVI are chosen according to the values in MISMIP+)
params = Params( accumulation_rate = 0.3);

# Now let's make our model! Note that we use the functional form of the bed (the array we plotted earlier has higher resolution than our model has)
model = Model(grid = grid,
            bed_elevation = mismip_plus_bed,
            initial_conditions = initial_conditions,
            params = params,
            solver_params = solver_params);

# ## Assembling the Simulation
# To get to steady state, we need to run our simulation for a good chunk of time, on the order of tens of 1000s of years. We'll run for 10000 years, with a timestep of half a year:
timestepping_params = TimesteppingParams(dt = 0.5,
                                        end_time = 100.,);
# NB!! We have to specify the end time as `10000.` (a float number) rather than `10000` (an integer) because WAVI.jl expects the same type for the timestep `dt` and the end time `end_time`.

# We'll output the solution along the way, and use this to convince ourselves later than we are in steady state. First let's make a directory to store the output
folder = joinpath(@__DIR__, "overdeepened_bed");
isdir(folder) && rm(folder, force = true, recursive = true);
mkdir(folder) ;

# Then define our output parameters. We'll output the thickness and grounded fraction every 200 years, and set the zip_format keyword argument to zip the output files.
output_params = OutputParams(outputs = (h = model.fields.gh.h,grfrac = model.fields.gh.grounded_fraction),
                            output_freq = 10.,
                            output_path = folder,
                            zip_format = "nc");

# Now we assemble our simulation, taking in the model, output parameters and timestepping_params:
simulation = Simulation(model = model, timestepping_params = timestepping_params, output_params = output_params);

# ## Timestepping
# Now all that's left to do is run our simulation! This is a long simulation and might take a while (~10 mins on my laptop)
run_simulation!(simulation);

# ## Visualization
# Let's have a look at the steady state thickness:
Plots.heatmap(simulation.model.grid.xxh[:,1]/1e3, simulation.model.grid.yyh[1,:]/1e3, simulation.model.fields.gh.h',
                xlabel = "x (km)",
                ylabel = "y (km)",
                colorbar_title = "ice thickness (m)")

# And add the grounding line, which is where the grounded fraction transitions between 0 and 1 (grounded_fraction takes the value 1 at fully grounded grid points and 0 at fully floating grid points.) Our choice of 0.5 is somewhat arbitrary here -- any value between 0 and 1 will do!
Plots.contour!(simulation.model.grid.xxh[:,1]/1e3,
            simulation.model.grid.yyh[1,:]/1e3,
            simulation.model.fields.gh.grounded_fraction',
            fill = false,
            levels = [0.5,0.5],
            linecolor = :blue,
            linewidth = 2)
plot!(size = (1000,550))

# You can see, by comparing with the plot of the bed earlier, that the grounding line sits on an overdeepened section of the bed!

# Finally, let's check that it's in steady state, by looking at the evolution of the volume above floatation:
filename = joinpath(@__DIR__, "overdeepened_bed", "outfile.nc");
h = ncread(filename, "h");
grfrac = ncread(filename, "grfrac");
time = ncread(filename, "TIME");
#compute the volume above floatation
vaf = zeros(1,length(time))
for i = 1:length(time)
    vaf[i] = volume_above_floatation(h[:,:,i], simulation.model.fields.gh.b, Ref(simulation.model.params), simulation.model.grid )
end
Plots.plot(time, vaf[:]/1e9,
             marker = true,
             label = :none,
             xlabel = "time",
             ylabel = "volume above floatation (km^3)",
             framestyle = :box)

# The volume above floatation reaches a plateau, suggesting that we have reached a steady state.
# Finally, we clear up the files we just outputted
rm(folder, force = true, recursive = true);


text:
using SimpleGraphs, ImplicitGraphs

import SimpleGraphs: SimpleGraph, SimpleDigraph

"""
`SimpleGraph(G::ImplicitGraph,A)` returns the induced `SimpleGraph`
of `G` with vertices in the collection `A`.
"""
function SimpleGraph(G::ImplicitGraph{T}, A) where {T}
    H = SimpleGraph{T}()
    for v in A
        if has(G, v)
            add!(H, v)
        end
    end
    for v in A
        for w in A
            if v != w && G[v, w]
                add!(H, v, w)
            end
        end
    end
    return H
end

"""
`SimpleDigraph(G::ImplicitGraph,A)` returns the induced `SimpleDigraph`
of `G` with vertices in the collection `A`.
"""
function SimpleDigraph(G::ImplicitGraph{T}, A) where {T}
    D = SimpleDigraph{T}()
    for v in A
        if has(G, v)
            add!(D, v)
        end
    end

    for v in A
        for w in A
            if G[v, w]
                add!(D, v, w)
            end
        end
    end
    return D
end


text:
export
    Simulation,
    simulate

"""
    Simulation

Describe a PDSampling: information about the initial point, the time of the
simulation, the function to sample from an IPP, etc.
"""
struct Simulation
    x0::Vector{Real}       # Starting point
    v0::Vector{Real}       # Starting velocity
    T::Real                # Simulation time
    nextevent::Function    # Appropriate simulation for first arrival time
    gll::Function          # Gradient of Log Lik (potentially CV)
    nextboundary::Function # Where/When is the next boundary hit
    lambdaref::Real        # Refreshment rate
    algname::String        # BPS, ZZ, GBPS
    # derived
    dim::Int               # dimensionality
    # optional named arguments
    mass::Matrix{Real}     # mass matrix (preconditioner)
    blocksize::Int         # increment the storage by blocks
    maxsimtime::Real       # max. simulation time (s)
    maxsegments::Int       # max. num. of segments
    maxgradeval::Int       # max. num. grad. evals
    refresh!::Function     # refreshment function (TODO: examples)

    # constructor
    function Simulation(x0, v0, T, nextevent, gradloglik, nextboundary,
                lambdaref=1.0, algname="BPS";
                mass=diagm(0=>ones(0)), blocksize=1000, maxsimtime=4e3,
                maxsegments=1_000_000, maxgradeval=100_000_000,
                refresh! = refresh_global! )
        # check none of the default named arguments went through
        @assert !(x0 == :undefined || v0 == :undefined || T == :undefined ||
                  nextevent == :undefined || gradloglik == :undefined || nextboundary == :undefined) "Essential arguments undefined"
        # basic check to see that things are reasonable
        @assert length(x0) == length(v0) > 0 "Inconsistent arguments"
        @assert T > 0.0 "Simulation time must be positive"
        @assert lambdaref >= 0.0 "Refreshment rate must be >= 0"
        ALGNAME = uppercase(algname)
        @assert (ALGNAME ∈ ["BPS", "ZZ","GBPS"]) "Unknown algorithm <$algname>"
        new( x0, v0, T,
             nextevent, gradloglik, nextboundary,
             lambdaref, ALGNAME, length(x0),
             mass, blocksize, maxsimtime,
             maxsegments, maxgradeval, refresh! )
    end
end
# Constructor with named arguments
function Simulation(;
            x0 = :undefined,
            v0 = :undefined,
            T = :undefined,
            nextevent = :undefined,
            gradloglik = :undefined,
            nextboundary = :undefined,
            lambdaref = 1.0,
            algname = "BPS",
            mass = diagm(0=>ones(0)),
            blocksize = 1000,
            maxsimtime = 4e3,
            maxsegments = Int(1e6),
            maxgradeval = Int(1e8),
            refresh! = refresh_global! )
    # calling the unnamed constructor
    Simulation( x0, v0, T,
                nextevent,
                gradloglik,
                nextboundary,
                lambdaref,
                algname;
                mass = mass,
                blocksize = blocksize,
                maxsimtime = maxsimtime,
                maxsegments = maxsegments,
                maxgradeval = maxgradeval,
                refresh! = refresh! )
end

"""
    simulate(sim)

Launch a PD simulation defined in the sim variable. Return the corresponding
Path and a dictionary of indicators (clocktime, ...).
"""
function simulate(sim::Simulation)
    # keep track of how much time we've been going for
    start = time()
    # dimensionality
    d = sim.dim
    # initial states
    x, v = copy(sim.x0), copy(sim.v0)
    # time counter, and segment counter
    t, i = 0.0, 1
    # counters for the number of effective loops and
    # for the number of evaluations of the gradient
    # this will be higher than the number of segments i
    lcnt, gradeval = 0, 0

    # storing by blocks of blocksize nodes at the time
    blocksize = sim.blocksize

    # storing xs as a single column for efficient resizing
    xs, ts = zeros(d*blocksize), zeros(blocksize)
    # store initial point
    xs[1:d] = x
    # mass matrix?
    mass = copy(sim.mass) # store it here as we may want to adapt it

    # check if nextevent takes 2 or 3 parameters
    nevtakes2 = (length(methods(sim.nextevent).ms[1].sig.parameters)-1) == 2

    # compute current reference bounce time
    lambdaref = sim.lambdaref
    tauref = (lambdaref>0.0) ? randexp() / lambdaref : Inf
    # Compute time to next boundary + normal
    (taubd, normalbd) = sim.nextboundary(x, v)

    # keep track of how many refresh events
    nrefresh = 0
    # keep track of how many boundary hits
    nboundary = 0
    # keep track of how many standard bounces
    nbounce = 0

    while (t < sim.T) && (gradeval < sim.maxgradeval)
        # increment the counter to keep track of the number of effective loops
        lcnt += 1
        # simulate first arrival from IPP
        bounce = nevtakes2 ? sim.nextevent(x, v) : sim.nextevent(x, v, tauref)
        # find next event (unconstrained case taubd=NaN (ignored))
        tau = min(bounce.tau, taubd, tauref)
        # standard bounce
        if tau == bounce.tau
            # there will be an evaluation of the gradient
            gradeval += 1
            # updating position/time
            t += tau
            x += tau*v
            # exploiting the memoryless property
            tauref -= tau
            # ---- BOUNCE ----
            g = sim.gll(x)
            if bounce.dobounce(g, v) # e.g.: thinning, acc/rej
                # if accept
                nbounce += 1
                if sim.algname == "BPS"
                    # if a mass matrix is provided
                    if length(mass) > 0
                        v = reflect_bps!(g, v, mass)
                    # standard BPS bounce
                    else
                        v = reflect_bps!(g,v)
                    end
                elseif sim.algname == "GBPS"
                    v = reflect_gbps(g, v)
                elseif sim.algname == "ZZ"
                    v = reflect_zz!(bounce.flipindex, v)
                end
            else
                # move closer to the boundary/refreshment time
                taubd -= tau
                # we don't need to record when rejecting
                continue
            end
        # hard bounce against boundary
        elseif tau == taubd
            #
            nboundary += 1
            # Record point epsilon from boundary for numerical stability
            x += (tau - 1e-10) * v
            t += tau
            # exploiting the memoryless property
            tauref -= tau
            # ---- BOUNCE (boundary) ----
            if sim.algname ∈ ["BPS", "GBPS"]
                # Specular reflection (possibly with mass matrix)
                if length(mass) > 0
                    v = reflect_bps!(normalbd, v, mass)
                else
                    v = reflect_bps!(normalbd,v)
                end
            elseif sim.algname == "ZZ"
                v = reflect_zz!(find((v.*normalbd).<0.0), v)
            end
        # random refreshment
        else
            #= to be in this part, lambdaref should be greater than 0.0
            because if lambdaref=0 then tauref=Inf. There may be a weird
            corner case in which an infinity filters through which we would
            skip =#
            if !isinf(tau)
                #
                nrefresh += 1
                #
                x += tau*v
                t += tau
                # ---- REFRESH ----
                v = sim.refresh!(v)
                # update tauref
                tauref = randexp()/lambdaref
            end
        end
        # check when the next boundary hit will occur
        (taubd, normalbd) = sim.nextboundary(x, v)
        # increment the counter for the number of segments
        i += 1

        # Increase storage on a per-need basis.
        if mod(i,blocksize)==0
            resize!(xs, length(xs) + d * blocksize)
            resize!(ts, length(ts) + blocksize)
        end

        # Storing path times
        ts[i] = t
         # storing columns vertically, a vector is simpler/cheaper to resize
        xs[((i-1) * d + 1):(i * d)] = x

        # Safety checks to break long loops every 100 iterations
        if mod(lcnt, 100) == 0
            if  (time() - start) > sim.maxsimtime
                println("Max simulation time reached. Stopping")
                break
            end
            if  i > sim.maxsegments
                println("Too many segments generated. Stopping")
                break
            end
        end
    end # End of while loop

    details = Dict(
        "clocktime" => time()-start,
        "ngradeval" => gradeval,
        "nloops"    => lcnt,
        "nsegments" => i,
        "nbounce"   => nbounce,
        "nboundary" => nboundary,
        "nrefresh"  => nrefresh
    )

    (Path(reshape(xs[1:(i*d)], (d,i)), ts[1:i]), details)
end


text:
# TODO:
# Install mclust if necessary

function load_or_install_mclust()
  R"""
  if ("mclust" %in% installed.packages()) {
    library(mclust)
  } else {
    cat("Package `mclust` not found. Installing `mclust`...\n")
    install.packages("mclust")
  }
  """
end

function rangev(a, b; length::Int)
  return collect(range(a, b, length=length))
end

function gen_idx(N)
  I = length(N)
  upper_idx = cumsum(N)
  lower_idx = [1; upper_idx[1:end-1] .+ 1]
  return [lower_idx[i]:upper_idx[i] for i in 1:I]
end

function subsampleData(y::Matrix{T}, percentage) where {T <: Number}
  if 0 < percentage < 1
    N = size(y, 1)
    N_subsample = round(Int, N * percentage)
    idx_subsample = StatsBase.sample(1:N, N_subsample, replace=true)
    return y[idx_subsample, :]
  else
    println("Info: percentage ∉  (0, 1). Not subsampling.")
    return y
  end
end

function preimpute!(y::Matrix{T}, missMean::AbstractFloat,
                    missSD::AbstractFloat=0.2) where {T <: Number}
  num_missing = sum(isnan.(y))
  # y[isnan.(y)] .= missMean .+ randn(num_missing) * missSD
  y[isnan.(y)] .= rand(TruncatedNormal(missMean, missSD, -Inf, 0))
end

# https://rdrr.io/cran/mclust/man/mclustModelNames.html
function smartInit(c::Constants, d::Data; iterMax::Int=10,
                   modelNames::String="VVI", warn::Bool=true, seed=0,
                   allowZDupCols=true)
  if seed > 0
    Random.seed!(seed)
  end

  R"set.seed($(seed))"
  if modelNames != "kmeans"
    load_or_install_mclust()
    Mclust = R"mclust::Mclust"
    println("Initialize with Mclust.")
  else
    kmeans = R"kmeans"
    println("Initialize with kmeans.")
  end

  y_imputed = deepcopy(d.y)
  I = d.I
  N = d.N
  J = d.J
  K = c.K
  L = c.L
  idx = gen_idx(N)

  # Mean of imputed values
  missMean = [c.beta[2, i] / (-2 * c.beta[1, i]) for i in 1:I]

  # Preimpute values
  for i in 1:I
    preimpute!(y_imputed[i], missMean[i])
  end

  R"set.seed($(seed))"
  if modelNames == "kmeans"
    clus = kmeans(vcat(y_imputed...), centers=K, iter=iterMax)
  else
    clus = Mclust(vcat(y_imputed...), G=K, modelNames=modelNames, warn=warn)
  end

  # Get order of class labels
  if modelNames == "kmeans"
    lam = [Int.(clus[:cluster])[idx[i]] for i in 1:I]
  else
    # println(clus[:classification])
    lam = [Int.(clus[:classification])[idx[i]] for i in 1:I]
  end
  lam = [Int8.(lam[i]) for i in 1:I]
  ord = [sortperm(lam[i]) for i in 1:I]

  # Get Z
  group_sizes = [sum(vcat(lam...) .== k) for k in 1:K]
  clus_sums = [sum(y_imputed[i][lam[i] .== k, :], dims=1)
               for k in 1:K, i in 1:I]
  clus_means = [sum(clus_sums[k, :]) / group_sizes[k] for k in 1:K]
  Z = Matrix{Bool}(vcat(clus_means...)' .> 0)

  while length(unique(Z, dims=2)) < K && !allowZDupCols
    println("Z columns are not unique! Flipping a bit in the Z matrix.")
    # TODO: Can i do something more sophisticated?
    Z .= filp_bits(Z, 1)
  end

  # Get W
  W = [(sum(lam[i] .== k) + 1/K) / N[i] for i in 1:I, k in 1:K]

  # Get alpha
  alpha = (sum(Z) * K) / (J*K - sum(Z))

  # Get v
  # v = vec(mean(Z .+ (1.0 / J), dims=1))
  v = fill(alpha / (alpha + K), K)

  (m0, s0, _, _) = params(c.delta_prior[0])
  (m1, s1, _, _) = params(c.delta_prior[1])
  delta = Dict(false => rand(TruncatedNormal(m0, s0, 0, Inf), L[0]),
               true  => rand(TruncatedNormal(m1, s1, 0, Inf), L[1]))

  # Get gam
  gam = [zeros(Int8, N[i], J) for i in 1:I]
  for i in 1:I
    for j in 1:J
      for n in 1:N[i]
        z = Z[j, lam[i][n]]
        mus_z = cumsum(delta[z]) * (-1) ^ (1 - z)
        gam[i][n, j] = argmin(abs.(y_imputed[i][n, j] .- mus_z))
      end
    end
  end

  # Get eta
  eta = Dict(Bool(z) => zeros(I, J, L[z]) for z in 0:1)

  for i in 1:I
    for j in 1:J
      counts = Dict(z => ones(L[z]) for z in 0:1)
      for n in 1:N[i]
        z = Z[j, lam[i][n]]
        l = gam[i][n, j]
        counts[z][l] += 1
      end
      eta[0][i, j, :] .= counts[0] / sum(counts[0])
      eta[1][i, j, :] .= counts[1] / sum(counts[1])
    end
  end

  # Get sig2
  sig2 = zeros(I)
  for i in 1:I
    for j in 1:J
      for n in 1:N[i]
        z = Z[j, lam[i][n]]
        l = gam[i][n, j]
        mus_z = cumsum(delta[z]) * (-1) ^ (1 - z)
        sig2[i] += (y_imputed[i][n, j] - mus_z[l]) ^ 2
      end
    end
    sig2[i] /= (N[i] * J)
  end

  eps = mean.(c.eps_prior)

  # Create State object
  state = State{Float64}()

  # Initialize State
  state.Z = Z
  state.delta = delta
  state.alpha = alpha
  state.v = v
  state.W = W
  state.sig2 = sig2
  state.eta = eta
  state.lam = lam
  state.gam = gam
  state.y_imputed = y_imputed
  state.eps = eps

  return state
end


text:
# based on https://dlmf.nist.gov/11.6.i
# not valid for real z < 0

const order_large = 10

_K_large(ν, z) = (1 / pi) * sum(gamma(k + 0.5) * (0.5z)^(ν - 2k - 1) / gamma(ν + 0.5 - k) for k in 0:order_large)

_K0_large(z) = _K_large(0, z)

_M_large(ν, z) = (1 / pi) * sum((-1)^(k + 1) * gamma(k + 0.5) * (0.5z)^(ν - 2k - 1) / gamma(ν + 0.5 - k) for k in 0:order_large)

_M0_large(z) = _M_large(0, z)


text:
# Write function to automate plotting fitted versus residuals plot
# - Plots standard fitted versus residuals plot for assessing linear model fits

function plot_fit_res(model)

    # Extract fitted/predicted values from model object
    pred = GLM.predict(model)

    # Extract residuals from model object
    res = GLM.residuals(model)

    # Create plot
    p = Gadfly.plot(
                # Add points layer
                Gadfly.layer(x = pred,
                             y = res,
                             Gadfly.Geom.point),
                # Add horizontal y = 0 line
                Gadfly.layer(yintercept = [0],
                             Gadfly.Geom.hline),
                # Change plot aesthetics
                Gadfly.Guide.xlabel("Fitted values", orientation=:horizontal),
                Gadfly.Guide.ylabel("Residuals"),
                Gadfly.Guide.title("Residuals vs Fitted"))
end

max_stars_repo_name: guysutton/GLMdiagnostics

max_issues_repo_name: guysutton/GLMdiagnostics

max_forks_repo_name: guysutton/GLMdiagnostics


text:
names = ["Sprinkles", "PeanutButter", "Frosting", "Sugar"]
capacity = [5, -1, 0, -1]
durability = [-1, 3, -1, 0]
flavor = [0, 0, 4, 0]
texture = [0, 0, 0, 2]
calories = [5, 1, 6, 8]

maxscore = 0
for parts in partitions(100, 4)
	for perms in permutations(parts)
		if dot(perms, calories) != 500
			continue
		end
		score = max(0,dot(perms, capacity))*max(0,dot(perms, durability))*max(0,dot(perms, flavor))*max(0,dot(perms, texture))
		if score > maxscore
			maxscore = score
		end
	end
end
println(maxscore)

max_stars_repo_name: jstrube/adventofcode

max_issues_repo_name: jstrube/adventofcode

max_forks_repo_name: jstrube/adventofcode


text:
using GLAbstraction, Colors, GeometryTypes, GLVisualize, Reactive

if !isdefined(:runtests) # We use the examples to generate docs and tests
    window = glscreen()
end

description = """
This example draws multiple equally sized lines efficiently.
"""

vx = -10:0.1:10;
vy = -8:0.1:8;
x = [i for i in vx, j in vy];
y = [j for i in vx, j in vy];
xyz = map(x, y) do x, y
    z = (x * y * (x^2 - y^2) ./ (x^2 + y^2 + eps())) / 5.0
    Point3f0(x, y, z)
end

# ma
lines3d = visualize(
    xyz, :lines,
    thickness = 1f0
)
_view(lines3d, window, camera=:perspective)

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


text:
"""
    conv_encoder(xsize, zsize, kernelsizes, channels, scalings[; activation, densedims])

Constructs a convolutional encoder.

# Arguments
- `xsize`: size of input - (h,w,c)
- `zsize`: latent space dimension
- `kernelsizes`: kernelsizes for conv. layers (only odd numbers are allowed)
- `channels`: channel numbers
- `scalings`: scalings vector
- `activation`: default relu
- `densedims`: if set, more than one dense layers are used

# Example
```julia-repl
julia> encoder = conv_encoder((64, 48, 1), 2, (3, 5, 5), (2, 4, 8), (2, 2, 2), densedims = (256))
Chain(Conv((3, 3), 1=>2, relu), MaxPool((2, 2), pad = (0, 0, 0, 0), stride = (2, 2)), Conv((5, 5), 2=>4, relu), MaxPool((2, 2), pad = (0, 0, 0, 0), stride = (2, 2)), Conv((5, 5), 4=>8, relu), MaxPool((2, 2), pad = (0, 0, 0, 0), stride = (2, 2)), #44, Dense(384, 2))

julia> encoder(randn(Float32, 64, 48, 1, 2))
2×2 Array{Float32,2}:
  0.247844   0.133781
 -0.605763  -0.494911
```
"""
function conv_encoder(xsize::Union{Tuple, Vector}, zsize::Int, kernelsizes::Union{Tuple, Vector},
    channels::Union{Tuple, Vector}, scalings::Union{Tuple, Vector};
    activation = relu, densedims::Union{Tuple, Vector} = [])
    nconv = length(kernelsizes)
    (nconv == length(channels) == length(scalings)) ? nothing : error("incompatible input dimensions")
    (length(xsize) == 3) ? nothing : error("xsize must be (h, w, c)")
    # also check that kernelsizes are all odd numbers
    (all(kernelsizes .% 2 .== 1)) ? nothing : error("not implemented for even kernelsizes")

    # initialize some stuff
    cins = vcat(xsize[3], channels[1:end-1]...) # channels in
    couts = channels # channels out
    ho = xsize[1]/(reduce(*, scalings)) # height before reshaping
    wo = xsize[2]/(reduce(*, scalings)) # width before reshaping
    (ho == floor(Int, ho)) ? ho = floor(Int, ho) : error("your input size and scaling is not compatible")
    (wo == floor(Int, wo)) ? wo = floor(Int, wo) : error("your input size and scaling is not compatible")
    din = ho*wo*channels[end]

    # now build a vector of layers to be used later
    layers = Array{Any,1}()
    # first add the convolutional and maxpooling layers
    for (k, ci, co, s) in zip(kernelsizes, cins, couts, scalings)
        pad = Int((k-1)/2)
        # paddding so that input and output size are same
        push!(layers, Conv((k,k), ci=>co, activation; pad = (pad, pad)))
        push!(layers, MaxPool((s,s)))
    end

    # reshape
    push!(layers, x -> reshape(x, din, :))

    # and dense layers
    ndense = length(densedims)
    dins = vcat(din, densedims...)
    douts = vcat(densedims..., zsize)
    dacts = vcat([activation for _ in 1:ndense]..., identity)
    for (_di, _do, _da) in zip(dins, douts, dacts)
        push!(layers, Dense(_di, _do, _da))
    end

    Flux.Chain(layers...)
end

"""
    conv_decoder(xsize, zsize, kernelsizes, channels, scalings[; activation, densedims, vec_output,
        vec_output_dim])

Constructs a convolutional encoder.

# Arguments
- `xsize`: size of input - (h,w,c)
- `zsize`: latent space dimension
- `kernelsizes`: kernelsizes for conv. layers (only odd numbers are allowed)
- `channels`: channel numbers
- `scalings`: scalings vector
- `activation`: default relu
- `densedims`: if set, more than one dense layers are used
- `vec_output`: output is vectorized (default false)
- `vec_output_dim`: determine what the final size of the vectorized input should be (e.g. add extra dimensions for variance estimation)

# Example
```julia-repl
julia> decoder = conv_decoder((64, 48, 1), 2, (5, 5, 3), (8, 4, 2), (2, 2, 2); densedims = (256))
Chain(Dense(2, 256, relu), Dense(256, 384, relu), #19, ConvTranspose((2, 2), 8=>8, relu), Conv((5, 5), 8=>4, relu), ConvTranspose((2, 2), 4=>4, relu), Conv((5, 5), 4=>2, relu), ConvTranspose((2, 2), 2=>2, relu), Conv((3, 3), 2=>1))

julia> y = decoder(randn(Float32, 2, 2));

julia> size(y)
(64, 48, 1, 2)
```
"""
function conv_decoder(xsize::Union{Tuple, Vector}, zsize::Int, kernelsizes::Union{Tuple, Vector},
    channels::Union{Tuple, Vector}, scalings::Union{Tuple, Vector};
    activation = relu, densedims::Union{Tuple, Vector} = [],
    vec_output = false, vec_output_dim = nothing)
    nconv = length(kernelsizes)
    (nconv == length(channels) == length(scalings)) ? nothing : error("incompatible input dimensions")
    (length(xsize) == 3) ? nothing : error("xsize must be (h, w, c)")
    # also check that kernelsizes are all odd numbers
    (all(kernelsizes .% 2 .== 1)) ? nothing : error("not implemented for even kernelsizes")

    # initialize some stuff
    cins = channels # channels in
    couts = vcat(channels[2:end]..., xsize[3]) # channels out
    ho = xsize[1]/(reduce(*, scalings)) # height after reshaping
    wo = xsize[2]/(reduce(*, scalings)) # width after reshaping
    (ho == floor(Int, ho)) ? ho = floor(Int, ho) : error("your input size and scaling is not compatible")
    (wo == floor(Int, wo)) ? wo = floor(Int, wo) : error("your input size and scaling is not compatible")
    dout = ho*wo*channels[1]

    # now build a vector of layers to be used later
    layers = Array{Any,1}()

    # start with dense layers
    ndense = length(densedims)
    dins = vcat(zsize, densedims...)
    douts = vcat(densedims..., dout)
    for (_di, _do) in zip(dins, douts)
        push!(layers, Dense(_di, _do, activation))
    end

    # reshape
    push!(layers, x -> reshape(x, ho, wo, channels[1], :))

    # add the transpose and convolutional layers
    acts = vcat([activation for _ in 1:nconv-1]..., identity)
    for (k, ci, co, s, act) in zip(kernelsizes, cins, couts, scalings, acts)
        pad = Int((k-1)/2)
        # use convtranspose for upscaling - there are other posibilities, however this seems to be ok
        push!(layers, ConvTranspose((s,s), ci=>ci, activation, stride = (s,s)))
        push!(layers, Conv((k,k), ci=>co, act; pad = (pad, pad)))
    end

    # if you want vectorized output
    if vec_output
        push!(layers, x -> reshape(x, reduce(*, xsize), :))
        # if you want the final vector to be a different size than the output
        # (e.g. you want extra dimensions for variance estimation)
        (vec_output_dim == nothing) ? nothing : push!(layers, Dense(reduce(*, xsize), vec_output_dim))
    end

    Flux.Chain(layers...)
end


text:
## Testing an optimization on a small set of files
##

## Main loop to extract in e2e the DBSCAN parameters.
## The list of votable is selected directly from the directory

using DataFrames
using CSV, Glob, Dates
using Statistics
import DataFrames

rootdir =  ENV["GAIA_ROOT"]

push!(LOAD_PATH,"$rootdir/run/src")
using GaiaClustering

## directory
wdir    = "$rootdir/products"
votdir  = "$rootdir/products/votable"
plotdir = "$rootdir/products/test"
ocdir   = "$rootdir/products/octest"
sclist  = "$rootdir/e2e_products/"

## load a liist of votable and update the file if done
## add results
##

function readlist_votable(filelist::String)
    vot = DataFrames.copy(CSV.read(filelist))
    return(vot)
end

function getdata(filevot)
    voname = filevot

    data       = read_votable(voname)
    df         = filter_data(data)
    dfcart     = add_cartesian(df)
    blck       = [[1,2,3],[4,5], [6,7,8]]
    wghtblck   = [4.0,5.0,1.0]
    norm       = "identity"

    dfcartnorm , scale8 = normalization_PerBlock(dfcart, blck, wghtblck , norm, false)
    return(df, dfcart , dfcartnorm)
end

function mcmc_params()
    minQ    = 2.7
    minstars = 40
    forcedminstars = 30
##
    epsmean   = 2.5
    epsdisp   = 1.5
    min_nei   = 10
    min_cl    = 15
    ncoredisp = 10
    w3dmean   = 6.0
    w3ddisp   = 4.0
    wvelmean  = 6.0
    wveldisp  = 4.0
    whrdmean  = 2.0
    whrddisp  = 1.5
## MCMC parameters
    nburnout  = 2000
    niter     = 15000
##
    pinit = GaiaClustering.abcfull(minQ, minstars, forcedminstars, epsmean, epsdisp, min_nei, min_cl, ncoredisp, w3dmean, w3ddisp ,
    wvelmean, wveldisp, whrdmean, whrddisp, nburnout , niter)
    return(pinit)
end

function dbscanmcmcfull_updt!(ismcmc, fileres, mc ,votname)
    epsm = mean(mc.eps)
    epsd = std(mc.eps)
    mneim = mean(mc.mne)
    mneid = std(mc.mne)
    mclm = mean(mc.mcl)
    mcld = std(mc.mcl)
    qcm = mean(mc.qc)
    qnm = mean(mc.qn)
    qcd = std(mc.qc)
    qnd = std(mc.qn)
    w3dm = mean(mc.w3d)
    w3dd = std(mc.w3d)
    wvelm = mean(mc.wvel)
    wveld = std(mc.wvel)
    whrdm = mean(mc.whrd)
    whrdd = std(mc.whrd)

    println("## DBSCAN/MCMC stats:")
    println("### ϵ : ",epsm," +/- ", epsd)
    println("### min_nei  : ", mneim," +/- ", mneid)
    println("### min_clus  : ", mclm,"+/- ", mcld)
    println("### W3d  : ", w3dm,"+/- ", w3dd)
    println("### Wvel  : ", wvelm,"+/- ", wveld)
    println("### Whrd  : ", whrdm,"+/- ", whrdd)
    println("### Qn  : ",qnm," +/- ", qnd)
    println("### Qc  : ",qcm," +/- ", qcd)
    println("##")

    if !ismcmc
        res = DataFrame(votname=votname, epsm = epsm, epsd=epsd, mneim=mneim,mneid=mneid,mclm=mclm,mcld=mcld,
            qcm=qcm,qcd=qcd, qnm=qnm,qnd=qnd,
            w3dm=w3dm,w3dd=w3dd,wvelm=wvelm,wveld=wveld,whrdm=whrdm,whrdd=whrdd)
        CSV.write(fileres,res,delim=';')
        initmcmc = true
        return(res)
    else
        res = DataFrames.copy(CSV.read(fileres, delim=";"))
        newrow = DataFrame(votname=votname,epsm = epsm, epsd=epsd, mneim=mneim,mneid=mneid,mclm=mclm,mcld=mcld,
            qcm=qcm,qcd=qcd, qnm=qnm,qnd=qnd,
            w3dm=w3dm,w3dd=w3dd,wvelm=wvelm,wveld=wveld,whrdm=whrdm,whrdd=whrdd)
        println("### add DBSCAN/MCMC FULL results ...")

        append!(res,newrow)
        CSV.write(fileres,res,delim=';')
        return(newrow)
    end
end

function check_mcmc(votname, fileres)
    try
        res = DataFrames.copy(CSV.read(fileres, delim=";"))
        if votname in res.votname
            return(true, true)
        else
            return(false , true)
        end
    catch
        println("### $fileres will be created...")
        return(false, false)
    end
end

### update basic parameters of the extracted cluster
####
function SCparameters_updt(fileres,sc::GaiaClustering.SCproperties,votname)
    try
        res = DataFrames.copy(CSV.read(fileres, delim=";"))
        newrow = DataFrame(votname=votname)
        newrow = DataFrame(votname=votname,nstars=sc.nstars,distance=sc.distance,l=sc.l,b=sc.b,
            vl=sc.vl,vb=sc.vb,vrad=sc.vrad , xdisp=sc.xdisp,ydisp=sc.ydisp,zdisp=sc.zdisp,
            vldisp=sc.vldisp,vbdisp=sc.vbdisp , vraddisp=sc.vraddisp)
        println("### add SC basic properties...")
        append!(res,newrow)
        CSV.write(fileres,res,delim=';')
        return(newrow)
    catch
        println("## No $fileres file, it will be created...")
        res = DataFrame(votname=votname,nstars=sc.nstars,distance=sc.distance,l=sc.l,b=sc.b,
            vl=sc.vl,vb=sc.vb,vrad=sc.vrad , xdisp=sc.xdisp,ydisp=sc.ydisp,zdisp=sc.zdisp,
            vldisp=sc.vldisp,vbdisp=sc.vbdisp , vraddisp=sc.vraddisp)
        CSV.write(fileres,res,delim=';')
        initmcmc = true
        return(res)
    end
end

## check for blacklist
##
function read_blacklist(blackname)

    if isfile(blackname)
        df= DataFrames.copy(CSV.read(blackname, delim=";"))
        blacklist= df.votname
    else
        blacklist= [""]
    end

    return(blacklist)
end
#######################
## Main loop
##
function main(filelist,fileres, fileSCres)
    let
        println("## Starting main loop ..")
        println("## It can be very long but will be resumed to the last reduced file.")
        wd= pwd()
        println("## Working directory: $wd")
        nfile= size(filelist)[1]
        println("## $nfile files to analyze...")
        totalTime= 0.

        # read a possible votname blacklist
        blackname= "blacklist-test.csv"
        blacklist= read_blacklist(blackname)
        println("## blacklist read...")

        s= size(filelist)

        for i in 1:nfile
            votname = filelist[i]
            mcmcfound , ismcmcfile = check_mcmc(votname, fileres)

            ## test blacklist
            if votname in blacklist
                println("##")
                println("## $votname in blacklist...")
                println("##")
                mcmcfound= true
            end

            if !mcmcfound
                println("## Starting with $votname")
                tstart= now()
                println("## Starting at $tstart")

                df , dfcart , dfcartnorm = getdata(votdir*"/"*votname)

                ## MCMC optimization
                params = mcmc_params()
                mc = abc_mcmc_dbscan_full2(dfcart, params)
                plot_dbscanfull_mcmc(plotdir, votname, mc , false)
                res = dbscanmcmcfull_updt!(ismcmcfile, fileres, mc ,votname)

                ## get the cluster and plot it
                println("## Extracting the cluster using DBSCAN/WEIGHTING with:")
                eps = res.epsm[1]
                min_nei = trunc(Int,res.mneim[1] + 0.5)
                min_cl = trunc(Int,res.mclm[1] + 0.5)
                w3d = res.w3dm[1]
                wvel = res.wvelm[1]
                whrd = res.whrdm[1]
                println("### ϵ : $eps")
                println("### min_neighbor: $min_nei")
                println("### min_cluster : $min_cl")
                println("### W3d : $w3d")
                println("### Wvel : $wvel")
                println("### Whrd : $whrd")

                mres = GaiaClustering.modelfull(eps,min_nei,min_cl,w3d,wvel,whrd)
                dfcartnorm = getDfcartnorm(dfcart, mres)
                labels = clusters(dfcartnorm.data ,eps  , 20, min_nei, min_cl)
                #labelmax , nmax = find_cluster_label(labels)
                labelmax , nmax = find_cluster_label2(labels, df, dfcart)
                println("### Label solution: $labelmax")
                println("### N stars: $nmax")

                ### best solution in terms of Qc...
                println("## Extracting and writing the OC stars...")
                export_df("$votname.0", ocdir, df , dfcart, labels , labelmax)
                scproperties0 = get_properties_SC(labels[labelmax] , df, dfcart)
                println("### ",scproperties0)
                plot_cluster2(plotdir, "$votname.0", labels[labelmax], scproperties0,  dfcart , false)

                ### all solutions are plotted/exported
                index=1
                for ilab in labels
                    scproperties = get_properties_SC(ilab , df, dfcart)
                    export_df("$votname.$index", ocdir, df , dfcart, labels , index)
                    plot_cluster2(plotdir, "$votname.$index", ilab, scproperties,  dfcart , false)
                    index +=1
                end
                ######

                SCparameters_updt(fileSCres, scproperties0, votname)

                tend= now()
                println("## Ending at $tend")
                duration= Dates.value(tend-tstart) / (1000*3600)
                totalTime += duration
                meanTime= totalTime / i
                ETA= meanTime * (nfile-i) / 24
                println("## Duration: $duration hours")
                println("## ETA: $ETA days")
                println("##\n##")

            end

        end
    end
    print("## Main loop done.")
end

###############################################################################
println("# Test of the optimization on a subset of targets. Check the code...")

cd(votdir)
votlist= glob("Ber*.vot")
cd(wdir)

main(votlist,"votlist-test_mcmc_full.csv", "votlist-test_SCproperties_full.csv")

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

estimate a covariance matrix that satisfies some given prior information and
minimizes a given convex spectral function

p ∈ 𝕊ᵈ is the covariance variable
minimize    f(p)                    (note: enforces p ⪰ 0)
subject to  tr(p) = 1               (normalize)
            B p = b                 (prior info as equalities)
            C p ≤ c                 (prior info as inequalities)
where f is a convex spectral function
=#

struct CovarianceEstJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    d::Int
    ext::MatSpecExt # formulation specifier
end

function build(inst::CovarianceEstJuMP{T}) where {T <: Float64}
    d = inst.d
    @assert d >= 1
    @assert is_domain_pos(inst.ext)

    p0 = randn(T, d, d)
    p0 = p0 * p0' + I / 2
    p0 ./= tr(p0)
    vec_dim = Cones.svec_length(d)
    p0_vec = zeros(T, vec_dim)
    Cones.smat_to_svec!(p0_vec, p0, sqrt(T(2)))

    model = JuMP.Model()
    JuMP.@variable(model, p[1:d, 1:d], Symmetric)
    JuMP.@constraint(model, tr(p) == 1)
    p_vec = zeros(JuMP.AffExpr, vec_dim)
    Cones.smat_to_svec!(p_vec, one(T) * p, sqrt(T(2)))

    # convex objective
    JuMP.@variable(model, epi)
    JuMP.@objective(model, Min, epi)
    add_homog_spectral(inst.ext, d, vcat(1.0 * epi, p_vec), model)

    # linear prior constraints
    lin_dim = round(Int, sqrt(d - 1))
    B = randn(T, lin_dim, vec_dim)
    b = B * p0_vec
    JuMP.@constraint(model, B * p_vec .== b)
    C = randn(T, lin_dim, vec_dim)
    c = C * p0_vec
    JuMP.@constraint(model, C * p_vec .<= c)

    # save for use in tests
    model.ext[:p_var] = p

    return model
end

function test_extra(inst::CovarianceEstJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    (stat == MOI.OPTIMAL) || return

    # check objective
    tol = eps(T)^0.2
    p_opt = JuMP.value.(model.ext[:p_var])
    λ = eigvals(Symmetric(p_opt, :U))
    @test minimum(λ) >= -tol
    obj_result = get_val(pos_only(λ), inst.ext)
    @test JuMP.objective_value(model) ≈ obj_result atol=tol rtol=tol
    return
end