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

"""Mirror callback function

Function sets `+refs/*:refs/*` refspecs and `mirror` flag for remote reference.
"""
function mirror_callback(remote::Ptr{Ptr{Cvoid}}, repo_ptr::Ptr{Cvoid},
                         name::Cstring, url::Cstring, payload::Ptr{Cvoid})
    ensure_initialized()
    # Create the remote with a mirroring url
    fetch_spec = "+refs/*:refs/*"
    err = ccall((:git_remote_create_with_fetchspec, :libgit2), Cint,
                (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cstring, Cstring, Cstring),
                remote, repo_ptr, name, url, fetch_spec)
    err != 0 && return Cint(err)

    # And set the configuration option to true for the push command
    config = GitConfig(GitRepo(repo_ptr,false))
    name_str = unsafe_string(name)
    err= try set!(config, "remote.$name_str.mirror", true)
         catch; -1
         finally close(config)
         end
    err != 0 && return Cint(err)
    return Cint(0)
end

"""
    LibGit2.is_passphrase_required(private_key) -> Bool

Return `true` if the `private_key` file requires a passphrase, `false` otherwise.
"""
function is_passphrase_required(private_key::AbstractString)
    !isfile(private_key) && return false

    # In encrypted private keys, the second line is "Proc-Type: 4,ENCRYPTED"
    return open(private_key) do f
        readline(f)
        readline(f) == "Proc-Type: 4,ENCRYPTED"
    end
end

function user_abort()
    ensure_initialized()
    # Note: Potentially it could be better to just throw a Julia error.
    ccall((:giterr_set_str, :libgit2), Cvoid,
          (Cint, Cstring), Cint(Error.Callback),
          "Aborting, user cancelled credential request.")
    return Cint(Error.EUSER)
end

function prompt_limit()
    ensure_initialized()
    ccall((:giterr_set_str, :libgit2), Cvoid,
          (Cint, Cstring), Cint(Error.Callback),
          "Aborting, maximum number of prompts reached.")
    return Cint(Error.EAUTH)
end

function exhausted_abort()
    ensure_initialized()
    ccall((:giterr_set_str, :libgit2), Cvoid,
          (Cint, Cstring), Cint(Error.Callback),
          "All authentication methods have failed.")
    return Cint(Error.EAUTH)
end

function authenticate_ssh(libgit2credptr::Ptr{Ptr{Cvoid}}, p::CredentialPayload, username_ptr)
    ensure_initialized()
    cred = p.credential::SSHCredential
    revised = false

    # Use a filled credential as-is on the first pass. Reset password on successive calls.
    if p.first_pass && isfilled(cred)
        revised = true
    elseif !p.first_pass
        cred.pass = ""
    end

    # first try ssh-agent if credentials support its usage
    if p.use_ssh_agent && username_ptr != Cstring(C_NULL) && (!revised || !isfilled(cred))
        err = ccall((:git_cred_ssh_key_from_agent, :libgit2), Cint,
                    (Ptr{Ptr{Cvoid}}, Cstring), libgit2credptr, username_ptr)

        p.use_ssh_agent = false  # use ssh-agent only one time
        err == 0 && return Cint(0)
    end

    if p.use_env && (!revised || !isfilled(cred))
        if isempty(cred.user) && username_ptr != Cstring(C_NULL)
            cred.user = unsafe_string(username_ptr)
        end

        if haskey(ENV, "SSH_KEY_PATH")
            cred.prvkey = ENV["SSH_KEY_PATH"]
        elseif isempty(cred.prvkey)
            for keytype in ("rsa", "ecdsa")
                private_key_file = joinpath(homedir(), ".ssh", "id_$keytype")
                if isfile(private_key_file)
                    cred.prvkey = private_key_file
                    break
                end
            end
        end

        cred.pubkey = Base.get(ENV, "SSH_PUB_KEY_PATH") do
            default = cred.prvkey * ".pub"
            if isempty(cred.pubkey) && isfile(default)
                default
            else
                cred.pubkey
            end
        end

        cred.pass = Base.get(ENV, "SSH_KEY_PASS", cred.pass)

        revised = true
        p.use_env = false
    end

    if p.remaining_prompts > 0 && (!revised || !isfilled(cred))
        if isempty(cred.user) || username_ptr == Cstring(C_NULL)
            url = git_url(scheme=p.scheme, host=p.host)
            response = Base.prompt("Username for '$url'", default=cred.user)
            response === nothing && return user_abort()
            cred.user = response
        end

        url = git_url(scheme=p.scheme, host=p.host, username=cred.user)

        # For SSH we need a private key location
        last_private_key = cred.prvkey
        if !isfile(cred.prvkey) || !revised || !haskey(ENV, "SSH_KEY_PATH")
            response = Base.prompt("Private key location for '$url'", default=cred.prvkey)
            response === nothing && return user_abort()
            cred.prvkey = expanduser(response)

            # Only update the public key if the private key changed
            if cred.prvkey != last_private_key
                cred.pubkey = cred.prvkey * ".pub"
            end
        end

        # For SSH we need a public key location. Avoid asking about the public key as
        # typically this will just annoy users.
        stale = !p.first_pass && cred.prvkey == last_private_key && cred.pubkey != cred.prvkey * ".pub"
        if isfile(cred.prvkey) && (stale || !isfile(cred.pubkey))
            response = Base.prompt("Public key location for '$url'", default=cred.pubkey)
            response === nothing && return user_abort()
            cred.pubkey = expanduser(response)
        end

        # Ask for a passphrase when the private key exists and requires a passphrase
        if isempty(cred.pass) && is_passphrase_required(cred.prvkey)
            if Sys.iswindows()
                response = Base.winprompt(
                    "Your SSH Key requires a password, please enter it now:",
                    "Passphrase required", cred.prvkey; prompt_username=false)
                response === nothing && return user_abort()
                cred.pass = response[2]
            else
                response = Base.getpass("Passphrase for $(cred.prvkey)")
                response === nothing && return user_abort()
                cred.pass = response
                isempty(cred.pass) && return user_abort()  # Ambiguous if EOF or newline
            end
        end

        revised = true

        p.remaining_prompts -= 1
        p.remaining_prompts <= 0 && return prompt_limit()
    end

    if !revised
        return exhausted_abort()
    end
    return ccall((:git_cred_ssh_key_new, :libgit2), Cint,
                 (Ptr{Ptr{Cvoid}}, Cstring, Cstring, Cstring, Cstring),
                 libgit2credptr, cred.user, cred.pubkey, cred.prvkey, cred.pass)
end

function authenticate_userpass(libgit2credptr::Ptr{Ptr{Cvoid}}, p::CredentialPayload)
    ensure_initialized()
    cred = p.credential::UserPasswordCredential
    revised = false

    # Use a filled credential as-is on the first pass. Reset password on successive calls.
    if p.first_pass && isfilled(cred)
        revised = true
    elseif !p.first_pass
        cred.pass = ""
    end

    if p.use_git_helpers && (!revised || !isfilled(cred))
        git_cred = GitCredential(p.config, p.url)

         # Use `deepcopy` to ensure shredding the `git_cred` does not shred the `cred`s copy
        cred.user = something(git_cred.username, "")
        cred.pass = deepcopy(something(git_cred.password, ""))
        Base.shred!(git_cred)
        revised = true

        p.use_git_helpers = false
    end

    if p.remaining_prompts > 0 && (!revised || !isfilled(cred))
        url = git_url(scheme=p.scheme, host=p.host)
        username = isempty(cred.user) ? p.username : cred.user
        if Sys.iswindows()
            response = Base.winprompt(
                "Please enter your credentials for '$url'", "Credentials required",
                username; prompt_username=true)
            response === nothing && return user_abort()
            cred.user, cred.pass = response
        else
            response = Base.prompt("Username for '$url'", default=username)
            response === nothing && return user_abort()
            cred.user = response

            url = git_url(scheme=p.scheme, host=p.host, username=cred.user)
            response = Base.getpass("Password for '$url'")
            response === nothing && return user_abort()
            cred.pass = response
            isempty(cred.pass) && return user_abort()  # Ambiguous if EOF or newline
        end

        revised = true

        p.remaining_prompts -= 1
        p.remaining_prompts <= 0 && return prompt_limit()
    end

    if !revised
        return exhausted_abort()
    end

    return ccall((:git_cred_userpass_plaintext_new, :libgit2), Cint,
                 (Ptr{Ptr{Cvoid}}, Cstring, Cstring),
                 libgit2credptr, cred.user, cred.pass)
end

"""
    credential_callback(...) -> Cint

A LibGit2 credential callback function which provides different credential acquisition
functionality w.r.t. a connection protocol. The `payload_ptr` is required to contain a
`LibGit2.CredentialPayload` object which will keep track of state and settings.

The `allowed_types` contains a bitmask of `LibGit2.Consts.GIT_CREDTYPE` values specifying
which authentication methods should be attempted.

Credential authentication is done in the following order (if supported):
- SSH agent
- SSH private/public key pair
- Username/password plain text

If a user is presented with a credential prompt they can abort the prompt by typing `^D`
(pressing the control key together with the `d` key).

**Note**: Due to the specifics of the `libgit2` authentication procedure, when
authentication fails, this function is called again without any indication whether
authentication was successful or not. To avoid an infinite loop from repeatedly
using the same faulty credentials, we will keep track of state using the payload.

For addition details see the LibGit2 guide on
[authenticating against a server](https://libgit2.org/docs/guides/authentication/).
"""
function credentials_callback(libgit2credptr::Ptr{Ptr{Cvoid}}, url_ptr::Cstring,
                              username_ptr::Cstring, allowed_types::Cuint,
                              p::CredentialPayload)
    err = Cint(0)

    # Parse URL only during the first call to this function. Future calls will use the
    # information cached inside the payload.
    if isempty(p.url)
        p.url = unsafe_string(url_ptr)
        m = match(URL_REGEX, p.url)

        p.scheme = something(m[:scheme], SubString(""))
        p.username = something(m[:user], SubString(""))
        p.host = m[:host]

        # When an explicit credential is supplied we will make sure to use the given
        # credential during the first callback by modifying the allowed types. The
        # modification only is in effect for the first callback since `allowed_types` cannot
        # be mutated.
        if p.explicit !== nothing
            cred = p.explicit

            # Copy explicit credentials to avoid mutating approved credentials.
            # invalidation fix from cred being non-inferrable
            p.credential = Base.invokelatest(deepcopy, cred)

            if isa(cred, SSHCredential)
                allowed_types &= Cuint(Consts.CREDTYPE_SSH_KEY)
            elseif isa(cred, UserPasswordCredential)
                allowed_types &= Cuint(Consts.CREDTYPE_USERPASS_PLAINTEXT)
            else
                allowed_types &= Cuint(0)  # Unhandled credential type
            end
        elseif p.cache !== nothing
            cred_id = credential_identifier(p.scheme, p.host)

            # Perform a deepcopy as we do not want to mutate approved cached credentials
            if haskey(p.cache, cred_id)
                # invalidation fix from p.cache[cred_id] being non-inferrable
                p.credential = Base.invokelatest(deepcopy, p.cache[cred_id])
            end
        end

        p.first_pass = true
    else
        p.first_pass = false
    end

    # use ssh key or ssh-agent
    if isset(allowed_types, Cuint(Consts.CREDTYPE_SSH_KEY))
        if p.credential === nothing || !isa(p.credential, SSHCredential)
            p.credential = SSHCredential(p.username)
        end
        err = authenticate_ssh(libgit2credptr, p, username_ptr)
        err == 0 && return err
    end

    if isset(allowed_types, Cuint(Consts.CREDTYPE_USERPASS_PLAINTEXT))
        if p.credential === nothing || !isa(p.credential, UserPasswordCredential)
            p.credential = UserPasswordCredential(p.username)
        end
        err = authenticate_userpass(libgit2credptr, p)
        err == 0 && return err
    end

    # No authentication method we support succeeded. The most likely cause is
    # that explicit credentials were passed in, but said credentials are incompatible
    # with the requested authentication method.
    if err == 0
        if p.explicit !== nothing
            ensure_initialized()
            ccall((:giterr_set_str, :libgit2), Cvoid, (Cint, Cstring), Cint(Error.Callback),
                  "The explicitly provided credential is incompatible with the requested " *
                  "authentication methods.")
        end
        err = Cint(Error.EAUTH)
    end
    return err
end

function credentials_callback(libgit2credptr::Ptr{Ptr{Cvoid}}, url_ptr::Cstring,
                              username_ptr::Cstring, allowed_types::Cuint,
                              payloads::Dict{Symbol, Any})
    p = payloads[:credentials]
    return credentials_callback(libgit2credptr, url_ptr, username_ptr, allowed_types, p)
end

function fetchhead_foreach_callback(ref_name::Cstring, remote_url::Cstring,
                                    oid_ptr::Ptr{GitHash}, is_merge::Cuint, payload::Any)
    fhead_vec = payload::Vector{FetchHead}
    Base.push!(fhead_vec, FetchHead(unsafe_string(ref_name), unsafe_string(remote_url),
        unsafe_load(oid_ptr), is_merge == 1))
    return Cint(0)
end

struct CertHostKey
    parent  :: Cint
    mask    :: Cint
    md5     :: NTuple{16,UInt8}
    sha1    :: NTuple{20,UInt8}
    sha256  :: NTuple{32,UInt8}
    type    :: Cint
    hostkey :: Ptr{Cchar}
    len     :: Csize_t
end

function verify_host_error(message::AbstractString)
    printstyled(stderr, "$message\n", color = :cyan, bold = true)
end

