text:
#=
This is the place, where I put functions, which are so annoying in OpenGL, that I felt the need to wrap them and make them more "Julian"
Its also to do some more complex error handling, not handled by the debug callback
=#
getnames(check_function::Function) = filter(check_function, uint32(0:65534))
# gets all the names currently boundo to programs
getProgramNames() = getnames(glIsProgram)
getShaderNames() = getnames(glIsShader)
getVertexArrayNames() = getnames(glIsVertexArray)

function glGetShaderiv(shaderID::GLuint, variable::GLenum)
    result = Ref{GLint}(-1)
    glGetShaderiv(shaderID, variable, result)
    result[]
end
function glShaderSource(shaderID::GLuint, shadercode::Vector{UInt8})
    shader_code_ptrs = Ptr{UInt8}[pointer(shadercode)]
    len = Ref{GLint}(length(shadercode))
    glShaderSource(shaderID, 1, shader_code_ptrs, len)
end
glShaderSource(shaderID::GLuint, shadercode::String) = glShaderSource(shaderID, Vector{UInt8}(shadercode))
function glGetAttachedShaders(program::GLuint)
    shader_count   = glGetProgramiv(program, GL_ATTACHED_SHADERS)
    length_written = GLsizei[0]
    shaders        = zeros(GLuint, shader_count)

    glGetAttachedShaders(program, shader_count, length_written, shaders)
    shaders[1:first(length_written)]
end

function ModernGL.glGetActiveUniformsiv(program::GLuint, index, var::GLenum)
    result = Ref{GLint}(-1)
    glGetActiveUniformsiv(program, 1, Ref{GLuint}(index), var, result)
    result[]
end

function glGetActiveUniform(programID::GLuint, index::Integer)
    actualLength   = GLsizei[1]
    uniformSize    = GLint[1]
    typ            = GLenum[1]
    maxcharsize    = glGetProgramiv(programID, GL_ACTIVE_UNIFORM_MAX_LENGTH)
    name           = zeros(GLchar, maxcharsize)

    glGetActiveUniform(programID, index, maxcharsize, actualLength, uniformSize, typ, name)

    actualLength[1] <= 0 &&  @error "No active uniform at given index. Index: $index"

    uname = unsafe_string(pointer(name), actualLength[1])
    uname = Symbol(replace(uname, r"\[\d*\]" => "")) # replace array brackets. This is not really a good solution.
    (uname, typ[1], uniformSize[1])
end
glGetActiveUniformName(program::GLuint, index::Integer) = glGetActiveUniform(program, index)[1]

function glGetActiveAttrib(programID::GLuint, index::Integer)
    actualLength   = GLsizei[1]
    attributeSize  = GLint[1]
    typ            = GLenum[1]
    maxcharsize    = glGetProgramiv(programID, GL_ACTIVE_ATTRIBUTE_MAX_LENGTH)
    name           = zeros(GLchar, maxcharsize)

    glGetActiveAttrib(programID, index, maxcharsize, actualLength, attributeSize, typ, name)

    actualLength[1] <= 0 && @error "No active uniform at given index. Index: $index"

    uname = unsafe_string(pointer(name), actualLength[1])
    uname = Symbol(replace(uname, r"\[\d*\]" => "")) # replace array brackets. This is not really a good solution.
    (uname, typ[1], attributeSize[1])
end
function glGetProgramiv(programID::GLuint, variable::GLenum)
    result = Ref{GLint}(-1)
    glGetProgramiv(programID, variable, result)
    result[]
end
function glGetIntegerv(variable::GLenum)
    result = Ref{GLint}(-1)
    glGetIntegerv(UInt32(variable), result)
    result[]
end

function glGenBuffers(n=1)
    result = GLuint[0]
    glGenBuffers(1, result)
    id = result[]
    if id <= 0
        @error "glGenBuffers returned invalid id. OpenGL Context active?"
    end
    id
end
function glGenVertexArrays()
    result = GLuint[0]
    glGenVertexArrays(1, result)
    id = result[1]
    if id <= 0
        @error "glGenVertexArrays returned invalid id. OpenGL Context active?"
    end
    id
end
function glGenTextures()
    result = GLuint[0]
    glGenTextures(1, result)
    id = result[1]
    if id <= 0
        @error "glGenTextures returned invalid id. OpenGL Context active?"
    end
    id
end
function glGenFramebuffers()
    result = GLuint[0]
    glGenFramebuffers(1, result)
    id = result[1]
    if id <= 0
        @error "glGenFramebuffers returned invalid id. OpenGL Context active?"
    end
    id
end

function glDeleteTextures(id::GLuint)
  arr = [id]
  glDeleteTextures(1, arr)
end
function glDeleteVertexArrays(id::GLuint)
  arr = [id]
  glDeleteVertexArrays(1, arr)
end
function glDeleteBuffers(id::GLuint)
  arr = [id]
  glDeleteBuffers(1, arr)
end

function glGetTexLevelParameteriv(target::GLenum, level, name::GLenum)
  result = GLint[0]
  glGetTexLevelParameteriv(target, level, name, result)
  result[1]
end

function glGenRenderbuffers(format::GLenum, attachment::GLenum, dimensions)
    renderbuffer = GLuint[0]
    glGenRenderbuffers(1, renderbuffer)
    glBindRenderbuffer(GL_RENDERBUFFER, renderbuffer[1])
    glRenderbufferStorage(GL_RENDERBUFFER, format, dimensions...)
    glFramebufferRenderbuffer(GL_FRAMEBUFFER, attachment, GL_RENDERBUFFER, renderbuffer[1])
    renderbuffer[1]
end

function glTexImage(ttype::GLenum, level::Integer, internalFormat::GLenum, w::Integer, h::Integer, d::Integer, border::Integer, format::GLenum, datatype::GLenum, data)
    glTexImage3D(GL_PROXY_TEXTURE_3D, level, internalFormat, w, h, d, border, format, datatype, C_NULL)
    for l in  0:level
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_3D, l, GL_TEXTURE_WIDTH)
        if result == 0
            @error "glTexImage 3D: width too large. Width: $w"
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_3D, l,GL_TEXTURE_HEIGHT)
        if result == 0
            @error "glTexImage 3D: height too large. height: $h"
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_3D, l, GL_TEXTURE_DEPTH)
        if result == 0
            @error "glTexImage 3D: depth too large. Depth: $d"
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_3D, l, GL_TEXTURE_INTERNAL_FORMAT)
        if result == 0
            @error "glTexImage 3D: internal format not valid. format: $(GLENUM(internalFormat).name)"
        end
    end
    glTexImage3D(ttype, level, internalFormat, w, h, d, border, format, datatype, data)
end
function glTexImage(ttype::GLenum, level::Integer, internalFormat::GLenum, w::Integer, h::Integer, border::Integer, format::GLenum, datatype::GLenum, data)
    glTexImage2D(GL_PROXY_TEXTURE_2D, level, internalFormat, w, h, border, format, datatype, C_NULL)
    for l in 0:level
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_2D, l, GL_TEXTURE_WIDTH)
        if result == 0
            @error "glTexImage 2D: width too large. Width: $w"
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_2D, l, GL_TEXTURE_HEIGHT)
        if result == 0
            @error "glTexImage 2D: height too large. height: $h"
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_2D, l, GL_TEXTURE_INTERNAL_FORMAT)
        if result == 0
            @error "glTexImage 2D: internal format not valid. format: $(GLENUM(internalFormat).name)"
        end
    end
    glTexImage2D(ttype, level, internalFormat, w, h, border, format, datatype, data)
end
function glTexImage(ttype::GLenum, level::Integer, internalFormat::GLenum, w::Integer, border::Integer, format::GLenum, datatype::GLenum, data)
    glTexImage1D(GL_PROXY_TEXTURE_1D, level, internalFormat, w, border, format, datatype, C_NULL)
    for l in 0:level
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_1D, l, GL_TEXTURE_WIDTH)
        if result == 0
            @error "glTexImage 1D: width too large. Width: $w"
        end
        result = glGetTexLevelParameteriv(GL_PROXY_TEXTURE_1D, l, GL_TEXTURE_INTERNAL_FORMAT)
        if result == 0
            @error "glTexImage 1D: internal format not valid. format: $(GLENUM(internalFormat).name)"
        end
    end
    glTexImage1D(ttype, level, internalFormat, w, border, format, datatype, data)
end

function compile_program(shaders::GLuint...)

    program = glCreateProgram()::GLuint
    glUseProgram(program)
    #attach new ones
    foreach(shaders) do shader
        glAttachShader(program, shader.id)
    end
    #link program
    glLinkProgram(program)
    if !islinked(program)
        for shader in shaders
            write(stdout, shader.source)
            println("---------------------------")
        end
        @error "program $program not linked. Error in: \n $(join(map(x-> string(x.name), shaders))), or, \n $(getinfolog(program))"
    end
    program
end

# function glsl_version_string()
#     glsl = split(unsafe_string(glGetString(GL_SHADING_LANGUAGE_VERSION)), ['.', ' '])
#     if length(glsl) >= 2
#         glsl = VersionNumber(parse(Int, glsl[1]), parse(Int, glsl[2]))
#         glsl.major == 1 && glsl.minor <= 2 && (@error "OpenGL shading Language version too low. Try updating graphic driver!")
#         glsl_version = string(glsl.major) * rpad(string(glsl.minor),2,"0")
#         return "#version $(glsl_version)\n"
#     else
#         @error "could not parse GLSL version: $glsl"
#     end
# end


text:
import MLJBase
import MonteCarloMeasurements
import Statistics

const rms_distribution = RMSDistribution()
const rms_expected = RMSExpected()
const rms_median = RMSMedian()

function (::RMSDistribution)(ŷ::MonteCarloMeasurements.MvParticles, y::Vector{<:Real})
    return MLJBase.rms(ŷ, y)
end

function (::RMSExpected)(ŷ::MonteCarloMeasurements.MvParticles, y::Vector{<:Real})
    return Statistics.mean(rms_distribution(ŷ, y))
end

function (::RMSMedian)(ŷ::MonteCarloMeasurements.MvParticles, y::Vector{<:Real})
    return Statistics.median(rms_distribution(ŷ, y))
end

ext: jl


text:
struct WeightingPotential{T, N, S} <: AbstractArray{T, N}
    data::Array{T, N}
    grid::Grid{T, N, S}
end

