text:
issettable(a::AbstractArray) = error("It is not known if $(typeof(a)) is settable")
issettable(::Array) = true
issettable(::BitArray) = true
issettable(::AbstractRange) = false
issettable(::LinearIndices) = false
issettable(::CartesianIndices) = false
issettable(a::SubArray) = issettable(parent(a))
issettable(a::PermutedDimsArray) = issettable(parent(a))
issettable(a::Base.ReinterpretArray) = issettable(parent(a))
issettable(a::Base.ReshapedArray) = issettable(parent(a))
issettable(::SparseVector) = true
issettable(::SparseMatrixCSC) = true

struct FrozenArray{T, N, A <: AbstractArray{T, N}} <: AbstractArray{T, N}
    a::A
end

Base.parent(a::FrozenArray) = a.a
Base.length(a::FrozenArray) = length(parent(a))
Base.size(a::FrozenArray) = size(parent(a))
Base.axes(a::FrozenArray) = axes(parent(a))
Base.@propagate_inbounds Base.getindex(a::FrozenArray, i...) = getindex(parent(a), i...)
issettable(a::FrozenArray) = false

function freeze(a::AbstractArray)
    if issettable(a)
        return FrozenArray{eltype(a), ndims(a), typeof(a)}(a)
    else
        return a
    end
end

function thaw(a::AbstractArray)
    if issettable(a)
        return a
    else
        out = similar(a)
        copy!(out, a)
        return out
    end
end

thaw(a::FrozenArray) = parent(a) # DECISION: Perhaps we want the defensive copy? Wait for the compiler to reason about this for `Array`?


text:
using state, breeder, evaluator, logger, replacement, criterion, initialiser,
      population

"""
Builder for algorithm.genetic.
"""
type GeneticAlgorithmDefinition
  output::AbstractString
  replacement::Replacement
  loggers::Vector{LoggerDefinition}
  termination::Vector{Criterion}
  population::PopulationDefinition
  evaluator::EvaluatorDefinition

  GeneticAlgorithmDefinition() =
    new("output", replacement.generational(), [], [])
end

type GeneticAlgorithm <: Algorithm
  state::State
  evaluator::Evaluator
  replacement::Replacement
  termination::Vector{Criterion}
  loggers::Vector{Logger}
  initialiser::Initialiser
  output::AbstractString

  GeneticAlgorithm() = new()
end

"""
Composes an `algorithm.genetic` instance from information provided within a
supplied builder.
"""
function compose!(def::GeneticAlgorithmDefinition)
  alg = GeneticAlgorithm()
  pop = population.compose!(def.population)
  alg.state = State(pop)
  alg.evaluator = compose!(def.evaluator,  pop)
  alg.replacement = def.replacement
  alg.termination = def.termination
  alg.initialiser = initialiser.DefaultInitialiser()
  alg.output = abspath(def.output)
  alg.loggers = Logger[compose!(l) for l in def.loggers]
  alg
end

"""
TODO: Explain algorithm.genetic
"""
function genetic(def::Function)
  d = GeneticAlgorithmDefinition()
  def(d)
  d
end

function run!(a::GeneticAlgorithm)
  reset!(a.state)
  prepare!(a.loggers, a.output)
  initialise!(a.initialiser, a.state.population)
  evaluate!(a.evaluator, a.state)
  scale!(a.state.population)

  # Record the best individual from the population.
  #pbest = gbest = best(a.state.population)

  #prepare!(a.state.population)

  while !any(c -> is_satisfied(c, a.state), a.termination)
    breed!(a.state.population)
    evaluate!(a.evaluator, a.state)
    scale!(a.state.population)

    # Record the best individual from the population.
    # --- What does this mean for co-evolution?
    #pbest = best(a.state.population)
    #gbest = best([pbest, gbest])

    # Debugging, for now.

    replace!(a.replacement, a.state)
    scale!(a.state.population)
    logger.call!(a.loggers, a.state)
    a.state.iterations += 1
  end
  logger.close!(a.loggers)
end

ext: jl


text:
gauss_nodes_weights(p::Type{P}, n) where {P <: Chebyshev} =
    FastGaussQuadrature.gausschebyshev(n)

gauss_nodes_weights(p::Type{P}, n) where {α, β, P <: Jacobi{α, β}} =
    FastGaussQuadrature.gaussjacobi(n, α, β)

include("Orthogonal/FastLegendre.jl")

eval_basis(::Type{P}, n, x::Float64) where {P <: Legendre} = FastLegendre.fastlegendre(n,x)
eval_basis(::Type{P}, n, x::Union{Int, Int32, Float16, Float32}) where {P <: Legendre} =
    FastLegendre.fastlegendre(n,float(x))
gauss_nodes_weights(p::Type{P}, n) where {P <: Legendre} =
    FastGaussQuadrature.gausslegendre(n)


text:
using SparseArrays

@inline function flat2_i(
    i :: Int64,
    j :: Int64,
    Nx :: Int64,
    Ny :: Int64,
)
    return i + (j-1) * Nx
end

@inline function flat3_i(
    k :: Int64,
    i :: Int64,
    j :: Int64,
    Nz :: Int64,
    Nx :: Int64,
    Ny :: Int64,
)
    return k + (i-1) * Nz + (j-1) * Nz * Nx
end

@inline function cyc(i::Int64, N::Int64)
    return mod(i-1, N) + 1
end

mutable struct AdvectionSpeedUpMatrix

    mtx_interp_U :: AbstractArray{Float64, 2}
    mtx_interp_V :: AbstractArray{Float64, 2}
    mtx_DIV_X    :: AbstractArray{Float64, 2}
    mtx_DIV_Y    :: AbstractArray{Float64, 2}
    mtx_DIV_Z    :: AbstractArray{Float64, 2}
    mtx_GRAD_X   :: AbstractArray{Float64, 2}
    mtx_GRAD_Y   :: AbstractArray{Float64, 2}
    mtx_GRAD_Z   :: AbstractArray{Float64, 2}
    mtx_CURV_X   :: AbstractArray{Float64, 2}
    mtx_CURV_Y   :: AbstractArray{Float64, 2}
    mtx_CURV_Z   :: AbstractArray{Float64, 2}

    function AdvectionSpeedUpMatrix(;
        gi             :: DisplacedPoleCoordinate.GridInfo,
        Nx             :: Int64,
        Ny             :: Int64,
        Nz_bone        :: Int64,
        Nz             :: AbstractArray{Int64, 2},
        mask3          :: AbstractArray{Float64, 3},
        noflux_x_mask3 :: AbstractArray{Float64, 3},
        noflux_y_mask3 :: AbstractArray{Float64, 3},
        Δzs            :: AbstractArray{Float64, 3},
        hs             :: AbstractArray{Float64, 3},
    )
   # mtx_GRAD_x = spzeros(Float64, Nx * Ny * Nz_bone, (Nx+1) * Ny * Nz_bone)
   # mtx_GRAD_y = spzeros(Float64, Nx * Ny * Nz_bone, Nx * (Ny+1) * Nz_bone)
   # mtx_GRAD_z = spzeros(Float64, Nx * Ny * Nz_bone, Nx * Ny * (Nz_bone+1))

   # mtx_CURV :: AbstractArray{Float64, 2}

        mtx_interp_U = spzeros(Float64, (Nx+1) * Ny    , Nx     * Ny    )
        mtx_interp_V = spzeros(Float64, Nx     * (Ny+1), Nx     * Ny    )
        mtx_DIV_X    = spzeros(Float64, Nx     * Ny    , (Nx+1) * Ny    )
        mtx_DIV_Y    = spzeros(Float64, Nx     * Ny    , Nx     * (Ny+1))
        mtx_DIV_Z    = spzeros(Float64, Nx     * Ny    , Nx     * Ny    )
        mtx_GRAD_X   = spzeros(Float64, (Nx+1) * Ny    , Nx     * Ny    )
        mtx_GRAD_Y   = spzeros(Float64, Nx     * (Ny+1), Nx     * Ny    )
        mtx_GRAD_Z   = spzeros(Float64, (Nz_bone+1)    , Nz_bone        )
        mtx_CURV_X   = spzeros(Float64, Nx     * Ny    , (Nx+1) * Ny    )
        mtx_CURV_Y   = spzeros(Float64, Nx     * Ny    , Nx     * (Ny+1))
        mtx_CURV_Z   = spzeros(Float64, Nz_bone        , (Nz_bone+1)    )

        println("Making Interp Matrix")
        # ===== [BEGIN] Making interp matrix =====
        # x
        for i=1:Nx+1, j=1:Ny  # iterate through bounds
            ib   = flat2_i(i           , j, Nx+1, Ny)
            ic_e = flat2_i(cyc(i  ,Nx) , j, Nx  , Ny)
            ic_w = flat2_i(cyc(i-1,Nx) , j, Nx  , Ny)

            #u_bnd[k, i, j] = u[k, i-1, j] * (1.0 - weight_e[i, j]) + u[k, i, j] * weight_e[i, j]
            mtx_interp_U[ib, ic_w] = 1.0 - gi.weight_e[i, j]
            mtx_interp_U[ib, ic_e] = gi.weight_e[i, j]
        end

        # y
        for i=1:Nx, j=2:Ny   # iterate through bounds
                    ib   = flat2_i(i, j  , Nx, Ny+1)
                    ic_n = flat2_i(i, j  , Nx, Ny  )
                    ic_s = flat2_i(i, j-1, Nx, Ny  )

                    #v_bnd[k, i, j] = v[k, i, j-1] * (1.0 - weight_n[i, j]) + v[k, i, j] * weight_n[i, j]
                    mtx_interp_V[ib, ic_s] = 1.0 - gi.weight_n[i, j]
                    mtx_interp_V[ib, ic_n] = gi.weight_n[i, j]
        end
        # ===== [END] Making interp matrix =====

        println("Making Divergence Matrix")
        # ===== [BEGIN] Making divergent matrix =====
        # x and y
        for i=1:Nx, j=1:Ny  # iterate through face centers
            ic = flat2_i(i, j, Nx  , Ny)

            # X direction
            ib_e   = flat2_i(i+1, j, Nx+1, Ny)
            ib_w   = flat2_i(i  , j, Nx+1, Ny)

            mtx_DIV_X[ic, ib_e] =   gi.DY[i+1, j] / gi.dσ[i, j]
            mtx_DIV_X[ic, ib_w] = - gi.DY[i  , j] / gi.dσ[i, j]

            # Y direction
            ib_n   = flat2_i(i, j+1, Nx, Ny+1)
            ib_s   = flat2_i(i, j  , Nx, Ny+1)

            mtx_DIV_Y[ic, ib_n] =   gi.DX[i, j+1] / gi.dσ[i, j]
            mtx_DIV_Y[ic, ib_s] = - gi.DX[i, j  ] / gi.dσ[i, j]
        end

        for k=1:Nz_bone
            ic = k

            # Z direction
            ib_t   = k
            ib_b   = k+1

            mtx_DIV_Z[ic, ib_t] =   1.0 / hs[k, i, j]
            mtx_DIV_Z[ic, ib_b] = - 1.0 / hs[k, i, j]

        end
        # ===== [END] Making divergent matrix =====

        println("Making GRAD Matrix")
        # ===== [BEGIN] Making GRAD matrix =====
        # x
        for i=1:Nx+1, j=1:Ny  # iterate through bounds
            for k=1:Nz[cyc(i, Nx), j]  # Bounds Nx+1 is the same as the bound 1
                if noflux_x_mask3[k, i, j] != 0.0
                    ib   = flat_i(k, i           , j, Nz_bone, Nx+1, Ny)
                    ic_e = flat_i(k, cyc(i  ,Nx) , j, Nz_bone, Nx  , Ny)
                    ic_w = flat_i(k, cyc(i-1,Nx) , j, Nz_bone, Nx  , Ny)

                    # ( qs[k, i, j] - qs[k, i-1, j] ) / gi.dx_w[i, j]
                    mtx_GRAD_X[ib, ic_e] =   1.0 / gi.dx_w[cyc(i, Nx), j]
                    mtx_GRAD_X[ib, ic_w] = - 1.0 / gi.dx_w[cyc(i, Nx), j]
                end
            end
        end

        # y
        for i=1:Nx, j=2:Ny   # iterate through bounds
            for k=1:Nz[i, j]
               if noflux_y_mask3[k, i, j] != 0.0
                    ib   = flat_i(k, i, j  , Nz_bone, Nx, Ny+1)
                    ic_n = flat_i(k, i, j  , Nz_bone, Nx, Ny  )
                    ic_s = flat_i(k, i, j-1, Nz_bone, Nx, Ny  )

                    # ( qs[k, i, j] - qs[k, i, j-1] ) / gi.dy_s[i, j]
                    mtx_GRAD_Y[ib, ic_n] =   1.0 / gi.dy_s[i, j]
                    mtx_GRAD_Y[ib, ic_s] = - 1.0 / gi.dy_s[i, j]
                end
            end
        end

        # z
        for i=1:Nx, j=1:Ny   # iterate through bounds

            if mask3[1, i, j] == 0.0
                continue
            end

            _Nz = Nz[i, j]

            # The frist layer of w is zero -- means do not assign any value to this row

            # Assign from the second row
            for k=2:_Nz
                ib   = flat_i(k  , i, j, Nz_bone+1, Nx, Ny)
                ic_t = flat_i(k-1, i, j, Nz_bone  , Nx, Ny)
                ic_b = flat_i(k  , i, j, Nz_bone  , Nx, Ny)

                # ( qs[k-1, i, j] - qs[k, i, j] ) / Δzs[k-1, i, j]
                mtx_GRAD_Z[ib, ic_t] =   1.0 / Δzs[k-1, i, j]
                mtx_GRAD_Z[ib, ic_b] = - 1.0 / Δzs[k-1, i, j]

                # Bottom is the same as the layer right above it.
                # This is necessary because the thickness of last layer might be
                # very thin due to topography to create instability during
                # doing adveciton.
                if k == _Nz
                    ib_b = flat_i(k+1, i, j, Nz_bone+1, Nx, Ny)
                    mtx_GRAD_Z[ib_b, ic_t] = mtx_GRAD_Z[ib, ic_t]
                    mtx_GRAD_Z[ib_b, ic_b] = mtx_GRAD_Z[ib, ic_b]
                end
            end

        end
        # ===== [END] Making GRAD matrix =====

        println("Making CURV Matrix")
        # ===== [BEGIN] Making CURV matrix =====
        for i=1:Nx, j=1:Ny
            for k=1:Nz[i, j]
                if mask3[k, i, j] == 0.0
                    break
                end

                ic = flat_i(k, i, j, Nz_bone, Nx  , Ny)

                # X direction
                # CURV_x[k, i, j] = ( GRAD_bnd_x[k, i+1, j  ] - GRAD_bnd_x[k  , i, j] ) / gi.dx_c[i, j]
                ib_e   = flat_i(k, i+1, j, Nz_bone, Nx+1, Ny)
                ib_w   = flat_i(k, i  , j, Nz_bone, Nx+1, Ny)

                mtx_CURV_X[ic, ib_e] =   1.0 / gi.dx_c[i, j]
                mtx_CURV_X[ic, ib_w] = - 1.0 / gi.dx_c[i, j]

                # Y direction
                # CURV_y[k, i, j] = ( GRAD_bnd_y[k, i  , j+1] - GRAD_bnd_y[k  , i, j] ) / gi.dy_c[i, j]
                ib_n   = flat_i(k, i, j+1, Nz_bone, Nx, Ny+1)
                ib_s   = flat_i(k, i, j  , Nz_bone, Nx, Ny+1)

                mtx_CURV_Y[ic, ib_n] =   1.0 / gi.dy_c[i, j]
                mtx_CURV_Y[ic, ib_s] = - 1.0 / gi.dy_c[i, j]

                # Z direction
                # CURV_z[k, i, j] = ( GRAD_bnd_z[k, i  , j  ] - GRAD_bnd_z[k+1, i, j] ) / hs[k, i, j]
                ib_t   = flat_i(k  , i, j, Nz_bone+1, Nx, Ny)
                ib_b   = flat_i(k+1, i, j, Nz_bone+1, Nx, Ny)

                mtx_CURV_Z[ic, ib_t] =   1.0 / hs[k, i, j]
                mtx_CURV_Z[ic, ib_b] = - 1.0 / hs[k, i, j]
            end
        end
        # ===== [END] Making CURV matrix =====

        return new(
            mtx_interp_U,
            mtx_interp_V,
            mtx_DIV_X,
            mtx_DIV_Y,
            mtx_DIV_Z,
            mtx_GRAD_X,
            mtx_GRAD_Y,
            mtx_GRAD_Z,
            mtx_CURV_X,
            mtx_CURV_Y,
            mtx_CURV_Z,
        )
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
import LibGit2
import Pkg