function certificate_callback(
    cert_p :: Ptr{CertHostKey},
    valid  :: Cint,
    host_p :: Ptr{Cchar},
    data_p :: Ptr{Cvoid},
)::Cint
    valid != 0 && return Consts.CERT_ACCEPT
    host = unsafe_string(host_p)
    cert_type = unsafe_load(convert(Ptr{Cint}, cert_p))
    transport = cert_type == Consts.CERT_TYPE_TLS ? "TLS" :
                cert_type == Consts.CERT_TYPE_SSH ? "SSH" : nothing
    if !NetworkOptions.verify_host(host, transport)
        # user has opted out of host verification
        return Consts.CERT_ACCEPT
    end
    if transport == "TLS"
        # TLS verification is done before the callback and indicated with the
        # incoming `valid` flag, so if we get here then host verification failed
        verify_host_error("TLS host verification: the identity of the server `$host` could not be verified. Someone could be trying to man-in-the-middle your connection. It is also possible that the correct server is using an invalid certificate or that your system's certificate authority root store is misconfigured.")
        return Consts.CERT_REJECT
    elseif transport == "SSH"
        # SSH verification has to be done here
        files = NetworkOptions.ssh_known_hosts_files()
        cert = unsafe_load(cert_p)
        check = ssh_knownhost_check(files, host, cert)
        valid = false
        if check == Consts.LIBSSH2_KNOWNHOST_CHECK_MATCH
            valid = true
        elseif check == Consts.LIBSSH2_KNOWNHOST_CHECK_NOTFOUND
            if Sys.which("ssh-keyscan") !== nothing
                msg = "Please run `ssh-keyscan $host >> $(files[1])` in order to add the server to your known hosts file and then try again."
            else
                msg = "Please connect once using `ssh $host` in order to add the server to your known hosts file and then try again. You may not be allowed to log in (wrong user and/or no login allowed), but ssh will prompt you to add a host key for the server which will allow libgit2 to verify the server."
            end
            verify_host_error("SSH host verification: the server `$host` is not a known host. $msg")
        elseif check == Consts.LIBSSH2_KNOWNHOST_CHECK_MISMATCH
            verify_host_error("SSH host verification: the identity of the server `$host` does not match its known hosts record. Someone could be trying to man-in-the-middle your connection. It is also possible that the server has changed its key, in which case you should check with the server administrator and if they confirm that the key has been changed, update your known hosts file.")
        else
            @error("unexpected SSH known host check result", check)
        end
        return valid ? Consts.CERT_ACCEPT : Consts.CERT_REJECT
    end
    @error("unexpected transport encountered, refusing to validate", cert_type)
    return Consts.CERT_REJECT
end

struct KnownHost
    magic :: Cuint
    node  :: Ptr{Cvoid}
    name  :: Ptr{Cchar}
    key   :: Ptr{Cchar}
    type  :: Cint
end

function ssh_knownhost_check(
    files :: AbstractVector{<:AbstractString},
    host  :: AbstractString,
    cert  :: CertHostKey,
)
    key = unsafe_wrap(Array, cert.hostkey, cert.len)
    return ssh_knownhost_check(files, host, key)
end

function ssh_knownhost_check(
    files :: AbstractVector{<:AbstractString},
    host  :: AbstractString,
    key   :: Vector{Cchar},
)
    if (m = match(r"^(.+):(\d+)$", host)) !== nothing
        host = m.captures[1]
        port = parse(Int, m.captures[2])
    else
        port = 22 # default SSH port
    end
    len = length(key)
    mask = Consts.LIBSSH2_KNOWNHOST_TYPE_PLAIN |
           Consts.LIBSSH2_KNOWNHOST_KEYENC_RAW
    session = @ccall "libssh2".libssh2_session_init_ex(
        C_NULL :: Ptr{Cvoid},
        C_NULL :: Ptr{Cvoid},
        C_NULL :: Ptr{Cvoid},
        C_NULL :: Ptr{Cvoid},
    ) :: Ptr{Cvoid}
    for file in files
        ispath(file) || continue
        hosts = @ccall "libssh2".libssh2_knownhost_init(
            session :: Ptr{Cvoid},
        ) :: Ptr{Cvoid}
        count = @ccall "libssh2".libssh2_knownhost_readfile(
            hosts :: Ptr{Cvoid},
            file  :: Cstring,
            1     :: Cint, # standard OpenSSH format
        ) :: Cint
        if count < 0
            @warn("Error parsing SSH known hosts file `$file`")
            @ccall "libssh2".libssh2_knownhost_free(hosts::Ptr{Cvoid})::Cvoid
            continue
        end
        check = @ccall "libssh2".libssh2_knownhost_checkp(
            hosts  :: Ptr{Cvoid},
            host   :: Cstring,
            port   :: Cint,
            key    :: Ptr{Cchar},
            len    :: Csize_t,
            mask   :: Cint,
            C_NULL :: Ptr{Ptr{KnownHost}},
        ) :: Cint
        if check == Consts.LIBSSH2_KNOWNHOST_CHECK_MATCH ||
            check == Consts.LIBSSH2_KNOWNHOST_CHECK_MISMATCH
            @ccall "libssh2".libssh2_knownhost_free(hosts::Ptr{Cvoid})::Cvoid
            @assert 0 == @ccall "libssh2".libssh2_session_free(session::Ptr{Cvoid})::Cint
            return check
        else
            @ccall "libssh2".libssh2_knownhost_free(hosts::Ptr{Cvoid})::Cvoid
            if check == Consts.LIBSSH2_KNOWNHOST_CHECK_FAILURE
                @warn("Error searching SSH known hosts file `$file`")
            end
            continue
        end
    end
    # name not found in any known hosts files
    @assert 0 == @ccall "libssh2".libssh2_session_free(session::Ptr{Cvoid})::Cint
    return Consts.LIBSSH2_KNOWNHOST_CHECK_NOTFOUND
end

"C function pointer for `mirror_callback`"
mirror_cb() = @cfunction(mirror_callback, Cint, (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cstring, Cstring, Ptr{Cvoid}))
"C function pointer for `credentials_callback`"
credentials_cb() = @cfunction(credentials_callback, Cint, (Ptr{Ptr{Cvoid}}, Cstring, Cstring, Cuint, Any))
"C function pointer for `fetchhead_foreach_callback`"
fetchhead_foreach_cb() = @cfunction(fetchhead_foreach_callback, Cint, (Cstring, Cstring, Ptr{GitHash}, Cuint, Any))
"C function pointer for `certificate_callback`"
certificate_cb() = @cfunction(certificate_callback, Cint, (Ptr{CertHostKey}, Cint, Ptr{Cchar}, Ptr{Cvoid}))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
# Energy calculation

export
    total_energy,
    kinetic_energy,
    temperature,
    potential_energy

"""
    total_energy(s, neighbors=nothing)

Calculate the total energy of the system.
If the interactions use neighbor lists, the neighbors should be computed
first and passed to the function.
"""
total_energy(s, neighbors=nothing) = kinetic_energy(s) + potential_energy(s, neighbors)

kinetic_energy_noconvert(s) = sum(mass.(s.atoms) .* sum.(abs2, s.velocities)) / 2

"""
    kinetic_energy(s)

Calculate the kinetic energy of the system.
"""
function kinetic_energy(s::System{D, G, T}) where {D, G, T}
    ke = kinetic_energy_noconvert(s)
    # Convert energy to per mol if required
    if dimension(s.energy_units) == u"𝐋^2 * 𝐌 * 𝐍^-1 * 𝐓^-2"
        return T(uconvert(s.energy_units, ke * Unitful.Na))
    else
        return T(uconvert(s.energy_units, ke))
    end
end

const temp_conversion_factor = ustrip(u"nm^2 * u * K^-1 * ps^-2", Unitful.k)

"""
    temperature(system)

Calculate the temperature of a system from the kinetic energy of the atoms.
"""
function temperature(s::System{D, G, T}) where {D, G, T}
    ke = kinetic_energy_noconvert(s)
    df = 3 * length(s) - 3
    if unit(ke) == NoUnits
        k = T(temp_conversion_factor)
    else
        k = T(uconvert(u"K^-1" * unit(ke), Unitful.k))
    end
    return 2 * ke / (df * k)
end

function check_energy_units(E, energy_units)
    if unit(E) != energy_units
        error("System energy units are ", energy_units, " but encountered energy units ",
                unit(E))
    end
end

@inline @inbounds function potential_energy_nounit(inters, coord_i, coord_j, atom_i, atom_j,
                                        box_size, energy_units, weight_14::Bool=false)
    dr = vector(coord_i, coord_j, box_size)
    sum(inters) do inter
        if weight_14
            E = potential_energy(inter, dr, coord_i, coord_j, atom_i, atom_j,
                                    box_size, true)
        else
            E = potential_energy(inter, dr, coord_i, coord_j, atom_i, atom_j,
                                    box_size)
        end
        check_energy_units(E, energy_units)
        return ustrip(E)
    end
end

@views function potential_energy_inters(inters, coords, atoms, neighbors, box_size,
                                        energy_units, weights_14)
    coords_i, atoms_i = getindices_i(coords, neighbors), getindices_i(atoms, neighbors)
    coords_j, atoms_j = getindices_j(coords, neighbors), getindices_j(atoms, neighbors)
    @inbounds energies = potential_energy_nounit.((inters,), coords_i, coords_j,
                                atoms_i, atoms_j, (box_size,), energy_units, weights_14)
    return sum(energies) * energy_units
end

"""
    potential_energy(s, neighbors=nothing)

Calculate the potential energy of the system using the pairwise, specific and
general interactions.
If the interactions use neighbor lists, the neighbors should be computed
first and passed to the function.

    potential_energy(inter::PairwiseInteraction, vec_ij, coord_i, coord_j,
                     atom_i, atom_j, box_size)
    potential_energy(inter::SpecificInteraction, coords_i, coords_j,
                     box_size)
    potential_energy(inter::SpecificInteraction, coords_i, coords_j,
                     coords_k, box_size)
    potential_energy(inter::SpecificInteraction, coords_i, coords_j,
                     coords_k, coords_l, box_size)
    potential_energy(inter, system, neighbors=nothing)

Calculate the potential energy due to a given interation type.
Custom interaction types should implement this function.
"""
function potential_energy(s::System{D, false, T}, neighbors=nothing) where {D, T}
    n_atoms = length(s)
    potential = zero(T) * s.energy_units

    for inter in values(s.pairwise_inters)
        if inter.nl_only
            if isnothing(neighbors)
                error("An interaction uses the neighbor list but neighbors is nothing")
            end
            @inbounds for ni in 1:neighbors.n
                i, j, weight_14 = neighbors.list[ni]
                dr = vector(s.coords[i], s.coords[j], s.box_size)
                if weight_14
                    potential += potential_energy(inter, dr, s.coords[i], s.coords[j], s.atoms[i],
                                                    s.atoms[j], s.box_size, true)
                else
                    potential += potential_energy(inter, dr, s.coords[i], s.coords[j], s.atoms[i],
                                                    s.atoms[j], s.box_size)
                end
            end
        else
            for i in 1:n_atoms
                for j in (i + 1):n_atoms
                    dr = vector(s.coords[i], s.coords[j], s.box_size)
                    potential += potential_energy(inter, dr, s.coords[i], s.coords[j], s.atoms[i],
                                                    s.atoms[j], s.box_size)
                end
            end
        end
    end

    for inter_list in values(s.specific_inter_lists)
        potential += potential_energy(inter_list, s.coords, s.box_size)
    end

    for inter in values(s.general_inters)
        potential += potential_energy(inter, s, neighbors)
    end

    return uconvert(s.energy_units, potential)
end

function potential_energy(s::System{D, true, T}, neighbors=nothing) where {D, T}
    potential = zero(T) * s.energy_units

    pairwise_inters_nonl = filter(inter -> !inter.nl_only, values(s.pairwise_inters))
    if length(pairwise_inters_nonl) > 0
        potential += potential_energy_inters(pairwise_inters_nonl, s.coords, s.atoms,
                        neighbors.all, s.box_size, s.energy_units, false)
    end

    pairwise_inters_nl = filter(inter -> inter.nl_only, values(s.pairwise_inters))
    if length(pairwise_inters_nl) > 0 && length(neighbors.close.nbsi) > 0
        potential += potential_energy_inters(pairwise_inters_nl, s.coords, s.atoms,
                        neighbors.close, s.box_size, s.energy_units, neighbors.close.weights_14)
    end

    for inter_list in values(s.specific_inter_lists)
        potential += potential_energy(inter_list, s.coords, s.box_size)
    end

    for inter in values(s.general_inters)
        potential += potential_energy(inter, s, neighbors)
    end

    return uconvert(s.energy_units, potential)
end

@views function potential_energy(inter_list::InteractionList2Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    (box_size,)))
end

@views function potential_energy(inter_list::InteractionList3Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    coords[inter_list.ks], (box_size,)))
end

@views function potential_energy(inter_list::InteractionList4Atoms, coords, box_size)
    return sum(potential_energy.(inter_list.inters, coords[inter_list.is], coords[inter_list.js],
                                    coords[inter_list.ks], coords[inter_list.ls], (box_size,)))
end

ext: jl


text:
module Parab
    export len_curve

    function len_curve(n)
        # Parameters
        f = x -> x^2
        a = 0
        b = 1
        # Calculation
        sum = 0
        h = (b-a)/n
        h² = h^2
        x0 = a
        f0 = f(x0)
        for i in 1:n
            x1 = a*(1-i/n) + b*(i/n) # no rounding error by added many h
            f1 = f(x1)
            sum += sqrt((f1-f0)^2+h²)
            x0, f0 = x1, f1
        end
        sum
    end

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# DTM is used as an abbreviation for DecisionTreeModels