@inline size(wp::WeightingPotential{T, N, S}) where {T, N, S} = size(wp.data)
@inline length(wp::WeightingPotential{T, N, S}) where {T, N, S} = length(wp.data)
@inline getindex(wp::WeightingPotential{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(wp.data, I...)
@inline getindex(wp::WeightingPotential{T, N, S}, i::Int) where {T, N, S} = getindex(wp.data, i)
@inline getindex(wp::WeightingPotential{T, N, S}, s::Symbol) where {T, N, S} = getindex(wp.grid, s)

function WeightingPotential(fss::PotentialSimulationSetup{T, 3, :cylindrical}; kwargs...)::WeightingPotential{T, 3, :cylindrical} where {T <: SSDFloat}
    return get_2π_potential(WeightingPotential{T, 3, :cylindrical}(fss.potential, fss.grid); kwargs...)
end
function WeightingPotential(fss::PotentialSimulationSetup{T, 3, :cartesian})::WeightingPotential{T, 3, :cartesian} where {T <: SSDFloat}
    return WeightingPotential{T, 3, :cartesian}(fss.potential, fss.grid)
end

@recipe function f( wp::WeightingPotential{T, 3, :cylindrical};
                    r = missing,
                    φ = missing,
                    z = missing ) where {T <: SSDFloat}
    g::Grid{T, 3, :cylindrical} = wp.grid

    seriescolor --> :viridis
    st --> :heatmap
    aspect_ratio --> 1
    foreground_color_border --> nothing
    tick_direction --> :out

    cross_section::Symbol, idx::Int = if ismissing(φ) && ismissing(r) && ismissing(z)
        :φ, 1
    elseif !ismissing(φ) && ismissing(r) && ismissing(z)
        φ_rad::T = T(deg2rad(φ))
        while !(g.φ.interval.left <= φ_rad <= g.φ.interval.right)
            if φ_rad > g.φ.interval.right
                φ_rad -= g.φ.interval.right - g.φ.interval.left
            elseif φ_rad < g.φ.interval.left
                φ_rad += g.φ.interval.right - g.φ.interval.left
            end
        end
        :φ, searchsortednearest(g.φ, φ_rad)
    elseif ismissing(φ) && !ismissing(r) && ismissing(z)
        :r, searchsortednearest(g.r, T(r))
    elseif ismissing(φ) && ismissing(r) && !ismissing(z)
        :z, searchsortednearest(g.z, T(z))
    else
        error(ArgumentError, ": Only one of the keywords `r, φ, z` is allowed.")
    end
    value::T = if cross_section == :φ
        g.φ[idx]
    elseif cross_section == :r
        g.r[idx]
    elseif cross_section == :z
        g.z[idx]
    end

    @series begin
        clims --> (0, 1)
        if cross_section == :φ
            title --> "Weighting Potential @$(cross_section) = $(round(rad2deg(value), sigdigits = 2))"
            xlabel --> "r / m"
            ylabel --> "z / m"
            size --> ( 400, 350 / (g.r[end] - g.r[1]) * (g.z[end] - g.z[1]) )
            g.r, g.z, wp.data[:, idx,:]'
        elseif cross_section == :r
            g.φ, g.z, wp.data[idx,:,:]'
        elseif cross_section == :z
            title --> "Weighting Potential @$(cross_section) = $(round(value, sigdigits = 2))"
            proj --> :polar
            g.φ, g.r, wp.data[:,:,idx]
        end
    end
end

const ScalarPotential{T, N, S} = Union{ElectricPotential{T, N, S}, WeightingPotential{T, N, S}, PointTypes{T, N, S}, ChargeDensity{T, N, S}}

function get_2π_potential(wp::ScalarPotential{T, 3, :cylindrical}, axφ::DiscreteAxis{T, :periodic, :periodic}, int::Interval{:closed, :open, T}) where {T}
    @assert int.right != 0 "Right boundary of φ interval is not allowed to be 0"
    Potential::Type = typeof(wp)
    l::Int = length( axφ )
    Δφ::T = int.right - int.left
    new_int::Interval{:closed, :open, T} = Interval{:closed, :open, T}(0, 2π)
    n::Int = Int(round(T(2π) / Δφ, sigdigits = 6))
    new_ticks::Vector{T} = Vector{T}(undef, l * n)
    new_pot::Array{T, 3} = Array{T, 3}(undef, size(wp, 1), l * n, size(wp, 3))
    for idx_n in 1:n
        for il in 1:l
            idx::Int = il + (idx_n - 1) * l
            new_ticks[idx] = (idx_n - 1) * Δφ + axφ.ticks[il]
            new_pot[:, idx, :] = wp[:, il, :]
        end
    end
    new_axφ::DiscreteAxis{T, :periodic, :periodic} = DiscreteAxis{T, :periodic, :periodic}( new_int, new_ticks )
    new_grid::Grid{T, 3, :cylindrical} = Grid{T, 3, :cylindrical}( (wp.grid[1], new_axφ, wp.grid[3]) )
    new_wp::Potential = Potential( new_pot, new_grid )
    return new_wp
end

function get_2π_potential(wp::ScalarPotential{T, 3, :cylindrical}, axφ::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval{:closed, :closed, T}) where {T}
    @assert int.right != 0 "Right boundary of φ interval is not allowed to be 0"
    Potential::Type = typeof(wp)
    l::Int = 2 * length( axφ ) - 2
    Δφ::T = int.right - int.left
    new_int::Interval{:closed, :open, T} = Interval{:closed, :open, T}(int.left, int.right + Δφ)
    new_ticks::Vector{T} = Vector{T}(undef, l)
    new_ticks[1:length( axφ )] = collect(axφ.ticks)
    new_pot::Array{T, 3} = Array{T, 3}(undef, size(wp, 1), l, size(wp, 3))
    new_pot[:, 1:length( axφ ), :] = wp.data[:, :, :]
    for i in 1:(length(axφ) - 2)
        new_ticks[length(axφ) + i] = 2 * int.right - axφ.ticks[length(axφ) - i]
        new_pot[:, length(axφ) + i, :] = wp[:, length(axφ) - i, :]
    end
    new_axφ::DiscreteAxis{T, :periodic, :periodic} = DiscreteAxis{T, :periodic, :periodic}( new_int, new_ticks )
    new_grid::Grid{T, 3, :cylindrical} = Grid{T, 3, :cylindrical}( (wp.grid[1], new_axφ, wp.grid[3]) )
    new_wp::Potential = Potential( new_pot, new_grid )
    return get_2π_potential(new_wp, new_axφ, new_int)
    # return new_wp
end

function extend_2D_to_3D_by_n_points(wp::ScalarPotential{T, 3, :cylindrical}, axφ::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval{:closed, :closed, T},
        n_points_in_φ::Int ) where {T}
    Potential::Type = typeof(wp)
    new_int::Interval{:closed, :open, T} = Interval{:closed, :open, T}(0, 2π)
    new_ticks::Vector{T} = Vector{T}(undef, n_points_in_φ)
    new_pot::Array{T, 3} = Array{T, 3}(undef, size(wp, 1), n_points_in_φ, size(wp, 3))
    Δφ::T = 2π / n_points_in_φ
    for i in 1:n_points_in_φ
        new_ticks[i] = (i - 1) * Δφ
        new_pot[:, i, :] = wp[:, 1, :]
    end
    new_axφ::DiscreteAxis{T, :periodic, :periodic} = DiscreteAxis{T, :periodic, :periodic}( new_int, new_ticks )
    new_grid::Grid{T, 3, :cylindrical} = Grid{T, 3, :cylindrical}( (wp.grid[1], new_axφ, wp.grid[3]) )
    new_wp::Potential = Potential( new_pot, new_grid )
    return new_wp
end

function get_2π_potential(wp::ScalarPotential{T, 3, :cylindrical}; n_points_in_φ::Union{Missing, Int} = missing) where {T}
    axφ::DiscreteAxis{T} = wp.grid[2]
    int::Interval = axφ.interval
    if int.right == 0 && length(axφ) == 1 # 2D
        if ismissing(n_points_in_φ)
            error(ArgumentError, ": First Argument is a 2D potential (only 1 point in φ). User has to set the keyword `n_points_in_φ::Int` (e.g. `n_points_in_φ = 18`) in order to get a 3D potential.")
        else
            return extend_2D_to_3D_by_n_points(wp, axφ, int, n_points_in_φ)
        end
    else
        return get_2π_potential(wp, axφ, int)
    end
end

function WeightingPotential(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(ustrip(nt.values[1]))
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    WeightingPotential{T, N, S}( nt.values, grid)
end
Base.convert(T::Type{WeightingPotential}, x::NamedTuple) = T(x)

function NamedTuple(ep::WeightingPotential{T, 3}) where {T}
    return (
        grid = NamedTuple(ep.grid),
        values = ep.data,
    )
end
Base.convert(T::Type{NamedTuple}, x::WeightingPotential) = T(x)

"""
    PointTypes(setup::PotentialSimulationSetup{T, 3, :cylindrical} ; kwargs...)::PointTypes{T, 3, :cylindrical}

Extracts the electric potential from `setup` and extrapolate it to an 2π grid.

For 2D grids (r and z) the user has to set the keyword `n_points_in_φ::Int`, e.g.: `n_points_in_φ = 36`.
"""
function PointTypes(setup::PotentialSimulationSetup{T, 3, :cylindrical} ; kwargs...)::PointTypes{T, 3, :cylindrical} where {T}
    return get_2π_potential(PointTypes{T, 3, :cylindrical}(setup.pointtypes, setup.grid); kwargs...)
end

@recipe function f( wp::WeightingPotential{T, 3, :cartesian};
                    # dim = missing, dimvalue = missing,
                    x = missing,
                    y = missing,
                    z = missing,
                    contours_equal_potential=false ) where {T}
    g::Grid{T, 3, :cartesian} = wp.grid

    seriescolor --> :viridis
    st --> :heatmap
    aspect_ratio --> 1
    foreground_color_border --> nothing
    tick_direction --> :out

    cross_section::Symbol, idx::Int = if ismissing(x) && ismissing(y) && ismissing(z)
        :x, 1
    elseif !ismissing(x) && ismissing(y) && ismissing(z)
        :x, searchsortednearest(g[:x], T(x))
    elseif ismissing(x) && !ismissing(y) && ismissing(z)
        :y, searchsortednearest(g[:y], T(y))
    elseif ismissing(x) && ismissing(y) && !ismissing(z)
        :z, searchsortednearest(g.z, T(z))
    else
        error(ArgumentError, ": Only one of the keywords `r, y, z` is allowed.")
    end
    value::T = if cross_section == :x
        g[:x][idx]
    elseif cross_section == :y
        g[:y][idx]
    elseif cross_section == :z
        g.z[idx]
    end

    @series begin
        title --> "Weighting Potential @$(cross_section) = $(round(value, sigdigits = 2))"
        clims --> (0, 1)
        if cross_section == :x
            xlabel --> "y / m"
            ylabel --> "z / m"
            g[:y], g.z, wp.data[idx, :, :]'
        elseif cross_section == :y
            xlabel --> "x / m"
            ylabel --> "z / m"
            g[:x], g.z, wp.data[:, idx, :]'
        elseif cross_section == :z
            xlabel --> "x / m"
            ylabel --> "y / m"
            g[:x], g[:y], wp.data[:,:,idx]'
        end
    end
    if contours_equal_potential
        @series begin
            seriescolor := :thermal
            st := :contours
            if cross_section == :x
                g[:y], g.z, wp.data[idx, :, :]'
            elseif cross_section == :y
                g[:x], g.z, wp.data[:, idx, :]'
            elseif cross_section == :z
                g[:x], g[:y], wp.data[:,:,idx]
            end
        end
    end
end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

@testset "Indices" begin

    @testset "EENS" begin
        # Example from: Analytical Model for Availability Assessment of Large-
        # Scale Offshore Wind Farms including their Collector System (p. 6)
        ugf = UGF(:power,   [0.0u"MW", 2.0u"MW", 4.0u"MW", 6.0u"MW", 8.0u"MW"],
                            [0.307000, 0.012600, 0.119070, 0.102060, 0.459270])
        pcc = _MSS.PropDict(:ugf => ugf)
        @test isapprox(_MSS.EENS(pcc), 28137.12u"MWh")
    end

    @testset "GRA" begin
        # Example from: Analytical Model for Availability Assessment of Large-
        # Scale Offshore Wind Farms including their Collector System (p. 6)
        ugf = UGF(:power,   [0.0u"MW", 1.0u"MW", 2.0u"MW", 3.0u"MW", 4.0u"MW"],
                            [0.010000, 0.018000, 0.170100, 0.145800, 0.656100])
        ntw = _MSS.PropDict(:ugf => ugf)
        @test isapprox(_MSS.GRO(ntw, 0.7), 2.8u"MW")
        @test isapprox(_MSS.GRA(ntw, 0.7), 0.8019)
    end

end


text:
push!(LOAD_PATH,"../src/")

using Documenter, LocalCongruence
LC = LocalCongruence

Documenter.makedocs(
	modules = [LocalCongruence],
    format = Documenter.HTML(
        assets = ["assets/LocalCongruence.css"],
        highlights = ["yaml"],
	),
	clean = true,
	sitename = "LocalCongruence.jl",
	authors = "Gianmaria Del Monte, Elia Onofri, Giorgio Scorzelli and Alberto Paoluzzi.",
	pages = [
		"1 - Home" => "index.md",
		"2 - Complexes and Chain Operators" => "theory.md",
		"3 - GraphBLAS Introduction" => "graph_blas.md",
		"Cell Congruence Enabling" => [
			"4.1 - Vertices Congruence" => "verticesCongruence.md",
			"4.2 - Chain Complex Congruence" => "chainComplexCongruence.md",
		],
		"Examples" => [
			"5.1 - Cube Grids" => "example_1.md",
			"5.2 - ..." => "example_2.md",
			"5.3 - ..." => "example_3.md",
		],
		"A - About the Authors" => "authors.md",
		"B - Bibliography" => "bibliography.md"
	]
)


text:
function sdl_colors(c::Colorant)
  sdl_colors(
      convert(ARGB{Colors.FixedPointNumbers.Normed{UInt8,8}}, c)
     )
end

sdl_colors(c::ARGB) = Int.(reinterpret.((red(c), green(c), blue(c), alpha(c))))


text:
mutable struct TextBoxEntry
    value::String
end

mutable struct SliderEntry
    value::Float64
end

mutable struct Value
  name::String
  accepted::Bool
end

mutable struct Column
  name::String
  split::Bool
  values::Vector{Value}
  _values::ListModel
  Column(name, split, values) = new(name, split, values, ListModel(values))
end

Column(name, values::AbstractArray) = Column(name, false, Value.(collect(values), true))

mutable struct ComboBoxEntry
    value::String
end

mutable struct ComboBoxType
    name::String
    chosen_value::String
    options::Vector{ComboBoxEntry}
    _options::ListModel
    ask_info::Bool
    text_info::String
    ComboBoxType(name, options, ask_info) = new(name, options[1].value, options, ListModel(options), ask_info, "")
end

ComboBoxType(name, option_names::AbstractArray{T}, ask_info = false) where {T <: AbstractString} =
    ComboBoxType(name, ComboBoxEntry.(option_names), ask_info)

mutable struct SpinBoxEntry{T<:Real}
    selected_value::T
    min_value::T
    max_value::T
end

mutable struct SpinBoxType{T<:Real}
    name::String
    values::Vector{SpinBoxEntry{T}}
    _values::ListModel
    SpinBoxType{T}(name, options::Vector{SpinBoxEntry{T}}) where {T<:Real} =
        new(name, options, ListModel(options))
end

SpinBoxType(name, values::Vector{T}) where T<:Real =
    SpinBoxType{T}(name, SpinBoxEntry.(values, values[1], values[2]))


text:
using SHA

println(bytes2hex(sha256("test")))

max_stars_repo_name: marchunter/litcoin

max_issues_repo_name: marchunter/litcoin

max_forks_repo_name: marchunter/litcoin


text:
include("./Spherical_harmonics.jl")
using FastGaussQuadrature
using Plots

mutable struct Distribution
    #### This structure will hold many arrays, it reduce the memory cost and time cost (through allocations)        ####
    #### But the user doesn't really have to deal with most of them. I think only Etot, Ltot, J, L, N are important ####
    #### V could be usefull to compute an entropy quickly, but one can easily do without                            ####

    # Orbital parameters
    J::Array{Float64, 3} # coupling coefficients
    L::Array{Float64, 1} # angular momentums

    # Quadrature parameters
    List_x::Array{Float64,1} # nodes of the quadratures
    List_w::Array{Float64, 1} # weigth of the quadratures
    List_Y::Array{Float64, 2} # spherical harmonics at the nodes

    # Microcanonical parameters
    N::Array{Float64, 1} # number of particles in each population
    Etot::Float64 # Total energy
    Ltot::Float64 # Total angular momentum

    # Order parameters
    β::Float64 # Lagrange multiplier for the energy
    γ::Float64 # Lagrange multiplier for the momentum
    M::Array{Float64, 1} # magnetisations -- one dimension but of size n_pop * lmax_half
    #### About the magnetization vector M : we store lexicographically first by the number of the population, then by the number of the mode ####
    #### i.e. M_{i * n_pop + l} corresponds to star i and mode l (not exactly, but that's the idea).                                         ####

    # pre-allocations
    U::Array{Float64, 1} # Will hold an array that is usefull to compute the cost function (upper integrals)
    V::Array{Float64, 1} # Will hold an array that is usefull to compute the cost function (lower integrals)
    dU::Array{Float64, 2} # Will hold an array that is usefull to compute the Jacobian
    dV::Array{Float64, 2} # Will hold an array that is usefull to compute the Jacobian

    cost::Array{Float64, 1} # Will hold the cost vector
    #### About the cost vector : first line will hold the consistency cost on energy, ####
    #### second line the cost on momentum, and next lines on magnetisations           ####
    jacobian::Array{Float64, 2} # Will hold the Jacobian
    #### About the Jacobian : Jᵢ,ⱼ = ∂Costᵢ / ∂xⱼ                                                                                ####
    #### i.e First line will contains the derivative of Etot, second line the derivatives of Ltot and next lines the derivatives of M ####
    #### and first column will contain the derivatives with respect to β, second column w.r.t. γ and next comumns w.r.t. M            ####

    # normalisations
    Norms::Array{Float64, 1} # to normalize the cost function
end

function Distribution(J::Array{Float64, 3}, L::Array{Float64, 1}, N::Array{Float64, 1}, Etot::Float64, Ltot::Float64;
                                        β = nothing, γ = nothing, M = nothing)
    #### Helps to create a distribution object when some parameters are unknown ####

    # getting the constants we will need
    n_pop = size(J)[1] # number of populations
    lmax_half = size(J)[3] # number of interesting modes
    lmax = 2 * lmax_half # number of modes
    size_mag = lmax_half * n_pop # size of the magnitudes vector
    size_jac = lmax_half * n_pop + 2 # size of the jacobian matrix
    n_quad = max(100, 2 * lmax) # number of quadratures
    N_stars = sum(N) # total number of stars

    # let's preallocate the arrays we will need
    List_x , List_w = gausslegendre(n_quad) # will hold the position and weight (respectively) of the quadrature nodes
    List_Y = zeros(lmax_half, n_quad) # Will hold the spherical harmonics
    U = zeros(size_jac) # Will hold an array that is usefull to compute the cost function
    V = zeros(n_pop) # Will hold an array that is usefull to compute the cost function
    dU = zeros(size_jac, size_jac) # Will hold an array that is usefull to compute the Jacobian
    dV = zeros(size_jac, size_jac) # Will hold an array that is usefull to compute the Jacobian
    cost = zeros(size_jac) # Will hold the cost vector
    jacobian = zeros(size_jac, size_jac) # Will hold the Jacobian of the cost function
    Norms = zeros(size_mag) # Will hold the vector used to renormalize every line of the cost function to 1

    # let's compute the spherical harmonics
    Pₗ = zeros(lmax_half) # this is temporary (allocates as early as possible an array used by Spherical_harmonics_a!)
    Yₗ = zeros(lmax_half) # this is temporary (allocates as early as possible an array used by Spherical_harmonics_a!)
    for i = 1:n_quad
        Spherical_harmonics_a!(List_x[i], lmax, Yₗ = Yₗ, Pₗ = Pₗ) # Yₗ now holds the spherical harmonics at this point x
        List_Y[:, i] = Yₗ
    end

    # let's compute the norms
    for star = 1:n_pop
        for l_half = 1:lmax_half
            l = 2 * l_half
            index = (star - 1) * lmax_half + l_half

            Norms[index] = (1 / N[star]) * sqrt((4 * π) / (2 * l + 1)) # || Mₗ[K] || <= N[K] * sqrt((2l + 1) / 4 π)
        end
    end

    # Let's get the initial point β, γ, M, by choice at random
    if isnothing(β)
        r = 2 * rand()
        β = r * (10 / (N_stars * J[1, 1, 1])) # Choice around an order of magnitude, from Takacs (fig. 4 left)
    end

    if isnothing(γ)
        r = 2 * rand()
        γ = r * (1.0 / L[1]) # Choice around an order of magnitude, from Roupas (fig. 7 right)
    end

    if isnothing(M)
        M = zeros(size_mag)

        for star = 1:n_pop # I will try to use a and b as indexes for all loops on populations
            for l_half = 1:lmax_half
                l = 2 * l_half
                index = (star - 1) * lmax_half + l_half

                r = 2 * rand() - 1 # random value between -1 and 1
                M[index] = N[star] * sqrt((2 * l + 1) / (4 * π)) * r
            end
        end
    else
        if (length(M) != lmax_half * n_pop)
            println("Error : lenght(M) = $(length(M)) whereas lmax = $(lmax) and n_pop = $(n_pop)")
        end
    end

    # Let's call the standard initializer for Distribution in order to create our distribution
    dist = Distribution(J, L, List_x, List_w, List_Y, N, Etot, Ltot, β, γ, M, U, V, dU, dV, cost, jacobian, Norms)

    return(dist)
end

function reset!(dist::Distribution)
    #### Resets V, cost, jacobian to the value they have when all integrals are 0.0 ####
    # getting the constants we will need
    n_pop = size(dist.J)[1] # number of populations
    lmax_half = size(dist.J)[3] # number of interesting modes
    lmax = 2 * lmax_half # number of modes
    size_mag = lmax_half * n_pop # size of the magnitudes vector
    size_jac = lmax_half * n_pop + 2 # size of the jacobian matrix

    # Access the dist's array we will need, using passage by reference
    Etot, Ltot, M, V = dist.Etot, dist.Ltot, dist.M, dist.V
    cost, jacobian = dist.cost, dist.jacobian

    # reset V to 0
    fill!(V, 0.0)

    # resest cost to [0, 0, M]
    cost[1] = Etot
    cost[2] = Ltot
    for a = 1:n_pop
        for l_half = 1:lmax_half
            line = (a - 1) * lmax_half + l_half + 2
            cost[line] = M[line - 2]
        end
    end

    # reset jacobian to [0, 0, ....; 0, 0, ....; 0, 0, Identity]
    for column= 1:(size_jac)
        for line = 1:(size_jac)
            jacobian[line, column] = 0.0

            if ((column == line) && (line > 2))
                jacobian[line, column] += 1.0
            end
        end
    end
end

function reset_v2!(dist::Distribution, star_1::Int64)
    #### Resets the lines of U, dU, dV we will use in this iteration over star_1 to 0.0 ####
    # getting the constants we will need
    n_pop = size(dist.J)[1] # number of stellar populations
    lmax_half = size(dist.J)[3] # number of interesting modes
    size_jac = lmax_half * n_pop + 2 # size of the jacobian matrix

    # Access the dist's array we will need, using passage by reference
    U, dU, dV = dist.U, dist.dU, dist.dV

    # Reset the lines of U
    U[1] = 0.0
    U[2] = 0.0
    for l_half = 1:lmax_half
        line = (star_1 - 1) * lmax_half + l_half + 2
        U[line] = 0.0
    end

    # Reset the lines of dU
    for column = 1:size_jac
        dU[1, column] = 0.0
        dU[2, column] = 0.0
        for l_half = 1:lmax_half
            line = (star_1 - 1) * lmax_half + l_half + 2
            dU[line, column] = 0.0
        end
    end

    # Reset the lines of dV
    for column = 1:size_jac
        dV[column] = 0.0
    end
end

function Cost_and_Jacobian!(dist::Distribution)
    ####  Update the cost function and its jacobian, which are 2 of the dist's attributes  ####
    #### NOTE : in all our integrals we forget the factor 2π coming from the integral on ϕ ####
    ####        because it is both on the upper and lower integral.                        ####
    # getting the constants we will need
    n_pop = size(dist.J)[1] # number of stellar populations
    lmax_half = size(dist.J)[3] # number of interesting modes
    n_quad = length(dist.List_x) # number of quadrature nodes
    size_mag = lmax_half * n_pop # size of the magnitudes vector
    size_jac = lmax_half * n_pop + 2 # size of the jacobian matrix

    # Access the dist's array we will need, using passage by reference
    J, L, List_x, List_w, List_Y, N, Etot, Ltot, β, γ, M, U, V, dU, dV, Norms = dist.J, dist.L, dist.List_x, dist.List_w, dist.List_Y, dist.N, dist.Etot, dist.Ltot, dist.β, dist.γ, dist.M, dist.U, dist.V, dist.dU, dist.dV, dist.Norms
    cost, jacobian = dist.cost, dist.jacobian

    # Reseting dist.V, dist.cost, and dist.jacobian -- It does modify V, cost and jacobian here because they were passed by reference
    reset!(dist)

    # Now we can start to comput the integrals
    for star_1 = 1:n_pop
        # reset U[K], dU[K], dV[K] to 0
        # V[K] = dist.V[star_1] has aready been reset to 0
        reset_v2!(dist, star_1)

        # First let's compute U, V, dU, dV for this population
        #
        # U is the upper integral of the cost function : ∫dL (N[K]  * [ϵ; L; Y] * exp)
        # V is the lower integral of the cost function : ∫dL (exp)
        # the cost function looks something like C = ∫dK (U / V)
        #
        # dU is the jacobian of U
        # dV is the gradient of V
        # the Jacobian looks something like J = ∫dK (dU / V - U dVᵀ / V²)

        for i = 1:n_quad # We use a Gauss-Legendre quadrature to to integrate over L

            ##### computing the local parameters (wᵢ, ϵ, momentum)
            # local node of the quadrature
            cos_θ = List_x[i] # position
            weight = List_w[i] # weight

            # ϵ the local one-particle energy
            ϵ = 0.0
            for l_half = 1:lmax_half
                spherical_harmonic = List_Y[l_half, i] # local spherical harmonics

                for star_2 = 1:n_pop
                    index = (star_2 - 1) * lmax_half + l_half
                    ϵ -= spherical_harmonic * M[index] * J[star_2, star_1, l_half]
                end
            end

            # local one-particle momentum
            momentum = L[star_1] * cos_θ

            # local unnormalized distribution function
            exponential = exp(- β * ϵ + γ * momentum)

            #### precomputing some prefactors
            Ne_dθ = N[star_1] * weight * exponential # this prefactor is f * dL
            e_dθ = weight * exponential # this prefactor is f * dL / N[k]

            ###### updating all the upper integrals (U), lower integrals (V), and derivatives (dU, dV)

            ### update U[K]
            U[1] += 0.5 * Ne_dθ * ϵ
            U[2] += Ne_dθ * momentum

            for l_half = 1:lmax_half
                # here, notice that Mₗ[K] = ∫dL (N[K] Yₗ exp) does not need an integration over K.
                # that why each iteration on star_1 will fill diferrent cells of U, but there is no need for a second loop on star_2
                line = (star_1 - 1) * lmax_half + l_half + 2
                spherical_harmonic = List_Y[l_half, i]

                U[line] += Ne_dθ * spherical_harmonic
            end

            ### update V[K]
            V[star_1] += e_dθ

            ### update dU[K]
            # some prefactors we will need
            specific_pref_dM_dβ = -1.0 * Ne_dθ * ϵ
            specific_pref_dM_dγ = Ne_dθ * momentum
            specific_pref_dE_dM = -0.5 * Ne_dθ * (1 - β * ϵ)
            specific_pref_dL_dM = Ne_dθ * β * momentum
            specific_pref_dM_dM = Ne_dθ * β

            # upper left square of the matrix
            dU[1, 1] += -0.5 * Ne_dθ * ϵ^(2) # derivative of C_E with respect to β
            dU[1, 2] += 0.5 * Ne_dθ * momentum * ϵ # derivative of C_E with respect to γ
            dU[2, 1] += -1.0 * Ne_dθ * ϵ * momentum # derivative of C_L with respect to β
            dU[2, 2] += Ne_dθ * momentum^(2) # derivative of C_L with respect to γ

            # upper lines of the matrix
            for l_half = 1:lmax_half
                spherical_harmonic = List_Y[l_half, i]

                for star_2 = 1:n_pop
                    # here, E[K] and L[K] can be derived with respect to Mₗ[K'], that's why we need a second loop on star_2
                    column = (star_2 - 1) * lmax_half + l_half + 2

                    coupling_coefficient = J[star_2, star_1, l_half]
                    dU[1, column] += specific_pref_dE_dM * coupling_coefficient * spherical_harmonic # derivative of C_E with respect to Mₗ[K]
                    dU[2, column] += specific_pref_dL_dM * coupling_coefficient * spherical_harmonic # derivative of C_L with respect to Mₗ[K]
                end
            end

            # left column of the matrix
            for l_half = 1:lmax_half
                # here, the derivatives of Mₗ[K] WITH RESPECT TO β AND γ do not need a second integration over star_2
                line = (star_1 - 1) * lmax_half + l_half + 2
                spherical_harmonic = List_Y[l_half, i]

                dU[line, 1] += specific_pref_dM_dβ * spherical_harmonic # derivative of C_Mₗ[K] with respect to β
                dU[line, 2] += specific_pref_dM_dγ * spherical_harmonic # derivative of C_Mₗ[K] with respect to γ
            end

            # lower right square of the matrix
            for l_half = 1:lmax_half
                spherical_harmonic_column = List_Y[l_half, i]

                for star_2 = 1:n_pop
                    # here we have to derivate Mₗ[K] with respect to Mₗ[K'] so we need a second loop on star_2
                    # but the derivative of Mₗ[K] with respect to Mₗ[K'] does not need a third integration over star_3
                    column = (star_2 - 1) * lmax_half + l_half + 2
                    coupling_coefficient = J[star_2, star_1, l_half]

                    @inbounds @simd for ll_half = 1:lmax_half
                        line = (star_1 - 1) * lmax_half + ll_half + 2
                        spherical_harmonic_line = List_Y[ll_half, i]

                        dU[line, column] += specific_pref_dM_dM * coupling_coefficient * spherical_harmonic_column * spherical_harmonic_line
                    end
                end
            end

            ### update dV[K]
            # some prefactors we will need
            specific_pref_dM = e_dθ * β

            # leftmost 2 coefficients
            dV[1] += -1.0 * e_dθ * ϵ
            dV[2] += e_dθ * momentum

            # rightmost lmax_half * n_pop coefficients
            for l_half = 1:lmax_half
                spherical_harmonic = List_Y[l_half, i]

                for star_2 = 1:n_pop
                    # here, V[K] can be derived with respect to any Mₗ[K'] so we need a second loop on star_2
                    column = (star_2 - 1) * lmax_half + l_half +2

                    dV[column] += specific_pref_dM * J[star_2, star_1, l_half] * spherical_harmonic
                end
            end

        end ##### end of the integrals over L

        # Now that U[K], V[K], dU[K] and dV[K] are known, we can integrate them over K to update cost and jacobian
        # of course, we are inside a loop on star_1, which represents the integration over K

        V_of_K = V[star_1]

        ### update cost
        # Precisely, Cost = [Etot, Ltot, M] - ∫dK (U[K] / V[K]), and remember that cost has been reinitialized to [Etot; Ltot; M]
        cost[1] -= U[1] / V_of_K
        cost[2] -= U[2] / V_of_K
        for l_half = 1:lmax_half
            line = (star_1 - 1) * lmax_half + l_half + 2

            cost[line] -= U[line] / V_of_K
        end

        #### update Jacobian
        # Precisely, Jacobian = [0, 0, ...; 0, 0, ...; 0, 0, I] - ∫dK (dU'[K] / V[K] - U[K] dV[K]ᵀ / V[K]²),
        # and remember that Jacobian has already been reinitialized to [0, 0, ...; 0, 0, ...; 0, 0, I]

        ## part U' / V

        # upper left square of the martrix
        jacobian[1, 1] -= (dU[1, 1] / V_of_K)
        jacobian[1, 2] -= (dU[1, 2] / V_of_K)
        jacobian[2, 1] -= (dU[2, 1] / V_of_K)
        jacobian[2, 2] -= (dU[2, 2] / V_of_K)

        # upper lines of the matrix
        for star_2 = 1:n_pop
            for l_half = 1:lmax_half
                column = (star_2 - 1) * lmax_half + l_half + 2

                jacobian[1, column] -= (dU[1, column] / V_of_K)
                jacobian[2, column] -= (dU[2, column] / V_of_K)
            end
        end

        # left columns of the matrix
        for l_half = 1:lmax_half
            line = (star_1 - 1) * lmax_half + l_half + 2

            jacobian[line, 1] -= (dU[line, 1] / V_of_K)
            jacobian[line, 2] -= (dU[line, 2] / V_of_K)
        end

        # lower right square or the matrix
        for star_2 = 1:n_pop
            for l_half = 1:lmax_half
                column = (star_2 - 1) * lmax_half + l_half + 2

                @inbounds @simd for ll_half = 1:lmax_half
                    line = (star_1 - 1) * lmax_half + ll_half + 2

                    jacobian[line, column] -= (dU[line, column] / V_of_K)
                end
            end
        end

        ### part U dVᵀ / V²

        # upper two lines of the matrix
        for column = 1:(size_jac)
            jacobian[1, column] += (U[1] * dV[column] / (V_of_K^(2)))
            jacobian[2, column] += (U[2] * dV[column] / (V_of_K^(2)))
        end

        # lower n_pop * lmax_half lines of the matrix
        for column = 1:(size_jac)
            for l_half = 1:lmax_half
                line = (star_1 - 1) * lmax_half + l_half + 2

                jacobian[line, column] += (U[line] * dV[column] / (V_of_K^(2)))
            end
        end

        ##### End of the treatment of population a
    end

    ##### Normalisations
    # We want all dimensions of the cost vector to be normalized to 1
    # To do that, we will make the renormalizations {Etot -> 1; Ltot -> 1; ∀ l, ∀ K,  N[K] sqrt((2l + 1) / 4 π) -> 1}

    # for the cost
    cost[1] /= Etot
    cost[2] /= Ltot

    for star_1 = 1:n_pop
        for l_half = 1:lmax_half
            line = (star_1 - 1) * lmax_half + l_half + 2

            if (isinf(Norms[line - 2]) && (cost[line] == 0.0))
                cost[line] = 0.0
            else
                cost[line] *= Norms[line - 2]
            end
        end
    end

    # for the jacobian
    for column = 1:(size_jac)
        jacobian[1, column] /= Etot
        jacobian[2, column] /= Ltot

        for star_1 = 1:n_pop
            for l_half = 1:lmax_half
                line = (star_1 - 1) * lmax_half + l_half + 2

                if (isinf(Norms[line - 2]) && (jacobian[line, column] == 0.0))
                    jacobian[line, column] = 0.0
                else
                    jacobian[line, column] *= Norms[line - 2]
                end
            end
        end
    end

    # I don't need to update dist.U, dist.V, dist.dU, dist.dV, dist.cost, dist.jacobian because they were passed by reference
end

function get_f(dist::Distribution, star::Int64, i::Int64)
    #### Returns a Float. It is the distribution function for this distribution, for a given population (star), at a given quadrature (index i) ####
    #### This shouldn't be used on an unoptimized distribution                                                                                ####
    # Passing by reference the dist's fields we will need
    J, L, List_x, List_Y, N, β, γ, M, V = dist.J, dist.L, dist.List_x, dist.List_Y, dist.N, dist.β, dist.γ, dist.M, dist.V

    # getting the constants we will need
    n_pop = size(J)[1]
    lmax_half = size(J)[3]

    # getting ϵ : the local one-particle energy
    ϵ = 0.0
    for l_half = 1:lmax_half
        spherical_harmonic = List_Y[l_half, i] # local spherical harmonics

        for star_2 = 1:n_pop
            index = (star_2 - 1) * lmax_half + l_half
            ϵ -= spherical_harmonic * M[index] * J[star, star_2, l_half]
        end
    end

    # getting the local one-particle momentum
    momentum = L[star] * List_x[i] # local momentum

    # getting the local exponential
    exponential = exp(- β * ϵ + γ * momentum)

    # return
    f = N[star] * exponential / V[star]
    return(f) # f(K, L) = N[K] * exp / ∫dL exp
end

function plot_that_dist(x_axis::String,
                        Etot_norm::Float64, s::Float64, lmax::Int64,
                        dist::Distribution,
                        Bins_m::Array{Tuple{Float64,Float64}, 1}, n_pop_m::Int64,
                        Bins_a::Array{Tuple{Float64,Float64}, 1}, n_pop_a::Int64,
                        Bins_e::Array{Tuple{Float64,Float64}, 1}, n_pop_e::Int64)
    ####  Returns a plot. It is a heatmap of the distribution function with y = cos(θ) and x = m, a or e ####
    #### Options for input x_axis are "m", "a" and "e"
    # Getting the constants we will need
    n_pop = n_pop_m * n_pop_a * n_pop_e

    # Let's get our x axis
    if (x_axis == "m")
        List_m = zeros(n_pop_m)

        for k_m = 1:n_pop_m
            m_down, m_up = Bins_m[k_m]
            List_m[k_m] = 10^((log10(m_up) + log10(m_down)) / 2)
        end
    elseif (x_axis == "a")
        List_a = zeros(n_pop_a)

        for k_a = 1:n_pop_a
            a_down, a_up = Bins_a[k_a]
            List_a[k_a] = 10^((log10(a_up) + log10(a_down)) / 2)
        end
    elseif (x_axis == "e")
        List_e = zeros(n_pop_e)

        for k_e = 1:n_pop_e
            e_down, e_up = Bins_e[k_e]
            List_e[k_e] = (e_up + e_down) / 2
        end
    else
        return "Error : x_axis should be m, a or e"
    end

    # Let's get our y axis : the positions of the quadratures
    n_quad = max(100, 2 * lmax)
    List_x, List_w = gausslegendre(n_quad)

    # Create the array that will become a heatmap
    if (x_axis == "m")
        List_f = zeros(n_quad, n_pop_m)

        for k_m = 1:n_pop_m

            #first let's get the raw vector of densities
            i_true = 1
            for i = 1:n_quad
                x = List_x[i]

                # first we find the index of the distribution's quadrature closest to our x
                while (x > dist.List_x[i_true])
                    i_true +=1
                end

                # then we read the distribution function there
                for k_a = 1:n_pop_a
                    for k_e = 1:n_pop_e
                        index = k_e + n_pop_e * (k_a - 1) + n_pop_e * n_pop_a * (k_m - 1)
                        f = get_f(dist, index, i_true)

                        List_f[i, k_m] += f # To get the distribution in mass, we have to integrate over a and e
                    end
                end
            end

            # now let's normalize with max
            max_f = -Inf
            for i = 1:n_quad
                f = List_f[i, k_m]

                if (f > max_f)
                    max_f = f
                end
            end

            for i = 1:n_quad
                List_f[i, k_m] /= max_f
            end
        end

    elseif (x_axis == "a")
        List_f = zeros(n_quad, n_pop_a)

        for k_a = 1:n_pop_a

            #first let's get the raw vector of densities
            i_true = 1
            for i = 1:n_quad
                x = List_x[i]

                # first we find the index of the distribution's quadrature closest to our x
                while (x > dist.List_x[i_true])
                    i_true +=1
                end

                # then we read the distribution function there
                for k_m = 1:n_pop_m
                    for k_e = 1:n_pop_e
                        index = k_e + n_pop_e * (k_a - 1) + n_pop_e * n_pop_a * (k_m - 1)
                        f = get_f(dist, index, i_true)

                        List_f[i, k_a] += f # To get the distribution in sma, we have to integrate over m and e
                    end
                end
            end

            # now let's normalize with max
            max_f = -Inf
            for i = 1:n_quad
                f = List_f[i, k_a]

                if (f > max_f)
                    max_f = f
                end
            end

            for i = 1:n_quad
                List_f[i, k_a] /= max_f
            end
        end

    elseif (x_axis == "e")
        List_f = zeros(n_quad, n_pop_e)

        for k_e = 1:n_pop_e

            #first let's get the raw vector of densities
            i_true = 1
            for i = 1:n_quad
                x = List_x[i]

                # first we find the index of the distribution's quadrature closest to our x
                while (x > dist.List_x[i_true])
                    i_true +=1
                end

                # then we read the distribution function there
                for k_m = 1:n_pop_m
                    for k_a = 1:n_pop_a
                        index = k_e + n_pop_e * (k_a - 1) + n_pop_e * n_pop_a * (k_m - 1)
                        f = get_f(dist, index, i_true)

                        List_f[i, k_e] += f # To get the distribution in eccentricity, we have to integrate over m and a
                    end
                end
            end

            # now let's normalize with max
            max_f = -Inf
            for i = 1:n_quad
                f = List_f[i, k_e]

                if (f > max_f)
                    max_f = f
                end
            end

            for i = 1:n_quad
                List_f[i, k_e] /= max_f
            end
        end

    else
        return "Error : x_axis should be m, a or e"
    end

    # Let's make the plot
    plotlyjs()

    if (x_axis == "m")
        my_plot = heatmap(List_m, List_x, List_f,
                        c = :thermal, colorbar_title = "f_eq",
                        xscale = :log10,
                        xlabel = "m", ylabel = "cos(θ)",
                        title = "E = $(round(Etot_norm, digits = 3)), s = $(round(s, digits = 2))")

    elseif (x_axis == "a")
        my_plot = heatmap(List_a, List_x, List_f,
                        c = :thermal, colorbar_title = "f_eq",
                        xscale = :log10,
                        xlabel = "a", ylabel = "cos(θ)",
                        title = "E = $(round(Etot_norm, digits = 3)), s = $(round(s, digits = 2))")

    elseif (x_axis == "e")
        my_plot = heatmap(List_e, List_x, List_f,
                        c = :thermal, colorbar_title = "f_eq",
                        xlabel = "e", ylabel = "cos(θ)",
                        title = "E = $(round(Etot_norm, digits = 3)), s = $(round(s, digits = 2))")

    else
        return "Error : x_axis should be m, a or e"
    end

    # Return
    return(my_plot)
end

max_stars_repo_name: NathanMagnan/VRROOMe

max_issues_repo_name: NathanMagnan/VRROOMe

max_forks_repo_name: NathanMagnan/VRROOMe


text:
module Flight

using Reexport
@reexport using BenchmarkTools

#general
include("utils.jl")
include("misc.jl")
include("modeling.jl")
include("plotting.jl")

#math
include("quaternions.jl")
include("attitude.jl")

#environment & dynamics
include("geodesy.jl")
include("terrain.jl")
include("atmosphere.jl")

include("kinematics.jl")
include("dynamics.jl")

#aircraft components
include("components.jl")
include("airdata.jl")
include("propulsion.jl")
include("landinggear.jl")
include("aerodynamics.jl")

#IO interfaces
include("input.jl")
include("output.jl")

# #aircraft
include("aircraft.jl")
include("c172.jl")

@reexport using .Utils
@reexport using .Misc
@reexport using .Modeling
@reexport using .Plotting

@reexport using .Quaternions
@reexport using .Attitude

@reexport using .Geodesy
@reexport using .Terrain
@reexport using .Atmosphere

@reexport using .Kinematics
@reexport using .Dynamics

@reexport using .Components
@reexport using .Airdata
@reexport using .Propulsion
@reexport using .LandingGear
@reexport using .Aerodynamics

@reexport using .Input
@reexport using .Output

@reexport using .Aircraft
@reexport using .C172

export ftest

ftest() = println("Welcome to Flight")

end


text:
using PyPlot
using LinearAlgebra
using Statistics
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCellDG
include("../../../../test/useful_routines.jl")
include("../analytical-solver.jl")
include("../transformation-elasticity-solver.jl")
PS = PlaneStrainSolver
TES = TransformationElasticitySolver

function lame_lambda(k, m)
    return k - 2m / 3
end

function angular_position(points)
    cpoints = points[1, :] + im * points[2, :]
    return rad2deg.(angle.(cpoints))
end

function normalized_maxnorm_error(pd, exactpd)
    return maximum(abs.(pd .- exactpd)) / abs(exactpd)
end

function exact_normalized_potential_difference(solver, V01, V02, ΔG0)
    return (ΔG0 + 1e9PS.interface_potential_difference(solver, V01, V02)) /
           abs(ΔG0)
end

polyorder = 2
nelmts = 17
penaltyfactor = 1e3

width = 1.0              # mm
K1, K2 = 247.0, 192.0    # GPa
mu1, mu2 = 126.0, 87.0   # GPa

lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)
stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)

rho1 = 3.93e3           # Kg/m^3
rho2 = 3.68e3           # Kg/m^3
V01 = 1.0 / rho1
V02 = 1.0 / rho2

ΔG0Jmol = -14351.0  # J/mol
molarmass = 0.147
ΔG0 = ΔG0Jmol / molarmass

theta0 = -0.067
transfstress =
    CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

meshwidth = [width, width]
numqp = required_quadrature_order(polyorder) + 2

solverbasis = LagrangeTensorProductBasis(2, polyorder)
interfacecenter = [0.5, 0.5]
interfaceradius = 0.4
outerradius = 1.2
CFL = 0.25

basispts = interpolation_points(solverbasis)
dim,numpts = size(basispts)
cgmesh = CutCellDG.CGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], numpts)
dgmesh = CutCellDG.DGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], basispts)
levelset = CutCellDG.LevelSet(
    x -> -circle_distance_function(x, interfacecenter, interfaceradius),
    cgmesh,
    solverbasis,
)
coeffs0 = copy(CutCellDG.coefficients(levelset))
paddedmesh = CutCellDG.BoundaryPaddedMesh(cgmesh, 1)