function all_versions(repo_url::AbstractString)
    _temp_dir::String = mktempdir()
    atexit(() -> rm(_temp_dir; force = true, recursive = true))
    _repo_path = joinpath(_temp_dir, "GITREPO")
    LibGit2.clone(repo_url, _repo_path)
    _repo = LibGit2.GitRepo(_repo_path)
    _ref_list = LibGit2.ref_list(_repo)
    _all_versions = Vector{VersionNumber}(undef, 0)
    for _ref in _ref_list
        if occursin(r"^refs\/tags\/v(\d*).(\d*).(\d*)$", _ref)
            _m = match(r"^refs\/tags\/v(\d*).(\d*).(\d*)$", _ref)
            _v = VersionNumber("$(_m[1]).$(_m[2]).$(_m[3])")
            push!(_all_versions, _v)
        end
    end
    rm(_temp_dir; force = true, recursive = true)
    unique!(_all_versions)
    sort!(_all_versions)
    return _all_versions
end

latest_semver_version(repo_url) = maximum(all_versions(repo_url))

function main()
    dir_containing_this_file::String = joinpath(splitpath(@__DIR__)...)

    package_name_to_url = Pkg.TOML.parsefile(
        joinpath(dir_containing_this_file, "package_url_list.toml")
        )["package_url_list"]
    PREDICTMD_URL = package_name_to_url["PredictMD"]
    PREDICTMDAPI_URL = package_name_to_url["PredictMDAPI"]
    PREDICTMDEXTRA_URL = package_name_to_url["PredictMDExtra"]
    PREDICTMDFULL_URL = package_name_to_url["PredictMDFull"]
    PREDICTMD_LATEST_SEMVER_VERSION = string(latest_semver_version(PREDICTMD_URL))
    PREDICTMDAPI_LATEST_SEMVER_VERSION = string(latest_semver_version(PREDICTMDAPI_URL))
    PREDICTMDEXTRA_LATEST_SEMVER_VERSION = string(latest_semver_version(PREDICTMDEXTRA_URL))
    PREDICTMDFULL_LATEST_SEMVER_VERSION = string(latest_semver_version(PREDICTMDFULL_URL))

    offline_toml_in_path = joinpath(dir_containing_this_file, "offline.toml.template")
    offline_toml_out_path = joinpath(dir_containing_this_file, "offline.toml")

    offline_toml = read(offline_toml_in_path, String)

    offline_toml = replace(offline_toml, "{{PREDICTMD_LATEST_SEMVER_VERSION}}" => PREDICTMD_LATEST_SEMVER_VERSION)
    offline_toml = replace(offline_toml, "{{PREDICTMDAPI_LATEST_SEMVER_VERSION}}" => PREDICTMDAPI_LATEST_SEMVER_VERSION)
    offline_toml = replace(offline_toml, "{{PREDICTMDEXTRA_LATEST_SEMVER_VERSION}}" => PREDICTMDEXTRA_LATEST_SEMVER_VERSION)
    offline_toml = replace(offline_toml, "{{PREDICTMDFULL_LATEST_SEMVER_VERSION}}" => PREDICTMDFULL_LATEST_SEMVER_VERSION)

    write(offline_toml_out_path, offline_toml)

    @info("Wrote file $(offline_toml_out_path)")
end

main()


text:
using Test
using DelimitedFiles
using Base.Iterators
using FastGaussQuadrature

# Those corresponding verifiable data are obtained from orginal matlab functions at
# https://bitbucket.org/sbarbot/bssa-2016237/src/master/
# https://bitbucket.org/sbarbot/bssa-2018058/src/default/
@testset "Volume Hex8" begin
    epsv11 = 11e-6
    epsv12 = 5e-6
    epsv13 = 6e-6
    epsv22 = 7e-6
    epsv23 = 8e-6
    epsv33 = 13e-6
    L = 40e3
    W = 20e3
    T = 10e3
    q1 = -10e3
    q2 = 10e3
    q3 = 1e3
    theta = 40.0
    G = 1.0
    nu = 0.25

    x1s = range(-10000.0, stop=10000.0, step=1000.0)
    x2s = range(-10000.0, stop=10000.0, step=1000.0)
    x3s = range(-30.0, stop=-20.0, step=5.0) # in real case should not be negative
    xxs = product(x1s, x2s, x3s)

    @testset "Displacement" begin
        u_truth = readdlm(joinpath(@__DIR__, "data/test_disp_vol_hex8.dat"), ' ', Float64)
        for (i, x) in enumerate(xxs)
            u = _disp_vol_hex8(x..., q1, q2, q3, L, T, W, theta, epsv11, epsv12, epsv13, epsv22, epsv23, epsv33, G, nu)
            @test u ≈ u_truth[i,:]
        end
    end

    @testset "Stress" begin
        u_truth = readdlm(joinpath(@__DIR__, "data/test_stress_vol_hex8.dat"), ' ', Float64)
        for (i, x) in enumerate(xxs)
            u = _stress_vol_hex8(x..., q1, q2, q3, L, T, W, theta, epsv11, epsv12, epsv13, epsv22, epsv23, epsv33, G, nu)
            @test ifelse(all(map(isnan, u)) && all(map(isnan, u_truth[i,:])), true, u ≈ u_truth[i,:])
        end
    end
end

@testset "Volume Tet4" begin
    epsv11 = 11e0
    epsv12 = 5e0
    epsv13 = 6e0
    epsv22 = 7e0
    epsv23 = 8e0
    epsv33 = 13e0
    G = 1.0
    nu = 0.25
    A = [1.0, 2.0, 3.0]
    B = [2.0, 5.0, 4.0]
    C = [11.0, 13.0, 12.0]
    D = [6.0, 5.0, 4.0]

    x1s = range(-100.0, stop=100.0, step=20.0)
    x2s = range(-100.0, stop=100.0, step=20.0)
    x3s = range(20.0, stop=30.0, step=5.0)
    xxs = product(x1s, x2s, x3s)
    qd = gausslegendre(15)

    @testset "Displacement" begin
        u_truth = readdlm(joinpath(@__DIR__, "data/test_disp_vol_tet4.dat"), ' ', Float64)
        for (i, x) in enumerate(xxs)
            u = _disp_vol_tet4(qd, x..., A, B, C, D, epsv11, epsv12, epsv13, epsv22, epsv23, epsv33, nu)
            @test u ≈ u_truth[i, :]
        end
    end

    @testset "Stress" begin
        u_truth = readdlm(joinpath(@__DIR__, "data/test_stress_vol_tet4.dat"), ' ', Float64)
        for (i, x) in enumerate(xxs)
            u = _stress_vol_tet4(qd, x..., A, B, C, D, epsv11, epsv12, epsv13, epsv22, epsv23, epsv33, G, nu)
            @test u ≈ u_truth[i, :]
        end
    end
end

@testset "Volume Inplane rect" begin
    epsv22 = 7e-6
    epsv23 = 8e-6
    epsv33 = 9e-6
    T = 10e3
    W = 20e3
    q2 = 15e3
    q3 = 1e-3
    phi = 0.0
    G = 1.0
    nu = 0.25

    x3 = range(-10000.0, stop=-1000.0, step=1000)
    x2 = range(-10000.0, stop=10000.0, step=2000)
    xxs = product(x2, x3)

    @testset "Stress" begin
        u_truth = readdlm(joinpath(@__DIR__, "data", "test_stress_volinplane_rect.dat"), ' ', Float64)
        for (i, x) in enumerate(xxs)
            u = _stress_volinplane_rect(x..., q2, q3, T, W, phi, epsv22, epsv23, epsv33, G, nu)
            @test u ≈ u_truth[i,:]
        end
    end
end

@testset "Coordinates: ENU <-> NED" begin
    x, y, z, qx, qy, qz, L, T, W, θ, ϵxx, ϵxy, ϵxz, ϵyy, ϵyz, ϵzz, G = rand(17)
    ν = 0.2
    u1 = disp_vol_hex8(x, y, -z, qx, qy, -qz, T, L, W, θ, ϵxx, ϵxy, ϵxz, ϵyy, ϵyz, ϵzz, G, ν)
    u2 = _disp_vol_hex8(y, x, z, qy, qx, qz, L, T, W, θ, ϵyy, ϵxy, -ϵyz, ϵxx, -ϵxz, ϵzz, G, ν)
    @test u1[1] ≈ u2[2] && u1[2] ≈ u2[1] && u1[3] ≈ -u2[3]

    u1 = stress_vol_hex8(x, y, -z, qx, qy, -qz, T, L, W, θ, ϵxx, ϵxy, ϵxz, ϵyy, ϵyz, ϵzz, G, ν)
    u2 = _stress_vol_hex8(y, x, z, qy, qx, qz, L, T, W, θ, ϵyy, ϵxy, -ϵyz, ϵxx, -ϵxz, ϵzz, G, ν)
    @test u1[1] ≈ u2[4] && u1[2] ≈ u2[2] && u1[3] ≈ -u2[5] && u1[4] ≈ u2[1] && u1[5] ≈ -u2[3] && u1[6] ≈ u2[6]

    u1 = stress_volinplane_rect(x, -z, qx, -qz, T, W, θ, ϵxx, ϵxz, ϵzz, G, ν)
    u2 = _stress_volinplane_rect(x, z, qx, qz, T, W, θ, ϵxx, -ϵxz, ϵzz, G, ν)
    @test u1[1] ≈ u2[1] && u1[2] ≈ -u2[2] && u1[3] ≈ u2[3]

    A, B, C, D = [rand(3) for _ in 1: 4]
    A2, B2, C2, D2 = map(x -> [x[2], x[1], -x[3]], [A, B, C, D])
    qd = gausslegendre(15)

    u1 = disp_vol_tet4(qd, x, y, -z, A2, B2, C2, D2, ϵxx, ϵxy, ϵxz, ϵyy, ϵyz, ϵzz, ν)
    u2 = _disp_vol_tet4(qd, y, x, z, A, B, C, D, ϵyy, ϵxy, -ϵyz, ϵxx, -ϵxz, ϵzz, ν)
    @test u1[1] ≈ u2[2] && u1[2] ≈ u2[1] && u1[3] ≈ -u2[3]

    u1 = stress_vol_tet4(qd, x, y, -z, A2, B2, C2, D2, ϵxx, ϵxy, ϵxz, ϵyy, ϵyz, ϵzz, G, ν)
    u2 = _stress_vol_tet4(qd, y, x, z, A, B, C, D, ϵyy, ϵxy, -ϵyz, ϵxx, -ϵxz, ϵzz, G, ν)
    @test u1[1] ≈ u2[4] && u1[2] ≈ u2[2] && u1[3] ≈ -u2[5] && u1[4] ≈ u2[1] && u1[5] ≈ -u2[3] && u1[6] ≈ u2[6]
end