#=
    @info "DTM: BK Possibly add trnidx and validx to the resulting ensemble. This is relevant in case of a CV sampling which is performed. Otherwise it is not possible to reconstruct the Excel statistics after the model has run."
    @info "DTM: BK 'Remove' variables which are not used by a model from the SAS/VBA/CSharp code (e.g. dim command in VBA)"
    @info "DTM: BK Add 'time_finished' to model result (and possibly the time needed for the modelling)"
    @info "DTM: BK Need to ensure that the code runs smoothly even when no split is found (e.g. minWeight too big)"
    @info "DTM: BK Possibly introduce an option not to calculate certain intermediate metrics (scoring, sorting, ...) -> faster runtime"
    @info "DTM: BK To check if sortperm! performance has improved in newer Julia versions"
    @info "DTM: BK To fix this: the code should not rely on dtmtable.mappings and dtmtable.candmatwomax anymore (if possible). The pools have the equivalent information"
    @info "DTM: BK To add purity improvement for each split! Stop splitting if it is > some_specified_threshold"
    @info "DTM: BK Do not use df_name_vector anymore, but the symbol in the subset or simply propertynames(features)"
    @info "DTM: BK Add yes/no to the top of each printed DOT graph (in the PDF output) indicating that 'yes goes left'"
    @info "DTM: BK Try to avoid branching in critical functions; use ifelse or similar constructs (this might improve the performance)"
=#

module DecisionTrees

export run_model,run_model_actual,define_eltypevector,prepare_dataframe_for_dtm!,prep_data_from_df_modernized,run_legacy,prep_data_from_df

# include PooledArraysDTM module
include("PooledArraysDTM.jl")
import .PooledArraysDTM: PooledArray

import Core
import Dates
import Statistics
import Random
using Distributed
import DelimitedFiles
import PyCall
import DataFrames
import DataFrames: DataFrame
import JLD2
import OnlineStats
import ProgressMeter
import StatsBase

import Base: length,eltype,==,hash,iterate,append!,isless,resize!,convert
import Statistics.mean

Random.seed!(1234)

include("types.jl")
include("cross_validation_samplers.jl") # needs to be included early on as types are in here
include("sorting.jl")
include(joinpath("splittingmeasures", "splittingmeasures.jl"))
include("helper_functions.jl")
include("apply_tree_fn.jl")
include("build_tree.jl")
include("boosting.jl")
include("bagging.jl")
include("runmodel.jl")
include("runmodel_multi.jl")
include("charts.jl")
include("dot_graphs.jl")
include("roc.jl")
include("cross_validation.jl")
include("show.jl")
include("partialDependence.jl")
include("apply_tree_fn_unseen.jl")
include("unpack.jl")

global const pyModPandas = PyCall.PyNULL()
global const pyModxlsxwriter = PyCall.PyNULL()

global const emptyMeanVarSeries = OnlineStats.Series(OnlineStats.Mean(), OnlineStats.Variance())
global const emptyCustomVariance = CustomVariance()
global const globalConstAllowableMethodsForDefineCandidates = ("equalWeight", "basedOnWeightVector")
global const UInt8VECTORemptySplitDef = Vector{Splitdef{UInt8}}(undef, 0)
global const UInt16VECTORemptySplitDef = Vector{Splitdef{UInt16}}(undef, 0)
global const UInt8emptySplitDef = Splitdef(0, 0, Symbol(), Vector{UInt8}(), -Inf, 0.0, 0.0) # formerly const_default_splitdef
global const UInt16emptySplitDef = Splitdef(0, 0, Symbol(), Vector{UInt16}(), -Inf, 0.0, 0.0)
global const defaultModelName = "dtmresult"
global const defaultModelNameWtihCSVext = string(defaultModelName, ".csv")
global nLevelsThreshold = 2002
global const global_number_of_num_f_warning_mandatory_field = "Abort. Settings do not contain any value for number_of_num_features (which is mandatory)."
global const stars = "*******************************************************************************"
global const UInt8emptyRulepath = Array{Rulepath{UInt8}}(undef, 0)
global const UInt16emptyRulepath = Array{Rulepath{UInt16}}(undef, 0)
global const UInt8emptyLeaf = Leaf(0, NaN, NaN, NaN, -1, UInt8emptyRulepath, NaN, NaN, 0)
global const UInt16emptyLeaf = Leaf(0, NaN, NaN, NaN, -1, UInt16emptyRulepath, NaN, NaN, 0)
global const UInt8emptyNode = Node(0, 0, Array{UInt8}(undef, 0), deepcopy(UInt8emptyLeaf), deepcopy(UInt8emptyLeaf), deepcopy(UInt8emptyRulepath))
global const UInt16emptyNode = Node(0, 0, Array{UInt16}(undef, 0), deepcopy(UInt16emptyLeaf), deepcopy(UInt16emptyLeaf), deepcopy(UInt16emptyRulepath))
global const global_const_shift_cols = 5
global const global_pldamod_valid_types = [String,Float64]

global const global_nameOfSettingsSheet = "ModelSettings"
global const global_nameOfModelStatisticsSheet = "ModelStatistics"
global const global_nameOfTWOWayValidationSheet = "2Way_ValidationGraphs"

global const globalValidfitForStatsAndCharts = ["rawRelativities","unsmoothedPerScore","smoothedPerScore"]
global const global_statsperiter_header = ["Iteration" "Correlation" "Std Dev Trn" "Std Dev Val" "Std Dev Ratio" "Lift Trn" "Lift Val" "Reversals Trn" "Reversals Val" "RelDiff TrnObs vs ValObs" "RelDiff ValObs vs TrnFitted" "Min Rel Trn" "Max Rel Trn" "Min Observed Ratio Trn" "Max Observed Ratio Trn"  "Min Rel Val" "Max Rel Val" "Min Observed Ratio Val" "Max Observed Ratio Val" "Gini Num (Two Arg Fn) Unweighted Trn" "Gini Num (Two Arg Fn) Unweighted Val" "Gini Num Trn (Single Agrument Fn)" "Gini Num Val (Single Argument Fn)" "RSquared of Ratio Trn" "RSquared of Ratio Val" "RSquared of Numerator Trn" "RSquared of Numerator Val" "RSS of Ratio Trn" "RSS of Ratio Val" "RSS of Numerator Trn" "RSS of Numerator Val" "Average Poisson Error Trn" "Average Poisson Error Val"]
global const GLOBALperfMeasursesWhereMaxIsBest = ["Correlation","Lift Trn","Lift Val","Max Rel Trn","Max Observed Ratio Trn","Max Rel Val","Max Observed Ratio Val","Gini Num (Two Arg Fn) Unweighted Trn","Gini Num (Two Arg Fn) Unweighted Val","Gini Num Trn (Single Agrument Fn)","Gini Num Val (Single Argument Fn)","RSquared of Ratio Trn","RSquared of Ratio Val","RSquared of Numerator Trn","RSquared of Numerator Val","RSS of Ratio Trn","RSS of Ratio Val","RSS of Numerator Trn","RSS of Numerator Val"]

global const global_debug_csharp = false
global const global_byte_order_mark = '\ufeff'
global const CSHARP_VALID_CHARS = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ_123456789äëöüïÿâêîôûàèìòùáéíóúý" # the list is probably much longer, but this will do for now

function __init__()
	# the following lines may trigger the installation of the respective python packages
    copy!(pyModPandas, PyCall.pyimport_conda("pandas", "pandas"))
    copy!(pyModxlsxwriter, PyCall.pyimport_conda("xlsxwriter", "xlsxwriter"))
end

function get_sha1()
	dtpath = dirname(@__FILE__)
	current_path = pwd()
	cd(dtpath)
	sha1_of_dt = "n/a"
	try
		cmd = `git rev-parse HEAD`
		sha1_of_dt = read(cmd, String)[1:end - 1] # need to discard \n at the end
	catch
		sha1_of_dt = "unable to determine version. Git invocation possibly failed"
	end
	cd(current_path)
	return sha1_of_dt
end

global const DoubleQuote = '\"' # note, this is here at the end only because Notepad++ messes up the context highlighting after this definition

end # end Module DTM


text:
#**************************************************************************************
# test_Pythical_Model.jl
# ====================== part of the GeoEfficiency.jl package.
#
#
#
#**************************************************************************************

using GeoEfficiency
using Base.Test


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()
btns = ["btn-$i" for i in 1:6]
app.layout = html_div() do
    html_div([html_button(btn, id = btn) for btn in btns]),
    html_div(id = "output")
end

callback!(app,
    Output("output", "children"),
    [Input(btn, "n_clicks") for btn in btns]
) do args...
    isempty(callback_context().triggered) && throw(PreventUpdate())
    trigger = callback_context().triggered[1]
    name = split(trigger.prop_id, ".")[1]
    return "Just clicked $(name) for the $(trigger.value) time!"
end

run_server(app)


text:
module SlideBlock

export  solveSlideBlockProblemRBFS,
        solveSlideBlockProblemGSU,
        solveSlideBlockProblemGSBF,
        solveSlideBlockProblemGSBF2,
        solveSlideBlockProblemGSAS,
        solveSlideBlockProblemGSAS2

using AIMACore

import AIMACore: search, goal_test, step_cost, actions, result, heuristic
import Base: ==, show, hash

const E_INVALID_INPUT_LENGTH = "The input parameter is of a wrong length"
const E_INVALID_INPUT = "The input parameter is invalid"
const E_INVALID_STATE = "The state is invalid"
const E_INVALID_ACTION = "The action is invalid"

XY(n, N) = mod(n, N) + 1, mod(div(n,N), N) + 1

dist(t1, t2)    = abs(t1[1] - t2[1]) + abs(t1[2] - t2[2])
dist(n1, n2, N) = dist(XY(n1, N), XY(n2, N))

struct Grid{T} <: State
    loc::Vector{Tuple{Int, Int}}
    h::UInt
end

function Grid{T}(v::Vector{Int}) where T
    N = T*T
    @assert length(v) == N E_INVALID_INPUT
    loc = Vector{Tuple{Int, Int}}(undef, N)
    sum = UInt(0)
    for i = 1:N
        n = v[i]
        @assert n < N E_INVALID_INPUT
        loc[n+1] = XY(i-1, T)
        sum += UInt(n)
        sum *= UInt(N)
    end
    return Grid{T}(loc, sum)
end

Grid(T::Int) = Grid{T}([x for x = 0:T*T-1])

dist(g1::Grid{T}, g2::Grid{T}) where T = sum(dist.(g1.loc, g2.loc))
h2(g1::Grid{T}, g2::Grid{T})   where T = count(g1.loc .!= g2.loc)

==(g1::Grid, g2::Grid) = hash(g1) == hash(g2)

hash(g::Grid, h::UInt) = xor(g.h, h)

function convert(::Type{Matrix{Int}}, g::Grid{T}) where T
    N = T*T
    m = Matrix{Int}(undef, T, T)
    for i = 1:N
        t = g.loc[i]
        m[t[2], t[1]] = i - 1
    end
    return m
end

Matrix{Int}(g::Grid{T}) where T = convert(Matrix{Int}, g)

function show(io::IO, g::Grid{T}) where T
    m = Matrix{Int}(g)
    for j = 1:T-1
        show(io, m[j,:])
        println(io, "")
    end
    show(io, m[T,:])
end

struct Move <: Action
    d::Symbol
end
Move(d::AbstractString) = Move(Symbol(d))

mutable struct SlideBlockProblem{SA <: SearchAlgorithm, T} <: Problem
    initial_state::Grid{T}
    goal_state::Grid{T}
    search_algorithm::SA
    h::Function
end

SlideBlockProblem(init::Grid{T},
                  sa::SA, fh::Function) where {SA <: SearchAlgorithm, T} =
    SlideBlockProblem(init, Grid(T), sa, fh)

function result(problem::SlideBlockProblem{SA,T},
                state::Grid{T}, action::Move) where {SA <: SearchAlgorithm, T}
    N = UInt(T*T)
    loc = deepcopy(state.loc)
    x1, y1 = loc[1]
    xn, yn = action.d == :left  ? (x1 - 1, y1) :
             action.d == :right ? (x1 + 1, y1) :
             action.d == :up    ? (x1, y1 - 1) :
             action.d == :down  ? (x1, y1 + 1) :
             throw(ErrorException(E_INVALID_ACTION))

    i = findfirst(x-> x == (xn,yn), loc)
    if i != 0 && 1 <= xn <= T && 1 <= yn <= T
        loc[1] = (xn, yn)
        loc[i] = (x1, y1)
        #Readjusting the hash
        hn = UInt(i-1)*(N ^ (N-UInt(xn-1+T*(yn-1))))
        h1 = UInt(i-1)*(N ^ (N-UInt(x1-1+T*(y1-1))))
        h = state.h - hn + h1
        return Grid{T}(loc, h)
    else
        throw(ErrorException(E_INVALID_STATE))
    end
end

step_cost(problem::SlideBlockProblem, state::Grid, action::Move) = 1

heuristic(problem::SlideBlockProblem, state::Grid) = problem.h(state, problem.goal_state)

goal_test(problem::SlideBlockProblem, state::Grid) = (state == problem.goal_state)

function actions(problem::SlideBlockProblem, state::Grid{T}) where T
    ret=[]
    x1, y1 = state.loc[1]
    x1 != 1 && push!(ret, Move(:left))
    x1 != T && push!(ret, Move(:right))
    y1 != 1 && push!(ret, Move(:up))
    y1 != T && push!(ret, Move(:down))
    return ret
end

function solveSlideBlockProblem(obj::SA, h::Function) where SA
    vinit = [7, 2, 4, 5, 0, 6, 8, 3, 1]
    problem = SlideBlockProblem(Grid{3}(vinit), obj, h)
    path = search(problem)
    path isa Symbol && return path
    ret=[]
    count = 0
    for iter in path
        push!(ret, iter.state)
    end
    return ret
end

GSD = GraphSearchDepth(Grid(3))
GSU = GraphSearchUniformCost(Grid(3))
GSBF = GraphSearchBestFirst(Grid(3))
GSAS = GraphSearchAStar(Grid(3))
RBFS = RecursiveBestFirstSearch(Grid(3))