nodalcoordinates = CutCellDG.nodal_coordinates(cgmesh)
elementsize = CutCellDG.element_size(cgmesh)
minelmtsize = minimum(elementsize)
maxelmtsize = maximum(elementsize)
penalty = penaltyfactor / minelmtsize * 0.5 * (lambda1 + mu1 + lambda2 + mu2)
tol = minelmtsize^(polyorder + 1)
boundingradius = 1.5 * maxelmtsize

angularposition = angular_position(nodalcoordinates .- interfacecenter)
sortidx = sortperm(angularposition)
angularposition = angularposition[sortidx]

cutmesh = CutCellDG.CutMesh(dgmesh, levelset)
refseedpoints, seedcellids = CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
spatialseedpoints =
    CutCellDG.map_to_spatial(refseedpoints, seedcellids, cutmesh)
spatialseedradius =
    vec(mapslices(norm, spatialseedpoints .- interfacecenter, dims = 1))
R0 = mean(spatialseedradius)
devR0 = std(spatialseedradius)
analyticalsolution = PS.CylindricalSolver(
    R0,
    outerradius,
    interfacecenter,
    lambda1,
    mu1,
    lambda2,
    mu2,
    theta0,
)
pd0 = TES.potential_difference_at_nodal_coordinates(
    cutmesh,
    solverbasis,
    levelset,
    stiffness,
    theta0,
    analyticalsolution,
    numqp,
    penalty,
    spatialseedpoints,
    seedcellids,
    V01,
    V02,
    ΔG0,
    tol,
    boundingradius,
)
exactpd0 =
    exact_normalized_potential_difference(analyticalsolution, V01, V02, ΔG0)