text:
"""
    FunctionField{X, Y, Z, C, F, G} <: AbstractField{X, Y, Z, F, G}

An `AbstractField` that returns a function evaluated at location `(X, Y, Z)` (and time, if
`C` is not `Nothing`) when indexed at `i, j, k`.
"""
struct FunctionField{X, Y, Z, C, F, G} <: AbstractField{X, Y, Z, F, G}
     func :: F
     grid :: G
    clock :: C

    """
        FunctionField{X, Y, Z}(func, grid; clock=nothing) where {X, Y, Z}

    Returns a `FunctionField` on `grid` and at location `X, Y, Z`.

    If `clock` is not specified, then `func` must be a function with signature
    `func(x, y, z)`. If clock is specified, `func` must be a function with signature
    `func(x, y, z, t)`, where `t` is internally determined from `clock.time`.

    A FunctionField will return the result of `func(x, y, z [, t])` at `X, Y, Z` on
    `grid` when indexed at `i, j, k`.
    """
    function FunctionField{X, Y, Z}(func, grid; clock=nothing) where {X, Y, Z}
        return new{X, Y, Z, typeof(clock), typeof(func), typeof(grid)}(func, grid, clock)
    end
end

"""
    FunctionField(L::Tuple, func, grid)

Returns a stationary `FunctionField` on `grid` and at location `L = (X, Y, Z)`,
where `func` is callable with signature `func(x, y, z)`.
"""
FunctionField(L::Tuple, func, grid) = FunctionField{L[1], L[2], L[3]}(func, grid)

# Ordinary functions needed for fields
architecture(::FunctionField) = nothing
data(f::FunctionField) = f
Base.parent(f::FunctionField) = f

@inline Base.getindex(f::FunctionField{X, Y, Z, <:Nothing}, i, j, k) where {X, Y, Z} =
    f.func(xnode(X, i, f.grid), ynode(Y, j, f.grid), znode(Z, k, f.grid))

@inline Base.getindex(f::FunctionField{X, Y, Z}, i, j, k) where {X, Y, Z} =
    f.func(xnode(X, i, f.grid), ynode(Y, j, f.grid), znode(Z, k, f.grid), f.clock.time)

@inline (f::FunctionField)(x, y, z) = f.func(x, y, z, f.clock.time)
@inline (f::FunctionField{X, Y, Z, <:Nothing})(x, y, z) where {X, Y, Z} = f.func(x, y, z)

# set! for function fields
set!(u, f::FunctionField) = set!(u, (x, y, z) -> f.func(x, y, z, f.clock.time))
set!(u, f::FunctionField{X, Y, Z, <:Nothing}) where {X, Y, Z} = set!(u, f.func)

Adapt.adapt_structure(to, f::FunctionField{X, Y, Z}) where {X, Y, Z} =
    FunctionField{X, Y, Z}(adapt(to, func), grid, clock)


text:
export defΩ, SimpleΩ, SimpleΩ

using ..IDS, ..Util, ..Tagging
using Random: AbstractRNG
import ..Basis: replacetags
import ..Var

"Simplest, immutable Omega"
struct SimpleΩ{TAGS <: Tags, T, S} <: AbstractΩ
  data::T
  tags::TAGS
  subspace::S
end

Basis.subspace(ω::SimpleΩ) = ω.subspace

# AbstractDict interface
Base.keys(ω::SimpleΩ) = Base.keys(ω.data)
Base.values(ω::SimpleΩ) = Base.values(ω.data)
Base.pairs(ω::SimpleΩ) = Base.pairs(ω.data)
Base.getindex(ω::SimpleΩ, id) = ω.data[id]

idtype(::SimpleΩ{TAGS, Dict{T, V}}) where {TAGS, T, V} = T

SimpleΩ(data) = SimpleΩ(data, Tags())
replacetags(ω::SimpleΩ, tags) = SimpleΩ(ω.data, tags)
# (T::Type{<:Distribution})(π::SimpleΩ, args...) = ω.data[scope(ω)]

function Var.recurse(exo::Var.PrimRandVar, ω::SimpleΩ)
  newid = append(ω.subspace, exo.id)
  newexo = Var.Member(newid, exo.class)
  ω.data[newexo]::eltype(newexo.class)
end

traithastag(t::Type{SimpleΩ{TAGS, T}}, tag) where {TAGS, T} = traithastag(TAGS, tag)
traits(::Type{SimpleΩ{TAGS, T}}) where {TAGS, T} = traits(TAGS)

function resolve(dist, id, ω::SimpleΩ)
  # @assert false
  @show ω
  @show id
  @show ω.data[id]
  d, val = ω.data[id]
  val::eltype(dist)
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# This file was generated, do not modify it. # hide
yc = coerce(y, OrderedFactor);


text:
# This file is for functions that explictly mess with the dimensions of a NameDimsArray

"""
    rename(nda::NamedDimsArray, names)

Returns a new `NameDimsArray` with the given dimension `names`.
`rename` outright replaces the names; while still wrapping the same backing array.
Unlike the constructor, it does not require that new names are compatible
with the old names (though you do still need to match the number of dimensions).
"""
rename(nda::NamedDimsArray, names) = NamedDimsArray(parent(nda), names)

function Base.dropdims(nda::NamedDimsArray; dims)
    numerical_dims = dim(nda, dims)
    data = dropdims(parent(nda); dims=numerical_dims)
    L = remaining_dimnames_after_dropping(names(nda), numerical_dims)
    return NamedDimsArray{L}(data)
end

function Base.permutedims(nda::NamedDimsArray{L}, perm) where {L}
    numerical_perm = dim(nda, perm)
    new_names = permute_dimnames(L, numerical_perm)

    return NamedDimsArray{new_names}(permutedims(parent(nda), numerical_perm))
end

function Base.selectdim(nda::NamedDimsArray, s::Symbol, i)
    return selectdim(nda, dim(nda, s), i)
end

for f in (
    :(Base.transpose),
    :(Base.adjoint),
    :(Base.permutedims),
    :(LinearAlgebra.pinv)
)
    # Vector
    @eval function $f(nda::NamedDimsArray{L,T,1}) where {L,T}
        new_names = (:_, first(L))
        return NamedDimsArray{new_names}($f(parent(nda)))
    end

    # Vector Double Transpose
    if f != :(Base.permutedims)
        @eval function $f(nda::NamedDimsArray{L,T,2,A}) where {L,T,A<:CoVector}
            new_names = (last(L),)  # drop the name of the first dimensions
            return NamedDimsArray{new_names}($f(parent(nda)))
        end
    end

    # Matrix
    @eval function $f(nda::NamedDimsArray{L,T,2}) where {L,T}
        new_names = (last(L), first(L))
        return NamedDimsArray{new_names}($f(parent(nda)))
    end
end

# reshape
# For now we only implement the version that drops dimension names
Base.reshape(x::NamedDimsArray, d::Vararg{Union{Colon, Int}}) = reshape(parent(x), d)


text:
using POMDPs
using POMDPSimulators
using StructArrays

function discounted_avg_return(h::AbstractSimHistory)
    γ = discount(h)
    r_total = sum(r.*γ^(i-1) for (i, r) in enumerate(h[:r]))
    return sum(r_total)/length(r_total)
end

function undiscounted_avg_return(h::AbstractSimHistory)
    r_total = sum(r for r in h[:r])
    return sum(r_total) / length(r_total)
end

function online_evaluate(solvers, mdp, nevals, max_steps)
    sims = Sim[]
    for (sname, solver) in solvers
        for i in 1:nevals
            push!(sims, Sim(mdp, solve(solver, mdp);
                            rng=MersenneTwister(42 + i),
                            max_steps=max_steps,
                            metadata=Dict(:policy=>sname))
                  )
        end
    end

    # Note: some bug in progress_pmap
    df = run_parallel(sims; show_progress=false) do sim, hist
        result = (nsteps=length(hist), discret=discounted_avg_return(hist), undiscret=undiscounted_avg_return(hist))
        # Note: Because this is being done lazily on the backend
        if :info in keys(hist[1])
            return merge(result, map(sum, fieldarrays(StructArray(hist[:info]))))
        end
        return result
    end

    return df
end

function offline_evaluate(solvers, mdp, nsolves, nevals, max_steps)
    sims = Sim[]
    for (sname, solver) in solvers
        for j in 1:nsolves
            sol = solve(solver, mdp)
            for i in 1:nevals
                push!(sims, Sim(mdp, sol; rng=MersenneTwister(42 +i),
                                max_steps=max_steps,
                                metadata=Dict(:policy=>sname))
                      )
            end
        end
    end

    df = run_parallel(sims; show_progress=false) do sim, hist
        result = (nsteps=length(hist), discret=discounted_avg_return(hist), undiscret=undiscounted_avg_return(hist))
        # Note: Because this is being done lazily on the backend
        if :info in keys(hist[1])
            return merge(result, map(sum, fieldarrays(StructArray(hist[:info]))))
        end
        return result
    end

    return df
end


text:
psi = qbimps(rand, 6, 2, 1)
T = cmpo_ising(1.0)

B_canonical(T, psi)
A_canonical(T, psi)


text:
using LinearAlgebra
import Base: iterate

struct SumUp2K
    n :: Int
    k :: Int
end

function iterate(A::SumUp2K, prev::Vector{Int})
    a = prev[1]
    next = prev
    i_lasta = 1
    for i = 1:A.n
        if prev[i] == a
            i_lasta = i
        end
        if prev[i] > a+1
            next[i] -= 1
            next[i_lasta] += 1
            return next, next
        end
    end
    nothing
end

function iterate(A::SumUp2K)
    start = ones(Int, A.n)
    extras = A.k - A.n
    extras == 0 && return start, start
    bin = A.n
    while extras > 0
        if extras > (A.n-1)
            start[bin] = A.n
            bin -= 1
            extras -= A.n-1
        else
            start[bin] += extras
            return start, start
        end
    end
    nothing
end
iterate(A, nothing) = iterate(A)

function natural_latin_square(n::Int, k::Int)
    A = SumUp2K(n, k)
    for possible_diag in A
        ans = natural_latin_square(n, diagm(0 => possible_diag))
        if !isnothing(ans)
            return ans
        end
    end
    nothing
end

function natural_latin_square(n::Int, A::Matrix{Int})
    for j in 1:n
        for i in 1:n
            if A[i, j] == 0
                for candidate in 1:n
                    if ispossible(A, candidate, i, j)
                        newA = copy(A)
                        newA[i, j] = candidate
                        ans = natural_latin_square(n, newA)
                        if !isnothing(ans)
                            return ans
                        end
                    end
                end
                return nothing
            end
        end
    end
    A
end

function ispossible(A, a, i, j)
    !(a in A[i, :]) && !(a in A[:, j])
end

function printmat(A::Matrix)
    m, n = size(A)
    for i = 1:m
        for j = 1:n
            print(A[i,j], " ")
        end
        println()
    end
end

function solution()
    n_cases = parse(Int, readline())
    for x=1:n_cases
        n, k = [parse(Int, i) for i in split(readline())]
        A = natural_latin_square(n, k)
        if isnothing(A)
            println("Case #$x: IMPOSSIBLE")
        else
            println("Case #$x: POSSIBLE")
            printmat(A)
        end
    end
end

solution()

# input
# 2
# 3 6
# 2 3

# output
# Case #1: POSSIBLE
# 2 1 3
# 3 2 1
# 1 3 2
# Case #2: IMPOSSIBLE


text:
using Spark
using Base.Test

# write your own tests here
sc = SparkContext(master="local")
txt = parallelize(sc, ["hello", "world"])
rdd = map_partitions(txt, it -> map(s -> length(s), it))

@test count(rdd) == 2
@test reduce(rdd, +) == 10

close(sc)


text:
module MyPackagenew

using ForwardDiff

include("extrafile.jl")

export my_funct, der
end # module


text:
using Test
using YaoLang


text:
# Based on crying baby problem (Kochenderfer 2015)
# State: hungry = true; not hungry = false
# Action: feed = true; do nothing = false
# Observation: crying = true; not crying = false

mutable struct Baby3RPOMDP <: RPOMDP{Bool, Bool, Symbol}
    r_feed::Float64
    r_hungry::Float64
    discount::Float64
end

mutable struct PBox
    lower::Distribution
    upper::Distribution
end
PBox(d::Distribution) = PBox(d,d) # degenerate PBox

Baby3RPOMDP(r_feed, r_hungry) = Baby3RPOMDP(r_feed, r_hungry, 0.9)
Baby3RPOMDP() = Baby3RPOMDP(-5.0, -10.0, 0.9)

# updater(problem::Baby3RPOMDP) = DiscreteUpdater(problem)

initial_state_distribution(::Baby3RPOMDP) = BoolDistribution(0.0)

const states_const = [true, false] # hungry, full
const actions_const = [true, false] # feed, no action
const observations_const = [:quiet, :crying, :yelling]

observations(::Baby3RPOMDP) = observations_const

n_states(::Baby3RPOMDP) = 2
state_index(::Baby3RPOMDP, s::Bool) = s ? 1 : 2
action_index(::Baby3RPOMDP, a::Bool) = a ? 1 : 2
observation_index(::Baby3RPOMDP, z::Symbol) = z == :quiet ? 1 : z == :crying ? 2 : 3
obs_index(p::Baby3RPOMDP, z::Symbol) = observation_index(p::Baby3RPOMDP, z::Symbol)
n_actions(::Baby3RPOMDP) = 2
n_observations(::Baby3RPOMDP) = 3

# transition data
const t_array = [BoolDistribution(0.0) BoolDistribution(1.0);
    BoolDistribution(0.0) BoolDistribution(0.1)]

function transition(rpomdp::Baby3RPOMDP, s::Bool, a::Bool)
    return t_array[state_index(rpomdp, s), action_index(rpomdp, a)]
end

# pbox representation of obesrvation uncertainty
const o_pbox = [PBox(Categorical([0.15, 0.55, 0.3]), Categorical([0.35, 0.6, 0.05])) PBox(Categorical([0.8, 0.1, 0.1]));
    PBox(Categorical([0.15, 0.55, 0.3]), Categorical([0.35, 0.6, 0.05])) PBox(Categorical([0.8, 0.1, 0.1]))]