solveSlideBlockProblemGSU() = solveSlideBlockProblem(GSU, dist)
solveSlideBlockProblemGSBF() = solveSlideBlockProblem(GSBF, dist)
solveSlideBlockProblemGSAS() = solveSlideBlockProblem(GSAS, dist)
solveSlideBlockProblemRBFS() = solveSlideBlockProblem(RBFS, dist)
solveSlideBlockProblemGSAS2() = solveSlideBlockProblem(GSAS, h2)
solveSlideBlockProblemGSBF2() = solveSlideBlockProblem(GSBF, h2)

end


text:
using CSFML.LibCSFML

mutable struct CellularAutomata
    matrix
end

function gol_number_of_adjacent_cells(i, j, matrix)
    # faster implementation :

    #   memory estimate:  0 bytes
    #   allocs estimate:  0
    #   --------------
    #   minimum time:     34.048 ns (0.00% GC)
    #   median time:      35.196 ns (0.00% GC)
    #   mean time:        35.230 ns (0.00% GC)
    #   maximum time:     101.239 ns (0.00% GC)
    #   --------------
    #   samples:          10000
    #   evals/sample:     993

    h, w = size(matrix)

    sum = 0
    if i == 1
        if j == 1
            #  __..
            # |0x
            # |xx

            sum = matrix[ 2, 1 ] + matrix[ 1, 2 ] + matrix[ 2, 2 ]
        elseif j == w
            #         ..__
            #           xo|
            #           xx|

            sum = matrix[ i+1, j ] + matrix[ i,j-1 ] + matrix[ i+1, j-1 ]
        else
            #   ..___..
            #     xox
            #     xxx

            sum = matrix[ i+1, j-1 ] + matrix[ i+1, j ] + matrix[ i+1, j+1 ] +
                        matrix[ i, j-1 ] + matrix[ i, j+1 ]
        end
    elseif i == h
        if j == 1
            # :
            # |xx
            # |ox
            #  --

            sum = matrix[ i-1, j ] + matrix[ i-1, j+1 ] + matrix[ i, j+1 ]
        elseif j == w
            #               :
            #             xx|
            #             xo|
            #             --

            sum = matrix[ i, j-1 ] + matrix[ i-1, j ] + matrix[ i-1, j-1 ]
        else
            #
            #      xxx
            #      xox
            #    ..---..

            sum = matrix[ i-1, j-1 ] + matrix[ i-1, j ] + matrix[ i-1, j+1 ] +
                matrix[ i, j-1 ] + matrix[ i, j+1 ]
        end
    else # 1 < i < h
        if j == 1
            #  :
            #  |xx
            #  |ox
            #  |xx
            #  :

            sum = matrix[ i-1, j+1 ] + matrix[ i, j+1 ] + matrix[ i+1, j+1] +
                 matrix[ i-1, j ] + matrix[ i+1, j ]
        elseif j == w
            #          :
            #        xx|
            #        xo|
            #        xx|
            #          :

            sum = matrix[ i-1, j-1 ] + matrix[ i, j-1 ] + matrix[ i+1, j-1 ] +
                matrix[ i-1, j ] + matrix[ i+1, j ]
        else
            #     xxx
            #     xox
            #     xxx

            sum = matrix[ i-1, j-1 ] + matrix[ i-1, j ] + matrix[ i-1, j+1 ] +
                  matrix[ i,   j-1 ]                    + matrix[ i,   j+1 ] +
                  matrix[ i+1, j-1 ] + matrix[ i+1, j ] + matrix[ i+1, j+1 ]
        end
    end
    sum
end

# GOL rules
function gol_compute_next_state(i, j, previous_state_matrix)
	adjacent_live_cell_count = gol_number_of_adjacent_cells(i, j, previous_state_matrix)
    previous_state = previous_state_matrix[i,j]
    next_state = 0
    if previous_state == 1
        if adjacent_live_cell_count < 2
            next_state = 0
        elseif adjacent_live_cell_count <= 3
            next_state = 1
        elseif adjacent_live_cell_count > 3
            next_state = 0
        end
    elseif previous_state == 0
        if adjacent_live_cell_count == 3
            next_state = 1
        else
            next_state = 0
        end
    end
    next_state
end

function update_ca!(ca::CellularAutomata)
    h, w = size(ca.matrix)
    next_matrix = similar(ca.matrix)
    for i=1:w, j=1:h
        next_matrix[i,j] = gol_compute_next_state(i, j, ca.matrix)
    end
    ca.matrix = next_matrix
end

function draw_ca(window,ca,cell_size)
    matrix_h, matrix_w = size(ca.matrix)
    square = sfRectangleShape_create()
    sfRectangleShape_setSize(square, sfVector2f(cell_size,cell_size) )
    sfRectangleShape_setFillColor(square,sfColor_fromRGBA(255,155,105,255) )
    for i in 0:matrix_h-1, j in 0:matrix_w-1
        value = ca.matrix[i + 1, j + 1]
        if value == 1
            sfRectangleShape_setPosition(square, sfVector2f(cell_size*i, cell_size*j) )
            sfRenderWindow_drawRectangleShape(window, square, C_NULL)
        end
    end
end

println("""=========================================================================================

                                ███▄ ▄███▓ ▄▄▄       ██▓ ███▄     █
                                ▓██▒▀█▀ ██▒▒████▄    ▓██▒ ██ ▀█   █
                                ▓██    ▓██░▒██  ▀█▄  ▒██▒▓██  ▀█ ██▒
                                ▒██    ▒██ ░██▄▄▄▄██ ░██░▓██▒  ▐▌██▒
                                ▒██▒   ░██▒ ▓█   ▓██▒░██░▒██░   ▓██░
                                ░ ▒░   ░  ░ ▒▒   ▓▒█░░▓  ░ ▒░   ▒ ▒
                                ░  ░      ░  ▒   ▒▒ ░ ▒ ░░ ░░   ░ ▒░
                                ░      ░     ░   ▒    ▒ ░   ░   ░ ░
                                    ░         ░  ░ ░           ░

===============================================================================================""")

Config = (
    window_size         = 500    ,   # pixels
    ca_matrix_size      = 50      ,   # matrix size
    target_frequency    = 10           # CA updates per seconds.
)

# init window
window_height = Config.window_size
window_width = Config.window_size
mode = sfVideoMode(window_height, window_width, 32)
window = sfRenderWindow_create(mode, "Game Of Life", sfResize | sfClose, C_NULL)
@assert window != C_NULL
sfWindow_setVerticalSyncEnabled(window,true)
sfWindow_setActive(window, sfTrue)

# init font
font = sfFont_createFromFile(joinpath("/julia/playgrounds/resources/NotoSans-Black.ttf"))
@assert font != C_NULL

# init ca
s = Config.ca_matrix_size
w, h = s, s
matrix = [ rand([0,1]) for i = 1:h, j = 1:w ]
ca = CellularAutomata(matrix)

# init matrix rendering engine
cell_size = Config.window_size / s

# init clock
clock = sfClock_create()
# Let's express time in seconds.
function get_time(clock ::Ptr{Nothing})
    sfTime_asSeconds(sfClock_getElapsedTime(clock))
end
function restart(clock ::Ptr{Nothing})
     sfTime_asSeconds(sfClock_restart(clock))
end

event_ref = Ref{sfEvent}()
target_frequency = Config.target_frequency
frame_start_timestamp = -Inf
time_frame_computing = 0

running = true
while (running)

    while Bool(sfRenderWindow_isOpen(window)) && running

        previous_frame_timestamp = frame_start_timestamp

#       # process events
        frame_start_timestamp = get_time(clock)
        frame_time = frame_start_timestamp - previous_frame_timestamp
        actual_frequency = 1 / frame_time

        while Bool(sfRenderWindow_pollEvent(window, event_ref))
            if event_ref[].type == sfEvtClosed
                sfRenderWindow_close(window)
                running = false
                println("Render window closed.")
            elseif event_ref[].type == sfEvtResized
                glViewport(9, 9, event_ref[].size.width, event_ref[].size.height)
                println("Viewport updated.")
            elseif event_ref[].type == sfEvtKeyPressed
                if event_ref[].key.code == sfKeyQ
                    running = false
                    println("Quiting application.")
                end
            end
        end

#       # update CA
        ca_start_timestamp = get_time(clock)
        time_event_processing = ca_start_timestamp - frame_start_timestamp
        #println("Event processing took $time_event_processing seconds.")

        update_ca!(ca)

#       # draw scene
        drawca_start_timestamp = get_time(clock)
        time_update_ca = drawca_start_timestamp - ca_start_timestamp

        sfRenderWindow_clear(window, sfColor_fromRGBA(24,20,18,255))

        # ----------------------------------------------------------------------------
        draw_ca(window,ca,cell_size)
        #draw_scene(window,ca,cell_size)
        # ----------------------------------------------------------------------------

#       # draw hud
        drawhud_start_timestamp = get_time(clock)
        time_drawca = drawhud_start_timestamp - drawca_start_timestamp

        sf_text = sfText_create()
        sfText_setPosition(sf_text, sfVector2f(4,window_height-20))
        sfText_setString(sf_text, "frame computing time : $time_frame_computing seconds.")
        sfText_setFont(sf_text, font)
        sfText_setCharacterSize(sf_text, 14)
        sfRenderWindow_drawText(window, sf_text, C_NULL)

        sf_text2 = sfText_create()
        sfText_setPosition(sf_text2, sfVector2f(4,window_height-40))
        sfText_setString(sf_text2, "fps : $actual_frequency")
        sfText_setFont(sf_text2, font)
        sfText_setCharacterSize(sf_text2, 14)
        sfRenderWindow_drawText(window, sf_text2, C_NULL)

#       # render
        render_start_timestamp = get_time(clock)
        sfRenderWindow_display(window)
#       # computing end
        frame_end_timestamp = get_time(clock)

        time_render = frame_end_timestamp - render_start_timestamp

        time_drawhud = render_start_timestamp - drawhud_start_timestamp
        time_frame_computing = frame_end_timestamp - frame_start_timestamp

        time_to_sleep = 1 / target_frequency - time_frame_computing

        time_sum = time_event_processing +
             time_drawca +
             time_drawhud +
             time_render

        time_report = """
        **
        event processing         : $time_event_processing
        draw Cellular Automata   : $time_drawca
        draw hud                 : $time_drawhud
        render                   : $time_render
        time sum                 : $time_sum
        time computing frame     : $time_frame_computing

        """
        println(time_report)
        if time_to_sleep > 0.005
            sleep(time_to_sleep-0.005)
        end
    end

end

#sfText_destroy(text)
sfFont_destroy(font)
#sfSprite_destroy(sprite)
#sfTexture_destroy(texture)
sfRenderWindow_destroy(window)

max_stars_repo_name: leroy-walton/playgrounds

max_issues_repo_name: leroy-walton/playgrounds

max_forks_repo_name: leroy-walton/playgrounds


text:
"""
Returns if a vector is inside another vector
"""
function in_vector(sub_vector, vector)
    for sv in sub_vector
        if !in(sv, vector)
            return false
        end
    end
    return true
end

"""
Gets array column index by a name
"""
function get_column_index(name, names)
    return findfirst(isequal(name), names)
end

"""
Filter rawdata by empty values
"""
function filter_raw_data_by_empty_values(datatype, depvar_data, expvars_data, panel_data=nothing, time_data=nothing)
    keep_rows = Array{Bool}(undef, size(depvar_data, 1))
    keep_rows .= true
    keep_rows .&= map(b->!b, ismissing.(depvar_data))

    for i = 1:size(expvars_data, 2)
        keep_rows .&= map(b->!b, ismissing.(expvars_data[:, i]))
    end

    depvar_data = convert(Array{datatype}, depvar_data[keep_rows, 1])
    expvars_data = convert(Array{datatype}, expvars_data[keep_rows, :])

    if panel_data != nothing
        panel_data = panel_data[keep_rows, 1]
    end

    if time_data != nothing
        time_data = time_data[keep_rows, 1]
    end

    return depvar_data, expvars_data, panel_data, time_data
end

"""
Filter data by empty values
"""
function filter_data_by_empty_values(data)
    depvar_data, expvars_data, panel_data, time_data = filter_raw_data_by_empty_values(
        data.datatype,
        data.depvar_data,
        data.expvars_data,
        data.panel_data,
        data.time_data
    )

    data.depvar_data = depvar_data
    data.expvars_data = expvars_data
    data.panel_data = panel_data
    data.time_data = time_data
    data.nobs = size(data.depvar_data, 1)
    data.removemissings = true

    return data
end

"""
Convert column by data content
"""
function convert_column(datatype, column)
    if column != nothing
        has_missings = false

        if size(column, 2) == 1
            has_missings |= findfirst(x -> ismissing(x), column) != nothing
        else
            for i in 1:size(column, 2)
                has_missings |= findfirst(x -> ismissing(x), column[:,i]) != nothing
            end
        end

        if has_missings
            return convert(Array{Union{Missing, datatype}}, column)
        else
            return convert(Array{datatype}, column)
        end
    end
    return nothing
end

"""
Convert rawdata by data content
"""
function convert_raw_data(datatype, depvar_data, expvars_data, panel_data=nothing, time_data=nothing)
    depvar_data = convert_column(datatype, depvar_data)
    expvars_data = convert_column(datatype, expvars_data)
    panel_data = convert_column(datatype == Float64 ? Int64 : Int32, panel_data)
    time_data = convert_column(datatype, time_data)

    return depvar_data, expvars_data, panel_data, time_data
end

"""
Convert data by data
"""
function convert_data(data)
    depvar_data, expvars_data, panel_data, time_data = convert_raw_data(
        data.datatype,
        data.depvar_data,
        data.expvars_data,
        data.panel_data,
        data.time_data
    )
    data.depvar_data = depvar_data
    data.expvars_data = expvars_data
    data.panel_data = panel_data
    data.time_data = time_data

    return data