err0 = normalized_maxnorm_error(pd0, exactpd0)
sortedpd0 = pd0[sortidx]

yrange = 0.1*abs(exactpd0)
fig,ax = PyPlot.subplots()
ax.plot(angularposition,sortedpd0)
ax.plot(angularposition,exactpd0*ones(length(sortedpd0)),"--")
ax.set_ylim(exactpd0 - yrange, exactpd0 + yrange)
ax.grid()
fig

# coeffs1 = TES.step_levelset(
#     levelset,
#     pd0,
#     spatialseedpoints,
#     seedcellids,
#     paddedmesh,
#     tol,
#     boundingradius,
#     CFL,
# )
#
#
# CutCellDG.update_coefficients!(levelset, coeffs1)
# cutmesh = CutCellDG.CutMesh(dgmesh, levelset)
# refseedpoints, seedcellids = CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
# spatialseedpoints =
#     CutCellDG.map_to_spatial(refseedpoints, seedcellids, cutmesh)
# spatialseedradius =
#     vec(mapslices(norm, spatialseedpoints .- interfacecenter, dims = 1))
# R1 = TES.average(spatialseedradius)
# devR1 = std(spatialseedradius)
# analyticalsolution = PS.CylindricalSolver(
#     R1,
#     outerradius,
#     interfacecenter,
#     lambda1,
#     mu1,
#     lambda2,
#     mu2,
#     theta0,
# )
# pd1 = TES.potential_difference_at_nodal_coordinates(
#     cutmesh,
#     basis,
#     levelset,
#     stiffness,
#     theta0,
#     analyticalsolution,
#     numqp,
#     penalty,
#     spatialseedpoints,
#     seedcellids,
#     V01,
#     V02,
#     ΔG0,
#     tol,
#     boundingradius,
# )
# exactpd1 =
#     exact_normalized_potential_difference(analyticalsolution, V01, V02, ΔG0)
# err1 = normalized_maxnorm_error(pd1, exactpd1)
#
# sortedpd1 = pd1[sortidx]
#
#
#
#
# coeffs2 = TES.step_levelset(
#     levelset,
#     pd1,
#     spatialseedpoints,
#     seedcellids,
#     paddedmesh,
#     tol,
#     boundingradius,
#     CFL,
# )
#
# CutCellDG.update_coefficients!(levelset, coeffs2)
# cutmesh = CutCellDG.CutMesh(dgmesh, levelset)
# refseedpoints, seedcellids = CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
# spatialseedpoints =
#     CutCellDG.map_to_spatial(refseedpoints, seedcellids, cutmesh)
# spatialseedradius =
#     vec(mapslices(norm, spatialseedpoints .- interfacecenter, dims = 1))
# R2 = TES.average(spatialseedradius)
# devR2 = std(spatialseedradius)
# analyticalsolution = PS.CylindricalSolver(
#     R2,
#     outerradius,
#     interfacecenter,
#     lambda1,
#     mu1,
#     lambda2,
#     mu2,
#     theta0,
# )
# pd2 = TES.potential_difference_at_nodal_coordinates(
#     cutmesh,
#     basis,
#     levelset,
#     stiffness,
#     theta0,
#     analyticalsolution,
#     numqp,
#     penalty,
#     spatialseedpoints,
#     seedcellids,
#     V01,
#     V02,
#     ΔG0,
#     tol,
#     boundingradius,
# )
# exactpd2 =
#     exact_normalized_potential_difference(analyticalsolution, V01, V02, ΔG0)
# err2 = normalized_maxnorm_error(pd2, exactpd2)
# sortedpd2 = pd2[sortidx]
#
#
#
# CutCellDG.update_coefficients!(levelset,coeffs0)
# maxerridx = sortperm(abs.(pd1),rev=true)
# xq = nodalcoordinates[:, maxerridx[1:20]]
# xcp, xcellid = CutCellDG.closest_points_on_zero_levelset(
#     xq,
#     spatialseedpoints,
#     seedcellids,
#     levelset,
#     tol,
#     boundingradius,
# )