# returns pbox representation
function observation(rpomdp::Baby3RPOMDP, a::Bool, sp::Bool)
    o_pbox[action_index(rpomdp, a), state_index(rpomdp, sp)]::PBox
end
observation(rpomdp::Baby3RPOMDP, s::Bool, a::Bool, sp::Bool) = observation(rpomdp, a, sp)

# lower and upper probability representation of the pbox
function pzinterval(rpomdp::Baby3RPOMDP)
    nz = n_observations(rpomdp)
    na = n_actions(rpomdp)
    ns = n_states(rpomdp)
    pzl = Array{Float64}(nz,na,ns)
    pzu = Array{Float64}(nz,na,ns)
    for (aind, a) in enumerate(RPOMDPModels.ordered_actions(rpomdp)), (tind, t) in enumerate(RPOMDPModels.ordered_states(rpomdp))
        uncset = RPOMDPs.observation(rpomdp, a, t)
        cdfl = cumsum(pdf.(uncset.lower, support(uncset.lower)))
        cdfu = cumsum(pdf.(uncset.upper, support(uncset.upper)))
        pzl[1,aind,tind] = cdfl[1]
        pzu[1,aind,tind] = cdfu[1]
        for zind = 2:nz
            pzl[zind,aind,tind] = cdfl[zind] - cdfu[zind - 1]
            pzu[zind,aind,tind] = cdfu[zind] - cdfl[zind - 1]
        end
    end
    pzl, pzu
end

function pzainterval(rpomdp::Baby3RPOMDP, a)
    aind = action_index(rpomdp, a)
    nz = n_observations(rpomdp)
    ns = n_states(rpomdp)
    pzl = Array{Float64}(ns,nz)
    pzu = Array{Float64}(ns,nz)
    for (tind, t) in enumerate(RPOMDPModels.ordered_states(rpomdp))
        uncset = RPOMDPs.observation(rpomdp, a, t)
        cdfl = cumsum(pdf.(uncset.lower, support(uncset.lower)))
        cdfu = cumsum(pdf.(uncset.upper, support(uncset.upper)))
        pzl[tind,1] = cdfl[1]
        pzu[tind,1] = cdfu[1]
        for zind = 2:nz
            pzl[tind,zind] = cdfl[zind] - cdfu[zind - 1]
            pzu[tind,zind] = cdfu[zind] - cdfl[zind - 1]
        end
    end
    pzl, pzu
end

function pinterval(rpomdp::RPOMDP)
    ns = n_states(rpomdp)
    nz = n_observations(rpomdp)
    na = n_actions(rpomdp)
    p = Array{Float64}(ns,nz,ns,na)
    plower = Array{Float64}(ns,nz,ns,na)
    pupper = Array{Float64}(ns,nz,ns,na)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na
        t = RPOMDPModels.ordered_states(rpomdp)[tind]
        s = RPOMDPModels.ordered_states(rpomdp)[sind]
        a = RPOMDPModels.ordered_actions(rpomdp)[aind]
        z = RPOMDPModels.ordered_observations(rpomdp)[zind]
        dt = RPOMDPModels.transition(rpomdp, s, a)
        pt = pdf(dt, t)
        pzl, pzu = RPOMDPModels.pzainterval(rpomdp, a)
        plower[tind, zind, sind, aind] = pzl[tind, zind] * pt
        pupper[tind, zind, sind, aind] = pzu[tind, zind] * pt
    end
    plower, pupper
end

function reward(rpomdp::Baby3RPOMDP, s::Bool, a::Bool)
    r = 0.0
    if s # hungry
        r += rpomdp.r_hungry
    end
    if a # feed
        r += rpomdp.r_feed
    end
    return r
end

discount(p::Baby3RPOMDP) = p.discount

function generate_o(p::Baby3RPOMDP, s::Bool, rng::AbstractRNG, otype::String)
    d = observation(p, true, s, otype) # obs distribution not action dependant
    return observations[rand(rng, d)]
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -9.734460520574642 <= q <= 2.116911729869532)

                     add_NL_constraint(m, :(swish(-0.8046949419530325 + -0.2278494590955038*swish(0.1770285029406451 + -0.5067338360802731*swish(-0.7998422425711689 + 0.5883237869246916*swish(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*swish(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + -0.09889031570146667*swish(-0.7622628423532047 + 0.6518912481623378*swish(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*swish(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5])))) + -0.18111143847071443*swish(-0.12446319977672715 + 0.7018706953017659*swish(-0.7998422425711689 + 0.5883237869246916*swish(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*swish(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + 0.5275006439747312*swish(-0.7622628423532047 + 0.6518912481623378*swish(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*swish(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))) + swish(-0.4025778812055263 + 0.08853333939955377*swish(0.1770285029406451 + -0.5067338360802731*swish(-0.7998422425711689 + 0.5883237869246916*swish(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*swish(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + -0.09889031570146667*swish(-0.7622628423532047 + 0.6518912481623378*swish(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*swish(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5])))) + -0.4269192713141794*swish(-0.12446319977672715 + 0.7018706953017659*swish(-0.7998422425711689 + 0.5883237869246916*swish(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.47235910506966317*swish(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))) + 0.5275006439747312*swish(-0.7622628423532047 + 0.6518912481623378*swish(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5])) + 0.4006903724184485*swish(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module TestExamples1DMHD

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1d_dgsem")

@testset "MHD" begin

  @trixi_testset "elixir_mhd_alfven_wave.jl with initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [1.440611823425164e-15, 1.1373567770134494e-14, 3.024482376149653e-15, 2.0553143516814395e-15, 3.9938347410210535e-14, 3.984545392098788e-16, 2.4782402104201577e-15, 1.551737464879987e-15],
      linf = [1.9984014443252818e-15, 1.3405943022348765e-14, 3.3584246494910985e-15, 3.164135620181696e-15, 7.815970093361102e-14, 8.881784197001252e-16, 2.886579864025407e-15, 2.942091015256665e-15],
      initial_condition = initial_condition_constant,
      tspan = (0.0,1.0))
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [1.13950171e-05, 5.70552422e-07, 2.35934311e-05, 2.35934311e-05, 4.88745541e-06, 1.19632242e-16, 2.38901782e-05, 2.38901782e-05],
      linf = [3.00941311e-05, 8.20108595e-07, 1.19890130e-04, 1.19890130e-04, 1.47569139e-05, 2.22044605e-16, 1.22875440e-04, 1.22875440e-04])
  end

  @trixi_testset "elixir_mhd_alfven_wave.jl with flux_derigs_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [1.4717982305221483e-5, 1.182263100134546e-5, 2.357120395081263e-5, 2.357120395088376e-5, 1.8950292712848924e-6, 1.1963224165731992e-16, 2.3661336206751837e-5, 2.366133620672686e-5],
      linf = [6.156423569070313e-5, 3.474803933270032e-5, 0.00011882557236660152, 0.00011882557236685132, 6.8618872851589074e-6, 2.220446049250313e-16, 0.00012195115527005584, 0.00012195115526993094],
      volume_flux = flux_derigs_etal)
  end

  @trixi_testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [5.86009540e-02, 8.16048158e-02, 5.46791194e-02, 5.46791194e-02, 1.54509265e-01, 4.13046273e-17, 5.47637521e-02, 5.47637521e-02],
      linf = [1.10014999e-01, 1.81982581e-01, 9.13611439e-02, 9.13611439e-02, 4.23831370e-01, 1.11022302e-16, 9.93731761e-02, 9.93731761e-02])
  end

  @trixi_testset "elixir_mhd_briowu_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_briowu_shock_tube.jl"),
      l2   = [0.17764301067932906, 0.19693621875378622, 0.3635136528288642, 0.0, 0.3757321708837591, 8.593007507325741e-16, 0.36473438378159656, 0.0],
      linf = [0.5601530250396535, 0.43867368105486537, 1.0960903616351099, 0.0, 1.0551794137886303, 4.107825191113079e-15, 1.5374410890043144, 0.0])
  end

  @trixi_testset "elixir_mhd_torrilhon_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_torrilhon_shock_tube.jl"),
      l2   = [0.45702524713047815, 0.4792453322653316, 0.3406561762729464, 0.4478037841848423, 0.9204646319222438, 1.3216517820475193e-16, 0.2889772516628848, 0.2552071770333965],
      linf = [1.221507578527215, 0.8927062606782392, 0.8536666880170606, 0.9739151339925485, 1.6602343391033632, 2.220446049250313e-16, 0.686743631371773, 0.6555428379163893])
  end

  @trixi_testset "elixir_mhd_ryujones_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ryujones_shock_tube.jl"),
      l2   = [2.34809441e-01, 3.92255943e-01, 8.23575546e-02, 1.75599624e-01, 9.61613519e-01, 6.60825891e-17, 2.15346454e-01, 1.07006529e-01],
      linf = [6.40732148e-01, 9.44889516e-01, 3.54932707e-01, 8.54060243e-01, 2.07757711e+00, 1.11022302e-16, 4.92584725e-01, 2.49526561e-01],
      tspan = (0.0, 0.1))
  end

  @trixi_testset "elixir_mhd_shu_osher_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shu_osher_shock_tube.jl"),
      l2   = [1.01126210e+00, 8.27157902e+00, 1.30882545e+00, 0.00000000e+00, 5.21930435e+01, 6.56538824e-16, 1.01022340e+00, 0.00000000e+00],
      linf = [2.87172004e+00, 2.26438057e+01, 4.16672442e+00, 0.00000000e+00, 1.35152372e+02, 3.44169138e-15, 2.83556069e+00, 0.00000000e+00],
      tspan = (0.0, 0.2))
  end

  @trixi_testset "elixir_mhd_shu_osher_shock_tube.jl with flipped shock direction" begin
    # Include this elixir to make `initial_condition_shu_osher_shock_tube_flipped` available, which is used below
    trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shu_osher_shock_tube.jl"), tspan=(0.0, 0.0))
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shu_osher_shock_tube.jl"),
      l2   = [1.01539817e+00, 8.29625810e+00, 1.29548008e+00, 0.00000000e+00, 5.23565514e+01, 3.18641825e-16, 1.00485291e+00, 0.00000000e+00],
      linf = [2.92876280e+00, 2.28341581e+01, 4.11643561e+00, 0.00000000e+00, 1.36966213e+02, 1.55431223e-15, 2.80548864e+00, 0.00000000e+00],
      initial_condition = initial_condition_shu_osher_shock_tube_flipped,
      boundary_conditions=BoundaryConditionDirichlet(initial_condition_shu_osher_shock_tube_flipped),
      tspan = (0.0, 0.2))
  end
end

end # module


text:
include("../../../temperature_render.jl")

mass_vars = [10, 100, 10e6]
spin_vars = [0.0, 0.5, 0.998]
angle_vars = [15.0, 40.0, 85.0]

mass_fixed = 10
spin_fixed = 0.998
angle_fixed = 85.0

for mass in mass_vars
    hmap, cache, title, new_img = temperature_render(;
                                                    mass=mass,
                                                    spin=spin_fixed,
                                                    obs_angle=angle_fixed,
                                                    tolerance=1e-12,
                                                    dtmax=0.5,
                                                    size_multiplier=10,
                                                    resolution=3000
                                                    )
    title!(title)
    png(hmap, "render_mass=$mass.png")
end

for spin in spin_vars
    hmap, cache, title, new_img = temperature_render(;
                                                    spin=spin,
                                                    mass=mass_fixed,
                                                    obs_angle=angle_fixed,
                                                    tolerance=1e-12,
                                                    dtmax=0.5,
                                                    size_multiplier=10,
                                                    resolution=3000
                                                    )
    title!(title)
    png(hmap, "render_spin=$spin.png")
end

for angle in angle_vars
    hmap, cache, title, new_img = temperature_render(;
                                                    obs_angle=angle,
                                                    spin=spin_fixed,
                                                    mass=mass_fixed,
                                                    tolerance=1e-12,
                                                    dtmax=0.5,
                                                    size_multiplier=10,
                                                    resolution=3000
                                                    )
    title!(title)
    png(hmap, "render_angle=$angle.png")
end


text:
@testset "NullActor" begin
    a = NullActor()
    @test act!(a, nothing)
    @test Actors.init!(a) == nothing
end

@testset "StepLimitActor" begin
    n = 10
    a = StepLimitActor(n)
    @test a.max_step_limit == n
    @test ! a.hit_max_step_limit
    a.hit_max_step_limit = true
    @test Actors.init!(a) == nothing
    @test ! a.hit_max_step_limit
end


text:
export wrapexpr,
argspellings,
eclass2type,
wrapexpr_header,
epreamble

"""
    wrapexpr_header(path; kw...)

Wraps C++ header at path.
"""
function wrapexpr_header(path; kw...)
    topcu = cindex.parse_header(path, cplusplus=true; kw...)
    ret = topcu |> analyze |> wrapexpr
    ret
end

function wrapexpr(c::WrappedClass, wc::WrapperConfig=WrapperConfig())
    cxxspell = cxxspelling(c)
    jlspell = jlspelling(c, wc)
    jlspell_destroy = jlspelling(Val{:Destructor} , wc)
    blk = quote end
    append!(blk.args, eclass2type(cxxspell, jlspell, jlspell_destroy).args)
    append!(blk.args, map(x -> wrapexpr(x, wc), c.constructors))
    append!(blk.args, map(x -> wrapexpr(x, wc), c.methods))
    blk
end

function wrapexpr(c::WrappedConstructor, wc::WrapperConfig=WrapperConfig())
    args = argspellings(c)
    cxxargs = argspellings_jl2cxx(c)
    sconstructor = c |> spelling |> Symbol
    cxxspell = cxxspelling(c)
    jlspell = jlspelling(c, wc)
    body = Expr(:block,
        eequals(:ptr, ecxxnew(cxxspell, cxxargs)),
        Expr(:call, cxxspell, :ptr)
    )
    efunction(jlspell, args, body)
end