end

"""
Copy ModelSelectionData
"""
function copy_data(data::ModelSelectionData)
    new_data = ModelSelectionData(
        copy(data.equation),
        data.depvar,
        copy(data.expvars),
        data.panel,
        data.time,
        copy(data.depvar_data),
        copy(data.expvars_data),
        if (data.panel_data != nothing) copy(data.panel_data) else data.panel_data end,
        if (data.time_data != nothing) copy(data.time_data) else data.time_data end,
        data.intercept,
        data.datatype,
        data.removemissings,
        data.nobs
    )

    new_data.extras = data.extras
    new_data.options = copy(data.options)
    new_data.previous_data = copy(data.previous_data)
    new_data.results = copy(data.results)

    return new_data
end

"""
Copy ModelSelectionData to another data
"""
function copy_data!(from_data::ModelSelectionData, to_data::ModelSelectionData)
    to_data.equation = from_data.equation
    to_data.depvar = from_data.depvar
    to_data.expvars = from_data.expvars
    to_data.panel = from_data.panel
    to_data.time = from_data.time
    to_data.depvar_data = from_data.depvar_data
    to_data.expvars_data = from_data.expvars_data
    to_data.panel_data = from_data.panel_data
    to_data.time_data = from_data.time_data
    to_data.intercept = from_data.intercept
    to_data.datatype = from_data.datatype
    to_data.removemissings = from_data.removemissings
    to_data.nobs = from_data.nobs
    to_data.extras = from_data.extras
    to_data.options = from_data.options
    to_data.previous_data = from_data.previous_data
    to_data.results = from_data.results

    return to_data
end

"""
Generate extra key
"""
function generate_extra_key(extra_key, extras)
    if !(extra_key in keys(extras))
        return extra_key
    else
        posfix = 2
        while Symbol(string(extra_key, "_", posfix)) in keys(extras)
            posfix = posfix + 1
        end
        return Symbol(string(extra_key, "_", posfix))
    end
end

"""
Returns if feature extraction module was selected
"""
function featureextraction_enabled(fe_sqr, fe_log, fe_inv, fe_lag, interaction)
    return fe_sqr != nothing || fe_log != nothing || fe_inv != nothing || fe_lag != nothing || interaction != nothing
end

"""
Returns if preliminary selection was selected
"""
function preliminaryselection_enabled(preliminaryselection)
    return preliminaryselection != nothing
end

"""
Validates if preliminary selecttion method exists
"""
function validate_preliminaryselection(preliminaryselection)
    return preliminaryselection in VALID_PRELIMINARYSELECTION
end

"""
Add result
"""
function addresult!(data, result)
    push!(data.results, result)

    return data
end

"""
Creates an array with keys and array positions
"""
function create_datanames_index(datanames)
    header = Dict{Symbol,Int64}()
    for (index, name) in enumerate(datanames)
        header[name] = index
    end
    return header
end

"""
Returns selected appropiate covariates for each iteration
"""
function get_selected_variables(order, datanames, intercept; fixedvariables=nothing, num_jobs=nothing, num_job=nothing, iteration_num=nothing)
    cols = zeros(Int64, 0)
    binary = string(order, base = 2)
    k = 1

    for order = 1:length(binary)
        if binary[length(binary) - order + 1] == '1'
            push!(cols, k)
        end
        k = k + 1
    end
    if intercept
        push!(cols, ModelSelection.get_column_index(:_cons, datanames))
    end
    return cols
end

"""
Add intercept
"""
function add_intercept!(data)
    data.expvars_data = hcat(data.expvars_data, ones(data.nobs))
    push!(data.expvars, :_cons)
end

"""
Remove intercept
"""
function remove_intercept!(data)
    cons_index = get_column_index(:_cons, data.expvars)
    data.expvars_data = hcat(data.expvars_data[:, 1:cons_index-1], data.expvars_data[:, cons_index+1:end])
    data.expvars = vcat(data.expvars[1:cons_index-1], data.expvars[cons_index+1:end])
end


text:
"""
    proc(com)

Extract the content of a single-brace lx command. For instance `\\com{foo}`
would be extracted to `foo`.
"""
proc(com) = F.content(com.braces[1])

"""
    html(s)

Mark a string as HTML to be included in Franklin-markdown. Line spacing is
to reduce issues with `<p>`.
"""
html(s) = "\n~~~$s~~~\n"

"""
    lxargs(s)

Extract function-style arguments.
Expect (arg1, arg2, ..., name1=val1, name2=val2, ...)

Example:

    julia> a = ":section, 1, 3, title=\"hello\", name=\"bar\""
    julia> lxargs(a)
    (Any[:section, 1, 3], Any[:title => "hello", :name => "bar"])
"""
function lxargs(s, fname="")
    isempty(s) && return [], []
    s = "(" * s * ",)"
    args = nothing
    try
        args = Meta.parse(s).args
    catch
        error("A \\$fname{...} had improper specs:\n$s; verify.")
    end
    i = findfirst(e -> isa(e, Expr), args)
    if isnothing(i)
        cand_args = args
        cand_kwargs = []
    else
        cand_args = args[1:i-1]
        cand_kwargs = args[i:end]
    end
    proc_args   = []
    proc_kwargs = []
    for arg in cand_args
        if arg isa QuoteNode
            push!(proc_args, arg.value)
        else
            push!(proc_args, arg)
        end
    end
    for kwarg in cand_kwargs
        kwarg isa Expr || error("In \\$fname{...}, expected arguments " *
                                "followed by keyword arguments but got: " *
                                "$s; verify.")
        push!(proc_kwargs, kwarg.args[1] => kwarg.args[2])
    end
    return proc_args, proc_kwargs
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 3cd07b90-9e21-4a74-8d99-a6375beb7d1e
begin
	using Revise
	using Pkg
	Pkg.activate("../")

	using BayesianModels
	using DelimitedFiles
	using LinearAlgebra
	using Plots
	using Zygote

	include("./plotting.jl")
end

# ╔═╡ 6958bde2-c50f-4afd-8d08-7dfbfd97169a
md"""
# Learning and Modeling the Languages Phonetics

*[Lucas Ondel](https://lucasondel.github.io/), 8 March 2022*
"""

# ╔═╡ b83311ab-ca4f-432e-8627-5207517a81db
md"""## Lingua Libre Formants data

We use the [Lingua Libre Formants](https://github.com/lucasondel/LinguaLibreFormants) data for our analysis.

### Getting the data

We download the data and store it in a local file.
"""

# ╔═╡ d1588a1c-babe-4223-bc6d-f353314ed2c0
localfile = "formants.tsv"

# ╔═╡ 262e0395-814f-4315-b25e-21c18a7a9a5f
begin
	# Database URL
	local dataurl = "https://raw.githubusercontent.com/lucasondel/LinguaLibreFormants/main/data/formants.tsv"

	# Download the data if it is not done already
	if ! isfile(localfile)
		run(`wget $dataurl`)
	end
end

# ╔═╡ 4bbdd3d2-eb62-4889-9c8a-20141c0b140a
md"""
### Data preparation

In a nutshell:
1. we extract the list of languages, speakers and phones
2. organized the data as a $3 \times N$ matrix where 3 is for 3 first formants and $N$ is the number of data points
3. for each data point we extract the annotation (language, speaker, phone)
"""

# ╔═╡ ee1d1c45-70f5-4fe7-a8e2-9134d45a2b95
begin
	local rawdata = readdlm(localfile)

	speakers = Set()
	languages = Set()
	phones = Set()
	labels = []
	data = []
	for (lang, spk, ipaphone, f1, f2, f3, duration) in eachrow(rawdata)
		# Remove the duration symbol
		ipaphone = replace(ipaphone, "ː" => "")

		lang, spk, ipaphone = Symbol(lang), Symbol(spk), Symbol(ipaphone)
		push!(data, vcat(f1, f2, f3))

		phone = (lang, ipaphone)
		spk = (lang, spk)

		push!(labels, (lang, spk, phone))
		push!(speakers, spk)
		push!(languages, lang)
		push!(phones, phone)
	end
	data = hcat(data...)
	labels = vcat(labels...)

	# Mean-variance normalization of the data
	local x̄ = sum(data, dims=2) / size(data, 2)
	local x̄² = sum(data .^ 2, dims=2) / size(data, 2)
	local var = x̄² - (x̄ .^ 2)
	data = (data .- x̄) ./ sqrt.(var)

	# Transform the languages, speakers and phone sets to ordered collections
	languages = collect(languages)
	speakers = collect(speakers)
	phones = collect(phones)

	# Build the element <-> id mappings
	lang2id = Dict(lang => i for (i, lang) in enumerate(languages))
	id2lang = Dict(i => lang for (i, lang) in enumerate(languages))
	spk2id = Dict(spk => i for (i, spk) in enumerate(speakers))
	id2spk = Dict(i => spk for (i, spk) in enumerate(speakers))
	phone2id = Dict(phone => i for (i, phone) in enumerate(phones))
	id2phone = Dict(i => phone for (i, phone) in enumerate(phones))

	# Transform the labels to use the id rather than symbols
	labels = [(lang2id[lang], spk2id[spk], phone2id[phone])
			  for (lang, spk, phone) in labels]

	# Group the data for each triplet (lang, speaker, phone)
	groups = Dict()
	for (label, sample) in zip(labels, eachcol(data))
		list = get(groups, label, [])
		push!(list, sample)
		groups[label] = list
	end
	groups = tuple((k => hcat(groups[k]...) for k in keys(groups))...)
end

# ╔═╡ 474abeb1-82f8-4a74-a03c-be1cf55df840
md"Some important constants:"

# ╔═╡ 31e56fed-e4f0-4c71-abb2-5c2eff5a9a47
D = size(data, 1)

# ╔═╡ 61cc1ab2-a615-4a9b-91f9-bc5909e29297
N = size(data, 2)

# ╔═╡ e08cc384-49dc-4585-a2ab-fbec9fcfe9d2
L = length(languages)

# ╔═╡ b8e7f92f-0065-46fe-ac6b-5adf29b3379f
S = length(speakers)

# ╔═╡ ef71661f-338e-4d57-86ed-fc2c308a744f
P = length(phones)

# ╔═╡ c01de0cd-43f8-49cc-a8c9-ffbbc182d9e9
begin
	md"""
	### Data statistics

	| \# languages | \# speakers  | \# phones |
	|:------------:|:------------:|:---------:|
	| $L           | $S           | $P        |

	"""
end

# ╔═╡ 8d416ab7-3ae7-4984-a7e6-1f56692d1b7b
begin
	mdstr = "| language | # speakers | phones |\n"
	mdstr *= "|:---------|:-----------:|:-------|\n"
	for lang in languages
		mdstr *= "| $lang | "

		ns = 0
		for (l, spk) in speakers
			if l == lang ns += 1end
		end
		mdstr *= "$ns | "

		for (l, phone) in phones
			if l == lang
				mdstr *= "$phone "
			end
		end
		mdstr *= "|\n"
	end
	Markdown.parse(mdstr)
end

# ╔═╡ 8b535f17-7a55-42da-ba7a-c161bec9fd28
md"""
## Model

### Embeddings
1.  $\boldsymbol{\lambda} \sim p(\boldsymbol{\lambda}) \triangleq \mathcal{N}( \mathbf{0}, \mathbf{I})$, $\boldsymbol{\lambda}$ → language embedding of dimension $Q_\lambda$
2.  $\boldsymbol{\sigma} \sim p(\boldsymbol{\sigma}) \triangleq \mathcal{N}( \mathbf{0}, \mathbf{I})$, $\boldsymbol{\sigma}$ → speaker embedding of dimension $Q_\sigma$
3.  $\boldsymbol{\pi} \sim p(\boldsymbol{\pi}) \triangleq \mathcal{N}( \mathbf{0}, \mathbf{I})$ $\boldsymbol{\pi}$ → phone embedding of dimension $Q_\pi$
"""

# ╔═╡ 5175712e-7f27-459e-8e6a-9bed325d6b4e
md"""
### Subspace construction
1.  $\mathbf{L}_i \sim p(\mathbf{L}) \triangleq \mathcal{N}(\mathbf{0}, \mathbf{I})$, $\{\mathbf{L}_0, \dots, \mathbf{L}_{Q_\lambda} \}$ → Language hyper-subspace bases of dimension $(Q_\pi + 1) \times K$
2.  $\mathbf{S}_j \sim p(\mathbf{S}) \triangleq \mathcal{N}(\mathbf{0}, \mathbf{I})$, $\{\mathbf{S}_0, \dots, \mathbf{S}_{Q_\sigma} \}$ → Speaker hyper-subspace bases $(Q_\pi + 1) \times K$
3.  $\mathbf{P} = \sum_{i=0}^{Q_\lambda} \lambda_i \mathbf{L}_i + \sum_{j=0}^{Q_\sigma} \sigma_j \mathbf{S}_j$ where $\lambda_0 = \sigma_0 = 1$, → Phonetic subspace bases of dimension $(Q_\pi + 1) \times K$
"""

# ╔═╡ 763b5f95-db02-4ede-b632-848a646232d9
md"""
### Phone model
1.  $\boldsymbol{\psi} = \mathbf{P}^\top \begin{bmatrix} \boldsymbol{\pi} \\ 1 \end{bmatrix}$, $\boldsymbol{\psi} \in \mathbb{R}^K$
2.  $\boldsymbol{\eta} = f(\boldsymbol{\psi})$
3.  $\mathbf{x}_n \sim p(\mathbf{x} | \boldsymbol{\eta}) \triangleq \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})$, where $\boldsymbol{\eta}$ are the natural parameters of the multivariate Normal distribution
"""