# coeffs3 = TES.step_levelset(
#     levelset,
#     pd2,
#     spatialseedpoints,
#     seedcellids,
#     paddedmesh,
#     tol,
#     boundingradius,
#     CFL,
# )


text:
using QSimulator

include("speed_test.jl")

speed_test(5)


text:
# Most of this file was copied from the PProf.jl package, and then adapted to
# export a profile of the current memory usage for a given object, via
# `MemorySummarySize.summarysize()` from this package.
# This code is pretty hacky, and I could probably do a better job re-using
# logic from the PProf package, but :shrug:.

# Import the PProf generated protobuf types from the PProf package:
import PProf
import PProf.perftools.profiles: ValueType, Sample, Function, Location, Line, Label
const PProfile = PProf.perftools.profiles.Profile

using ProtoBuf
using OrderedCollections

using MemoryInspector.MemorySummarySize: FieldResult

struct FieldTraversalNode
    path::String
    parent_field::Union{Nothing,FieldResult}
    name::String
    field::FieldResult
end

# Resolves from `key` to the index (zero-based) in the dict.
# Useful for the Strings table
#
# NOTE: We must use Int64 throughout this package (regardless of system word-size) b/c the
# proto file specifies 64-bit integers.
function _enter!(dict::OrderedDict{T, Int64}, key::T) where T
    if haskey(dict, key)
        return dict[key]
    else
        l = Int64(length(dict))
        dict[key] = l
        return l
    end
end

using Base.StackTraces: StackFrame

# TODO:
# - Mappings

"""
    @mem_pprof [kwargs...] object
    @mem_pprof webport=62261 out="mem-inspect-profile" [...] object

Collects detailed breakdown of the recursive size information for `object`, via
the functionality provided by MemoryInspector.jl, and exports it to a profile in the
`pprof` format, and (optionally) opens a `pprof` web-server for interactively
viewing the results using the PProf.jl package.

The following flags are available, copied from the PProf.jl package:

If `web=true`, the web-server is opened in the background. Re-running `@pprof()` will refresh
the web-server to use the new output.

If you manually edit the output file or want to point at an existing profile object,
`PProf.refresh()` will refresh the server without overwriting the output file.
`PProf.kill()` will kill the server.

# Arguments:
- `object`: The object (or expression) whose size information you want to profile.

# Keyword Arguments
- `web::Bool`: Whether to launch the `go tool pprof` interactive webserver for viewing results.
- `webhost::AbstractString`: If using `web`, which host to launch the webserver on.
- `webport::Integer`: If using `web`, which port to launch the webserver on.
- `out::String`: Filename for output.
- `drop_frames`: frames with function_name fully matching regexp string will be dropped from the samples,
                 along with their successors.
- `keep_frames`: frames with function_name fully matching regexp string will be kept, even if it matches drop_functions.
- `ui_relative_percentages`: Passes `-relative_percentages` to pprof. Causes nodes
  ignored/hidden through the web UI to be ignored from totals when computing percentages.
"""
macro mem_pprof(exs...)
    args = exs[1:end-1]
    e = exs[end]
    :(_pprof($(string(e)), $(esc(e)); $((esc(a) for a in args)...)))
end