isoperator(s::String) = startswith(s, "operator")
isoperator(m::CXXMethod) = m |> spelling |> isoperator
isoperator(m::WrappedMethod) = m |> raw |> isoperator
jlspelling_operator(x::CXXMethod) = x |> spelling |> jlspelling_operator
jlspelling_operator(x::String) = x[9:end] |> Symbol

cxxspelling(x) = Symbol(spelling(x))
jlspelling(x, wc::WrapperConfig) = jlspelling(raw(x), wc)
function jlspelling(x::CXXMethod, wc::WrapperConfig)
    if isoperator(x)
        # should also import from Base for overloading!
        x |> jlspelling_operator
    else
        x |> spelling |> wc.rename.method |> Symbol
    end
end
function jlspelling(x::ClassDecl, wc::WrapperConfig)
    x |> spelling |> wc.rename.class |> Symbol
end
jlspelling(x::WrappedConstructor, wc::WrapperConfig) = jlspelling(x.parent, wc)
jlspelling(x::Type{Val{:Destructor}}, wc::WrapperConfig) = wc.rename.destructor |> Symbol

function wrapexpr(x::WrappedMethod, wc::WrapperConfig)
    if isoperator(x)
        # wrapexpr_operator2(x, wc)
        println("not wrapping operator $x")
        :()
    else
        wrapexpr_method(x, wc)
    end
end

function wrapexpr_method(m::WrappedMethod, wc::WrapperConfig=WrapperConfig())
    jlspell = jlspelling(m, wc)
    jlargs = argspellings(m)
    cxxspell = cxxspelling(m)
    cxxargs = argspellings_jl2cxx(m)
    obj = :obj
    T = jlspelling(m.parent, wc)
    objT = etyped(obj, T)
    body = quote
        ret = $(emethodcall(:(jl2cxx($obj)), cxxspell, cxxargs))
        cxx2jl(ret)
    end
    efunction(jlspell, [objT; jlargs], body)
end

function methodsymbols(x::WrappedClass, wc::WrapperConfig=WrapperConfig())
    Symbol[jlspelling(m, wc) for m in x.methods]
end
function exportsymbols(x::WrappedClass, wc::WrapperConfig)
    [jlspelling(x, wc), methodsymbols(x, wc)...]
end
function exportsymbols(x::WrappedHeader, wc::WrapperConfig)
    mapreduce(node -> exportsymbols(node, wc), vcat, [], x.classnodes)
end

function wrapexpr(x::WrappedHeader, wc::WrapperConfig=WrapperConfig())
    content = mapreduce(node -> wrapexpr(node, wc).args, vcat, [], x.classnodes)

    Expr(:block,
    eexport(exportsymbols(x, wc)),
    epreamble().args..., # TODO this should not be reproduced for every header, but only once per project
    content...
    )
end

# binary operator
# function wrapexpr_operator2(m::WrappedMethod, wc::WrapperConfig)
#     op = jlspelling(m, wc)
#     body = Expr(:block,
#         eequals(:ocxx1, :(jl2cxx(obj1))),
#         eequals(:ocxx2, :(jl2cxx(obj2))),
#         eicxx( "&(\$ocxx1) $op &(\$ocxx2);" )
#     )
#     f = Expr(Symbol("."), :Base, QuoteNode(op))
#     T = jlspelling(m.parent, wc)
#     args = (etyped(:obj1, T), etyped(:obj2, T))
#     efunction(f, args, body)
# end

argspellings(m::WrappedMethod) = argspellings(m.args)
argspellings(m::WrappedConstructor) = argspellings(m.args)
function argspellings(args::Vector)
    ret = Symbol[]
    for (i, arg) in enumerate(args)
        nm = arg |> spelling
        if nm == ""
            arg_i = Symbol("arg_$i")
        else
            arg_i = Symbol(nm)
        end
        push!(ret, arg_i)
    end
    ret
end
function argspellings_jl2cxx(m)
    args = argspellings(m)
    [:(jl2cxx($arg)) for arg in args]
end

function epreamble()
    quote
        using Cxx

        """
            jl2cxx(x)

        Convert julia object to cxx object. This method is called on each argument before doing a
        C++ method call.
        """
        jl2cxx(x) = x
        jl2cxx(s::AbstractString) = pointer(String(s))

        """
            cxx2jl(x)

        Convert C++ object to julia object. This method is called after each C++ method call.
        """
        cxx2jl(x) = x
        # cxx2jl(p::Ptr{UInt8}) = unsafe_string(p)?

    end
end

function edestroy(jlspell_class, jlspell_destroy)
    obj = "obj"
    body = eicxx("delete \$($obj.pointer);")
    obj_typed = etyped(Symbol(obj), jlspell_class)
    efunction(jlspell_destroy, (obj_typed,), body)
end

function eclass2type(cxxspell, jlspell, jlspell_destroy)
    sT = jlspell
    sPtrT = Symbol("Ptr", jlspell)
    blk = quote
        typealias $sPtrT $(ecxxtype(cxxspell))
        immutable $sT
            pointer::$sPtrT
            $sT(p::$sPtrT) = new(p)
        end
        function cxx2jl(ptr::$sPtrT)
            $sT(ptr)
        end
        function jl2cxx(obj::$sT)
            obj.pointer
        end
    end
    push!(blk.args, edestroy(jlspell, jlspell_destroy))
    blk
end


text:
module GreensFunctions

using Statistics
using Random
using LinearAlgebra
using FFTW
using Parameters

using ..Models: AbstractModel, mulMᵀ!, mulM!, update_model!
using ..Utilities: get_index, get_site, get_τ, reshaped

using ..KPMPreconditioners

export EstimateGreensFunction, update!, measure_GΔ0, measure_GΔ0_G0Δ, measure_GΔ0_GΔ0, measure_GΔΔ_G00, estimate

###########################################
## CODE FOR CALCULATING GREEN'S FUNCTION ##
###########################################

"""
A type for facilitating the stochastic estimation needed to make measurements
"""
mutable struct EstimateGreensFunction{T<:AbstractFloat,Tfft<:AbstractFFTs.Plan,Tifft<:AbstractFFTs.Plan}

    """
    Number of random vectors to be use.
    """
    nᵥ::Int

    """
    First in pair of random vectors to convolve together.
    """
    n₁::Int

    """
    Second in pair of random vectors to convolve together.
    """
    n₂::Int

    """
    Number of degrees of freedom.
    """
    NL::Int

    """
    Length of imaginary time axis.
    """
    L::Int

    """
    Number of sites in physical lattice.
    """
    N::Int

    """
    Size of lattice in direction of third lattice vector.
    """
    L₃::Int

    """
    Size of lattice in direction of second lattice vector.
    """
    L₂::Int

    """
    Size of lattice in direction of first lattice vector.
    """
    L₁::Int

    """
    Number of sites/orbitals per unit cell.
    """
    nₛ::Int

    """
    All the random vectors.
    """
    R::Matrix{T}

    """
    All solutions to linear systems associated with random vectors.
    """
    M⁻¹R::Matrix{T}

    """
    First random vector of length NL.
    """
    r₁::Vector{T}

    """
    Second random Vector of length NL.
    """
    r₂::Vector{T}

    """
    M⁻¹⋅r₁ of length NL.
    """
    M⁻¹r₁::Vector{T}

    """
    M⁻¹⋅r₁ of length NL.
    """
    M⁻¹r₂::Vector{T}

    """
    Time ordered Green's function G[Δ,0]=Gᵣ(τ)=⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ(0)⟩ where Δ=(τ,r)
    """
    GΔ0::Array{Complex{T},6}

    """
    G[Δ,Δ]⋅G[0,0]=⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ₊ᵣ(τ)⟩⋅⟨T⋅cᵢ(0)⋅cᵀᵢ(0)⟩ where Δ=(τ,r)
    """
    GΔΔ_G00::Array{Complex{T},6}

    """
    G[Δ,0]⋅G[Δ,0]=⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ(0)⟩⋅⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ(0)⟩ where Δ=(τ,r)
    """
    GΔ0_GΔ0::Array{Complex{T},6}

    """
    G[Δ,0]⋅G[0,Δ] =⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ(0)⟩⋅⟨T⋅cᵢ(0)⋅cᵀᵢ₊ᵣ(τ)⟩  where Δ=(τ,r)
    G[Δ,0]⋅G[-Δ,0]=⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ(0)⟩⋅⟨T⋅cᵢ₋ᵣ(-τ)⋅cᵀᵢ(0)⟩ where Δ=(τ,r)
    """
    GΔ0_G0Δ::Array{Complex{T},6}

    """
    Forward FFT plan.
    """
    pfft::Tfft

    """
    Inverse FFT plan.
    """
    pifft::Tifft

    """
    Array of dimension [2L,n,L₁,L₂,L₃] used for performing convolution.
    """
    a::Array{Complex{T},5}

    """
    Array of dimension [2L,n,L₁,L₂,L₃] used for preforming convolution.
    """
    b::Array{Complex{T},5}

    """
    Array of dimension [2L,n,n,L₁,L₂,L₃] used for preforming convolution.
    """
    ab′::Array{Complex{T},6}

    """
    Array of dimension [2L,n,n,L₁,L₂,L₃] used for preforming convolution.
    """
    ab″::Array{Complex{T},6}

    """
    Temporary storage array of dimension [2L,n,L₁,L₂,L₃].
    """
    z::Array{Complex{T},5}

    """
    Constructor for GreensFunction.
    """
    function EstimateGreensFunction(model::AbstractModel{T},nᵥ::Int=2) where {T<:AbstractFloat}

        nᵥ      = max(2,nᵥ)
        NL      = model.Ndim
        N       = model.Nsites
        L       = model.Lτ
        L₃      = model.lattice.L3
        L₂      = model.lattice.L2
        L₁      = model.lattice.L1
        nₛ      = model.lattice.unit_cell.norbits
        R       = zeros(T,NL,nᵥ)
        M⁻¹R    = zeros(T,NL,nᵥ)
        r₁      = zeros(T,NL)
        r₂      = zeros(T,NL)
        M⁻¹r₁   = zeros(T,NL)
        M⁻¹r₂   = zeros(T,NL)
        GΔ0     = zeros(Complex{T},2L,nₛ,nₛ,L₁,L₂,L₃)
        GΔΔ_G00 = zeros(Complex{T},2L,nₛ,nₛ,L₁,L₂,L₃)
        GΔ0_GΔ0 = zeros(Complex{T},2L,nₛ,nₛ,L₁,L₂,L₃)
        GΔ0_G0Δ = zeros(Complex{T},2L,nₛ,nₛ,L₁,L₂,L₃)
        a       = zeros(Complex{T},2L,nₛ,L₁,L₂,L₃)
        b       = zeros(Complex{T},2L,nₛ,L₁,L₂,L₃)
        ab′     = zeros(Complex{T},2L,nₛ,nₛ,L₁,L₂,L₃)
        ab″     = zeros(Complex{T},2L,nₛ,nₛ,L₁,L₂,L₃)
        z       = zeros(Complex{T},2L,nₛ,L₁,L₂,L₃)
        pfft    = plan_fft(a, (1,3,4,5), flags=FFTW.PATIENT)
        pifft   = plan_ifft(ab′, (1,4,5,6), flags=FFTW.PATIENT)
        Tfft    = typeof(pfft)
        Tifft   = typeof(pifft)

        return new{T,Tfft,Tifft}(nᵥ,1,2,NL,L,N,L₃,L₂,L₁,nₛ,R,M⁻¹R,r₁,r₂,M⁻¹r₁,M⁻¹r₂,GΔ0,GΔΔ_G00,GΔ0_GΔ0,GΔ0_G0Δ,pfft,pifft,a,b,ab′,ab″,z)
    end
end

"""
Update Green's functions with new random vectors.
"""
function update!(estimator::EstimateGreensFunction, model::T, preconditioner=I) where {T<:AbstractModel}

    @unpack r₁, M⁻¹r₁, R, M⁻¹R = estimator

    # update preconditioner
    KPMPreconditioners.setup!(preconditioner)

    # iterate over number or random vectors
    for i in 1:estimator.nᵥ

        # initialize random vector with new random number
        randn!(r₁)
        fill!(M⁻¹r₁,0.0)

        # solve linear system
        if model.mul_by_M
            model.transposed = false
            # solve M⋅x=r₁ ==> x=M⁻¹⋅r₁
            iters, err, flag = ldiv!(M⁻¹r₁, model, r₁, preconditioner)
        else
            model.transposed = false
            # solve MᵀM⋅x=Mᵀr₁ ==> x=[MᵀM]⁻¹⋅Mᵀr₁=M⁻¹r₁
            Mᵀr₁ = model.v″
            mulMᵀ!(Mᵀr₁, model, r₁)
            iters, err, flag = ldiv!(M⁻¹r₁, model, Mᵀr₁, preconditioner)
        end

        # store result
        @views @. R[:,i]    = r₁
        @views @. M⁻¹R[:,i] = M⁻¹r₁
    end

    return nothing
end