# ╔═╡ 841a3f0e-8975-49d2-959e-afe832ed41e4
md"""
### Parameter function mapping

The function $f$ transforms the parameters $\boldsymbol{\psi}$ from the real space into the natural parameter space. It is defined as follows:

```math
	\begin{align}
	\boldsymbol{\psi} &= \begin{bmatrix}
		 \boldsymbol{\psi}^{(\boldsymbol{\Lambda\mu})} \\
		 \boldsymbol{\psi}^{(\mathbf{L})}_1 \\
		 \boldsymbol{\psi}^{(\mathbf{L})}_2
	\end{bmatrix} \\
	\boldsymbol{\Lambda\mu} &= \boldsymbol{\psi}^{(\boldsymbol{\Lambda\mu})} \\
	\text{diag}(\mathbf{L}) &= \exp \{ \boldsymbol{\psi}^{(\mathbf{L})}_1 \} \\
	\text{vech}^{\dagger}(\mathbf{L}) &= \boldsymbol{\psi}^{(\mathbf{L})}_2 \\
	\boldsymbol{\Lambda} &= \mathbf{L} \mathbf{L}^\top
	\end{align}
```
where $\text{vech}^{\dagger}$ is the half-vectorization of a matrix **not including the diagonal**
"""

# ╔═╡ b0761304-49df-4626-a677-488295dc786d
function f(ψ::AbstractVector)
		L = length(ψ)
		D = Int((-(3/2) + sqrt((9/4) + 2L)) )

		Λμ = ψ[1:D]
		A = diagm(exp.(1e-3 .+ ψ[D+1:2*D]))
		L = A + CompressedLowerTriangular(D, 1, ψ[2*D+1:end])
		L⁻¹ = inv(L)

		Λ = L * L'
		μᵀΛμ = ((L⁻¹' * L⁻¹) * Λμ)' * Λμ
		logdetΛ = 2*sum(ψ[D+1:2*D])

		Λμ, μᵀΛμ, Λ, logdetΛ
	end

# ╔═╡ 4154a267-326d-495c-aa60-156560aa5f72
Qλ = 2

# ╔═╡ 4c69229d-a272-49eb-9dc1-484c53bad338
pλ = Normal(Qλ)

# ╔═╡ eebb75b8-8ece-4e26-972b-ed962ce1dd42
Qσ = 2

# ╔═╡ 3a303018-4560-4553-96f6-1600b0e1d6c6
pσ = Normal(Qσ)

# ╔═╡ b629f310-dbf7-4b6e-9c3d-199a31a1bfb5
Qπ = 2

# ╔═╡ 51d10c40-5dce-47e2-aeba-958beadda7b8
pπ = Normal(Qπ)

# ╔═╡ c352df19-b653-43f3-b490-046447cea85a
phonemodel = NormalModel((pλ=pλ, pσ=pσ, pπ=pπ))

# ╔═╡ c47dcaf4-049a-45e1-853f-d589283b2038
# The dimension of the parameter space is determined
# by the dimension of the observed space
K = Int(2D + (D^2 - D) / 2)

# ╔═╡ add6cb90-47bd-40e8-bbbb-34938714a79b
pL = JointNormalFixedCov(zeros(Qπ + 1, K))

# ╔═╡ 5415acd0-5e5f-40db-97eb-91b659a1ac5b
pS = JointNormalFixedCov(zeros(Qπ + 1, K))

# ╔═╡ b69c63d7-3f45-45d2-8278-05bc2fb5f3f1
md"""
### Hyper-parameters

| name          | type            | value|  description                  |
|:--------------|:----------------|:----:|:------------------------------|
| ``Q_\lambda`` | hyper-parameter | $Qλ  | language embedding dimension  |
| ``Q_\sigma``  | hyper-parameter | $Qσ  | speaker embedding dimension   |
| ``Q_\pi``     | hyper-parameter | $Qπ  | phone embedding dimension     |
| ``D``         | constant        | $D   | observed space dimension      |
| ``K``         | constant        | $K   | phone model natural parameter dimension|
| ``L`` 		| constant 		  | $L   | number of languages in the data |
| ``S``         | constant        | $S   | number of speakers in  the data |
| ``P``         | constant        | $P   | number of phones in the data  |
"""

# ╔═╡ 946af84f-a5f4-4780-a8c5-13d4acbb9a80
md"""
### Model parameters summary

```math
\begin{align}

\boldsymbol{\theta} &= \Big\{ \{\mathbf{L}_i\}, \{\mathbf{S}_j \}, \{\boldsymbol{\lambda}_l \}, \{\boldsymbol{\sigma}_s \}, \{ \boldsymbol{\pi}_p \} \Big\} \\
\{\mathbf{L}_i \}, \; &\forall i \in \{1, \dots, Q_\lambda \} \\
\{\mathbf{S}_j \}, \;  &\forall j \in \{1, \dots, Q_\sigma \} \\
\{\mathbf{\boldsymbol{\lambda}}_l \}, \; &\forall l \in \{1, \dots, L\} \\
\{\mathbf{\boldsymbol{\sigma}}_s \}, \; &\forall s \in \{1, \dots, S\} \\
\{\mathbf{\boldsymbol{\pi}_s} \}, \; &\forall p \in \{1, \dots, P\}
\end{align}
```
"""

# ╔═╡ f2f5e843-3505-421e-9c9c-71a564387bb2
# Complete model: hierarchical generalized subspace model.
hgsm = GSMPhonetic(f, pL, pS, phonemodel)

# ╔═╡ 1dd9bfa7-b7c9-4075-9c6c-2f33871d5e0a
md"""
### Variational posterior

```math
\begin{align}
	q(\boldsymbol{\theta}) &= \prod_{m}q(\mathbf{L}_m) \prod_{n} q(\mathbf{S}_n) \prod_{l} q(\boldsymbol{\lambda}_l) \prod_{s} q(\boldsymbol{\sigma}_s) \prod_{p} q(\boldsymbol{\pi}_p) \\
	q(\mathbf{L}_m) &= \delta(\mathbf{L}_m - \boldsymbol{\mu}^{L}_m) \\
	q(\mathbf{S}_n) &=  \delta(\mathbf{S}_n - \boldsymbol{\mu}^{S}_n) \\
	q(\boldsymbol{\lambda}_l) &= \mathcal{N}(\boldsymbol{\mu}^{\lambda}_l, \boldsymbol{\Sigma}^{\lambda}_l) \\
	q(\boldsymbol{\sigma}_s) &= \mathcal{N}(\boldsymbol{\mu}^{\sigma}_s, \boldsymbol{\Sigma}^{\sigma}_s) \\
		q(\boldsymbol{\pi}_p) &= \mathcal{N}(\boldsymbol{\mu}^{\pi}_p, \boldsymbol{\Sigma}^{\pi}_p)
\end{align}
```
"""

# ╔═╡ 6f1225eb-8df1-4991-8292-c0b5f2c45c10
q₀L = [JointNormalFixedCov(randn(Qπ + 1, K)) for i = 1:(Qλ+1)]

# ╔═╡ 349ef71a-55f8-4c9a-8216-d6fbce8d42e2
q₀S = [JointNormalFixedCov(randn(Qπ + 1, K)) for j = 1:(Qσ+1)]

# ╔═╡ b2875ebb-5b3b-4dad-b1e6-f5b38b231b03
L

# ╔═╡ 27fa667c-4a76-4fed-9afc-7482dc1339a0
q₀λ = [Normal(Qλ) for lang in languages]

# ╔═╡ 3a5cadd7-d15b-4da3-8a0b-9e7678514ff0
q₀σ = [Normal(Qσ) for spk in speakers]

# ╔═╡ e343d7f3-9a50-4952-b429-40b982b2b49f
q₀π = [Normal(Qπ) for phone in phones]

# ╔═╡ 42ccdf0d-8a5f-404d-8965-cccc7f85d164
q₀θ = (q₀L, q₀S, q₀λ, q₀σ, q₀π)

# ╔═╡ 85e89c27-0371-42a9-8b3c-4fc677d1f7a7
md"""
### Training

"""

# ╔═╡ cde6061d-57af-4aae-a964-4cafabfdd5f1
begin
	qθ = q₀θ

	loss = []
	for t in 1:300000
		μ = getstats(hgsm, qθ)
		(l, (∇μ,)) = withgradient(μ -> elbo(hgsm, groups, qθ, μ) / N, μ)
		qθ = newposterior(hgsm, qθ, ∇μ; lrate=1e-2)

		push!(loss, l)
	end
	loss
end

# ╔═╡ 7175c939-bc9d-4b60-b145-b09cac018b7f
plot(loss[10000:1:end], legend=false)

# ╔═╡ 4274c7f4-38d5-4bf0-87ae-417d715d9b3b
begin
	local qλ = qθ[3]
	local p = plot(legend=false)
	for (i, q) in enumerate(qλ)
		x, diagxxᵀ, hxxᵀ = unpack(q, μ(q))
		m = x
		Σ = diagm(diagxxᵀ) + CompressedSymmetric(2, 1, hxxᵀ) - x * x'
		annotate!(m[1], m[2], id2lang[i])
		plotnormal!(m, Σ, σ=2, fillalpha=0.4, linealpha=0, label=false)
	end
	p
end

# ╔═╡ ef349252-cf35-42bd-b562-7ec28bef7cdd
begin
	local qσ = qθ[4]
	local p = plot(legend=false)
	for (i, q) in enumerate(qσ)
		x, diagxxᵀ, hxxᵀ = unpack(q, μ(q))
		m = x
		Σ = diagm(diagxxᵀ) + CompressedSymmetric(2, 1, hxxᵀ) - x * x'
		plotnormal!(m, Σ, σ=2, fillalpha=0.4, linealpha=0, label=false)
	end
	p
end

# ╔═╡ b9bba539-ba31-420e-a718-a952db80d4c2
begin
	local qπ = qθ[5]
	local p = plot(legend=false)
	for (i, q) in enumerate(qπ)
		if phones[i][1] != :ita continue end
		x, diagxxᵀ, hxxᵀ = unpack(q, μ(q))
		m = x
		Σ = diagm(diagxxᵀ) + CompressedSymmetric(2, 1, hxxᵀ) - x * x'
		annotate!(m[1], m[2], id2phone[i][2])
		plotnormal!(m, Σ, σ=2, fillcolor=:blue, fillalpha=0.2, linealpha=0, label=false)
	end
	p
end

# ╔═╡ ab7167f2-41c4-4944-b769-d7d27c9bda12
id2phone

# ╔═╡ 8b939634-e3f2-4a29-91de-4ed56787bc78
function minimum_divergence(hgsm, qθ)
	qL, qS, qλ, qσ, qπ = qθ

	μπ = map(μ, qπ)
	stdπ = stdparams.(qπ, μπ)

	μπ₀ = sum(μπ) / length(qπ)
	m₀, Σ₀ = stdparams(hgsm.model.pθ.pπ, μπ₀)
	L₀ = cholesky(Σ₀).L
	L₀⁻¹ = L₀

	[Normal(L₀⁻¹ * (m - m₀), Symmetric(L₀⁻¹ * Σ * L₀⁻¹')) for (m, Σ) in stdπ]
end

# ╔═╡ b44032e6-bdff-4f43-a225-70cddc2d89cd
minimum_divergence(hgsm, qθ)

# ╔═╡ d045db5e-e646-4f04-9df3-be877cd505f9
m, Σ = stdparams(qλ[1], unpack(qλ[1], sum(μ, qλ) / length(qλ)))

# ╔═╡ b518b13e-f470-4c4b-8c2e-e3b0803beb62
qθ[1][1].M

# ╔═╡ ec85819b-1b0f-4f0e-9602-46c98402bfdf
Lm⁻¹ = inv(Lm)

# ╔═╡ Cell order:
# ╟─6958bde2-c50f-4afd-8d08-7dfbfd97169a
# ╠═3cd07b90-9e21-4a74-8d99-a6375beb7d1e
# ╟─b83311ab-ca4f-432e-8627-5207517a81db
# ╠═d1588a1c-babe-4223-bc6d-f353314ed2c0
# ╠═262e0395-814f-4315-b25e-21c18a7a9a5f
# ╟─4bbdd3d2-eb62-4889-9c8a-20141c0b140a
# ╠═ee1d1c45-70f5-4fe7-a8e2-9134d45a2b95
# ╟─474abeb1-82f8-4a74-a03c-be1cf55df840
# ╠═31e56fed-e4f0-4c71-abb2-5c2eff5a9a47
# ╠═61cc1ab2-a615-4a9b-91f9-bc5909e29297
# ╠═e08cc384-49dc-4585-a2ab-fbec9fcfe9d2
# ╠═b8e7f92f-0065-46fe-ac6b-5adf29b3379f
# ╠═ef71661f-338e-4d57-86ed-fc2c308a744f
# ╟─c01de0cd-43f8-49cc-a8c9-ffbbc182d9e9
# ╟─8d416ab7-3ae7-4984-a7e6-1f56692d1b7b
# ╟─8b535f17-7a55-42da-ba7a-c161bec9fd28
# ╠═4c69229d-a272-49eb-9dc1-484c53bad338
# ╠═3a303018-4560-4553-96f6-1600b0e1d6c6
# ╠═51d10c40-5dce-47e2-aeba-958beadda7b8
# ╟─5175712e-7f27-459e-8e6a-9bed325d6b4e
# ╠═add6cb90-47bd-40e8-bbbb-34938714a79b
# ╠═5415acd0-5e5f-40db-97eb-91b659a1ac5b
# ╟─763b5f95-db02-4ede-b632-848a646232d9
# ╠═c352df19-b653-43f3-b490-046447cea85a
# ╟─841a3f0e-8975-49d2-959e-afe832ed41e4
# ╠═b0761304-49df-4626-a677-488295dc786d
# ╟─b69c63d7-3f45-45d2-8278-05bc2fb5f3f1
# ╠═4154a267-326d-495c-aa60-156560aa5f72
# ╠═eebb75b8-8ece-4e26-972b-ed962ce1dd42
# ╠═b629f310-dbf7-4b6e-9c3d-199a31a1bfb5
# ╠═c47dcaf4-049a-45e1-853f-d589283b2038
# ╟─946af84f-a5f4-4780-a8c5-13d4acbb9a80
# ╠═f2f5e843-3505-421e-9c9c-71a564387bb2
# ╟─1dd9bfa7-b7c9-4075-9c6c-2f33871d5e0a
# ╠═6f1225eb-8df1-4991-8292-c0b5f2c45c10
# ╠═349ef71a-55f8-4c9a-8216-d6fbce8d42e2
# ╠═b2875ebb-5b3b-4dad-b1e6-f5b38b231b03
# ╠═27fa667c-4a76-4fed-9afc-7482dc1339a0
# ╠═3a5cadd7-d15b-4da3-8a0b-9e7678514ff0
# ╠═e343d7f3-9a50-4952-b429-40b982b2b49f
# ╠═42ccdf0d-8a5f-404d-8965-cccc7f85d164
# ╟─85e89c27-0371-42a9-8b3c-4fc677d1f7a7
# ╠═cde6061d-57af-4aae-a964-4cafabfdd5f1
# ╠═7175c939-bc9d-4b60-b145-b09cac018b7f
# ╠═4274c7f4-38d5-4bf0-87ae-417d715d9b3b
# ╠═ef349252-cf35-42bd-b562-7ec28bef7cdd
# ╠═b9bba539-ba31-420e-a718-a952db80d4c2
# ╠═ab7167f2-41c4-4944-b769-d7d27c9bda12
# ╠═8b939634-e3f2-4a29-91de-4ed56787bc78
# ╠═b44032e6-bdff-4f43-a225-70cddc2d89cd
# ╠═d045db5e-e646-4f04-9df3-be877cd505f9
# ╠═b518b13e-f470-4c4b-8c2e-e3b0803beb62
# ╠═ec85819b-1b0f-4f0e-9602-46c98402bfdf

max_stars_repo_path: examples/Learning and Modeling Languages Phonetics.jl

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
blockrowsupport(_, A, k) = blockaxes(A,2)
""""
    blockrowsupport(A, k)

gives an iterator containing the possible non-zero blocks in the k-th block-row of A.
"""
blockrowsupport(A, k) = blockrowsupport(MemoryLayout(A), A, k)
blockrowsupport(A) = blockrowsupport(A, blockaxes(A,1))

blockcolsupport(_, A, j) = blockaxes(A,1)

""""
    blockcolsupport(A, j)

gives an iterator containing the possible non-zero blocks in the j-th block-column of A.
"""
blockcolsupport(A, j) = blockcolsupport(MemoryLayout(A), A, j)
blockcolsupport(A) = blockcolsupport(A, blockaxes(A,2))

blockcolstart(A...) = first(blockcolsupport(A...))
blockcolstop(A...) = last(blockcolsupport(A...))
blockrowstart(A...) = first(blockrowsupport(A...))
blockrowstop(A...) = last(blockrowsupport(A...))

for Func in (:blockcolstart, :blockcolstop, :blockrowstart, :blockrowstop)
    @eval $Func(A, i::Block{1}) = $Func(A, Int(i))
end

abstract type AbstractBlockLayout <: MemoryLayout end
struct BlockLayout{LAY} <: AbstractBlockLayout end

similar(M::MulAdd{<:AbstractBlockLayout,<:AbstractBlockLayout}, ::Type{T}, axes) where {T,N} =
    similar(BlockArray{T}, axes)

MemoryLayout(::Type{<:PseudoBlockArray{T,N,R}}) where {T,N,R} = MemoryLayout(R)
MemoryLayout(::Type{<:BlockArray{T,N,R}}) where {T,N,R} = BlockLayout{typeof(MemoryLayout(R))}()

sublayout(::BlockLayout{LAY}, ::Type{NTuple{N,BlockSlice1}}) where {LAY,N} = LAY()
sublayout(BL::BlockLayout, ::Type{<:NTuple{N,BlockSlice}}) where N = BL

conjlayout(::Type{T}, ::BlockLayout{LAY}) where {T<:Complex,LAY} = BlockLayout{typeof(conjlayout(T,LAY()))}()
conjlayout(::Type{T}, ::BlockLayout{LAY}) where {T<:Real,LAY} = BlockLayout{LAY}()

transposelayout(::BlockLayout{LAY}) where LAY = BlockLayout{typeof(transposelayout(LAY()))}()

#############
# BLAS overrides
#############

function materialize!(M::MatMulVecAdd{<:AbstractBlockLayout,<:AbstractStridedLayout,<:AbstractStridedLayout})
    α, A, x_in, β, y_in = M.α, M.A, M.B, M.β, M.C
    if length(x_in) != size(A,2) || length(y_in) != size(A,1)
        throw(DimensionMismatch())
    end

    # impose block structure
    y = PseudoBlockArray(y_in, (axes(A,1),))
    x = PseudoBlockArray(x_in, (axes(A,2),))

    _fill_lmul!(β, y)

    for J = blockaxes(A,2)
        for K = blockcolsupport(A,J)
            muladd!(α, view(A,K,J), view(x,J), one(α), view(y,K))
        end
    end
    y_in
end

function _block_muladd!(α, A, X, β, Y)
    _fill_lmul!(β, Y)
    for J = blockaxes(X,2), N = blockcolsupport(X,J), K = blockcolsupport(A,N)
        muladd!(α, view(A,K,N), view(X,N,J), one(α), view(Y,K,J))
    end
    Y
end

mul_blockscompatible(A, B, C) = blockisequal(axes(A,2), axes(B,1)) &&
    blockisequal(axes(A,1), axes(C,1)) &&
    blockisequal(axes(B,2), axes(C,2))

function materialize!(M::MatMulMatAdd{<:AbstractBlockLayout,<:AbstractBlockLayout,<:AbstractBlockLayout})
    α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
    if mul_blockscompatible(A,B,C)
        _block_muladd!(α, A, B, β, C)
    else # use default
        materialize!(MulAdd{UnknownLayout,UnknownLayout,UnknownLayout}(α, A, B, β, C))
    end
end

function materialize!(M::MatMulMatAdd{<:AbstractBlockLayout,<:AbstractBlockLayout,<:AbstractColumnMajor})
    α, A, X, β, Y_in = M.α, M.A, M.B, M.β, M.C
    Y = PseudoBlockArray(Y_in, (axes(A,1), axes(Y_in,2)))
    _block_muladd!(α, A, X, β, Y)
    Y_in
end

function materialize!(M::MatMulMatAdd{<:AbstractBlockLayout,<:AbstractColumnMajor,<:AbstractColumnMajor})
    α, A, X_in, β, Y_in = M.α, M.A, M.B, M.β, M.C
    _fill_lmul!(β, Y_in)
    X = PseudoBlockArray(X_in, (axes(A,2), axes(X_in,2)))
    Y = PseudoBlockArray(Y_in, (axes(A,1), axes(Y_in,2)))
    _block_muladd!(α, A, X, β, Y)
    Y_in
end

function materialize!(M::MatMulMatAdd{<:AbstractColumnMajor,<:AbstractBlockLayout,<:AbstractColumnMajor})
    α, A_in, X, β, Y_in = M.α, M.A, M.B, M.β, M.C
    _fill_lmul!(β, Y_in)
    A = PseudoBlockArray(A_in, (axes(A_in,1),axes(X,1)))
    Y = PseudoBlockArray(Y_in, (axes(Y_in,1),axes(X,2)))
    _block_muladd!(α, A, X, β, Y)
    Y_in
end

####
# Triangular
####

@inline hasmatchingblocks(A) = blockisequal(axes(A)...)

triangularlayout(::Type{Tri}, ::ML) where {Tri,ML<:AbstractBlockLayout} = Tri{ML}()

_triangular_matrix(::Val{'U'}, ::Val{'N'}, A) = UpperTriangular(A)
_triangular_matrix(::Val{'L'}, ::Val{'N'}, A) = LowerTriangular(A)
_triangular_matrix(::Val{'U'}, ::Val{'U'}, A) = UnitUpperTriangular(A)
_triangular_matrix(::Val{'L'}, ::Val{'U'}, A) = UnitLowerTriangular(A)

function _matchingblocks_triangular_mul!(::Val{'U'}, UNIT, A::AbstractMatrix{T}, dest) where T
    # impose block structure
    b = PseudoBlockArray(dest, (axes(A,1),))

    for K = blockaxes(A,1)
        b_2 = view(b, K)
        Ũ = _triangular_matrix(Val('U'), UNIT, view(A, K,K))
        materialize!(Lmul(Ũ, b_2))
        JR = (K+1):last(blockrowsupport(A,K))
        if !isempty(JR)
            muladd!(one(T), view(A, K, JR), view(b,JR), one(T), b_2)
        end
    end
    dest
end

function _matchingblocks_triangular_mul!(::Val{'L'}, UNIT, A::AbstractMatrix{T}, dest) where T
    # impose block structure
    b = PseudoBlockArray(dest, (axes(A,1),))

    N = blocksize(A,1)

    for K = N:-1:1
        b_2 = view(b, Block(K))
        L̃ = _triangular_matrix(Val('L'), UNIT, view(A, Block(K,K)))
        materialize!(Lmul(L̃, b_2))
        JR = blockrowstart(A,K):Block(K-1)
        if !isempty(JR)
            muladd!(one(T), view(A, Block(K), JR), view(b,JR), one(T), b_2)
        end
    end

    dest
end

@inline function materialize!(M::MatLmulVec{<:TriangularLayout{UPLO,UNIT,<:AbstractBlockLayout},
                                   <:AbstractStridedLayout}) where {UPLO,UNIT}
    U,x = M.A,M.B
    T = eltype(M)
    @boundscheck size(U,1) == size(x,1) || throw(BoundsError())
    if hasmatchingblocks(U)
        _matchingblocks_triangular_mul!(Val(UPLO), Val(UNIT), triangulardata(U), x)
    else # use default
        x_1 = PseudoBlockArray(copy(x), (axes(U,2),))
        x_2 = PseudoBlockArray(x, (axes(U,1),))
        _block_muladd!(one(T), U, x_1, zero(T), x_2)
    end
end

for UNIT in ('U', 'N')
    @eval begin
        @inline function materialize!(M::MatLdivVec{<:TriangularLayout{'U',$UNIT,<:AbstractBlockLayout},
                                        <:AbstractStridedLayout})
            U,dest = M.A,M.B
            T = eltype(dest)

            A = triangulardata(U)
            if !hasmatchingblocks(A) # Use default for now
                return materialize!(Ldiv{TriangularLayout{'U',$UNIT,UnknownLayout},
                                         typeof(MemoryLayout(dest))}(U, dest))
            end

            @boundscheck size(A,1) == size(dest,1) || throw(BoundsError())

            # impose block structure
            b = PseudoBlockArray(dest, (axes(A,1),))

            N = blocksize(A,1)

            for K = N:-1:1
                b_2 = view(b, Block(K))
                Ũ = _triangular_matrix(Val('U'), Val($UNIT), view(A, Block(K,K)))
                materialize!(Ldiv(Ũ, b_2))

                if K ≥ 2
                    KR = blockcolstart(A, K):Block(K-1)
                    V_12 = view(A, KR, Block(K))
                    b̃_1 = view(b, KR)
                    muladd!(-one(T), V_12, b_2, one(T), b̃_1)
                end
            end

            dest
        end

        @inline function materialize!(M::MatLdivVec{<:TriangularLayout{'L',$UNIT,<:AbstractBlockLayout},
                                        <:AbstractStridedLayout})
            L,dest = M.A, M.B
            T = eltype(dest)
            A = triangulardata(L)
            if !hasmatchingblocks(A) # Use default for now
                return materialize!(Ldiv{TriangularLayout{'L',$UNIT,UnknownLayout},
                                         typeof(MemoryLayout(dest))}(L, dest))
            end

            @boundscheck size(A,1) == size(dest,1) || throw(BoundsError())

            # impose block structure
            b = PseudoBlockArray(dest, (axes(A,1),))

            N = blocksize(A,1)

            for K = 1:N
                b_2 = view(b, Block(K))
                L̃ = _triangular_matrix(Val('L'), Val($UNIT), view(A, Block(K,K)))
                materialize!(Ldiv(L̃, b_2))

                if K < N
                    KR = Block(K+1):blockcolstop(A, K)
                    V_12 = view(A, KR, Block(K))
                    b̃_1 = view(b, KR)
                    muladd!(-one(T), V_12, b_2, one(T), b̃_1)
                end
            end

            dest
        end
    end
end


text:
function tr6!(r::Robot)#главная функция
    side=Sud
    side_to_border=Ost
    side_of_movement=Nord
    vert=0
    hor=0
    vert_1=0
    hor_1=0
    while isborder(r,side)==false
    move!(r,side)
    vert=vert + 1
    end
    side=next_move(side)
    while isborder(r,side)==false
    move!(r,side)
    hor=hor + 1
    end
    side=Sud
    while isborder(r,side)==false
    move!(r,side)
    vert_1=vert_1 + 1
    end
    side=next_move(side)
    while isborder(r,side)==false
    move!(r,side)
    hor_1=hor_1 + 1
    end
    find_border!(r,side_to_border,side_of_movement)
    #Поиск перегородки
    mark_perimetr!(r,side_to_border,side_of_movement)
    #Маркировка перегородки
    side=West
    while isborder(r,side)==false
    move!(r,side)
    end
    side=next_move(side)
    while isborder(r,side)==false
    move!(r,side)
    end
    side=Sud
    while isborder(r,side)==false
    move!(r,side)
    end
    side=next_move(side)
    while isborder(r,side)==false
    move!(r,side)
    end
    #Робот снова в углу
    side=Ost
    movements!(r,side,hor_1)
    side=Nord
    movements!(r,side,vert_1)
    side=Ost
    movements!(r,side,hor)
    side=Nord
    movements!(r,side,vert)
    #Перемещение робота по запомненному пути

    end

    function find_border!(r::Robot,side_to_border::HorizonSide, side_of_movement::HorizonSide)
    while isborder(r,side_to_border)==false
    if isborder(r,side_of_movement)==false
    move!(r,side_of_movement)
    else
    move!(r,side_to_border)
    side_of_movement=inverse(side_of_movement)
    end
    end
    end

    function find_corner!(r::Robot,side::HorizonSide)
    vert=0
    hor=0
    vert_1=0
    hor_1=0
    while isborder(r,side)==false
    move!(r,side)
    end
    side=next_move(side)
    while isborder(r,side)==false
    move!(r,side)
    end
    side=Sud
    while isborder(r,side)==false
    move!(r,side)
    end
    side=next_move(side)
    while isborder(r,side)==false
    move!(r,side)
    end
    end

    function mark_perimetr!(r::Robot,side_to_border::HorizonSide, side_of_movement::HorizonSide)
    for i ∈ 1:4
    while isborder(r,side_to_border)==true
    putmarker!(r)
    move!(r,side_of_movement)
    end
    side_to_border=next_bord(side_to_border)
    side_of_movement=next_move(side_of_movement)
    putmarker!(r)
    movements!(r,side_of_movement,1)
    i=i+1
    end
    end

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

    next_bord(side_to_border::HorizonSide)=HorizonSide(mod(Int(side_to_border)+3,4))

    next_move(side_of_movement::HorizonSide)=HorizonSide(mod(Int(side_of_movement)+3,4))

    movements!(r::Robot,side::HorizonSide,vert::Int)=
    for i in 1:vert
    move!(r,side)
    end


text:
using ModelingToolkit
@variables t
D = Differential(t)

@variables xxx(t)
[xxx]

function InPin(;name)
  @parameters x
  defaults = Dict(x => 13.37)
  ODESystem(Equation[],t,Num[],[x]; name, defaults)
end

function OutPin(;name)
  vars = @variables x(t) xx(t)
  defaults = Dict(x => 0.0)
  ODESystem(Equation[D(x)~xx],t,[x,xx],Num[]; name, defaults)
end

function create_pins(in::Integer, out::Integer)
  inpins = [InPin(;name=Symbol("x$(i)_InPin")) for i in 1:in]
  outpins = [OutPin(;name=Symbol("x$(i)_OutPin")) for i in 1:out]
  inpins, outpins
end

function Network(;name)
    @variables x(t)=0.0 a(t)
    @parameters p=0.4
    eqs = [D(x) ~ a, a ~ 0.1p]
    ODESystem(eqs;name)
end

function Model(;name)
    op = OutPin(;name=:op)
    net = Network(;name=:net)
    eqs = [op.xx ~ net.x]
    # eqs = Equation[]
    systems = [op,net]
    ODESystem(eqs,t;name,systems)
end

# using Symbolics

model = Model(;name=:model)
sys = structural_simplify(model)
@nonamespace n = model.net
@nonamespace o = model.op
u0 = Dict(
    n.x => 0.0,
    o.x => 0.0,
)
p = Dict(
    n.p => 0.4,
)
defs = ModelingToolkit.get_defaults(sys)
prob = ODEProblem(sys, u0, (0.0,1.0), p)
sol = solve(prob, Tsit5())
plot(sol[o.x])
plot(sol[n.x])
sol
plot(sol(sol.t,Val{1})[1,:])
plot(sol(sol.t,Val{1})[2,:])
plot(sol[1,:])
plot(sol[2,:])

p = [0.1]
u0b = rand(1,10)
function prob_func(prob,i,repeat)
    u0 = @view u0b[:,i]
    remake(prob; u0)
end

function output_func(sol,i)
    sol[:, end], false
end

function loss(p,model)
    ensemble_prob = EnsembleProblem(prob; prob_func, output_func, safetycopy=false) # TODO: safetycopy ???
    sol = solve(ensemble_prob, Tsit5(), EnsembleThreads(), trajectories=size(u0b,2),
              sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)), save_everystep=false, save_start=false)
    sum(Array(sol))
end

cbg = function (p,l;doplot=false)
    display(l)
    return false
end

optfun = GalacticOptim.OptimizationFunction((θ,p) -> loss(θ,model), GalacticOptim.AutoZygote())
optfunc = GalacticOptim.instantiate_function(optfun, p, GalacticOptim.AutoZygote(), nothing)
optprob = GalacticOptim.OptimizationProblem(optfunc, p)
GalacticOptim.solve(optprob, ADAM(), cb = cbg, maxiters=100)

using DiffEqCallbacks, OrdinaryDiffEq, LinearAlgebra
prob = ODEProblem((du,u,p,t) -> du .= u, rand(4), (0.0,1.0))
saved_values = SavedValues(Float64, Vector{Float64})
cb = SavingCallback((u,t,integrator)->integrator(t,Val{1})[:,1], saved_values, saveat = 0.0:0.1:1.0)
sol = solve(prob, Tsit5(), saveat=0.1, callback=cb)
as = Array(sol)
plot(s)
sol[:]
saved_values.saveval
plot(saved_values.saveval)

using ModelingToolkit
using ModelingToolkit: get_defaults
using OrdinaryDiffEq

@variables t
D = Differential(t)

function SubSys(;name)
  @variables x(t)=0.0
  @parameters p=0.3
  ODESystem([x~p],t,[x],[p];name)
end
function Network(;name)
  @variables x(t)=0.0
  @named subsys = SubSys()
  ODESystem([D(x)~subsys.x],t,[x],[];name,systems=[subsys])
end

@named net = Network()
sys = structural_simplify(net)
prob = ODEProblem(sys,get_defaults(sys),(0.0,0.1))
sol = solve(prob,Tsit5())

@nonamespace x = net.x
@nonamespace subsys = net.subsys
sol[x]
sol[subsys.x]

sol[net.x]
sol[net.subsys.x]


text:
let
    warmup()
end


text:
using EnvNoise
using Test, BenchmarkTools, Base.Threads, DelimitedFiles
using JLD

#@time @threads for i=100:200
#    filename = "R:\\CMPH-Windfarm Field Study\\Hallett\\set1\\Recording-3.$i.pti"
#    ptiread(filename)
#end

#signal = rand(8192*600)

#Aweighting = weightingFilter("A-weighting", 8192)
#@btime yA = filt(Aweighting,signal)

dir = "R:\\CMPH-Windfarm Field Study\\Hallett\\set2\\"
allpti = filter(x->endswith(x, ".pti"), readdir(dir; join=true))

#filename = "R:\\CMPH-Windfarm Field Study\\Hallett\\set1\\Recording-3.1.pti"
for i=1#:100

    filename = allpti[i]
    signal, Fs, Date, Time = ptiread(filename)

    res = 0.1
    nfft = floor(Int,Fs/res)
    ovlap = floor(Int, 0.5*nfft)
    pxx = pwelch(signal, hanning(nfft),ovlap,nfft,Fs)[1]
    pxx = Float32.(pxx) # using Float32 precision

    # save results
    set_i = filename[38:41]
    name_i = filename[43:end-4]
    savedir = "R:\\CMPH-Windfarm Field Study\\Duc Phuc Nguyen\\3. Spectrum quantification\\Hallett_spectrum"
    #save("$savedir\\$set_i-$name_i.jld", "pxx", pxx, "date", Date,"time",Time) # save spectrum
    #save("$savedir\\info-$set_i-$name_i.jld", "date", Date,"time",Time) # save DateTime

end

#pxx2 = load("pxx2.jld")
#writedlm("pxx.txt", pxx)

#=
@time open("pxx.txt", "w") do file
    write(file, pxx)
end
=#

#@testset "EnvNoise.jl" begin
    # Write your tests here.
#end


text:
#=

Author: Benjamin Hemingway

This code returns the Gaussian quadrature nodes and weights associated
with a Truncated Normal distribution.

This code is based on the method outlined by John Burkardt in the following
presentation:
https://people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf
=#
#---------------------------------------------------------------------#
# Generate Moments
#---------------------------------------------------------------------#
function normalDistMoment(μ::Real,σ::Real, n::Int)
  _normalDistMoment(μ, σ, n, -Inf, Inf)
end
function normalDistMoment(μ::Real,σ::Real, n::Int, l::Float64, u::Float64)
  _normalDistMoment(μ, σ, n, l, u)
end
function normalDistMoment(n::Int)
  _normalDistMoment(n, -Inf, Inf)
end
function normalDistMoment(n::Int, l::Float64, u::Float64)
  _normalDistMoment(n, l, u)
end

function _normalDistMoment(μ::Real,σ::Real, n::Int, l::Real, u::Real)

  f0 = Distributions.Normal(0,1)
  mom = zeros(n+1)
  mom[1] = 1.0

  if l>-Inf && u<Inf # If the distribution is doubly truncated
    α = (l - μ)/σ
    β = (u - μ)/σ

    # Calculate the recursive part of moments
    L = Array{Float64,1}(undef,n+1)
    L[1] = 1.0
    L[2] = ( Distributions.pdf(f0,α)-Distributions.pdf(f0,β)
    ) /  ( Distributions.cdf(f0,β)-Distributions.cdf(f0,α) )

    for i = (3:n+1)
      L[i] = (( α^(i-2)*Distributions.pdf(f0,α)-β^(i-2)*
      Distributions.pdf(f0,β) ) / ( Distributions.cdf(f0,β)
      -Distributions.cdf(f0,α) ) + (i-2)*L[i-2])
    end

    for k = 1:n, i =0:k
      mom[k+1] = (binomial(k,i)*L[i+1])*(σ^i)*(μ^(k-i)) + mom[k+1]
    end

  elseif l>-Inf && u==Inf # If the distribution is upper truncated
    α = (l - μ)/σ
    L = Array{Float64,1}(undef,n+1)
    L[1] = 1.0
    L[2] = ( Distributions.pdf(f0,α) )/( Distributions.ccdf(f0,α) )
    for i = (3:n+1)
      L[i] = (((α^(i-2))*Distributions.pdf(f0,α))/(Distributions.ccdf(f0,α))+(i-2)*L[i-2])
    end

    for k = 1:n, i =0:k
      mom[k+1] = ((binomial(k,i)*L[i+1])*(σ^i)*(μ^(k-i)) + mom[k+1])
    end

  elseif l==-Inf && u<Inf # If the distribution is lower truncated
    L = Array{Float64,1}(undef,n+1)
    β = (u - μ)/σ

    L[1] = 1.0
    L[2] = ( -Distributions.pdf(f0,β) )/( Distributions.cdf(f0,β) )

    for i = (3:n+1)
      L[i] = (-(β^(i-2))*Distributions.pdf(f0,β) / Distributions.cdf(f0,β)  + (i-2)*L[i-2])
    end

    for k = 1:n, i =0:k
      mom[k+1] = ((binomial(k,i)*L[i+1])*(σ^i)*(μ^(k-i)) + mom[k+1])
    end

  else # If the distribution is non-truncated
    mom[2] = μ
    for i = 3:n+1
      ind = floor(Int,(i-1)/2)
      for j = 0:ind
        if 2j-1<0
          # Convention is (-1)!!=1 but not followed by Combinatorics.doublefactorial
          mom[i] += binomial(i-1,2j)*σ^(2j)*μ^(i-2j-1)
        else
          mom[i] += binomial(i-1,2j)*σ^(2j)*μ^(i-2j-1)*Combinatorics.doublefactorial(2j-1)
        end
      end
    end
  end
  # Return the vector of moments
  return mom
end

function _normalDistMoment(n::Int, α::Real, β::Real)
  # calculation for standard normal
  f0 = Distributions.Normal(0,1)
  mom = zeros(n+1)
  mom[1] = 1.0

  if α>-Inf && β<Inf # If the distribution is doubly truncated

    # Calculate the recursive part of moments
    L = Array{Float64,1}(undef,n+1)
    L[1] = 1.0
    L[2] = ( Distributions.pdf(f0,α)-Distributions.pdf(f0,β)
    ) /  ( Distributions.cdf(f0,β)-Distributions.cdf(f0,α) )

    for i = (3:n+1)
      L[i] = (( α^(i-2)*Distributions.pdf(f0,α)-β^(i-2)*
      Distributions.pdf(f0,β) ) / ( Distributions.cdf(f0,β)
      -Distributions.cdf(f0,α) ) + (i-2)*L[i-2])
    end

    for k = 1:n
      mom[k+1] += L[k+1]
    end

  elseif α>-Inf && β==Inf # If the distribution is upper truncated
    L = Array{Float64,1}(undef,n+1)
    L[1] = 1.0
    L[2] = ( Distributions.pdf(f0,α) )/( Distributions.ccdf(f0,α) )
    for i = (3:n+1)
      L[i] = (((α^(i-2))*Distributions.pdf(f0,α))/(Distributions.ccdf(f0,α))+(i-2)*L[i-2])
    end

    for k = 1:n
      mom[k+1] += L[k+1]
    end

  elseif α==-Inf && β<Inf # If the distribution is lower truncated
    L = Array{Float64,1}(undef,n+1)

    L[1] = 1.0
    L[2] = ( -Distributions.pdf(f0,β) )/( Distributions.cdf(f0,β) )

    for i = (3:n+1)
      L[i] = (-(β^(i-2))*Distributions.pdf(f0,β) / Distributions.cdf(f0,β)  + (i-2)*L[i-2])
    end

    for k = 1:n
      mom[k+1] += L[k+1]
    end

  else # If the distribution is non-truncated
    mom[2] = 0.0
    for i=3:n+1
      if rem(i-1,2)==0
        mom[i] += Combinatorics.doublefactorial(i-2)
      end
    end
  end

  # Return the vector of moments
  return mom

end