function _pprof(top_level_path::String, @nospecialize(obj),
               ;
               web::Bool = true,
               webhost::AbstractString = "localhost",
               webport::Integer = 62261,  # Use a different port than PProf (chosen via rand(33333:99999))
               out::AbstractString = "mem-inspect-profile.pb.gz",
               #from_c::Bool = true,
               drop_frames::Union{Nothing, AbstractString} = nothing,
               keep_frames::Union{Nothing, AbstractString} = nothing,
               ui_relative_percentages::Bool = true,
            )
    field_summary = MemorySummarySize.summarysize(obj, parentname=top_level_path)

    period = UInt64(0x1)

    @assert !isempty(basename(out)) "`out=` must specify a file path to write to. Got unexpected: '$out'"
    if !endswith(out, ".pb.gz")
        out = "$out.pb.gz"
        @info "Writing output to $out"
    end

    string_table = OrderedDict{AbstractString, Int64}()
    enter!(string) = _enter!(string_table, string)
    enter!(::Nothing) = _enter!(string_table, "nothing")
    ValueType!(_type, unit) = ValueType(_type = enter!(_type), unit = enter!(unit))

    # Setup:
    enter!("")  # NOTE: pprof requires first entry to be ""
    # Functions need a uid, we'll use the pointer for the method instance
    seen_funcs = Set{UInt64}()
    funcs = Dict{UInt64, Function}()

    seen_locs = Set{UInt64}()
    locs  = Dict{UInt64, Location}()

    sample_type = [
        ValueType!("instances", "count"), # Mandatory
        ValueType!("size", "bytes")
    ]

    prof = PProfile(
        sample = [], location = [], _function = [],
        mapping = [], string_table = [],
        sample_type = sample_type, default_sample_type = 2, # size
        period = period, period_type = ValueType!("heap", "bytes")
    )

    if drop_frames !== nothing
        prof.drop_frames = enter!(drop_frames)
    end
    if keep_frames !== nothing
        prof.keep_frames = enter!(keep_frames)
    end

    # Mappings to PProf:
    #  Type of the element => Function
    #  instance objectid => IP
    #  path to the element => Stack Trace

    # start decoding backtraces
    location_ids = Vector{UInt64}()
    nodes_stack = Vector{Union{Nothing,FieldTraversalNode}}()
    # start with the root node argument.
    push!(nodes_stack, FieldTraversalNode("", nothing, top_level_path, field_summary))

    while !isempty(nodes_stack)
        node = pop!(nodes_stack)
        if node === nothing
            # nothing means we've finished iterating all children of a node, and are popping
            # up to the parent.
            pop!(location_ids)
            continue
        end
        field = node.field

        # At each node in the tree, we create an entry for that element
        if node.parent_field === nothing  # Special-case the first node
            path = node.name
        else
            path = MemoryInspector._path_str(node.path, node.parent_field, node.name)
        end

        ip = field.objectid

        # A backtrace consists of a set of IP (Instruction Pointers), each IP points
        # a single line of code and `litrace` has the necessary information to decode
        # that IP to a specific frame (or set of frames, if inlining occured).

        # if we have already seen this IP avoid decoding it again
        #if ip in seen_locs
        #    push!(location_ids, ip)
        #    continue
        #end
        #push!(seen_locs, ip)

        # Decode the IP into information about this stack frame (or frames given inlining)
        location = Location(;line=[])

        # Record a new frame for this object instance
        begin
            frame = field.type

            if ip == 0
                ip = UInt64(frame.uid)
            end
            location.address = UInt64(frame.uid)
            location.id = ip
            # TODO: Line numbers? (A line number is required, so we use fake number here)
            push!(location.line, Line(function_id = ip, line = 1))

            # Store the function in our functions dict
            funcProto = Function()
            funcProto.id = ip
            funcProto.name = enter!(string(frame))
            # TODO: Line numbers?
            funcProto.start_line = convert(Int64, 1)
            # TODO: File names?
            # TODO: We could consider generating temp files with the `dump()` of the
            # type? :/ that's not ideal...
            # Or maybe we can find the file from the method table for the type?
            #file = Base.find_source_file(file)
            file = "<nofile>"
            funcProto.filename   = enter!(file)
            funcProto.system_name = funcProto.name
            funcs[ip] = funcProto
        end

        locs[ip] = location
        push!(location_ids, ip)

        begin
            # Record for this element in the struct tree.
            flat_size = field.size - sum([0, (c.size for (_,c) in field.children)...])
            value = [
                1,              # events  TODO: should we actually put the size here instead?
                flat_size,     # size of the element
            ]
            push!(prof.sample, Sample(;
                    label= [Label(key=enter!(""), str=enter!(path))],
                    location_id = reverse(location_ids),
                    value = value,
                ))
        end

        push!(nodes_stack, nothing)
        for (name, child_field) in field.children
            push!(nodes_stack, FieldTraversalNode(path, field, name, child_field))
        end
    end

    # Build Profile
    prof.string_table = collect(keys(string_table))
    # If from_c=false funcs and locs should NOT contain C functions
    prof._function = collect(values(funcs))
    prof.location  = collect(values(locs))

    # Write to disk
    open(out, "w") do io
        writeproto(io, prof)
    end

    if web
        refresh_pprof(webhost = webhost, webport = webport, file = out,
                      ui_relative_percentages = ui_relative_percentages,
        )
    end

    out
end

const proc = Ref{Union{Base.Process, Nothing}}(nothing)

"""
    refresh_pprof(; webhost = "localhost", webport = 57599, file = "profile.pb.gz",
                    ui_relative_percentages = true)
Start or restart the go pprof webserver.
- `webhost::AbstractString`: Which host to launch the webserver on.
- `webport::Integer`: Which port to launch the webserver on.
- `file::String`: Profile file to open.
- `ui_relative_percentages::Bool`: Passes `-relative_percentages` to pprof. Causes nodes
  ignored/hidden through the web UI to be ignored from totals when computing percentages.
"""
function refresh_pprof(; webhost::AbstractString = "localhost",
                         webport::Integer = 57599,
                         file::AbstractString = "profile.pb.gz",
                         ui_relative_percentages::Bool = true,
                      )

    if proc[] === nothing
        # The first time, register an atexit hook to kill the web server.
        atexit(kill_pprof)
    else
        # On subsequent calls, restart the pprof web server.
        Base.kill(proc[])
    end

    relative_percentages_flag = ui_relative_percentages ? "-relative_percentages" : ""

    proc[] = PProf.pprof_jll.pprof() do pprof_path
        open(pipeline(`$pprof_path -http=$webhost:$webport $relative_percentages_flag $file`))
    end
end

"""
    kill_pprof()
Kills the pprof server if running.
"""
function kill_pprof()
    if proc[] !== nothing
        Base.kill(proc[])
        proc[] = nothing
    end
end


text:
module MyUtils

export all_nacols, drop_all_nacols!
export floydgen
export boxplot

include("dataframes.jl")
include("misc.jl")
include("plotlyplots.jl")

end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

## Unary operators ##

"""
    conj!(A)

Transform an array to its complex conjugate in-place.

See also [`conj`](:func:`conj`).
"""
function conj!{T<:Number}(A::AbstractArray{T})
    for i in eachindex(A)
        A[i] = conj(A[i])
    end
    return A
end

for f in (:-, :~, :conj, :sign)
    @eval begin
        function ($f)(A::AbstractArray)
            F = similar(A)
            RF, RA = eachindex(F), eachindex(A)
            if RF == RA
                for i in RA
                    F[i] = ($f)(A[i])
                end
            else
                for (iF, iA) in zip(RF, RA)
                   F[iF] = ($f)(A[iA])
                end
            end
            return F
        end
    end
end

(-)(A::AbstractArray{Bool}) = reshape([ -A[i] for i in eachindex(A) ], size(A))

real(A::AbstractArray) = reshape([ real(x) for x in A ], size(A))
imag(A::AbstractArray) = reshape([ imag(x) for x in A ], size(A))

function !(A::AbstractArray{Bool})
    F = similar(A)
    RF, RA = eachindex(F), eachindex(A)
    if RF == RA
        for i in RA
            F[i] = !A[i]
        end
    else
        for (iF, iA) in zip(RF, RA)
            F[iF] = !A[iA]
        end
    end
    return F
end

## Binary arithmetic operators ##

promote_array_type(F, ::Type, ::Type, T::Type) = T
promote_array_type{S<:Real, A<:AbstractFloat}(F, ::Type{S}, ::Type{A}, ::Type) = A
promote_array_type{S<:Integer, A<:Integer}(F, ::Type{S}, ::Type{A}, ::Type) = A
promote_array_type{S<:Integer, A<:Integer}(::typeof(./), ::Type{S}, ::Type{A}, T::Type) = T
promote_array_type{S<:Integer, A<:Integer}(::typeof(.\), ::Type{S}, ::Type{A}, T::Type) = T
promote_array_type{S<:Integer}(::typeof(./), ::Type{S}, ::Type{Bool}, T::Type) = T
promote_array_type{S<:Integer}(::typeof(.\), ::Type{S}, ::Type{Bool}, T::Type) = T
promote_array_type{S<:Integer}(F, ::Type{S}, ::Type{Bool}, T::Type) = T

for f in (:+, :-, :div, :mod, :&, :|, :xor)
    @eval ($f)(A::AbstractArray, B::AbstractArray) =
        _elementwise($f, promote_eltype_op($f, A, B), A, B)
end
function _elementwise(op, ::Type{Any}, A::AbstractArray, B::AbstractArray)
    promote_shape(A, B) # check size compatibility
    return broadcast(op, A, B)
end
function _elementwise{T}(op, ::Type{T}, A::AbstractArray, B::AbstractArray)
    F = similar(A, T, promote_shape(A, B))
    RF, RA, RB  = eachindex(F), eachindex(A), eachindex(B)
    if RF == RA == RB
        for i in RA
            @inbounds F[i] = op(A[i], B[i])
        end
    else
        for (iF, iA, iB) in zip(RF, RA, RB)
            @inbounds F[iF] = op(A[iA], B[iB])
        end
    end
    return F
end

for f in (:.+, :.-, :.*, :./, :.\, :.^, :.÷, :.%, :.<<, :.>>, :div, :mod, :rem, :&, :|, :xor)
    @eval begin
        function ($f){T}(A::Number, B::AbstractArray{T})
            R = promote_op($f, typeof(A), T)
            S = promote_array_type($f, typeof(A), T, R)
            S === Any && return [($f)(A, b) for b in B]
            F = similar(B, S)
            RF, RB = eachindex(F), eachindex(B)
            if RF == RB
                for i in RB
                    @inbounds F[i] = ($f)(A, B[i])
                end
            else
                for (iF, iB) in zip(RF, RB)
                    @inbounds F[iF] = ($f)(A, B[iB])
                end
            end
            return F
        end
        function ($f){T}(A::AbstractArray{T}, B::Number)
            R = promote_op($f, T, typeof(B))
            S = promote_array_type($f, typeof(B), T, R)
            S === Any && return [($f)(a, B) for a in A]
            F = similar(A, S)
            RF, RA = eachindex(F), eachindex(A)
            if RF == RA
                for i in RA
                    @inbounds F[i] = ($f)(A[i], B)
                end
            else
                for (iF, iA) in zip(RF, RA)
                    @inbounds F[iF] = ($f)(A[iA], B)
                end
            end
            return F
        end
    end
end

# familiar aliases for broadcasting operations of array ± scalar (#7226):
(+)(A::AbstractArray{Bool},x::Bool) = A .+ x
(+)(x::Bool,A::AbstractArray{Bool}) = x .+ A
(-)(A::AbstractArray{Bool},x::Bool) = A .- x
(-)(x::Bool,A::AbstractArray{Bool}) = x .- A
(+)(A::AbstractArray,x::Number) = A .+ x
(+)(x::Number,A::AbstractArray) = x .+ A
(-)(A::AbstractArray,x::Number) = A .- x
(-)(x::Number,A::AbstractArray) = x .- A

## data movement ##

function flipdim{T}(A::Array{T}, d::Integer)
    nd = ndims(A)
    1 ≤ d ≤ nd || throw(ArgumentError("dimension $d is not 1 ≤ $d ≤ $nd"))
    sd = size(A, d)
    if sd == 1 || isempty(A)
        return copy(A)
    end

    B = similar(A)

    nnd = 0
    for i = 1:nd
        nnd += Int(size(A,i)==1 || i==d)
    end
    if nnd==nd
        # flip along the only non-singleton dimension
        for i = 1:sd
            B[i] = A[sd+1-i]
        end
        return B
    end

    d_in = size(A)
    leading = d_in[1:(d-1)]
    M = prod(leading)
    N = length(A)
    stride = M * sd

    if M==1
        for j = 0:stride:(N-stride)
            for i = 1:sd
                ri = sd+1-i
                B[j + ri] = A[j + i]
            end
        end
    else
        if isbits(T) && M>200
            for i = 1:sd
                ri = sd+1-i
                for j=0:stride:(N-stride)
                    offs = j + 1 + (i-1)*M
                    boffs = j + 1 + (ri-1)*M
                    copy!(B, boffs, A, offs, M)
                end
            end
        else
            for i = 1:sd
                ri = sd+1-i
                for j=0:stride:(N-stride)
                    offs = j + 1 + (i-1)*M
                    boffs = j + 1 + (ri-1)*M
                    for k=0:(M-1)
                        B[boffs + k] = A[offs + k]
                    end
                end
            end
        end
    end
    return B
end

"""
    rotl90(A)

Rotate matrix `A` left 90 degrees.

```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rotl90(a)
2×2 Array{Int64,2}:
 2  4
 1  3
```
"""
function rotl90(A::AbstractMatrix)
    ind1, ind2 = indices(A)
    B = similar(A, (ind2,ind1))
    n = first(ind2)+last(ind2)
    for i=indices(A,1), j=ind2
        B[n-j,i] = A[i,j]
    end
    return B
end

"""
    rotr90(A)

Rotate matrix `A` right 90 degrees.

```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rotr90(a)
2×2 Array{Int64,2}:
 3  1
 4  2
```
"""
function rotr90(A::AbstractMatrix)
    ind1, ind2 = indices(A)
    B = similar(A, (ind2,ind1))
    m = first(ind1)+last(ind1)
    for i=ind1, j=indices(A,2)
        B[j,m-i] = A[i,j]
    end
    return B