"""
Setup Green's function convolutions with specified random vectors.
"""
function setup!(estimator::EstimateGreensFunction, n₁::Int, n₂::Int) where {T<:AbstractModel}

    # unpack variables
    @unpack nᵥ, r₁, r₂, M⁻¹r₁, M⁻¹r₂, R, M⁻¹R, a, b, L = estimator

    # assign current random vectors and solutions
    estimator.n₁    = n₁
    estimator.n₂    = n₂
    @views @. M⁻¹r₁ = M⁻¹R[:,n₁]
    @views @. r₁    = R[:,n₁]
    @views @. M⁻¹r₂ = M⁻¹R[:,n₂]
    @views @. r₂    = R[:,n₂]

    # initialize measured values to zero
    fill!(estimator.GΔ0,0.0)
    fill!(estimator.GΔ0_GΔ0,0.0)
    fill!(estimator.GΔ0_G0Δ,0.0)
    fill!(estimator.GΔΔ_G00,0.0)

    # vector to be convolved together
    a = estimator.a
    b = estimator.b

    # calcualte G[Δ,0]
    z = estimator.z
    antiperiodic_copy!(a,M⁻¹r₁,L)
    antiperiodic_copy!(z,M⁻¹r₂,L)
    @. a = (a+z)/sqrt(2.0)
    antiperiodic_copy!(b,r₁,L)
    antiperiodic_copy!(z,r₂,L)
    @. b = (b+z)/sqrt(2.0)
    convolve!(estimator.GΔ0,a,b,estimator)

    # calculate G[Δ,0]⋅G[Δ,0]
    periodic_product!(a,M⁻¹r₁,M⁻¹r₂,L)
    periodic_product!(b,r₁,r₂,L)
    convolve!(estimator.GΔ0_GΔ0,a,b,estimator)

    # calculate G[Δ,Δ]⋅G[0,0]
    periodic_product!(a,M⁻¹r₂,r₂,L)
    periodic_product!(b,M⁻¹r₁,r₁,L)
    convolve!(estimator.GΔΔ_G00,a,b,estimator)

    # calculate G[Δ,0]⋅G[0,Δ]
    periodic_product!(a,M⁻¹r₁,r₂,L)
    periodic_product!(b,M⁻¹r₂,r₁,L)
    convolve!(estimator.GΔ0_G0Δ,a,b,estimator)

    return nothing
end

"""
Measure time ordered Green's function G[Δ,0]=Gᵣ(τ)=⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ(0)⟩ where Δ=(τ,r)
"""
function measure_GΔ0(estimator::EstimateGreensFunction,l₁::Int,l₂::Int,l₃::Int,o₁::Int,o₂::Int,τ::Int)

    L   = estimator.L
    GΔ0 = estimator.GΔ0
    return GΔ0[mod1(τ+1,2L),o₂,o₁,l₁+1,l₂+1,l₃+1]
end

"""
Measure G[Δ,0]⋅G[Δ,0]=⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ(0)⟩⋅⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ(0)⟩ where Δ=(τ,r)
"""
function measure_GΔ0_GΔ0(estimator::EstimateGreensFunction,l₁::Int,l₂::Int,l₃::Int,o₁::Int,o₂::Int,τ::Int)

    L       = estimator.L
    GΔ0_GΔ0 = estimator.GΔ0_GΔ0
    return GΔ0_GΔ0[mod1(τ+1,2L),o₂,o₁,l₁+1,l₂+1,l₃+1]
end

"""
Measure G[Δ,Δ]⋅G[0,0]=⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ₊ᵣ(τ)⟩⋅⟨T⋅cᵢ(0)⋅cᵀᵢ(0)⟩ where Δ=(τ,r)
"""
function measure_GΔΔ_G00(estimator::EstimateGreensFunction,l₁::Int,l₂::Int,l₃::Int,o₁::Int,o₂::Int,τ::Int)

    L       = estimator.L
    GΔΔ_G00 = estimator.GΔΔ_G00
    return GΔΔ_G00[mod1(τ+1,2L),o₂,o₁,l₁+1,l₂+1,l₃+1]
end

"""
Measure G[Δ,0]⋅G[0,Δ] =⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ(0)⟩⋅⟨T⋅cᵢ(0)⋅cᵀᵢ₊ᵣ(τ)⟩  where Δ=(τ,r)
        G[Δ,0]⋅G[-Δ,0]=⟨T⋅cᵢ₊ᵣ(τ)⋅cᵀᵢ(0)⟩⋅⟨T⋅cᵢ₋ᵣ(-τ)⋅cᵀᵢ(0)⟩ where Δ=(τ,r)
"""
function measure_GΔ0_G0Δ(estimator::EstimateGreensFunction,l₁::Int,l₂::Int,l₃::Int,o₁::Int,o₂::Int,τ::Int)

    L       = estimator.L
    GΔ0_G0Δ = estimator.GΔ0_G0Δ
    return GΔ0_G0Δ[mod1(τ+1,2L),o₂,o₁,l₁+1,l₂+1,l₃+1]
end

"""
Measure time ordered Green's function Gᵢ-ⱼ(τ₂-τ₁)=⟨T⋅cᵢ(τ₂)⋅cᵀⱼ(τ₁)⟩ where -β<(τ₂-τ₁)<β
"""
function estimate(estimator::EstimateGreensFunction,i::Int,j::Int,τ₂::Int,τ₁::Int,σ::Int)

    m = get_index(τ₁,j,estimator.L)
    n = get_index(τ₂,i,estimator.L)
    if σ==1
        Gᵢⱼτ₁τ₂ =  estimator.M⁻¹r₁[n] * estimator.r₁[m]
    elseif σ==2
        Gᵢⱼτ₁τ₂ =  estimator.M⁻¹r₂[n] * estimator.r₂[m]
    else
        throw(DomainError())
    end
    return Gᵢⱼτ₁τ₂
end

"""
Calculate convolution a⋆b.
"""
function convolve!(ab::AbstractArray,a::AbstractArray,b::AbstractArray,estimator::EstimateGreensFunction)

    b′  = estimator.z
    ab′ = estimator.ab′
    ab″ = estimator.ab″
    L   = estimator.L
    N   = estimator.N
    L₃  = estimator.L₃
    L₂  = estimator.L₂
    L₁  = estimator.L₁
    nₛ  = estimator.nₛ

    # forward FFT of a
    a′ = a
    mul!(b′,estimator.pfft,a)
    copyto!(a′,b′)

    # forward FFT of b
    mul!(b′,estimator.pfft,b)

    # normalization factor
    V = 2*L*N/nₛ

    # perform elementwise multiplication
    @fastmath @inbounds for k₃ in 1:L₃
        for k₂ in 1:L₂
            for k₁ in 1:L₁
                for s₁ in 1:nₛ
                    for s₂ in 1:nₛ
                        for ω in 1:2L
                            nk₃ = mod1(-k₃+2,L₃)
                            nk₂ = mod1(-k₂+2,L₂)
                            nk₁ = mod1(-k₁+2,L₁)
                            nω  = mod1(-ω+2,2L)
                            ab′[ω,s₂,s₁,k₁,k₂,k₃] = a′[ω,s₂,k₁,k₂,k₃] * b′[nω,s₁,nk₁,nk₂,nk₃] / V
                        end
                    end
                end
            end
        end
    end

    # perform inverse FFT
    mul!(ab″,estimator.pifft,ab′)

    # add result to output vector
    @. ab += ab″

    return nothing
end

"""
Copy so the vector x=[ x(1) , ... , x(τ) , ... , x(L) ] to a vector y so that
y=[ x(1) , ... , x(τ) , ... , x(L) , -x(1) , ... , -x(τ) , ... , -x(L) ]
"""
function antiperiodic_copy!(y::AbstractArray,x::AbstractArray,L::Int)

    N  = div(length(x),L)
    x′ = reshaped(x,L,N)
    y′ = reshaped(y,2L,N)
    @fastmath @inbounds for i in 1:N
        for τ in 1:L
            y′[τ,i]   =  x′[τ,i]
            y′[τ+L,i] = -x′[τ,i]
        end
    end
    return nothing
end

"""
Multiply so that for x=[ x(1) , ... , x(τ) , ... , x(L) ] and y=[ y(1) , ... , y(τ) , ... , y(L) ] then
z=[ x(1)⋅y(1) , x(2)⋅y(2) , ... , x(L)⋅y(L) , x(1)⋅y(1) , ... , x(L)⋅y(L) , x(L)⋅y(L) ]
"""
function periodic_product!(z::AbstractArray,y::AbstractArray,x::AbstractArray,L::Int)

    N  = div(length(x),L)
    x′ = reshaped(x,L,N)
    y′ = reshaped(y,L,N)
    z′ = reshaped(z,2L,N)
    @fastmath @inbounds for i in 1:N
        for τ in 1:L
            val       = y′[τ,i] * x′[τ,i]
            z′[τ,i]   = val
            z′[τ+L,i] = val
        end
    end

    return nothing
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
"""
    ReducedSpaceEvaluator{T} <: AbstractNLPEvaluator

Evaluator working in the reduced space corresponding to the
control variable `u`. Once a new point `u` is passed to the evaluator,
the user needs to call the method `update!` to find the corresponding
state `x(u)` satisfying the equilibrium equation `g(x(u), u) = 0`.

Taking as input a given `AbstractFormulation`, the reduced evaluator
builds the bounds corresponding to the control `u` and the state `x`,
and initiate an `AutoDiffFactory` tailored to the problem. The reduced evaluator
could be instantiated on the main memory, or on a specific device (currently,
only CUDA is supported).

"""
mutable struct ReducedSpaceEvaluator{T} <: AbstractNLPEvaluator
    model::AbstractFormulation
    x::AbstractVector{T}
    p::AbstractVector{T}
    λ::AbstractVector{T}

    x_min::AbstractVector{T}
    x_max::AbstractVector{T}
    u_min::AbstractVector{T}
    u_max::AbstractVector{T}

    constraints::Array{Function, 1}
    g_min::AbstractVector{T}
    g_max::AbstractVector{T}

    buffer::AbstractNetworkBuffer
    autodiff::AutoDiffFactory
    ∇gᵗ::AbstractMatrix
    linear_solver::LinearSolvers.AbstractLinearSolver
    ε_tol::Float64
end

function ReducedSpaceEvaluator(
    model, x, u, p;
    constraints=Function[state_constraint, power_constraints],
    linear_solver=DirectSolver(),
    ε_tol=1e-12,
    verbose_level=VERBOSE_LEVEL_NONE,
)
    # First, build up a network buffer
    buffer = get(model, PhysicalState())
    # Initiate adjoint
    λ = similar(x)
    # Build up AutoDiff factory
    jx, ju, adjoint_f = init_autodiff_factory(model, buffer)
    ad = AutoDiffFactory(jx, ju, adjoint_f)

    u_min, u_max = bounds(model, Control())
    x_min, x_max = bounds(model, State())

    MT = model.AT
    g_min = MT{eltype(x), 1}()
    g_max = MT{eltype(x), 1}()
    for cons in constraints
        cb, cu = bounds(model, cons)
        append!(g_min, cb)
        append!(g_max, cu)
    end

    return ReducedSpaceEvaluator(model, x, p, λ, x_min, x_max, u_min, u_max,
                                 constraints, g_min, g_max,
                                 buffer,
                                 ad, jx.J, linear_solver, ε_tol)
end
function ReducedSpaceEvaluator(
    datafile;
    device=CPU(),
    options...
)
    # Load problem.
    pf = PS.PowerNetwork(datafile)
    polar = PolarForm(pf, device)

    x0 = initial(polar, State())
    p = initial(polar, Parameters())
    uk = initial(polar, Control())

    return ReducedSpaceEvaluator(
        polar, x0, uk, p; options...
    )
end

n_variables(nlp::ReducedSpaceEvaluator) = length(nlp.u_min)
n_constraints(nlp::ReducedSpaceEvaluator) = length(nlp.g_min)

initial(nlp::ReducedSpaceEvaluator) = initial(nlp.model, Control())
bounds(nlp::ReducedSpaceEvaluator, ::Variables) = nlp.u_min, nlp.u_max
bounds(nlp::ReducedSpaceEvaluator, ::Constraints) = nlp.g_min, nlp.g_max

function update!(nlp::ReducedSpaceEvaluator, u; verbose_level=0)
    x₀ = nlp.x
    jac_x = nlp.autodiff.Jgₓ
    # Transfer x, u, p into the network cache
    transfer!(nlp.model, nlp.buffer, nlp.x, u, nlp.p)
    # Get corresponding point on the manifold
    conv = powerflow(nlp.model, jac_x, nlp.buffer, tol=nlp.ε_tol;
                     solver=nlp.linear_solver, verbose_level=verbose_level)
    if !conv.has_converged
        @warn("Newton-Raphson algorithm failed to converge ($(conv.norm_residuals))")
        return conv
    end
    ∇gₓ = nlp.autodiff.Jgₓ.J
    nlp.∇gᵗ = LinearSolvers.get_transpose(nlp.linear_solver, ∇gₓ)
    # Switch preconditioner to transpose mode
    if isa(nlp.linear_solver, LinearSolvers.AbstractIterativeLinearSolver)
        LinearSolvers.update!(nlp.linear_solver, nlp.∇gᵗ)
    end

    # Update value of nlp.x with new network state
    get!(nlp.model, State(), nlp.x, nlp.buffer)
    # Refresh value of the active power of the generators
    refresh!(nlp.model, PS.Generator(), PS.ActivePower(), nlp.buffer)
    return conv
end

function objective(nlp::ReducedSpaceEvaluator, u)
    # Take as input the current cache, updated previously in `update!`.
    cost = cost_production(nlp.model, nlp.buffer.pg)
    # TODO: determine if we should include λ' * g(x, u), even if ≈ 0
    return cost
end

# compute inplace reduced gradient (g = ∇fᵤ + (∇gᵤ')*λₖ)
# equivalent to: g = ∇fᵤ - (∇gᵤ')*λₖ_neg
# (take λₖ_neg to avoid computing an intermediate array)
function _reduced_gradient!(g, ∇fᵤ, ∇gᵤ, λₖ_neg)
    g .= ∇fᵤ
    mul!(g, transpose(∇gᵤ), λₖ_neg, -1.0, 1.0)
end

function gradient!(nlp::ReducedSpaceEvaluator, g, u)
    buffer = nlp.buffer
    xₖ = nlp.x
    ∇gₓ = nlp.autodiff.Jgₓ.J
    # Evaluate Jacobian of power flow equation on current u
    ∇gᵤ = jacobian(nlp.model, nlp.autodiff.Jgᵤ, buffer)
    # Evaluate adjoint of cost function and update inplace AdjointStackObjective
    ∂cost(nlp.model, nlp.autodiff.∇f, buffer)

    ∇fₓ, ∇fᵤ = nlp.autodiff.∇f.∇fₓ, nlp.autodiff.∇f.∇fᵤ
    # Update (negative) adjoint
    λₖ_neg = nlp.λ
    LinearSolvers.ldiv!(nlp.linear_solver, λₖ_neg, nlp.∇gᵗ, ∇fₓ)
    _reduced_gradient!(g, ∇fᵤ, ∇gᵤ, λₖ_neg)
    return nothing
end

function constraint!(nlp::ReducedSpaceEvaluator, g, u)
    xₖ = nlp.x
    ϕ = nlp.buffer
    # First: state constraint
    mf = 1
    mt = 0
    for cons in nlp.constraints
        m_ = size_constraint(nlp.model, cons)
        mt += m_
        cons_ = @view(g[mf:mt])
        cons(nlp.model, cons_, ϕ)
        mf += m_
    end
end

function jacobian_structure!(nlp::ReducedSpaceEvaluator, rows, cols)
    m, n = n_constraints(nlp), n_variables(nlp)
    idx = 1
    for c in 1:m #number of constraints
        for i in 1:n # number of variables
            rows[idx] = c ; cols[idx] = i
            idx += 1
        end
    end
end

function jacobian!(nlp::ReducedSpaceEvaluator, jac, u)
    model = nlp.model
    xₖ = nlp.x
    ∇gₓ = nlp.autodiff.Jgₓ.J
    ∇gᵤ = nlp.autodiff.Jgᵤ.J
    nₓ = length(xₖ)
    MT = nlp.model.AT
    μ = similar(nlp.λ)
    ∂obj = nlp.autodiff.∇f
    cnt = 1

    for cons in nlp.constraints
        mc_ = size_constraint(nlp.model, cons)
        for i_cons in 1:mc_
            jacobian(model, cons, i_cons, ∂obj, nlp.buffer)
            jx, ju = ∂obj.∇fₓ, ∂obj.∇fᵤ
            # Get adjoint
            LinearSolvers.ldiv!(nlp.linear_solver, μ, nlp.∇gᵗ, jx)
            jac[cnt, :] .= (ju .- ∇gᵤ' * μ)
            cnt += 1
        end
    end
end

function jtprod!(nlp::ReducedSpaceEvaluator, cons, jv, u, v; start=1)
    model = nlp.model
    xₖ = nlp.x
    ∇gₓ = nlp.autodiff.Jgₓ.J
    ∇gᵤ = nlp.autodiff.Jgᵤ.J
    nₓ = length(xₖ)
    μ = nlp.λ

    ∂obj = nlp.autodiff.∇f
    # Get adjoint
    jtprod(model, cons, ∂obj, nlp.buffer, v)
    jvx, jvu = ∂obj.∇fₓ, ∂obj.∇fᵤ
    # jv .+= (ju .- ∇gᵤ' * μ)
    LinearSolvers.ldiv!(nlp.linear_solver, μ, nlp.∇gᵗ, jvx)
    jv .+= jvu
    mul!(jv, transpose(∇gᵤ), μ, -1.0, 1.0)
end
function jtprod!(nlp::ReducedSpaceEvaluator, jv, u, v)
    μ = nlp.λ
    ∇gᵤ = nlp.autodiff.Jgᵤ.J
    ∂obj = nlp.autodiff.∇f
    jvx = ∂obj.jvₓ
    jvu = ∂obj.jvᵤ
    fill!(jvx, 0)
    fill!(jvu, 0)

    fr_ = 0
    for cons in nlp.constraints
        n = size_constraint(nlp.model, cons)
        mask = fr_+1:fr_+n
        vv = @view v[mask]
        # Compute jtprod of current constraint
        jtprod(nlp.model, cons, ∂obj, nlp.buffer, vv)
        jvx .+= ∂obj.∇fₓ
        jvu .+= ∂obj.∇fᵤ
        fr_ += n
    end
    LinearSolvers.ldiv!(nlp.linear_solver, μ, nlp.∇gᵗ, jvx)
    jv .+= jvu
    mul!(jv, transpose(∇gᵤ), μ, -1.0, 1.0)
    return
end

# Utils function
function primal_infeasibility!(nlp::ReducedSpaceEvaluator, cons, u)
    constraint!(nlp, cons, u) # Evaluate constraints
    (n_inf, err_inf, n_sup, err_sup) = _check(cons, nlp.g_min, nlp.g_max)
    return max(err_inf, err_sup)
end
function primal_infeasibility(nlp::ReducedSpaceEvaluator, u)
    cons = similar(nlp.g_min) ; fill!(cons, 0)
    return primal_infeasibility!(nlp, cons, u)
end

# Printing
function sanity_check(nlp::ReducedSpaceEvaluator, u, cons)
    println("Check violation of constraints")
    print("Control  \t")
    (n_inf, err_inf, n_sup, err_sup) = _check(u, nlp.u_min, nlp.u_max)
    @printf("UB: %.4e (%d)    LB: %.4e (%d)\n",
            err_sup, n_sup, err_inf, n_inf)
    print("State    \t")
    (n_inf, err_inf, n_sup, err_sup) = _check(nlp.x, nlp.x_min, nlp.x_max)
    @printf("UB: %.4e (%d)    LB: %.4e (%d)\n",
            err_sup, n_sup, err_inf, n_inf)
    print("Constraints\t")
    (n_inf, err_inf, n_sup, err_sup) = _check(cons, nlp.g_min, nlp.g_max)
    @printf("UB: %.4e (%d)    LB: %.4e (%d)\n",
            err_sup, n_sup, err_inf, n_inf)
end

function Base.show(io::IO, nlp::ReducedSpaceEvaluator)
    n = n_variables(nlp)
    m = n_constraints(nlp)
    println(io, "A ReducedSpaceEvaluator object")
    println(io, "    * device: ", nlp.model.device)
    println(io, "    * #vars: ", n)
    println(io, "    * #cons: ", m)
    println(io, "    * constraints:")
    for cons in nlp.constraints
        println(io, "        - ", cons)
    end
    print(io, "    * linear solver: ", nlp.linear_solver)
end

function reset!(nlp::ReducedSpaceEvaluator)
    # Reset adjoint
    fill!(nlp.λ, 0)
    # Reset initial state
    x0 = initial(nlp.model, State())
    copy!(nlp.x, x0)
    # Reset buffer
    nlp.buffer = get(nlp.model, PhysicalState())
end


text:
"""
requireents \\

For calculating the required entropies  \\
Arguments:  probs, dictionary of probabilities  \\
            calc, the required entropies  \\
Returns:    ents, dictionary calc => entropy  \\
"""

function requireents(probs::Dict{String, Array{Float64}}, calc::Array{String, 1}; base::Int=2)
  ents = Dict{String, Float64}()
  in("HX", calc) && (ents["HX"] = entropy(probs["PX"], base))
  in("HY", calc) && (ents["HY"] = entropy(probs["PY"], base))
  in("HXY", calc) && (ents["HXY"] = condentropy(probs["PY"], probs["PXY"], base=base))
  return ents
end


text:
include("./melspectrogram.jl")
using FFTW

"""
    mfcc(signal, sr, window=ones, windowsize=2048, overlap=512, numcep=26)

Compute the Mel-Frequency Cepstral Coefficients of the given signal and return
an Array where each column contains the coefficients of a given frame obtained
from the stft.

**Parameters:**

`signal`: The signal on which the coefficients will be computed.

`sr`: Sample rate of the signal.

`window`: The window function to be applied on the signal by the STFT.

`windowsize`: The size of the STFT window.

`overlap`: The overlap of two adjacent windows of the STFT in samples.

`numcep`: The amount of cepstral coefficients to be extracted from the signal.
"""
function mfcc(
	signal, sr, window=ones, windowsize=2048, overlap=512, numcep=26
)

	# Evaluate the melspectrogram
	melspec = melspectrogram(signal, sr, window, windowsize, overlap, numcep)

	#Take the Type-III DCT of the log for each coeficient
	return mapslices(frame->idct(log.(frame)), melspec .+ eps(Float64), dims=1)

end

max_stars_repo_name: gabriel-ss/music-autotagger

max_issues_repo_name: gabriel-ss/music-autotagger

max_forks_repo_name: gabriel-ss/music-autotagger


text:
import MaximumLikelihoodProblems

import CategoricalArrays
import DataFrames
import Econometrics
import Statistics
import Test

β_true = [1.0 2.0 3.0; 4.0 5.0 6.0]

# @show β_true
# @show β_hat

Test.@test typeof(β_hat) == typeof(β_true)
Test.@test ndims(β_hat) == ndims(β_true)
Test.@test size(β_hat) == size(β_true)
absolute_error = abs.(β_hat - β_true)
square_error = abs2.(β_hat - β_true)
absolute_error_proportional = absolute_error ./ abs.(β_true)
square_error_proportional = square_error ./ abs.(β_true)
Test.@test sum(absolute_error) < 0.60
Test.@test sum(square_error) < 0.06
Test.@test sum(absolute_error_proportional) < 0.35
Test.@test sum(square_error_proportional) < 0.03
Test.@test maximum(absolute_error) < 0.30
Test.@test maximum(square_error) < 0.05
Test.@test maximum(absolute_error_proportional) < 0.20
Test.@test maximum(square_error_proportional) < 0.03
Test.@test Statistics.mean(absolute_error) < 0.95
Test.@test Statistics.mean(square_error) < 0.0095
Test.@test Statistics.mean(absolute_error_proportional) < 0.05
Test.@test Statistics.mean(square_error_proportional) < 0.0055

gradient_vector_at_θ_hat = MaximumLikelihoodProblems.Internal.gradient_vector(transformed_gradient_problem,
                                                                              θ_hat)
hessian_matrix_at_θ_hat = MaximumLikelihoodProblems.Internal.hessian_matrix(transformed_gradient_problem,
                                                                            θ_hat)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_zero(gradient_vector_at_θ_hat)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_hermitian(hessian_matrix_at_θ_hat)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_negative_definite(hessian_matrix_at_θ_hat;
                                                                                fuzz_factor = 10)
Test.@test MaximumLikelihoodProblems.Internal._is_approximately_positive_definite(-hessian_matrix_at_θ_hat;
                                                                                fuzz_factor = 10)

external_df = DataFrames.DataFrame()
external_df[!, :X_2] = X[:, 2]
n = size(X, 1)
external_y_string = Vector{String}(undef, n)
for i = 1:n
    external_y_string[i] = string(argmax(y[i, :]))
end
external_y_string_categorical = CategoricalArrays.CategoricalArray(external_y_string;
                                                                   ordered = false)
external_df[!, :y] = external_y_string_categorical
external_formula = Econometrics.@formula(y ~ 1 + X_2)
# external_model_econometrics = Econometrics.fit(Econometrics.EconometricModel,
                                               # external_formula,
                                               # external_df)
# external_model_coef = Econometrics.coef(external_model_econometrics)
# beta_hat_external_model_econometrics = reshape(external_model_coef, (2, 3))
# Test.@test all(isapprox.(β_hat, beta_hat_external_model_econometrics; atol = 1e-2))
# external_model_absolute_error = abs.(β_hat - beta_hat_external_model_econometrics)
# external_model_square_error = abs2.(β_hat - beta_hat_external_model_econometrics)
# external_model_absolute_error_proportional = external_model_absolute_error ./ abs.(beta_hat_external_model_econometrics)
# external_model_square_error_proportional = external_model_square_error ./ abs.(beta_hat_external_model_econometrics)
# Test.@test sum(external_model_absolute_error) < 0.07
# Test.@test sum(external_model_square_error) < 0.001
# Test.@test sum(external_model_absolute_error_proportional) < 0.04
# Test.@test sum(external_model_square_error_proportional) < 1e-4
# Test.@test maximum(external_model_absolute_error) < 0.01
# Test.@test maximum(external_model_square_error) < 1e-4
# Test.@test maximum(absolute_error_proportional) < 0.1
# Test.@test maximum(external_model_square_error_proportional) < 1e-4
# Test.@test Statistics.mean(external_model_absolute_error) < 0.01
# Test.@test Statistics.mean(external_model_square_error) < 1e-4
# Test.@test Statistics.mean(external_model_absolute_error_proportional) < 0.01
# Test.@test Statistics.mean(external_model_square_error_proportional) < 1e-4

# coverage for the "failed to converge" code paths
β_hat_initial_guess = zeros(size_β)
θ_hat_initial = (; β = β_hat_initial_guess)
θ_hat = MaximumLikelihoodProblems.fit(transformed_gradient_problem,
                                      θ_hat_initial;
                                      max_iterations = 10,
                                      throw_convergence_exception = false)
@info("The previous warning message (\"Warning: Failed to converge after 10 iterations\") was expected. It is a normal part of the test suite.")
Test.@test_throws(MaximumLikelihoodProblems.ConvergenceException,
                  MaximumLikelihoodProblems.fit(transformed_gradient_problem,
                                                θ_hat_initial;
                                                max_iterations = 10,
                                                throw_convergence_exception = true))
@info("The previous error message (\"Error: Failed to converge after 10 iterations\") was expected. It is a normal part of the test suite.")
Test.@test_throws(MaximumLikelihoodProblems.ConvergenceException,
                  MaximumLikelihoodProblems.fit(transformed_gradient_problem,
                                                θ_hat_initial;
                                                max_iterations = 10))
@info("The previous error message (\"Error: Failed to converge after 10 iterations\") was expected. It is a normal part of the test suite.")


text:
# left.jl
# "special cases" related left vector multiplication like x = y'*A
# The subtlety is that y' is a AdjointAbsVec or a TransposeAbsVec
# which is a subtype of AbstractMatrix but it is really "vector like"
# so we want handle it like a vector (Issue#99)
# So this is an exception to the left multiplication rule by a AbstractMatrix
# that usually makes a WrappedMap.

# x = y'*A ⇐⇒ x' = (A'*y)
"""
    *(x::LinearAlgebra.AdjointAbsVec, A::LinearMap)::AdjointAbsVec

Compute the right-action of the linear map `A` on the adjoint vector `x`
and return an adjoint vector.

## Examples
```jldoctest; setup=(using LinearAlgebra, LinearMaps)
julia> A=LinearMap([1.0 2.0; 3.0 4.0]); x=[1.0, 1.0]; x'A
1×2 Adjoint{Float64,Array{Float64,1}}:
 4.0  6.0
```
"""
Base.:(*)(y::LinearAlgebra.AdjointAbsVec, A::LinearMap) = adjoint(A' * y')

"""
    *(x::LinearAlgebra.TransposeAbsVec, A::LinearMap)::TransposeAbsVec

Compute the right-action of the linear map `A` on the transpose vector `x`
and return a transpose vector.

## Examples
```jldoctest; setup=(using LinearAlgebra, LinearMaps)
julia> A=LinearMap([1.0 2.0; 3.0 4.0]); x=[1.0, 1.0]; transpose(x)*A
1×2 Transpose{Float64,Array{Float64,1}}:
 4.0  6.0
```
"""
Base.:(*)(y::LinearAlgebra.TransposeAbsVec, A::LinearMap) = transpose(transpose(A) * transpose(y))

# multiplication with vector/matrix
const AdjointAbsVecOrMat{T} = Adjoint{T,<:AbstractVecOrMat}
const TransposeAbsVecOrMat{T} = Transpose{T,<:AbstractVecOrMat}

function mul!(x::AbstractMatrix, y::AdjointAbsVecOrMat, A::LinearMap)
    check_dim_mul(x, y, A)
    _unsafe_mul!(x', A', y')
    return x
end

function mul!(x::AbstractMatrix, y::AdjointAbsVecOrMat, A::LinearMap, α::Number, β::Number)
    check_dim_mul(x, y, A)
    _unsafe_mul!(x', conj(α)*A', y', true, conj(β))
    return x
end

function mul!(x::AbstractMatrix, y::TransposeAbsVecOrMat, A::LinearMap)
    check_dim_mul(x, y, A)
    _unsafe_mul!(transpose(x), transpose(A), transpose(y))
    return x
end

function mul!(x::AbstractMatrix, y::TransposeAbsVecOrMat, A::LinearMap, α::Number, β::Number)
    check_dim_mul(x, y, A)
    _unsafe_mul!(transpose(x), α*transpose(A), transpose(y), true, β)
    return x
end


text:
mutable struct FedDCDServer{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Vector{Int64}, T5<:Vector{FedDCDClient}, T6<:Matrix{Float64}} <:AbstractServer
    Xtest::T3                       # testing data
    Ytest::T4                       # testing label
    num_classes::T1                 # number of classes
    num_clients::T1                 # number of clients
    participation_rate::T2          # participation rate
    clients::T5                     # set of clients
    W::T6                           # aggregation of uploaded model updates
    τ::T1                           # number of particiating clients
    selectedIndices::T4             # selected clients
    η::T2                           # learning rate
    function FedDCDServer(Xtest::SparseMatrixCSC{Float64, Int64}, Ytest::Vector{Int64}, config::Dict{String, Real})
        num_classes = config["num_classes"]
        num_clients = config["num_clients"]
        participation_rate = config["participation_rate"]
        η = config["learning_rate"]
        τ = floor(Int64, num_clients * participation_rate)
        _, d = size(Xtest)
        clients = Vector{FedDCDClient}(undef, num_clients)
        W = zeros(Float64, d, num_classes)
        selectedIndices = zeros(Int64, τ)
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Vector{Int64}, Vector{FedDCDClient}, Matrix{Float64}}(Xtest, Ytest, num_classes, num_clients, participation_rate, clients, W, τ, selectedIndices, η)
    end
end

# Select clients
function select!(
    server::FedDCDServer
    )
    server.selectedIndices .= randperm(server.num_clients)[1:server.τ]
    return nothing
end

# Send global model
function sendModel!(
    server::FedDCDServer
)
    # Only send model to selected clients
    Threads.@threads for i = 1:server.τ
        idx = server.selectedIndices[i]
        server.clients[idx].y .-= (server.η * server.clients[idx].λ) * (server.clients[idx].W - server.W)
    end
    return nothing
end

# Send global model to all clients
function sendModelToAllClients!(
    server::FedDCDServer
)
    # Only send model to selected clients
    for client in server.clients
        client.W = copy(server.W)
    end
    return nothing
end

# Aggregation
function aggregate!(
    server::FedDCDServer
)
    # Take the average of the selected clients' model
    fill!(server.W, 0.0)
    for i = 1:server.τ
        idx = server.selectedIndices[i]
        server.W += server.clients[idx].W
    end
    server.W ./= server.τ
    return nothing
end

########################################################################################################
mutable struct FedDCDServerNN{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Vector{Int64}, T5<:Vector{FedDCDClientNN}, T6<:Flux.Chain, T7<:Flux.OneHotArray} <:AbstractServer
    Xtest::T3                       # testing data
    XtestT::T3                      # copy of transpose
    Ytest::T7                       # testing label
    num_classes::T1                 # number of classes
    num_clients::T1                 # number of clients
    participation_rate::T2          # participation rate
    clients::T5                     # set of clients
    W::T6                           # aggregation of uploaded model updates
    τ::T1                           # number of particiating clients
    selectedIndices::T4             # selected clients
    η::T2                           # learning rate
    decay::T2                       # decay rate
    function FedDCDServerNN(Xtest::SparseMatrixCSC{Float64, Int64}, Ytest::Vector{Int64}, dim::Int64, config::Dict{String, Real})
        num_classes = config["num_classes"]
        num_clients = config["num_clients"]
        participation_rate = config["participation_rate"]
        η = config["learning_rate"]
        decay = config["decay_rate"]
        W = Chain( Dense(780, dim, relu, bias=false), Dense(dim, 10, bias=false), NNlib.softmax)
        τ = floor(Int64, num_clients * participation_rate)
        XtestT = copy(Xtest')
        Ytest = Flux.onehotbatch(Ytest, 1:num_classes)
        clients = Vector{FedDCDClientNN}(undef, num_clients)
        selectedIndices = zeros(Int64, τ)
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Vector{Int64}, Vector{FedDCDClientNN}, Flux.Chain, Flux.OneHotArray}(Xtest, XtestT, Ytest, num_classes, num_clients, participation_rate, clients, W, τ, selectedIndices, η, decay)
    end
end

# Select clients
function select!(
    server::FedDCDServerNN
    )
    server.selectedIndices .= randperm(server.num_clients)[1:server.τ]
    return nothing
end

# Send global model
function sendModel!(
    server::FedDCDServerNN
)
    # Only send model to selected clients
    Threads.@threads for i = 1:server.τ
        idx = server.selectedIndices[i]
        c = server.clients[idx]
        for j = 1:length(c.y)
            c.y[j] .-= server.η * ( params(c.W)[j] - params(server.W)[j] )
        end
    end
    return nothing
end

# Aggregation
function aggregate!(
    server::FedDCDServerNN
)
    # Take the average of the selected clients' model
    l = length(params(server.W))
    for j = 1:l
        fill!(params(server.W)[j], 0.0)
    end
    for i = 1:server.τ
        idx = server.selectedIndices[i]
        c = server.clients[idx]
        for j = 1:l
            params(server.W)[j] .+= (1/server.τ) * params(c.W)[j]
        end
    end
    return nothing
end

# synchronization
function synchronize!(server::FedDCDServerNN)
    for i in server.selectedIndices
        c = server.clients[i]
        for j = 1:length(params(server.W))
            params(c.W)[j] .= params(server.W)[j]
        end
    end
end

# update learning rate
function decay_lr!(server::FedDCDServerNN)
    server.η *= (1 - server.decay)
end

########################################################################################################
mutable struct AccFedDCDServer{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Vector{Int64}, T5<:Vector{AccFedDCDClient}, T6<:Matrix{Float64}} <:AbstractServer
    Xtest::T3                       # testing data
    Ytest::T4                       # testing label
    num_classes::T1                 # number of classes
    num_clients::T1                 # number of clients
    participation_rate::T2          # participation rate
    clients::T5                     # set of clients
    W::T6                           # aggregation of uploaded model updates
    τ::T1                           # number of particiating clients
    selectedIndices::T4             # selected clients
    η::T2                           # learning rate
    function AccFedDCDServer(Xtest::SparseMatrixCSC{Float64, Int64}, Ytest::Vector{Int64}, config::Dict{String, Real})
        num_classes = config["num_classes"]
        num_clients = config["num_clients"]
        participation_rate = config["participation_rate"]
        η = config["learning_rate"]
        τ = floor(Int64, num_clients * participation_rate)
        _, d = size(Xtest)
        clients = Vector{AccFedDCDClient}(undef, num_clients)
        W = zeros(Float64, d, num_classes)
        selectedIndices = zeros(Int64, τ)
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Vector{Int64}, Vector{AccFedDCDClient}, Matrix{Float64}}(Xtest, Ytest, num_classes, num_clients, participation_rate, clients, W, τ, selectedIndices, η)
    end
end

# Select clients
function select!(server::AccFedDCDServer)
    server.selectedIndices .= randperm(server.num_clients)[1:server.τ]
    return nothing
end

# Send global model
function sendModel!(server::AccFedDCDServer, round::Int64)
    # Only send model to selected clients
    Threads.@threads for i = 1:server.τ
        idx = server.selectedIndices[i]
        if round == 1
            server.clients[idx].y .-= (server.η * server.clients[idx].λ) * (server.clients[idx].W - server.W)
        else
            a = server.clients[idx].a; b = server.clients[idx].b; r = server.clients[idx].r
            θ = a*r / (a^2 + b)
            server.clients[idx].z .-= (server.η * server.clients[idx].λ * θ) * (server.clients[idx].W - server.W)
        end
    end
    return nothing
end

# Aggregation
function aggregate!(server::AccFedDCDServer)
    # Take the average of the selected clients' model
    fill!(server.W, 0.0)
    for i = 1:server.τ
        idx = server.selectedIndices[i]
        server.W += server.clients[idx].W
    end
    server.W ./= server.τ
    return nothing
end

# Calculate objective value
function getObjValue(
    server::Union{FedDCDServer, AccFedDCDServer}
)
    objValue = 0.0
    totalNumData = 0
    for i = 1:server.num_clients
        n = size(server.clients[i].Xtrain, 1)
        totalNumData += n
        objValue += n * obj(
            server.clients[i].Xtrain,
            server.clients[i].Ytrain,
            server.W,
            server.clients[i].λ
         )
        # objValue += n * getObjValue(server.clients[i])
    end
    return objValue / totalNumData
end

########################################################################################################
mutable struct AccFedDCDServerNN{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Vector{Int64}, T5<:Vector{AccFedDCDClientNN}, T6<:Flux.Chain, T7<:Flux.OneHotArray} <:AbstractServer
    Xtest::T3                       # testing data
    XtestT::T3                      # copy of transpose
    Ytest::T7                       # testing label
    num_classes::T1                 # number of classes
    num_clients::T1                 # number of clients
    participation_rate::T2          # participation rate
    clients::T5                     # set of clients
    W::T6                           # aggregation of uploaded model updates
    τ::T1                           # number of particiating clients
    selectedIndices::T4             # selected clients
    η::T2                           # learning rate
    function AccFedDCDServerNN(Xtest::SparseMatrixCSC{Float64, Int64}, Ytest::Vector{Int64}, dim::Int64, config::Dict{String, Real})
        num_classes = config["num_classes"]
        num_clients = config["num_clients"]
        participation_rate = config["participation_rate"]
        η = config["learning_rate"]
        W = Chain( Dense(780, dim, relu, bias=false), Dense(dim, 10, bias=false), NNlib.softmax)
        τ = floor(Int64, num_clients * participation_rate)
        XtestT = copy(Xtest')
        Ytest = Flux.onehotbatch(Ytest, 1:num_classes)
        clients = Vector{AccFedDCDClientNN}(undef, num_clients)
        selectedIndices = zeros(Int64, τ)
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Vector{Int64}, Vector{AccFedDCDClientNN}, Flux.Chain, Flux.OneHotArray}(Xtest, XtestT, Ytest, num_classes, num_clients, participation_rate, clients, W, τ, selectedIndices, η)
    end
end

# Select clients
function select!(server::AccFedDCDServerNN)
    server.selectedIndices .= randperm(server.num_clients)[1:server.τ]
    return nothing
end

# Send global model
function sendModel!(server::AccFedDCDServerNN, round::Int64)
    # Only send model to selected clients
    Threads.@threads for i = 1:server.τ
        idx = server.selectedIndices[i]
        client = server.clients[idx]
        if round == 1
            for j = 1:length(client.y)
                client.y[j] .-= server.η * (params(client.W)[j] - params(server.W)[j])
            end
        else
            a = client.a; b = client.b; r = client.r
            θ = a*r / (a^2 + b)
            for j = 1:length(client.z)
                client.z[j] .-= server.η * θ * (params(client.W)[j] - params(server.W)[j])
            end
        end
    end
    return nothing
end

# synchronization
function synchronize!(server::AccFedDCDServerNN)
    for i in server.selectedIndices
        c = server.clients[i]
        for j = 1:length(params(server.W))
            params(c.W)[j] .= params(server.W)[j]
        end
    end
end

# Aggregation
function aggregate!(
    server::AccFedDCDServerNN
)
    # Take the average of the selected clients' model
    l = length(params(server.W))
    for j = 1:l
        fill!(params(server.W)[j], 0.0)
    end
    for i = 1:server.τ
        idx = server.selectedIndices[i]
        c = server.clients[idx]
        for j = 1:l
            params(server.W)[j] .+= (1/server.τ) * params(c.W)[j]
        end
    end
    return nothing
end

# Calculate objective value
function getObjValue(
    server::Union{FedDCDServerNN, AccFedDCDServerNN}
)
    objValue = 0.0
    totalNumData = 0
    for i = 1:server.num_clients
        c = server.clients[i]
        n = size(c.Xtrain, 1)
        totalNumData += n
        objValue += n * obj( c.XtrainT, c.Ytrain, server.W, c.λ)
    end
    return objValue / totalNumData
end