end
"""
    rot180(A)

Rotate matrix `A` 180 degrees.

```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rot180(a)
2×2 Array{Int64,2}:
 4  3
 2  1
```
"""
function rot180(A::AbstractMatrix)
    B = similar(A)
    ind1, ind2 = indices(A,1), indices(A,2)
    m, n = first(ind1)+last(ind1), first(ind2)+last(ind2)
    for j=ind2, i=ind1
        B[m-i,n-j] = A[i,j]
    end
    return B
end
"""
    rotl90(A, k)

Rotate matrix `A` left 90 degrees an integer `k` number of times.
If `k` is zero or a multiple of four, this is equivalent to a `copy`.

```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rotl90(a,1)
2×2 Array{Int64,2}:
 2  4
 1  3

julia> rotl90(a,2)
2×2 Array{Int64,2}:
 4  3
 2  1

julia> rotl90(a,3)
2×2 Array{Int64,2}:
 3  1
 4  2

julia> rotl90(a,4)
2×2 Array{Int64,2}:
 1  2
 3  4
```
"""
function rotl90(A::AbstractMatrix, k::Integer)
    k = mod(k, 4)
    k == 1 ? rotl90(A) :
    k == 2 ? rot180(A) :
    k == 3 ? rotr90(A) : copy(A)
end
"""
    rotr90(A, k)

Rotate matrix `A` right 90 degrees an integer `k` number of times. If `k` is zero or a
multiple of four, this is equivalent to a `copy`.

```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rotr90(a,1)
2×2 Array{Int64,2}:
 3  1
 4  2

julia> rotr90(a,2)
2×2 Array{Int64,2}:
 4  3
 2  1

julia> rotr90(a,3)
2×2 Array{Int64,2}:
 2  4
 1  3

julia> rotr90(a,4)
2×2 Array{Int64,2}:
 1  2
 3  4
```
"""
rotr90(A::AbstractMatrix, k::Integer) = rotl90(A,-k)
"""
    rot180(A, k)

Rotate matrix `A` 180 degrees an integer `k` number of times.
If `k` is even, this is equivalent to a `copy`.

```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rot180(a,1)
2×2 Array{Int64,2}:
 4  3
 2  1

julia> rot180(a,2)
2×2 Array{Int64,2}:
 1  2
 3  4
```
"""
rot180(A::AbstractMatrix, k::Integer) = mod(k, 2) == 1 ? rot180(A) : copy(A)

## Transpose ##

"""
    transpose!(dest,src)

Transpose array `src` and store the result in the preallocated array `dest`, which should
have a size corresponding to `(size(src,2),size(src,1))`. No in-place transposition is
supported and unexpected results will happen if `src` and `dest` have overlapping memory
regions.
"""
transpose!(B::AbstractMatrix, A::AbstractMatrix) = transpose_f!(transpose, B, A)

"""
    ctranspose!(dest,src)

Conjugate transpose array `src` and store the result in the preallocated array `dest`, which
should have a size corresponding to `(size(src,2),size(src,1))`. No in-place transposition
is supported and unexpected results will happen if `src` and `dest` have overlapping memory
regions.
"""
ctranspose!(B::AbstractMatrix, A::AbstractMatrix) = transpose_f!(ctranspose, B, A)
function transpose!(B::AbstractVector, A::AbstractMatrix)
    indices(B,1) == indices(A,2) && indices(A,1) == 1:1 || throw(DimensionMismatch("transpose"))
    copy!(B, A)
end
function transpose!(B::AbstractMatrix, A::AbstractVector)
    indices(B,2) == indices(A,1) && indices(B,1) == 1:1 || throw(DimensionMismatch("transpose"))
    copy!(B, A)
end
function ctranspose!(B::AbstractVector, A::AbstractMatrix)
    indices(B,1) == indices(A,2) && indices(A,1) == 1:1 || throw(DimensionMismatch("transpose"))
    ccopy!(B, A)
end
function ctranspose!(B::AbstractMatrix, A::AbstractVector)
    indices(B,2) == indices(A,1) && indices(B,1) == 1:1 || throw(DimensionMismatch("transpose"))
    ccopy!(B, A)
end

const transposebaselength=64
function transpose_f!(f,B::AbstractMatrix,A::AbstractMatrix)
    inds = indices(A)
    indices(B,1) == inds[2] && indices(B,2) == inds[1] || throw(DimensionMismatch(string(f)))

    m, n = length(inds[1]), length(inds[2])
    if m*n<=4*transposebaselength
        @inbounds begin
            for j = inds[2]
                for i = inds[1]
                    B[j,i] = f(A[i,j])
                end
            end
        end
    else
        transposeblock!(f,B,A,m,n,first(inds[1])-1,first(inds[2])-1)
    end
    return B
end
function transposeblock!(f,B::AbstractMatrix,A::AbstractMatrix,m::Int,n::Int,offseti::Int,offsetj::Int)
    if m*n<=transposebaselength
        @inbounds begin
            for j = offsetj+(1:n)
                for i = offseti+(1:m)
                    B[j,i] = f(A[i,j])
                end
            end
        end
    elseif m>n
        newm=m>>1
        transposeblock!(f,B,A,newm,n,offseti,offsetj)
        transposeblock!(f,B,A,m-newm,n,offseti+newm,offsetj)
    else
        newn=n>>1
        transposeblock!(f,B,A,m,newn,offseti,offsetj)
        transposeblock!(f,B,A,m,n-newn,offseti,offsetj+newn)
    end
    return B
end

function ccopy!(B, A)
    RB, RA = eachindex(B), eachindex(A)
    if RB == RA
        for i = RB
            B[i] = ctranspose(A[i])
        end
    else
        for (i,j) = zip(RB, RA)
            B[i] = ctranspose(A[j])
        end
    end
end

"""
    transpose(A)

The transposition operator (`.'`).
"""
function transpose(A::AbstractMatrix)
    ind1, ind2 = indices(A)
    B = similar(A, (ind2, ind1))
    transpose!(B, A)
end
function ctranspose(A::AbstractMatrix)
    ind1, ind2 = indices(A)
    B = similar(A, (ind2, ind1))
    ctranspose!(B, A)
end
ctranspose{T<:Real}(A::AbstractVecOrMat{T}) = transpose(A)

transpose(x::AbstractVector) = [ transpose(v) for i=of_indices(x, OneTo(1)), v in x ]
ctranspose{T}(x::AbstractVector{T}) = T[ ctranspose(v) for i=of_indices(x, OneTo(1)), v in x ]

# see discussion in #18364 ... we try not to widen type of the resulting array
# from cumsum or cumprod, but in some cases (+, Bool) we may not have a choice.
rcum_promote_type{T<:Number}(op, ::Type{T}) = promote_op(op, T)
rcum_promote_type{T}(op, ::Type{T}) = T

# handle sums of Vector{Bool} and similar.   it would be nice to handle
# any AbstractArray here, but it's not clear how that would be possible
rcum_promote_type{T,N}(op, ::Type{Array{T,N}}) = Array{rcum_promote_type(op,T), N}

for (f, f!, fp, op) = ((:cumsum, :cumsum!, :cumsum_pairwise!, :+),
                       (:cumprod, :cumprod!, :cumprod_pairwise!, :*) )
    # in-place cumsum of c = s+v[range(i1,n)], using pairwise summation
    @eval function ($fp){T}(v::AbstractVector, c::AbstractVector{T}, s, i1, n)
        local s_::T # for sum(v[range(i1,n)]), i.e. sum without s
        if n < 128
            @inbounds s_ = v[i1]
            @inbounds c[i1] = ($op)(s, s_)
            for i = i1+1:i1+n-1
                @inbounds s_ = $(op)(s_, v[i])
                @inbounds c[i] = $(op)(s, s_)
            end
        else
            n2 = n >> 1
            s_ = ($fp)(v, c, s, i1, n2)
            s_ = $(op)(s_, ($fp)(v, c, ($op)(s, s_), i1+n2, n-n2))
        end
        return s_
    end

    @eval function ($f!)(result::AbstractVector, v::AbstractVector)
        li = linearindices(v)
        li != linearindices(result) && throw(DimensionMismatch("input and output array sizes and indices must match"))
        n = length(li)
        if n == 0; return result; end
        i1 = first(li)
        @inbounds result[i1] = v1 = v[i1]
        n == 1 && return result
        ($fp)(v, result, v1, i1+1, n-1)
        return result
    end

    @eval function ($f){T}(v::AbstractVector{T})
        return ($f!)(similar(v, rcum_promote_type($op, T)), v)
    end
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
#
# Move coordinates to minimum perioddic images
#

function image(x :: Float64, input :: MDInput)
  image = x%input.side
  if image <= -input.side/2.0
    image = image + input.side
  elseif image > input.side/2.0
    image = image - input.side
  end
  return image
end


text:
using GeometryTypes
using ComputerVision
using Test

digital_camera = DigitalCamera()

camera_model = get_model(digital_camera)
pictures = get_pictures(digital_camera)
@test typeof(camera_model) <: AbstractCameraModel
@test isnothing(pictures)

intrinsics = get_intrinsics(camera_model)
@test typeof(intrinsics) <: AbstractIntrinsicParameters
f = get_focal_length(intrinsics)
@test f == 50
@test @inferred get_focal_length(intrinsics) == 50


text:
# Example on how to plot Y_N for a bi-objective discrete optimization case
# July 2017

using PyPlot

title(L"$Y_N$")
xlabel(L"$z_1$")
ylabel(L"$z_2$")
axis([0, 40, 0, 40])
plot( z1, z2, color="green", linestyle="", marker="o", label="YN" )

max_stars_repo_path: examples/programs and data/plotYn.jl

max_stars_repo_name: faithyyyy/vOptSolver

max_issues_repo_path: examples/programs and data/plotYn.jl

max_issues_repo_name: faithyyyy/vOptSolver

max_forks_repo_path: examples/programs and data/plotYn.jl

max_forks_repo_name: faithyyyy/vOptSolver


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
using GenerativeAD.Models
using BSON
using StatsBase: fit!, predict, sample

using Flux
using MLDataPattern

s = ArgParseSettings()
@add_arg_table! s begin
   "max_seed"
		required = true
		arg_type = Int
		help = "seed"
	"dataset"
		required = true
		arg_type = String
		help = "dataset"
	"anomaly_classes"
		arg_type = Int
		default = 10
		help = "number of anomaly classes"
	"method"
		arg_type = String
		default = "leave-one-out"
		help = "method for data creation -> \"leave-one-out\" or \"leave-one-in\" "
    "contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, max_seed, anomaly_classes, method, contamination = parsed_args

modelname ="DeepSVDD"

function sample_params()
	# first sample the number of layers
	nlayers = rand(2:4)
	kernelsizes = reverse((3,5,7,9)[1:nlayers])
	channels = reverse((16,32,64,128)[1:nlayers])
	scalings = reverse((1,2,2,2)[1:nlayers])

	par_vec = (
		2 .^(3:8),
		10f0 .^(-4:-3),
		10f0 .^(-4:-3),
		[true, false],
		2 .^ (6:7),
		["relu", "swish", "tanh"],
		["soft-boundary", "one-class"],
		[0.01f0, 0.1f0, 0.5f0, 0.99f0], # paper 0.1
		[1e-6], # from paper
		1:Int(1e8)
	)
	argnames = (
		:zdim,
		:lr_ae,
		:lr_svdd,
		:batchnorm,
		:batch_size,
		:activation,
		:objective,
		:nu,
		:decay,
		:init_seed
	)
	parameters = (;zip(argnames, map(x->sample(x, 1)[1], par_vec))...)
	return merge(parameters, (nlayers=nlayers, kernelsizes=kernelsizes,
		channels=channels, scalings=scalings))
end

function fit(data, parameters)
	# define models (Generator, Discriminator)
	all_parameters = merge(
		parameters,
		(
			idim=size(data[1][1])[1:3],
			iters=5000,
			check_every=30,
			ae_check_every = 200,
			patience=10,
			ae_iters=10000
		)
	)
	svdd = GenerativeAD.Models.conv_ae_constructor(;all_parameters...) |> gpu

	# define optimiser
	try
		global info, fit_t, _, _, _ = @timed fit!(svdd, data, all_parameters)
	catch e
		println("Error caught => $(e).")
		return (fit_t = NaN, model = nothing, history = nothing, n_parameters = NaN), []
	end

	training_info = (
		fit_t = fit_t,
		model = info[2]|>cpu,
		history_svdd = info[1][1], # losses through time
		history_ae = info[1][2],
		npars = info[3], # number of parameters
		iters = info[4] # optim iterations of model
		)

	return training_info, [(x -> GenerativeAD.Models.anomaly_score_gpu(svdd, x), parameters)]
end

#_________________________________________________________________________________________________

try_counter = 0
max_tries = 10*max_seed
cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
while try_counter < max_tries
	parameters = sample_params()

	for seed in 1:max_seed
		for i in 1:anomaly_classes
			savepath = datadir("experiments/images_$(method)$cont_string/$(modelname)/$(dataset)/ac=$(i)/seed=$(seed)")
			mkpath(savepath)

			data = GenerativeAD.load_data(dataset, seed=seed, anomaly_class_ind=i, method=method, contamination=contamination)

			@info "Trying to fit $modelname on $dataset with parameters $(parameters)..."
			@info "Train/validation/test splits: $(size(data[1][1], 4)) | $(size(data[2][1], 4)) | $(size(data[3][1], 4))"
			@info "Number of features: $(size(data[1][1])[1:3])"

			# here, check if a model with the same parameters was already tested
			if GenerativeAD.check_params(savepath, parameters)
				# fit model
				training_info, results = fit(data, parameters)
				# saving model separately
				if training_info.model !== nothing
					tagsave(joinpath(savepath, savename("model", parameters, "bson", digits=5)), Dict("model"=>training_info.model), safe = true)
					training_info = merge(training_info, (model = nothing,))
				end
				save_entries = merge(training_info, (modelname = modelname, seed = seed, dataset = dataset, anomaly_class = i,
					contamination=contamination))
				# now loop over all anomaly score funs
				for result in results
					GenerativeAD.experiment(result..., data, savepath; save_entries...)
				end
				global try_counter = max_tries + 1
			else
				@info "Model already present, sampling new hyperparameters..."
				global try_counter += 1
			end
		end
	end
end
(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 37b018c4-34b0-11eb-234d-4fbeaa683873
using LinearAlgebra, SparseArrays, Plots, DifferentialEquations

# ╔═╡ baed3ace-34ad-11eb-0cf6-19357238b1de
md"## Heat equation with Dirichlet boundary conditions

The domain is given by the one-dimensional rod $\Omega = (0,L)$ and the time by $t \in (0,T).$ The rod has the physical properties
- thermal conductivity $\lambda$,
- specific heat capacity $c$,
- mass density $\rho$.

### Heat Equation

The diffusivity constant is given by $\alpha = \frac{\lambda}{c \rho}$. The heat equation is noted as

$\dot{\vartheta}(t,x) = \alpha \frac{\partial^2}{\partial x^2} \vartheta(t,x)$

for $(t,x) \in (0,T) \times \Omega$ with initial condition

$\vartheta(0,x) = \vartheta_{0}(x)$

for $x \in \overline{\Omega}$ and boundary conditions  $\qquad \vartheta(\cdot,0) = \vartheta(\cdot,L) = 0$.
"

# ╔═╡ 8f5d25dc-34ae-11eb-2aa7-39e601f648a1
md"#### Code: implementing physical constants"

# ╔═╡ 8f493cca-34ae-11eb-3e60-31ab68674568
# Length of rod
L = 0.5

# ╔═╡ 8f149f6a-34ae-11eb-01f2-5983e256a7de
# Thermal conductivity
λ = 45.0

# ╔═╡ 8efa0286-34ae-11eb-1952-e790d141e19b
# Specific heat capacity
c = 480.0

# ╔═╡ 8ecb6bba-34ae-11eb-2c16-819675412725
# Mass density
ρ = 7800.0

# ╔═╡ e6433430-34af-11eb-318d-2199461e5ec9
# Thermal diffusivity
α = λ/(c*ρ)

# ╔═╡ 37c91fae-34b0-11eb-0f82-639173602ac0
md"### Spatial approximation

The one-dimensional rod is discretized as a one-dimensional grid with N points and the finite discretization is noted by

$\Delta x = \frac{L}{N-1}.$

The resulting grid points have the position $x^{0} = 0$, $x^{n} = n~\Delta x$ and $x^{N-1} = L$. The second order derivative is approximated using the [Taylor series](https://en.wikipedia.org/wiki/Taylor_series) as

$\left. \frac{\partial^2 f(x)}{\partial x^2} \right\rvert_{x = \tilde{x}} \approx \frac{1}{\Delta x^2} \left[ f(\tilde{x} - \Delta x) - 2 f(\tilde{x}) + f(\tilde{x} + \Delta x)  \right].$

This [Finite Difference](https://en.wikipedia.org/wiki/Finite_difference_method) scheme is written in matrix-vector form as

$\frac{\partial^2 \vartheta(t,x)}{\partial x^2} \approx \frac{1}{\Delta x^2}
\begin{pmatrix}
-2 & 1 & 0 & \cdots &  & 0 \\
1 & -2 & 1 & 0  & \cdots & \vdots \\
0 & 1 & -2 & 1 & \ddots &  \\
 & & \ddots & \ddots & \ddots &  & \\
& & & 1 & -2 & 1 \\
& & &  & 1 & -2 \\
\end{pmatrix}
~
\begin{pmatrix}
\vartheta(t, x^{0}) \\
\vartheta(t, x^{1}) \\
\vdots \\
\vartheta(t, x^{N-1}) \\
\end{pmatrix}.$
"

# ╔═╡ 65049b1a-3549-11eb-33ee-9dfae609164d
# Number of grid elements
N = 101

# ╔═╡ a36fd440-3549-11eb-0e93-0d36244400ef
# Finite discretization
Δx = L/(N-1)

# ╔═╡ c18a976a-3549-11eb-2991-799f558e4493
# 1-dimensional grid
xspan = 0 : Δx : L

# ╔═╡ d7d5be46-3549-11eb-3a12-afab5ab4bce5
# Diffusion matrix
M = spdiagm(-1 => ones(N-1), 0 => -2*ones(N), 1 => ones(N-1));

# ╔═╡ 045066e2-354a-11eb-164b-09fc7095de47
Matrix(M)[1:5,1:5]

# ╔═╡ 37822ec8-34b0-11eb-175d-7745382c2955
md"### Boundary conditions

The Dirichlet boundary condition states that the value at both sides is fixed as

$\qquad \vartheta(\cdot,0) = \vartheta(\cdot, x^{0}) = 0$

on the left boundary and

$\vartheta(\cdot,L) = \vartheta(\cdot, x^{N-1})  = 0$

on the right boundary. This means, there is no dynamical (e.g. diffusive) behaviour at both boundaries and thus the first and the last row of the diffusion matrix M is set to 0. "

# ╔═╡ 3768b7e0-34b0-11eb-0772-f5f4c95cf8f5
# First row
M[1,1] = M[1,2] = 0;

# ╔═╡ 37529500-34b0-11eb-1ee4-b5837bea23c9
# Last row
M[end,end] = M[end,end-1] = 0;

# ╔═╡ 373652f0-34b0-11eb-2fea-4f14b3ee9140
md"### Heat equation as ODE

After the spatial approximation the heat equation has a form of

$\dot{\theta}(t) = \frac{\alpha}{\Delta x^2} ~ M ~ \theta(t)$

with M as diffusion matrix and $\quad \theta(t) = \left( \vartheta(t, x^{0}), \cdots, \vartheta(t, x^{N-1}) \right)^{\top}$. This is an ordinary differential equation (ODE) that can be solved with common solvers like forward Euler method or Runge-Kutta scheme."

# ╔═╡ 371c6d52-34b0-11eb-12ab-a51e10fe33c5
# Heat Equation as ODE
function heat_eq(dθ, θ, p, t)
    return dθ .= (1/Δx^2) * α * M * θ
end

# ╔═╡ 37013e12-34b0-11eb-3d23-9d43ca3f06a3
md"### Initial Conditions

The initial data of the original heat equation is assumed as $\quad \vartheta_{0}(x) := m \left[ L~x - x^2 \right]~$ with  $m > 0$. This function is approximated to gain the initial conditions of the ODE

$\theta(0) = \left( \vartheta_{0}(x^{0}), \cdots , \vartheta_{0}(x^{N-1}) \right)^{\top}.$
"

# ╔═╡ 78b7fa20-354c-11eb-3754-d5d515d0a2d5
#Initial conditions: θ(0) = m * (L * x - x^2)
function initial_data(x, params)
    # α = params[1] # unused here
    L = params[2]
    m = params[3]

    return m*( L * x - x^2)
end

# ╔═╡ 77906a42-354c-11eb-397e-9dca847ab964
# Initial heat distribution
θ₀ = zeros(N);

# ╔═╡ 93696cf4-354d-11eb-0af3-9b27a8186fbd
# Amplification of initial data
@bind amp html"<input type='range' min='0' max='100' step='1'>"

# ╔═╡ a91c8b80-354d-11eb-03ec-8f4386a69708
m = amp

# ╔═╡ 86f0f21c-354d-11eb-2a32-299183f9118b
# Parameters
param = [α, L, amp]

# ╔═╡ 666a8878-354d-11eb-22dc-e94a72ede22d
let

	for i = 1 : length(xspan)
		θ₀[i] = initial_data(xspan[i], param)
	end

	plot(xspan, θ₀, xlabel="Position x", ylabel="Temperature", legend=false)
end

# ╔═╡ cfc0e592-354d-11eb-2b06-458d198c38e5
md"### Numerical stability

Before the approximated heat equation can be simulated, the [numerical stability](https://en.wikipedia.org/wiki/Numerical_stability) has to be proven with the [von Neumann stability analysis](https://en.wikipedia.org/wiki/Von_Neumann_stability_analysis).

In case of the forward Euler method one finds the approach

$\vartheta(t^{k+1}, x^{n}) = \vartheta(t^{k}, x^{n}) + \alpha \frac{\Delta t}{\Delta x^2} \left( \vartheta(t^{k}, x^{n-1}) - 2 ~ \vartheta(t^{k}, x^{n}) + \vartheta(t^{k}, x^{n+1}) \right)$

or equivalent

$\vartheta(t^{k+1}, x^{n}) = \left(1 - 2 \alpha \frac{\Delta t}{\Delta x^2} \right) \vartheta(t^{k}, x^{n}) + \alpha \frac{\Delta t}{\Delta x^2}  \left( \vartheta(t^{k}, x^{n-1}) + \vartheta(t^{k}, x^{n+1}) \right).$

The sampling time $\Delta t$ has to be chosen to guarantee

$\left(1 - 2 \alpha \frac{\Delta t}{\Delta x^2} \right) > 0.$

Therefore, the upper limit of the sampling time is given by

$\Delta t < \frac{1}{2} \frac{\Delta x^2}{\alpha}.$
"

# ╔═╡ 25789968-3552-11eb-1e25-35c0d1b1c461
# Upper limit of sampling time
ul = 0.5 * Δx^2 / α

# ╔═╡ 56daee6e-3554-11eb-37cd-49ad4ac56805
# Sampling time
Δt = 0.8

# ╔═╡ 56bb225a-3554-11eb-0062-6d03b8ea772f
# time discretization
tspan = (0.0, 1000.0)

# ╔═╡ 569a20f0-3554-11eb-0479-f1e583620682
md"### Code: Solving the ODE"

# ╔═╡ 56790ed8-3554-11eb-3c51-b9b31a4c340e
# Solving the ODE
let
	prob = ODEProblem( heat_eq, θ₀, tspan ) # ODE Problem

	sol = solve(prob,Euler(),dt=Δt,progress=true, save_everystep=false,
				save_start=true) # Solving the ODE

	plot(xspan, sol.u[2], xlabel = "Position x", ylabel="Temperature", legend=false)
end

# ╔═╡ Cell order:
# ╟─baed3ace-34ad-11eb-0cf6-19357238b1de
# ╟─8f5d25dc-34ae-11eb-2aa7-39e601f648a1
# ╠═8f493cca-34ae-11eb-3e60-31ab68674568
# ╠═8f149f6a-34ae-11eb-01f2-5983e256a7de
# ╠═8efa0286-34ae-11eb-1952-e790d141e19b
# ╠═8ecb6bba-34ae-11eb-2c16-819675412725
# ╠═e6433430-34af-11eb-318d-2199461e5ec9
# ╟─37c91fae-34b0-11eb-0f82-639173602ac0
# ╠═37b018c4-34b0-11eb-234d-4fbeaa683873
# ╠═65049b1a-3549-11eb-33ee-9dfae609164d
# ╠═a36fd440-3549-11eb-0e93-0d36244400ef
# ╠═c18a976a-3549-11eb-2991-799f558e4493
# ╠═d7d5be46-3549-11eb-3a12-afab5ab4bce5
# ╠═045066e2-354a-11eb-164b-09fc7095de47
# ╟─37822ec8-34b0-11eb-175d-7745382c2955
# ╠═3768b7e0-34b0-11eb-0772-f5f4c95cf8f5
# ╠═37529500-34b0-11eb-1ee4-b5837bea23c9
# ╟─373652f0-34b0-11eb-2fea-4f14b3ee9140
# ╠═371c6d52-34b0-11eb-12ab-a51e10fe33c5
# ╟─37013e12-34b0-11eb-3d23-9d43ca3f06a3
# ╠═78b7fa20-354c-11eb-3754-d5d515d0a2d5
# ╠═77906a42-354c-11eb-397e-9dca847ab964
# ╠═93696cf4-354d-11eb-0af3-9b27a8186fbd
# ╠═a91c8b80-354d-11eb-03ec-8f4386a69708
# ╠═86f0f21c-354d-11eb-2a32-299183f9118b
# ╠═666a8878-354d-11eb-22dc-e94a72ede22d
# ╟─cfc0e592-354d-11eb-2b06-458d198c38e5
# ╠═25789968-3552-11eb-1e25-35c0d1b1c461
# ╠═56daee6e-3554-11eb-37cd-49ad4ac56805
# ╠═56bb225a-3554-11eb-0062-6d03b8ea772f
# ╟─569a20f0-3554-11eb-0479-f1e583620682
# ╠═56790ed8-3554-11eb-3c51-b9b31a4c340e