text:
@testset "Models" begin
    @info "Testing models..."

    @testset "Model constructor errors" begin
        grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
        @test_throws TypeError NonhydrostaticModel(architecture=CPU, grid=grid)
        @test_throws TypeError NonhydrostaticModel(architecture=GPU, grid=grid)
        @test_throws TypeError NonhydrostaticModel(grid=grid, boundary_conditions=1)
        @test_throws TypeError NonhydrostaticModel(grid=grid, forcing=2)
        @test_throws TypeError NonhydrostaticModel(grid=grid, background_fields=3)
    end

    topos = ((Periodic, Periodic, Periodic),
             (Periodic, Periodic,  Bounded),
             (Periodic,  Bounded,  Bounded),
             (Bounded,   Bounded,  Bounded))

    for topo in topos
        @testset "$topo model construction" begin
            @info "  Testing $topo model construction..."
            for arch in archs, FT in float_types
		        arch isa GPU && topo == (Bounded, Bounded, Bounded) && continue

                grid = RegularRectilinearGrid(FT, topology=topo, size=(16, 16, 2), extent=(1, 2, 3))
                model = NonhydrostaticModel(grid=grid, architecture=arch)

                @test model isa NonhydrostaticModel
            end
        end
    end

    @testset "Adjustment of halos in NonhydrostaticModel constructor" begin
        @info "  Testing adjustment of halos in NonhydrostaticModel constructor..."

        default_grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 2, 3))
        funny_grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 2, 3), halo=(1, 3, 4))

        # Model ensures that halos are at least of size 1
        model = NonhydrostaticModel(grid=default_grid)
        @test model.grid.Hx == 1 && model.grid.Hy == 1 && model.grid.Hz == 1

        model = NonhydrostaticModel(grid=funny_grid)
        @test model.grid.Hx == 1 && model.grid.Hy == 3 && model.grid.Hz == 4

        # Model ensures that halos are at least of size 2
        for scheme in (CenteredFourthOrder(), UpwindBiasedThirdOrder())
            model = NonhydrostaticModel(advection=scheme, grid=default_grid)
            @test model.grid.Hx == 2 && model.grid.Hy == 2 && model.grid.Hz == 2

            model = NonhydrostaticModel(advection=scheme, grid=funny_grid)
            @test model.grid.Hx == 2 && model.grid.Hy == 3 && model.grid.Hz == 4
        end

        # Model ensures that halos are at least of size 3
        for scheme in (WENO5(), UpwindBiasedFifthOrder())
            model = NonhydrostaticModel(advection=scheme, grid=default_grid)
            @test model.grid.Hx == 3 && model.grid.Hy == 3 && model.grid.Hz == 3

            model = NonhydrostaticModel(advection=scheme, grid=funny_grid)
            @test model.grid.Hx == 3 && model.grid.Hy == 3 && model.grid.Hz == 4
        end

        # Model ensures that halos are at least of size 2 with AnisotropicBiharmonicDiffusivity
        model = NonhydrostaticModel(closure=AnisotropicBiharmonicDiffusivity(), grid=default_grid)
        @test model.grid.Hx == 2 && model.grid.Hy == 2 && model.grid.Hz == 2

        model = NonhydrostaticModel(closure=AnisotropicBiharmonicDiffusivity(), grid=funny_grid)
        @test model.grid.Hx == 2 && model.grid.Hy == 3 && model.grid.Hz == 4
    end

    @testset "Model construction with single tracer and nothing tracer" begin
        @info "  Testing model construction with single tracer and nothing tracer..."
        for arch in archs
            grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 2, 3))

            model = NonhydrostaticModel(grid=grid, architecture=arch, tracers=:c, buoyancy=nothing)
            @test model isa NonhydrostaticModel

            model = NonhydrostaticModel(grid=grid, architecture=arch, tracers=nothing, buoyancy=nothing)
            @test model isa NonhydrostaticModel
        end
    end

    @testset "Setting model fields" begin
        @info "  Testing setting model fields..."
        for arch in archs, FT in float_types
            N = (4, 4, 4)
            L = (2π, 3π, 5π)

            grid = RegularRectilinearGrid(FT, size=N, extent=L)
            model = NonhydrostaticModel(architecture=arch, grid=grid)

            u, v, w = model.velocities
            T, S = model.tracers

            # Test setting an array
            T₀ = rand(FT, size(grid)...)
            T_answer = deepcopy(T₀)

            set!(model; enforce_incompressibility=false, T=T₀)

            @test Array(interior(T)) ≈ T_answer

            # Test setting functions
            u₀(x, y, z) = 1 + x + y + z
            v₀(x, y, z) = 2 + sin(x * y * z)
            w₀(x, y, z) = 3 + y * z
            T₀(x, y, z) = 4 + tanh(x + y - z)
            S₀(x, y, z) = 5

            set!(model, enforce_incompressibility=false, u=u₀, v=v₀, w=w₀, T=T₀, S=S₀)

            xC, yC, zC = nodes((Center, Center, Center), model.grid; reshape=true)
            xF, yF, zF = nodes((Face, Face, Face), model.grid; reshape=true)

            # Form solution arrays
            u_answer = u₀.(xF, yC, zC)
            v_answer = v₀.(xC, yF, zC)
            w_answer = w₀.(xC, yC, zF)
            T_answer = T₀.(xC, yC, zC)
            S_answer = S₀.(xC, yC, zC)

            Nx, Ny, Nz = size(model.grid)

            u_cpu = XFaceField(CPU(), grid)
            v_cpu = YFaceField(CPU(), grid)
            w_cpu = ZFaceField(CPU(), grid)
            T_cpu = CenterField(CPU(), grid)
            S_cpu = CenterField(CPU(), grid)

            set!(u_cpu, u)
            set!(v_cpu, v)
            set!(w_cpu, w)
            set!(T_cpu, T)
            set!(S_cpu, S)

            values_match = [
                            all(u_answer .≈ interior(u_cpu)),
                            all(v_answer .≈ interior(v_cpu)),
                            all(w_answer[:, :, 2:Nz] .≈ interior(w_cpu)[:, :, 2:Nz]),
                            all(T_answer .≈ interior(T_cpu)),
                            all(S_answer .≈ interior(S_cpu)),
                           ]

            @test all(values_match)

            # Test that update_state! works via u boundary conditions
            @test u_cpu[1, 1, 1] == u_cpu[Nx+1, 1, 1]  # x-periodicity
            @test u_cpu[1, 1, 1] == u_cpu[1, Ny+1, 1]  # y-periodicity
            @test all(u_cpu[1:Nx, 1:Ny, 1] .== u_cpu[1:Nx, 1:Ny, 0])     # free slip at bottom
            @test all(u_cpu[1:Nx, 1:Ny, Nz] .== u_cpu[1:Nx, 1:Ny, Nz+1]) # free slip at top

            # Test that enforce_incompressibility works
            set!(model, u=0, v=0, w=1, T=0, S=0)
            ϵ = 10 * eps(FT)
            set!(w_cpu, w)
            @test all(abs.(interior(w_cpu)) .< ϵ)

            # Test setting the background_fields to a Field
            U_field = XFaceField(arch, grid)
            U_field .= 1
            model = NonhydrostaticModel(grid=grid, architecture=arch, background_fields = (u=U_field,))
            @test model.background_fields.velocities.u isa Field

	    U_field = CenterField(arch, grid)
	    @test_throws ArgumentError NonhydrostaticModel(grid=grid, architecture=arch, background_fields = (u=U_field,))
        end
    end
end


text:
# Discern connectivity from the coordinates in a Geometry

using ArgMacros
using Logging
using ProtoBuf
using TFRecord

append!(LOAD_PATH, ["smu/jlout/", "smu/topology_from_geometry/", "smu/geometry"])

using BondLengthDistributions
using dataset_pb2
using TopologyFromGeometry

function get_topology_from_geometry(bond_length_distributions::AllBondLengthDistributions,
                bond_topology::BondTopology,
                geometry::Geometry,
                output_stream::IOStream)::Bool
  matching_parameters = MatchingParameters()
  result = bond_topologies_from_geom(bond_length_distributions, bond_topology,
                                     geometry, matching_parameters)
  hasproperty(result, :bond_topology) || return false
  number_bond_topologies = length(result.bond_topology)
  to_write = string(number_bond_topologies)
  for bt in result.bond_topology
    to_write *= ",$(round(bt.score, digits=3))"
    number_bond_topologies == 1 && break
    to_write *= bt.is_starting_topology ? ",t" : ",f"
  end
  println(output_stream, to_write)
  true
end

function main()
  @inlinearguments begin
    @argumentrequired String input_fname "-i" "--input"
    @argumentrequired String output_fname "-o" "--output"
    @argumentrequired String bonds "-b" "--bonds"
    @argumentoptional Int64 nprocess "-N" "--nprocess"
    @argumentflag exclude_non_bonded "-x" "--xcldnonbond"
    @argumentflag debug "-debug" "--debug"
    @argumentflag verbose "-v" "--verbose"
  end
  if debug
    logger=Logging.SimpleLogger(stderr,Logging.Debug)
    global_logger(logger)
  end
# global_logger(Logging.SimpleLogger(stderr, Logging.Error))

  @info("input $(input_fname) bonds $(bonds) output $(output_fname) $nprocess")
  flush(stdout)
  nprocess === nothing && (nprocess = typemax(Int64))

  bond_length_distributions = AllBondLengthDistributions()
  exclude_non_bonded !== nothing && (bond_length_distributions.include_non_bonded = false)

  add_from_files!(bonds, bond_length_distributions) || @error("Cannot build bond length distribution $bonds")

  molecules_read = 0
  molecules_processed = 0

  open(output_fname, "w") do output_stream
    for conformer in TFRecord.read([input_fname], record_type=Conformer)
      molecules_read += 1
#     println(conformer)
      if ! hasproperty(conformer, :optimized_geometry)
        @warn("No optimized geometry, skipping")
        continue
      end

      if ! hasproperty(conformer, :bond_topologies)
        @warn("No bond_topology, skipping")
        continue
      end
      @info("processing $(conformer.conformer_id)")
      if get_topology_from_geometry(bond_length_distributions, conformer.bond_topologies[1],
                                    conformer.optimized_geometry, output_stream)
        molecules_processed += 1
      end
      molecules_read > nprocess && break
    end
  end
  verbose && println("Read $(molecules_read) molecules, processed $(molecules_processed)")
  0
end

main()

ext: jl

max_stars_repo_name: IanAWatson/smu

max_issues_repo_name: IanAWatson/smu

max_forks_repo_name: IanAWatson/smu


text:
"""
    Aqua.test_deps_compat(package)

Test that `Project.toml` of `package` list all `compat` for `deps`.

# Arguments
- `packages`: a top-level `Module`, a `Base.PkgId`, or a collection of
  them.
"""
test_deps_compat
function test_deps_compat(packages)
    @testset "$(result.label)" for result in analyze_deps_compat(packages)
        @debug result.label result
        @test result ⊜ true
    end
end

analyze_deps_compat(packages) = [_analyze_deps_compat_1(pkg) for pkg in aspkgids(packages)]

function _analyze_deps_compat_1(pkg::PkgId)
    result = root_project_or_failed_lazytest(pkg)
    result isa LazyTestResult && return result
    root_project_path = result
    return _analyze_deps_compat_2(pkg, root_project_path, TOML.parsefile(root_project_path))
end

# For supporting Julia 1.8-DEV and above which give us a tuple instead of a string
_unwrap_name(x::Tuple) = first(x)
_unwrap_name(x::String) = x
_unwrap_name(x::Nothing) = x
function _analyze_deps_compat_2(pkg::PkgId, root_project_path, prj)
    label = "$pkg"

    deps = get(prj, "deps", nothing)
    if deps === nothing
        return LazyTestResult(label, "`$root_project_path` does not have `deps`", true)
    end
    compat = get(prj, "compat", nothing)
    if compat === nothing
        return LazyTestResult(label, "`$root_project_path` does not have `compat`", false)
    end

    stdlib_name_from_uuid = stdlibs()
    stdlib_deps = filter!(
        !isnothing,
        [_unwrap_name(get(stdlib_name_from_uuid, UUID(uuid), nothing)) for (_, uuid) in deps],
    )
    missing_compat = setdiff(setdiff(keys(deps), keys(compat)), stdlib_deps)
    if !isempty(missing_compat)
        msg = join(
            [
                "`$root_project_path` does not specify `compat` for:"
                sort!(collect(missing_compat))
            ],
            "\n",
        )
        return LazyTestResult(label, msg, false)
    end

    return LazyTestResult(
        label,
        "`$root_project_path` specifies `compat` for all `deps`",
        true,
    )
end


text:
using Test

include("vehicle-purchase.jl")

@testset "Affordability" begin
    @test canibuy("1974 Ford Pinto", 516, 100) == "Yes! I'm getting a 1974 Ford Pinto."
    @test canibuy("2014 Bugatti Veyron", 562_500, 5000) == "Damn! No 2014 Bugatti Veyron for me."
    @test canibuy("2020 Gazelle Medeo", 3000, 50) == "I'll have to be frugal if I want a 2020 Gazelle Medeo."
end

@testset "Licence" begin
    @test licence("2014 Bugatti Veyron", "car") == "The 2014 Bugatti Veyron requires a licence to operate."
    @test licence("2020 Gazelle Medeo", "bike") == "The 2020 Gazelle Medeo requires no licence to operate."
end

@testset "Registration Fees" begin
    @test fee(  562_500,  6,  "car") == 2250
    @test fee(   25_000,  3,  "car") ==  175
    @test fee(   34_000, 30,  "car") ==   25
    @test fee(    3_000,  0, "bike") ==    0
end


text:
struct DiagonallyDominantBridge{T, F, G} <: MOIB.Constraint.AbstractBridge
    # |Qij| variables
    abs_vars::Vector{MOI.VariableIndex}
    # |Qij| ≥ +Qij
    abs_plus::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},
                                         MOI.GreaterThan{T}}}
    # |Qij| ≥ -Qij
    abs_minus::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},
                                          MOI.GreaterThan{T}}}
    # inequalities Qjj ≥ sum_{i ≠ j} |Qij|
    dominance::Vector{MOI.ConstraintIndex{F, MOI.GreaterThan{T}}}
end

function MOIB.Constraint.bridge_constraint(
    ::Type{DiagonallyDominantBridge{T, F, G}},
    model::MOI.ModelLike, f::MOI.AbstractVectorFunction,
    s::SOS.DiagonallyDominantConeTriangle) where {T, F, G}

    @assert MOI.output_dimension(f) == MOI.dimension(s)
    n = s.side_dimension
    g = F[zero(F) for i in 1:n]
    fs = MOIU.eachscalar(f)
    num_off_diag = MOI.dimension(s) - n
    CI = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}
    abs_vars = Vector{MOI.VariableIndex}(undef, num_off_diag)
    abs_plus = Vector{CI}(undef, num_off_diag)
    abs_minus = Vector{CI}(undef, num_off_diag)
    k = 0
    koff = 0
    for j in 1:n
        # Qjj ≥ sum_{i ≠ j} |Qij|
        for i in 1:(j-1)
            k += 1
            koff += 1
            # abs ≥ |Qij|
            abs_vars[koff] = MOI.add_variable(model)
            fabs = MOI.SingleVariable(abs_vars[koff])
            MOIU.operate!(-, T, g[j], fabs)
            MOIU.operate!(-, T, g[i], fabs)
            abs_plus[koff] = MOI.add_constraint(
                model, MOIU.operate(+, T, fabs, fs[k]), MOI.GreaterThan(0.0))
            abs_minus[koff] = MOI.add_constraint(
                model, MOIU.operate(-, T, fabs, fs[k]), MOI.GreaterThan(0.0))
        end
        k += 1
        MOIU.operate!(+, T, g[j], fs[k])
    end
    dominance = map(f -> MOI.add_constraint(model, f, MOI.GreaterThan(0.0)), g)
    return DiagonallyDominantBridge{T, F, G}(abs_vars, abs_plus, abs_minus,
                                          dominance)
end

function MOI.supports_constraint(::Type{<:DiagonallyDominantBridge},
                                 ::Type{<:MOI.AbstractVectorFunction},
                                 ::Type{<:SOS.DiagonallyDominantConeTriangle})
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:DiagonallyDominantBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:DiagonallyDominantBridge{T, F}}) where {T, F}
    added = [(F, MOI.GreaterThan{T})]
    if F != MOI.ScalarAffineFunction{T}
        push!(added, (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}))
    end
    return added
end
function MOIB.Constraint.concrete_bridge_type(
    ::Type{<:DiagonallyDominantBridge{T}},
    G::Type{<:MOI.AbstractVectorFunction},
    ::Type{SOS.DiagonallyDominantConeTriangle}) where T

    S = MOIU.scalar_type(G)
    F = MOIU.promote_operation(-, T, S, MOI.SingleVariable)
    return DiagonallyDominantBridge{T, F, G}
end

# Attributes, Bridge acting as an model
function MOI.get(bridge::DiagonallyDominantBridge, ::MOI.NumberOfVariables)
    return length(bridge.abs_vars)
end
function MOI.get(bridge::DiagonallyDominantBridge, ::MOI.ListOfVariableIndices)
    return bridge.abs_vars
end
function MOI.get(bridge::DiagonallyDominantBridge{T, MOI.ScalarAffineFunction{T}},
                 ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},
                                           MOI.GreaterThan{T}}) where T
    return length(bridge.abs_plus) + length(bridge.abs_minus) + length(bridge.dominance)
end
function MOI.get(bridge::DiagonallyDominantBridge{T},
                 ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},
                                           MOI.GreaterThan{T}}) where T
    return length(bridge.abs_plus) + length(bridge.abs_minus)
end
function MOI.get(bridge::DiagonallyDominantBridge{T, F},
                 ::MOI.NumberOfConstraints{F, MOI.GreaterThan{T}}) where {T, F}
    return length(bridge.dominance)
end
function MOI.get(bridge::DiagonallyDominantBridge{T, MOI.ScalarAffineFunction{T}},
                 ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},
                                               MOI.GreaterThan{T}}) where T
    return vcat(bridge.abs_plus, bridge.abs_minus, bridge.dominance)
end
function MOI.get(bridge::DiagonallyDominantBridge{T},
                 ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},
                                               MOI.GreaterThan{T}}) where T
    return vcat(bridge.abs_plus, bridge.abs_minus)
end
function MOI.get(bridge::DiagonallyDominantBridge{T, F},
                 ::MOI.ListOfConstraintIndices{F, MOI.GreaterThan{T}}) where {T, F}
    return bridge.dominance
end

# Indices
function MOI.delete(model::MOI.ModelLike, bridge::DiagonallyDominantBridge)
    for ci in bridge.dominance
        MOI.delete(model, ci)
    end
    for ci in bridge.abs_plus
        MOI.delete(model, ci)
    end
    for ci in bridge.abs_minus
        MOI.delete(model, ci)
    end
    for vi in bridge.abs_vars
        MOI.delete(model, vi)
    end
end

# Attributes, Bridge acting as a constraint
function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::DiagonallyDominantBridge)
    return SOS.DiagonallyDominantConeTriangle(length(bridge.dominance))
end
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::DiagonallyDominantBridge{T, F, G}) where {T, F, G}
    set = MOI.get(model, MOI.ConstraintSet(), bridge)
    H = MOIU.scalar_type(G)
    g = Vector{H}(undef, MOI.dimension(set))
    k = 0
    koff = 0
    for j in 1:MOI.side_dimension(set)
        for i in 1:(j - 1)
            k += 1
            koff += 1
            func = MOI.get(model, attr, bridge.abs_plus[koff])
            g[k] = MOIU.convert_approx(H, MOIU.remove_variable(
                func, bridge.abs_vars))
        end
        k += 1
        func = MOI.get(model, attr, bridge.dominance[j])
        g[k] = MOIU.convert_approx(H, MOIU.remove_variable(
            func, bridge.abs_vars))
    end
    return MOIU.vectorize(g)
end

# TODO ConstraintPrimal

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::DiagonallyDominantBridge{T}) where T
    dominance_dual = MOI.get(model, attr, bridge.dominance)
    side_dim = length(dominance_dual)
    dim = MOI.dimension(SOS.DiagonallyDominantConeTriangle(side_dim))
    dual = Array{T}(undef, dim)
    k = 0
    for j in 1:side_dim
        for i in 1:(j-1)
            k += 1
            # Need to divide by 2 because of the custom scalar product for this
            # cone
            dual[k] = (- dominance_dual[i] - dominance_dual[j]) / 2
        end
        k += 1
        dual[k] = dominance_dual[j]
    end
    @assert k == dim
    return dual
end


text:
for dtype in ("", "Dist")
    for (elty, ext) in ((:Float32, :s),
                        (:Float64, :d),
                        (:Complex64, :c),
                        (:Complex128, :z))
        mat = Symbol(dtype, "Matrix")
        @eval begin
            function inverse!(A::$mat{$elty})
                err = ccall(($(string("ElInverse", dtype, "_", ext)), libEl), Cuint,
                    (Ptr{Void},),
                    A.obj,)
                err == 0 || throw(ElError(err))
                return A
            end
        end
    end
end


text:
function kdotp(lgir::LGIrrep{D}, αβγ=nothing; timereversal::Bool=true) where D
    lg = group(lgir) # TODO: maybe e.g., skip identity?
    hs = Hermitian.(gellmann(Crystalline.irdim(lgir); skip_identity=false, norm_identity=true))
    hs = hs
    N = length(hs)
    G = length(lg)

    Γs = lgir(αβγ)
    kv = position(lgir)

    # assembling matrices
    # TODO: there is some minimal necessary set of operators to include here; possibly just
    #       the generators?
    A = zeros(Float64, N*D*G, N*D)
    for g in eachindex(lg)
        Γ  = Γs[g]
        Γ⁻¹ = inv(Γ)
        for m in 1:N
            for n in 1:N
                tr_hₘΓhₙΓ⁻¹ = tr(hs[m]'*Γ*hs[n]*Γ⁻¹)
                if abs(imag(tr_hₘΓhₙΓ⁻¹)) > ATOL_DEFAULT
                    error("element of real matrix had non-negligible imaginary part $tr_hₘΓhₙΓ⁻¹")
                end
                tr_hₘΓhₙΓ⁻¹_re = real(tr_hₘΓhₙΓ⁻¹)
                abs(tr_hₘΓhₙΓ⁻¹_re) > ATOL_DEFAULT || continue
                for i in 1:D # δⱼᵢ
                    A[(g-1)*(N*D) + (i-1)*N + m, (i-1)*N + n] = tr_hₘΓhₙΓ⁻¹_re
                end
            end
        end
    end

    γ2 = zeros(Float64, N*D*G, N*D)
    for (g, op) in enumerate(lg)
        R = rotation(op) # apparently, there should _not_ be a transpose here¹⁾
        # ¹⁾ g is imagined as acting on a parameter "k", rather than a "k-vector" per se.
        # R = inv(R) # R'
        for i in 1:D
            for j in 1:D
                Rⱼᵢ = R[j,i]
                iszero(Rⱼᵢ) && continue
                for m in 1:N # δₙₘ
                    γ2[(g-1)*(N*D) + (i-1)*N + m, (j-1)*N + m] = 2Rⱼᵢ
                end
            end
        end
    end

    if Crystalline.iscorep(lgir) && !timereversal
        error("provided irrep is a corep (i.e., is modified by and assumes time-reversal symmetry), but the timereversal keyword is set to `false`.")
    end

    # if assumed present, check if time-reversal has an impact
    if timereversal # TODO: also check if already a corep; if so, TR _must_ be assumed
        g₀ = kv_to_negative_kv_operation(lgir)
        if !isnothing(g₀)
            # Find a representation of 𝒯g₀ [Γ(𝒯g₀)K, with K denoting complex conjugation]
            # and then apply the constraint
            #                    Γ(𝒯g₀)H(k)Γ(𝒯g₀)⁻¹ = H*(−g₀k)
            Γ𝒯g₀ = find_antiunitary_irrep(lgir, g₀)
            Γ𝒯g₀⁻¹ = inv(Γ𝒯g₀) # TODO: could just be Γ𝒯g₀†, right?

            A𝒯 = zeros(Float64, N*D, N*D)
            for m in 1:N
                for n in 1:N
                    tr_hₘΓhₙΓ⁻¹ = tr(hs[m]'*Γ𝒯g₀*hs[n]*Γ𝒯g₀⁻¹)
                    if abs(imag(tr_hₘΓhₙΓ⁻¹)) > ATOL_DEFAULT
                        error("element of real matrix had non-negligible imaginary part $tr_hₘΓhₙΓ⁻¹")
                    end
                    tr_hₘΓhₙΓ⁻¹_re = real(tr_hₘΓhₙΓ⁻¹)
                    abs(tr_hₘΓhₙΓ⁻¹_re) > ATOL_DEFAULT || continue
                    for i in 1:D # δⱼᵢ
                        A𝒯[(i-1)*N + m, (i-1)*N + n] = tr_hₘΓhₙΓ⁻¹_re
                    end
                end
            end

            γ2𝒯 = zeros(Float64, N*D, N*D)
            R₀ = rotation(g₀) # see ¹⁾ regarding transposition/no transposition
            conj_prefs = [real(tr(conj(h)'*h)) for h in hs] # ⟨hₘ*|hₘ⟩_F = ±2
            for i in 1:D
                for j in 1:D
                    R₀ⱼᵢ = R₀[j,i]
                    iszero(R₀ⱼᵢ) && continue
                    for m in 1:N # δₙₘ
                        γ2𝒯[(i-1)*N + m, (j-1)*N + m] = -conj_prefs[m]*R₀ⱼᵢ
                    end
                end
            end
            A  = vcat(A, A𝒯)   # TODO: do this without allocating twice?
            γ2 = vcat(γ2, γ2𝒯)
        end
    end

    # solving nullspace (c = [c_11, …, c_N1, c_12, …, cN2, …, c_1D, …, c_ND] = cₙᵈ)
    _cs = nullspace(A - γ2, atol=ATOL_DEFAULT)
    cs = sparsify_columns(_cs, atol=ATOL_DEFAULT)

    # drop very small parts
    map!(x->abs(x) > ATOL_DEFAULT ? x : 0.0, cs, cs)

    # split coefficients into meaningful structure
    cs′ = [ntuple(d->cs[(d-1)*N+1:d*N,a], Val(D)) for a in 1:size(cs, 2)]

    # remove/prune coefficients/matrices that do not appear in the allowed expansions
    redundant_hs = Int[]
    for n in 1:N
        redundant = true
        for cₐ in cs′
            for cₐᵈ in cₐ
                !iszero(cₐᵈ[n]) && (redundant=false; break)
            end
            redundant || break
        end
        redundant && push!(redundant_hs, n)
    end
    deleteat!(hs, redundant_hs)
    for cₐ in cs′
        for cₐᵈ in cₐ
            deleteat!(cₐᵈ, redundant_hs)
        end
    end

    # normalize coefficients for pretty-printing
    for cₐ in cs′
        scale = zero(Float64)
        for cₐᵈ in cₐ
            v, idx = findmax(abs, cₐᵈ)
            scale = abs(scale) < v ? cₐᵈ[idx] : scale
        end
        for cₐᵈ in cₐ
            cₐᵈ ./= scale
        end
    end

    return KPHamiltonian{D}(lgir, hs, cs′)
end

# matrix sparsification (adapted from https://math.stackexchange.com/a/4227787)
function sparsify_columns(A; atol=ATOL_DEFAULT)
    n = size(A,2)
    if n ≤ 1
        return A
    else
        # by assumption, the n columns of A are linearly independent
        rank(A) ≠ n && error("columns of A must be independent")

        # row rank == column rank: so, we can find n linearly independent rows in A
        indep_rows = Vector{eltype(A)}[]
        for row in eachrow(A)
            # if row is linearly independent with all existing ones, add it to the list
            if isempty(nullspace(hcat(indep_rows..., row); atol=atol))
                push!(indep_rows, row)
            end
        end
        return A * inv(hcat(indep_rows...))'
    end
end


text:
using MPI
using StaticArrays
using ClimateMachine
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.MPIStateArrays
using ClimateMachine.DGmethods
using ClimateMachine.DGmethods.NumericalFluxes
using Printf
using LinearAlgebra
using Logging
using GPUifyLoops

import ClimateMachine.DGmethods:
    BalanceLaw,
    vars_state_auxiliary,
    vars_state_conservative,
    vars_state_gradient,
    vars_state_gradient_flux,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    LocalGeometry,
    boundary_state!,
    init_state_auxiliary!,
    init_state_conservative!,
    init_ode_state,
    update_auxiliary_state!,
    vars_integrals,
    vars_reverse_integrals,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

struct IntegralTestModel{dim} <: BalanceLaw end

vars_reverse_integrals(::IntegralTestModel, T) = @vars(a::T, b::T)
vars_integrals(::IntegralTestModel, T) = @vars(a::T, b::T)
vars_state_auxiliary(m::IntegralTestModel, T) = @vars(
    int::vars_integrals(m, T),
    rev_int::vars_reverse_integrals(m, T),
    coord::SVector{3, T},
    a::T,
    b::T,
    rev_a::T,
    rev_b::T
)

vars_state_conservative(::IntegralTestModel, T) = @vars()
vars_state_gradient_flux(::IntegralTestModel, T) = @vars()

flux_first_order!(::IntegralTestModel, _...) = nothing
flux_second_order!(::IntegralTestModel, _...) = nothing
source!(::IntegralTestModel, _...) = nothing
boundary_state!(_, ::IntegralTestModel, _...) = nothing
init_state_conservative!(::IntegralTestModel, _...) = nothing
wavespeed(::IntegralTestModel, _...) = 1

function init_state_auxiliary!(
    ::IntegralTestModel{dim},
    aux::Vars,
    g::LocalGeometry,
) where {dim}
    x, y, z = aux.coord = g.coord
    if dim == 2
        aux.a = x * y + z * y
        aux.b = 2 * x * y + sin(x) * y^2 / 2 - (z - 1)^2 * y^3 / 3
        y_top = 3
        a_top = x * y_top + z * y_top
        b_top = 2 * x * y_top + sin(x) * y_top^2 / 2 - (z - 1)^2 * y_top^3 / 3
        aux.rev_a = a_top - aux.a
        aux.rev_b = b_top - aux.b
    else
        aux.a = x * z + z^2 / 2
        aux.b = 2 * x * z + sin(x) * y * z - (1 + (z - 1)^3) * y^2 / 3
        zz_top = 3
        a_top = x * zz_top + zz_top^2 / 2
        b_top =
            2 * x * zz_top + sin(x) * y * zz_top -
            (1 + (zz_top - 1)^3) * y^2 / 3
        aux.rev_a = a_top - aux.a
        aux.rev_b = b_top - aux.b
    end
end

function update_auxiliary_state!(
    dg::DGModel,
    m::IntegralTestModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    indefinite_stack_integral!(dg, m, Q, dg.state_auxiliary, t, elems)
    reverse_indefinite_stack_integral!(dg, m, Q, dg.state_auxiliary, t, elems)

    return true
end

@inline function integral_load_auxiliary_state!(
    m::IntegralTestModel,
    integrand::Vars,
    state::Vars,
    aux::Vars,
)
    x, y, z = aux.coord
    integrand.a = x + z
    integrand.b = 2 * x + sin(x) * y - (z - 1)^2 * y^2
end

@inline function integral_set_auxiliary_state!(
    m::IntegralTestModel,
    aux::Vars,
    integral::Vars,
)
    aux.int.a = integral.a
    aux.int.b = integral.b
end

@inline function reverse_integral_load_auxiliary_state!(
    m::IntegralTestModel,
    integral::Vars,
    state::Vars,
    aux::Vars,
)
    integral.a = aux.int.a
    integral.b = aux.int.b
end

@inline function reverse_integral_set_auxiliary_state!(
    m::IntegralTestModel,
    aux::Vars,
    integral::Vars,
)
    aux.rev_int.a = integral.a
    aux.rev_int.b = integral.b
end

using Test
function run(mpicomm, dim, Ne, N, FT, ArrayType)

    brickrange = ntuple(j -> range(FT(0); length = Ne[j] + 1, stop = 3), dim)
    topl = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = ntuple(j -> true, dim),
    )

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )
    dg = DGModel(
        IntegralTestModel{dim}(),
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    Q = init_ode_state(dg, FT(0))
    dQdt = similar(Q)

    dg(dQdt, Q, nothing, 0.0)

    # Wrapping in Array ensure both GPU and CPU code use same approx
    @test Array(dg.state_auxiliary.data[:, 1, :]) ≈
          Array(dg.state_auxiliary.data[:, 8, :])
    @test Array(dg.state_auxiliary.data[:, 2, :]) ≈
          Array(dg.state_auxiliary.data[:, 9, :])
    @test Array(dg.state_auxiliary.data[:, 3, :]) ≈
          Array(dg.state_auxiliary.data[:, 10, :])
    @test Array(dg.state_auxiliary.data[:, 4, :]) ≈
          Array(dg.state_auxiliary.data[:, 11, :])
end

let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    numelem = (5, 5, 5)
    lvls = 1

    polynomialorder = 4

    for FT in (Float64,) #Float32)
        for dim in 2:3
            err = zeros(FT, lvls)
            for l in 1:lvls
                @info (ArrayType, FT, dim)
                run(
                    mpicomm,
                    dim,
                    ntuple(j -> 2^(l - 1) * numelem[j], dim),
                    polynomialorder,
                    FT,
                    ArrayType,
                )
            end
        end
    end
end

nothing

max_stars_repo_name: sandreza/CLIMA

max_issues_repo_name: sandreza/CLIMA

max_forks_repo_name: sandreza/CLIMA


text:
export Halfspace

struct Halfspace{T <: Real} <: ConvexSet
    a::Vector{T}
    b::T
end

function project!(s::Halfspace{T}, v::Vector{T}, y::Vector{T}) where {T <: Real}
   c = transpose(s.a)* y
  if c > s.b
    v = y - ((c - s.b) / sum(abs2, s.a)) * s.a
  else
    v = y
  end
end

function project(s::Halfspace{<:Real}, y::Vector{<:Real})
    project!(s, similar(y), y)
end


text:
using Test
using PrettyTables
using DataFrames

println("Text backend")
println("============")
println()
include("./text_backend.jl")

println("HTML backend")
println("============")
println()
include("./html_backend.jl")


text:
################################################################################
# [3] SUPPORTING FUNCTIONS
################################################################################

# ------------------------------------------------------
# [3.1] Data for VAR
# ------------------------------------------------------
function VARData(y::Array, p; Cons = true, trend = false, τ = 1)
    T, n = size(y)
    T = T - p
    r = p + 1
    Y = y[r:end, :]
    X = Array{Float64}(undef, T, n * p + Cons + trend * τ)
    for i = 1:p
        X[:, n*(i-1)+1:i*n] = y[r-i:end-i, :]
    end
    if Cons
        X[:, end-trend*τ] = ones(T)
    end
    if trend
        for i = 1:τ
            X[:, end-i+1] = (1:T) .^ (τ - i + 1)
        end
    end
    return Y, X
end

# ------------------------------------------------------
# [3.2] OLS estimation
# ------------------------------------------------------
function OLSbetas(Y::Array, X::Array, m, p; Xblock = false, Nx = 0)
    B = (X' * X) \ (X' * Y)
    if Xblock
        Bres = ROLS(X, B, m, p, Nx)
        B = Bres
    end
    f = B'
    f = f[:, 1:m*p]
    Φ = [f; I(m * (p - 1)) zeros(m * (p - 1), m)]

    if ~isa(Φ, Array)
        Φ = Array(Φ) # Just to apply eigvals later
    end
    return B, Φ
end

# ------------------------------------------------------
# [3.3] Restricted OLS
# ------------------------------------------------------
function ROLS(x, β, m, p, nx)
    row, cols = size(β)
    nres = (m - nx) * nx * p
    RR = zeros(row, cols)

    for i = 1:p
        RR[m * (i - 1)+nx+1:m*i, 1:nx] = ones(m - nx, nx)
    end

    R = zeros(nres, length(β))
    index = findall(x -> x == 1, RR)

    for i = 1:nres
        R[i, (index[i][2]-1)*row+index[i][1]] = 1
    end
    # Objectives
    r = zeros(nres, 1)
    α = vec(β)

    # Restricted Estimation
    X = kron(I(m), x)
    XX = inv(X' * X)
    αR = α + (XX * R') * (inv((R * XX * R')) * (r - R * α))
    αR[vec(RR).==1] .= 0  # Jusdt to eliminate the round difference
    return reshape(αR, (row, cols))
end

# --------------------------------------------------------
# [3.4] Simulating β
# --------------------------------------------------------
# from a diffuse inverse wishart prior
function BetaDraw(disti, rc::Int, cc::Int, m::Int, p::Int)
    status = 1
    β_draw = Array{Float64,2}(undef, rc, cc)
    Φ_draw = Array{Float64,2}(undef, m * p, m * p)
    while status >= 1
        b_draw = rand(disti)
        β_draw = reshape(b_draw, (rc, cc))
        f = β_draw'
        f = f[:, 1:m*p]
        Φ_draw = [f; I(m * (p - 1)) zeros(m * (p - 1), m)]
        if ~isa(Φ_draw, Array)
            Φ_draw = Array(Φ_draw)
        end
        status = maximum(abs.(eigvals(Φ_draw)))
    end
    return β_draw, Φ_draw
end
# ------------------------------------------------------
# [3.4] Calculating matrix of cumulative FEVD
# ------------------------------------------------------
function Λmatrix(IRF, nvar, m; Lτ = 1, Uτ = 5, Xblock = false, nx = 0)
    Λ = zeros(m, m)
    for i = 1:Uτ
        R_tilde = IRF[nvar, :, i]
        Λ = Λ + ((Uτ - 1) + 1 - max(Lτ - 1, i - 1)) * (R_tilde * R_tilde')
    end
    if Xblock
        Λ = Λ[1:nx, 1:nx]
    end
    return Λ
end

# ------------------------------------------------------
# [3.5] IRF and FEVD functions
# ------------------------------------------------------
function irf_fevd(Θ, Γ, h, m)
    # --------------------------------------------------
    # [3.5.1] Impulse Response functions
    n = size(Θ)[1]
    IRFred = Array{Float64,3}(undef, n, n, h)
    IRFred[:, :, 1] = Θ
    IRF = Array{Float64,3}(undef, m, m, h)
    IRF[:, :, 1] = Θ[1:m, :1:m] * Γ

    for i = 2:h
        IRFred[:, :, i] = Θ * IRFred[:, :, i-1]
        IRF[:, :, i] = IRFred[1:m, 1:m, i] * Γ
    end

    # --------------------------------------------------
    # [3.5.2] Forecast Error variance Decomposition
    MSPEkj = cumsum(IRF .^ 2, dims = 3)
    MSPEk = sum(MSPEkj, dims = 2)
    FEVD = MSPEkj ./ MSPEk
    return IRF, FEVD
end

# ------------------------------------------------------
# [3.6] Get quantiles
# ------------------------------------------------------
function Qntls(x, nmodls, quintls, m, h)
    xvec = [vec(x[:, :, i]) for i = 1:nmodls]
    xvec = hcat(xvec...)
    sort!(xvec, dims = 2)
    indxs = convert(Array{Int64}, floor.(quintls .* nmodls))
    xqnt = [reshape(xvec[:, indxs[i]], (m, h)) for i = 1:length(indxs)]
    xtup = Tuple(x for x in xqnt)
    perci = "perct" .* string.(Int.(floor.(quintls * 100)))
    perci = Symbol.(perci)
    keys = (:perct1, :perct2, :perct3)
    xNtup = (; zip(perci, xtup)...)
    return xNtup
end

max_stars_repo_name: ChrisVeJa/GlobalShock

max_issues_repo_name: ChrisVeJa/GlobalShock

max_forks_repo_name: ChrisVeJa/GlobalShock


text:
function Base.resize!(A::GBMatrix, nrows::Integer, ncols::Integer)
    libgb.GrB_Matrix_resize(A, libgb.GrB_Index(nrows), libgb.GrB_Index(ncols))
    return A
end

function Base.resize!(v::GBVector, nrows::Integer)
    libgb.GrB_Vector_resize(v, libgb.GrB_Index(nrows))
    return v
end


text:
using Mocking

compare_faces_1_patch = @patch function compare_faces(a...; b...)
    return Dict{String, Any}("SourceImageFace" => Dict{String, Any}("Confidence" => 99.9991226196289, "BoundingBox" => Dict{String, Any}("Width" => 0.25, "Left" => 0.31888890266418457, "Height" => 0.33481481671333313, "Top" => 0.4933333396911621)), "FaceMatches" => Any[Dict{String, Any}("Similarity" => 100, "Face" => Dict{String, Any}("Confidence" => 99.9991226196289, "BoundingBox" => Dict{String, Any}("Width" => 0.25, "Left" => 0.31888890266418457, "Height" => 0.33481481671333313, "Top" => 0.4933333396911621)))])
end

create_collection_1_patch = @patch function create_collection(a...; b...)
    return Dict{String, Any}("CollectionArn" => "aws:rekognition:us-west-2:123456789012:collection/myphotos", "StatusCode" => 200)
end

delete_collection_1_patch = @patch function delete_collection(a...; b...)
    return Dict{String, Any}("StatusCode" => 200)
end

delete_faces_1_patch = @patch function delete_faces(a...; b...)
    return Dict{String, Any}("DeletedFaces" => Any["ff43d742-0c13-5d16-a3e8-03d3f58e980b"])
end

detect_faces_1_patch = @patch function detect_faces(a...; b...)
    return Dict{String, Any}("OrientationCorrection" => "ROTATE_0", "FaceDetails" => Any[Dict{String, Any}("Quality" => Dict{String, Any}("Brightness" => 37.60169982910156, "Sharpness" => 80), "Landmarks" => Any[Dict{String, Any}("Y" => 0.40819624066352844, "X" => 0.6394737362861633, "Type" => "eyeLeft"), Dict{String, Any}("Y" => 0.41039225459098816, "X" => 0.7266660928726196, "Type" => "eyeRight"), Dict{String, Any}("Y" => 0.44240960478782654, "X" => 0.6912462115287781, "Type" => "eyeRight"), Dict{String, Any}("Y" => 0.46700039505958557, "X" => 0.6306198239326477, "Type" => "mouthDown"), Dict{String, Any}("Y" => 0.47114261984825134, "X" => 0.7215608954429626, "Type" => "mouthUp")], "Pose" => Dict{String, Any}("Roll" => 0.9950747489929199, "Yaw" => 13.693790435791016, "Pitch" => 4.050806522369385), "Confidence" => 100, "BoundingBox" => Dict{String, Any}("Width" => 0.23999999463558197, "Left" => 0.5555555820465088, "Height" => 0.18000000715255737, "Top" => 0.33666667342185974))])
end

detect_labels_1_patch = @patch function detect_labels(a...; b...)
    return Dict{String, Any}("Labels" => Any[Dict{String, Any}("Confidence" => 99.25072479248047, "Name" => "People"), Dict{String, Any}("Confidence" => 99.25074005126953, "Name" => "Person")])
end

index_faces_1_patch = @patch function index_faces(a...; b...)
    return Dict{String, Any}("OrientationCorrection" => "ROTATE_0", "FaceRecords" => Any[Dict{String, Any}("FaceDetail" => Dict{String, Any}("Quality" => Dict{String, Any}("Brightness" => 29.23358917236328, "Sharpness" => 80), "Landmarks" => Any[Dict{String, Any}("Y" => 0.6248345971107483, "X" => 0.3976764678955078, "Type" => "eyeLeft"), Dict{String, Any}("Y" => 0.6317117214202881, "X" => 0.4810936450958252, "Type" => "eyeRight"), Dict{String, Any}("Y" => 0.7111940383911133, "X" => 0.41986238956451416, "Type" => "noseLeft"), Dict{String, Any}("Y" => 0.7497701048851013, "X" => 0.40525302290916443, "Type" => "mouthDown"), Dict{String, Any}("Y" => 0.7558549642562866, "X" => 0.4753248989582062, "Type" => "mouthUp")], "Pose" => Dict{String, Any}("Roll" => 4.707281112670898, "Yaw" => -24.438663482666016, "Pitch" => -9.713645935058594), "Confidence" => 99.9991226196289, "BoundingBox" => Dict{String, Any}("Width" => 0.25, "Left" => 0.31888890266418457, "Height" => 0.33481481671333313, "Top" => 0.4933333396911621)), "Face" => Dict{String, Any}("FaceId" => "ff43d742-0c13-5d16-a3e8-03d3f58e980b", "Confidence" => 99.9991226196289, "ImageId" => "465f4e93-763e-51d0-b030-b9667a2d94b1", "BoundingBox" => Dict{String, Any}("Width" => 0.25, "Left" => 0.31888890266418457, "Height" => 0.33481481671333313, "Top" => 0.4933333396911621))), Dict{String, Any}("FaceDetail" => Dict{String, Any}("Quality" => Dict{String, Any}("Brightness" => 34.951698303222656, "Sharpness" => 160), "Landmarks" => Any[Dict{String, Any}("Y" => 0.290842205286026, "X" => 0.6006892323493958, "Type" => "eyeLeft"), Dict{String, Any}("Y" => 0.29609042406082153, "X" => 0.6808141469955444, "Type" => "eyeRight"), Dict{String, Any}("Y" => 0.3522595763206482, "X" => 0.6395332217216492, "Type" => "noseLeft"), Dict{String, Any}("Y" => 0.38689887523651123, "X" => 0.5892083048820496, "Type" => "mouthDown"), Dict{String, Any}("Y" => 0.394125759601593, "X" => 0.674560010433197, "Type" => "mouthUp")], "Pose" => Dict{String, Any}("Roll" => 2.1029529571533203, "Yaw" => 6.716655254364014, "Pitch" => -4.683138370513916), "Confidence" => 99.99950408935547, "BoundingBox" => Dict{String, Any}("Width" => 0.24444444477558136, "Left" => 0.5144444704055786, "Height" => 0.32592591643333435, "Top" => 0.15111111104488373)), "Face" => Dict{String, Any}("FaceId" => "8be04dba-4e58-520d-850e-9eae4af70eb2", "Confidence" => 99.99950408935547, "ImageId" => "465f4e93-763e-51d0-b030-b9667a2d94b1", "BoundingBox" => Dict{String, Any}("Width" => 0.24444444477558136, "Left" => 0.5144444704055786, "Height" => 0.32592591643333435, "Top" => 0.15111111104488373)))])
end

list_collections_1_patch = @patch function list_collections(a...; b...)
    return Dict{String, Any}("CollectionIds" => Any["myphotos"])
end

list_faces_1_patch = @patch function list_faces(a...; b...)
    return Dict{String, Any}("Faces" => Any[Dict{String, Any}("FaceId" => "1c62e8b5-69a7-5b7d-b3cd-db4338a8a7e7", "Confidence" => 100, "ImageId" => "147fdf82-7a71-52cf-819b-e786c7b9746e", "BoundingBox" => Dict{String, Any}("Width" => 0.23999999463558197, "Left" => 0.5555559992790222, "Height" => 0.18000000715255737, "Top" => 0.336667001247406)), Dict{String, Any}("FaceId" => "29a75abe-397b-5101-ba4f-706783b2246c", "Confidence" => 100, "ImageId" => "147fdf82-7a71-52cf-819b-e786c7b9746e", "BoundingBox" => Dict{String, Any}("Width" => 0.22074100375175476, "Left" => 0.30963000655174255, "Height" => 0.16555599868297577, "Top" => 0.7066670060157776)), Dict{String, Any}("FaceId" => "38271d79-7bc2-5efb-b752-398a8d575b85", "Confidence" => 99.99829864501953, "ImageId" => "d5631190-d039-54e4-b267-abd22c8647c5", "BoundingBox" => Dict{String, Any}("Width" => 0.24222199618816376, "Left" => 0.3233329951763153, "Height" => 0.3234420120716095, "Top" => 0.5)), Dict{String, Any}("FaceId" => "3b01bef0-c883-5654-ba42-d5ad28b720b3", "Confidence" => 99.99210357666016, "ImageId" => "812d9f04-86f9-54fc-9275-8d0dcbcb6784", "BoundingBox" => Dict{String, Any}("Width" => 0.04747769981622696, "Left" => 0.37388700246810913, "Height" => 0.03555560111999512, "Top" => 0.2477779984474182)), Dict{String, Any}("FaceId" => "4839a608-49d0-566c-8301-509d71b534d1", "Confidence" => 99.99919891357422, "ImageId" => "812d9f04-86f9-54fc-9275-8d0dcbcb6784", "BoundingBox" => Dict{String, Any}("Width" => 0.07121659815311432, "Left" => 0.2937690019607544, "Height" => 0.05333330109715462, "Top" => 0.35666701197624207)), Dict{String, Any}("FaceId" => "70008e50-75e4-55d0-8e80-363fb73b3a14", "Confidence" => 99.99949645996094, "ImageId" => "d5631190-d039-54e4-b267-abd22c8647c5", "BoundingBox" => Dict{String, Any}("Width" => 0.24333299696445465, "Left" => 0.5155559778213501, "Height" => 0.3249259889125824, "Top" => 0.1513350009918213)), Dict{String, Any}("FaceId" => "7f5f88ed-d684-5a88-b0df-01e4a521552b", "Confidence" => 99.92639923095703, "ImageId" => "812d9f04-86f9-54fc-9275-8d0dcbcb6784", "BoundingBox" => Dict{String, Any}("Width" => 0.05044509842991829, "Left" => 0.7002969980239868, "Height" => 0.03777780011296272, "Top" => 0.18777799606323242)), Dict{String, Any}("FaceId" => "895b4e2c-81de-5902-a4bd-d1792bda00b2", "Confidence" => 99.99469757080078, "ImageId" => "812d9f04-86f9-54fc-9275-8d0dcbcb6784", "BoundingBox" => Dict{String, Any}("Width" => 0.07270029932260513, "Left" => 0.13946600258350372, "Height" => 0.05555560067296028, "Top" => 0.46333301067352295)), Dict{String, Any}("FaceId" => "8be04dba-4e58-520d-850e-9eae4af70eb2", "Confidence" => 99.99949645996094, "ImageId" => "465f4e93-763e-51d0-b030-b9667a2d94b1", "BoundingBox" => Dict{String, Any}("Width" => 0.24444399774074554, "Left" => 0.5144439935684204, "Height" => 0.3259260058403015, "Top" => 0.15111100673675537)), Dict{String, Any}("FaceId" => "908544ad-edc3-59df-8faf-6a87cc256cf5", "Confidence" => 99.9999008178711, "ImageId" => "3c731605-d772-541a-a5e7-0375dbc68a07", "BoundingBox" => Dict{String, Any}("Width" => 0.25222599506378174, "Left" => 0.3783380091190338, "Height" => 0.18888899683952332, "Top" => 0.2355560064315796)), Dict{String, Any}("FaceId" => "ff43d742-0c13-5d16-a3e8-03d3f58e980b", "Confidence" => 99.99909973144531, "ImageId" => "465f4e93-763e-51d0-b030-b9667a2d94b1", "BoundingBox" => Dict{String, Any}("Width" => 0.25, "Left" => 0.31888899207115173, "Height" => 0.33481499552726746, "Top" => 0.49333301186561584))])
end

search_faces_1_patch = @patch function search_faces(a...; b...)
    return Dict{String, Any}("FaceMatches" => Any[Dict{String, Any}("Similarity" => 99.97222137451172, "Face" => Dict{String, Any}("FaceId" => "8be04dba-4e58-520d-850e-9eae4af70eb2", "Confidence" => 99.99949645996094, "ImageId" => "465f4e93-763e-51d0-b030-b9667a2d94b1", "BoundingBox" => Dict{String, Any}("Width" => 0.24444399774074554, "Left" => 0.5144439935684204, "Height" => 0.3259260058403015, "Top" => 0.15111100673675537))), Dict{String, Any}("Similarity" => 97.04154968261719, "Face" => Dict{String, Any}("FaceId" => "29a75abe-397b-5101-ba4f-706783b2246c", "Confidence" => 100, "ImageId" => "147fdf82-7a71-52cf-819b-e786c7b9746e", "BoundingBox" => Dict{String, Any}("Width" => 0.22074100375175476, "Left" => 0.30963000655174255, "Height" => 0.16555599868297577, "Top" => 0.7066670060157776))), Dict{String, Any}("Similarity" => 95.94520568847656, "Face" => Dict{String, Any}("FaceId" => "908544ad-edc3-59df-8faf-6a87cc256cf5", "Confidence" => 99.9999008178711, "ImageId" => "3c731605-d772-541a-a5e7-0375dbc68a07", "BoundingBox" => Dict{String, Any}("Width" => 0.25222599506378174, "Left" => 0.3783380091190338, "Height" => 0.18888899683952332, "Top" => 0.2355560064315796)))], "SearchedFaceId" => "70008e50-75e4-55d0-8e80-363fb73b3a14")
end

search_faces_by_image_1_patch = @patch function search_faces_by_image(a...; b...)
    return Dict{String, Any}("SearchedFaceBoundingBox" => Dict{String, Any}("Width" => 0.25, "Left" => 0.31888890266418457, "Height" => 0.33481481671333313, "Top" => 0.4933333396911621), "SearchedFaceConfidence" => 99.9991226196289, "FaceMatches" => Any[Dict{String, Any}("Similarity" => 99.97036743164062, "Face" => Dict{String, Any}("FaceId" => "38271d79-7bc2-5efb-b752-398a8d575b85", "Confidence" => 99.99829864501953, "ImageId" => "d5631190-d039-54e4-b267-abd22c8647c5", "BoundingBox" => Dict{String, Any}("Width" => 0.24222199618816376, "Left" => 0.3233329951763153, "Height" => 0.3234420120716095, "Top" => 0.5)))])
end


text:
function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--nv"
            help = "Number of vintages"
            arg_type = Int
            default = 2
        "--nsrc"
            help = "Number of sources per vintage"
            arg_type = Int
            default = 8
        "--vm"
            help = "Type of virtual machine"
            arg_type = String
            default = "Standard_E8s_v3"
        "--nth"
            help = "Number of threads on a single node"
            arg_type = Int
            default = 4
        "--niter"
            help = "JRM iterations"
            arg_type = Int
            default = 8
        "--bs"
            help = "batchsize in JRM iterations"
            arg_type = Int
            default = 4
        "--snr"
            help = "SNR of noisy data"
            arg_type = Float64
            default = 0.0
        "--gamma"
            help = "Weighting on common component"
            arg_type = Float64
            default = 1.0
    end
    return parse_args(s)
end


text:
module CubicSplineNR

using Compat
using ForwardDiff

include("wrap_spline1.jl")
include("multi-splines.jl")

export CubicSpline, MutableCubicSpline, ImmutableCubicSpline,
    interp_level, update!
export MultiCubicSpline

# package code goes here

end # module


text:
# Test to rewrite the reconstructions functions with inputs the data instead of some noise level.

push!(LOAD_PATH, "../models")
push!(LOAD_PATH, "../")
using SparseInverseProblems
using SuperResModels
using TestCases
using Utils

using PyCall
@pyimport numpy as np
@pyimport matplotlib.pyplot as plt

# Containers to handle parameters
type StaticParameters
    # allowed frequencies
    f_c::Int64
    # spatial max val.
    x_max::Float64
    # ???
    num_x::Int64
    # Weights on the allowed frequencies.
    filter_x::Array{Float64}
end
type DynamicParameters
    # 2K+1 is the number of time samples.
    K::Int64
    # sampling rate
    tau::Float64
    # maxium allowed speed
    v_max::Float64
    # ????
    num_v::Int64
end

### Static parameters and container for use.
f_c = 20
x_max = 1.0
num_x = 10*f_c
filter_x =  ones(2*f_c+1);
static_parameters = StaticParameters(f_c, x_max, num_x, filter_x)

model_static = SuperResModels.Fourier1d(static_parameters.f_c, static_parameters.x_max, static_parameters.filter_x, static_parameters.num_x)
### Dynamic parameters and container for use.
K = 2
tau = 1.0/(K*1.0)
# the maximum velocity on the domain is given by x_max/2 + K tau v_max = 1
v_max = 0.5
num_v = 10
dynamic_parameters = DynamicParameters(K, tau, v_max, num_v)
model_dynamic = SuperResModels.DynamicFourier1d(model_static, dynamic_parameters.v_max, dynamic_parameters.tau, dynamic_parameters.K, dynamic_parameters.num_v*dynamic_parameters.K)

# We check the dimensions of each case.
(thetas, weights) = TestCases.cloud_1d_full(x_max, v_max, 0.1,10,K, tau, 1)
println(" thetas = ", thetas)
println(" weights = ", weights)

# Size of the dynamic reconstruction
target_dynamic = Utils.generate_target(model_dynamic, thetas, weights)
println(" Size of the dynamic reconstruction is: ", size(target_dynamic) )

# Size of a single static reconstruction
target_static = Utils.target_to_static(target_dynamic, model_dynamic.times)
#for k = 1: length(model_dynamic.times)
#thetas_t = Utils.to_static(thetas, model_dynamic.times[k], model_static.x_max)
#push!(target_static, Utils.generate_target(model_static,thetas_t, weights))
#println(" Size of a single static reconstruction is: ", size(target_static[k]) )
#end

res = zeros(size(target_dynamic))
for k = 1:length(model_dynamic.times)
    for j = 1:length(target_static[1])
	index = (k-1)*length(target_static[1]) + j;
        res[index] = target_dynamic[index] - target_static[k][j]
	println("res = ", res[(k-1)*length(target_static[1]) + j], "target static = ", target_static[k][j], "target dynamic = ", target_dynamic[(k-1)*length(target_static[1]) + j] )
    end
end
println(" The difference is: ", norm(res)/norm(target_dynamic))
#println(" Target dynamic = ", target_dynamic)
#println(" Target static = ", target_static)


text:
# # Basics: elevation data

# In this example, we will look at elevation data from the worldclim 2 data,
# crop it for Western Europe, and then change the resolution to aggregate the
# data. The first step is to get the worldclim layer for elevation:

using SimpleSDMLayers
using StatsPlots
import Statistics

#-

elevation = convert(Float32, SimpleSDMPredictor(WorldClim, Elevation))

# Thanks to the integration with Plots and StatsPlots, we can very rapidly
# visualize these data:

heatmap(elevation, c=:cividis, frame=:box)
xaxis!("Longitude")
yaxis!("Latitude")

# Let's also have a look at the density while we're at it:

density(elevation, frame=:zerolines, c=:grey, fill=(0, :grey, 0.2), leg=false)
xaxis!("Elevation")

# The next step is to clip the data to the region of interest. This requires a
# the coordinates of the bounding box as two tuples (for longitude and latitude)
# -- we can also make a quick heatmap to see what the region looks like:

elevation_europe = clip(elevation; left=-11.0, right=31.5, bottom=29.0, top=71.5)

#-

heatmap(elevation_europe, c=:cividis, aspectratio=1, frame=:box)

# The next step will be to coarsen these data, which requires to give the number
# of cells to merge alongside each dimension. This number of cells must be a
# divider of the grid size, which we can view with:

size(elevation_europe)

# In an ideal world, we could want to find a number of cells that is the same both
# for latitude and longitude, and one approach is to finagle our way into a
# correct grid by changing the clipping region.

# In this case, we will use a coarsening scale of `(5,5)`, which gives us a total
# of 25 cells in the aggregated result. Our aggregation function will be `mean`
# (so we report the average elevation across these cells):

elevation_europe_coarse = coarsen(elevation_europe, Statistics.mean, (5, 5))

# Once again, we can plot these data:

heatmap(elevation_europe_coarse, aspectratio=1, c=:cividis, frame=:box)

# Finally, we can compare our different clipping and approximations to the overall
# dataset:

density(elevation, frame=:zerolines, c=:grey, fill=(0, :grey, 0.5), lab="")
density!(elevation_europe, c=:black, lab="Raw data")
density!(elevation_europe_coarse, c=:darkgrey, lab="Average")
xaxis!("Elevation")


text:
"""
    module BandTools

Utility functions for manipulating band structures. These need to be defined early in the code so other Modules have acccess to them.
"""
module BandTools

using LinearAlgebra
using SpecialFunctions

function gaussian(de, smearing=0.01)
    return 0.5 * (erfc.(de/smearing))
end

function gaussian_derivative(de, smearing=0.01)
    return (-0.5 / sqrt(pi)) * exp.(-1.0*(de/smearing).^2)
end

"""
    function calc_fermi(eigs, weights, nelec, smearing = 0.01)

calculate fermi energy using bisection
"""
function calc_fermi(eigs, weights, nelec, smearing = 0.01)

#    println("calc_fermi eigs ", eigs)

    efermi_max = maximum(eigs[:])+0.1
    efermi_min = minimum(eigs[:])-0.1

    efermi = (efermi_max+efermi_min)/2.0

    norm = sum(weights)

    n2= Int64(round(nelec/2))

#    println(size(eigs), " n2 $n2")

    if  abs(nelec - n2*2) < 1e-7
        if length(size(eigs)) == 2 && size(eigs)[2] > 1
            gap  = minimum(eigs[:,n2+1]) - maximum(eigs[:,n2])   #large gaps become numerically unstable with small smearing
        else
            gap  = minimum(eigs[n2+1]) - maximum(eigs[n2])   #large gaps become numerically unstable with small smearing
        end

        if gap > 0.05  #if gap is positive/large and we have an integer number of electrons
            smearing = max(smearing, gap / 10.0)
        end
        #    println("gap $gap smearing $smearing n2 $n2")

    end

    for iter = 1:35

        efermi = (efermi_max+efermi_min)/2.0

        occ = gaussian.(eigs.-efermi, smearing)

        n = sum(sum(occ, dims=2).*weights)/norm
#        println(iter, "  ", n,"  " ,  efermi)
        if n > (nelec/2.0)
            efermi_max = efermi
        else
            efermi_min = efermi
        end
    end

    return efermi

end

"""
    function band_energy(eigs, weights, nelec, smearing = 0.01; returnk=false, returnocc=false, returnef=false, returnboth=false)

Calculate band energy. Has options for additional return variables. Calculates fermi energy internally.
"""
function band_energy(eigs, weights, nelec, smearing = 0.01; returnk=false, returnocc=false, returnef=false, returnboth=false)

    efermi = calc_fermi(eigs, weights, nelec, smearing)
    norm = sum(weights)

    occ = gaussian.(eigs.-efermi, smearing)

    ek = sum(eigs.*occ, dims=2).*weights * 2.0 / norm
    #    energy = sum(sum(eigs.*occ, dims=2).*weights) * 2.0 / norm
    energy=sum(ek)

    if returnk
        return energy,ek
    elseif returnocc
        return energy, occ
    elseif returnef
        return energy, efermi
    elseif returnboth
        return energy, efermi, occ
    else
        return energy
    end

end

"""
    function smearing_energy(eigs, weights, efermi, smearing = 0.01)

Smearing contribution to total energy. If you don't take this int account in metals or
small gap semiconductors, your energy is not variational
in the correct way, and things like forces and stresses become wrong.
"""
function smearing_energy(eigs, weights, efermi, smearing = 0.01)

#    efermi = calc_fermi(eigs, weights, nelec, smearing)
    norm = sum(weights)

    dgauss = gaussian_derivative(eigs.-efermi, smearing)

    ek = smearing * sum(dgauss, dims=2).*weights * 2.0 / norm

    energy=sum(ek)

    return energy

end

end


text:
# # [Basic Interaction](@id basic_interaction)
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/example_placeholder.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/example_placeholder.ipynb)

# In this example page I will show how to work with basic function that are meant for working with this registration algorithms.

# ## First how to generate some images that we can register
# We can do this with the [`gen_init`](@ref) for example.

## using the package
using LAP_julia;

## default arguments:
img, imgw, flow = gen_init();

## get a chess board image with a uniform flow:
img_chess, imgw_chess, flow_uniform = gen_init(:chess, :uniform, flow_args=[1 + 1im]);

# See also: [`gen_chess`](@ref), [`gen_quad_flow`](@ref), [`gen_quad_flow`](@ref), [`gen_lena`](@ref), [`gen_chess`](@ref), [`gen_anhir`](@ref)

# ## Next we would like to display these generated images and flows
# For that we will use the [`imgshow`](@ref), [`showflow`](@ref) and [`imgoverlay`](@ref) functions.

# See a single image:
imgshow(img)
#-----------
# See the default random flow:
showflow(flow)
#-----------
# See the uniform flow:
showflow(flow_uniform)
#-----------
# See the differences between 2 images:
imgoverlay(img, imgw)

# See also: [`imgshowflow`](@ref), [`warp_imgshowflow`](@ref), [`gen_quad_flow`](@ref), [`gen_lena`](@ref), [`gen_chess`](@ref), [`gen_anhir`](@ref),


text:
#=
P13 (**) Run-length encoding of a list (direct solution).
Implement the so-called run-length encoding data compression method directly. I.e. don't explicitly create the sublists containing the duplicates, as in problem P09, but only count them. As in problem P11, simplify the result list by replacing the singleton lists (1 X) by X.

Example:
* (encode-direct '(a a a a b c c a a d e e e e))
((4 A) B (2 C) (2 A) D (4 E))
=#

function _encode_direct(coll, current_pair, pairs)
    if isempty(coll)
        vcat(pairs, Vector{}[current_pair])
    else
        count, item = current_pair
        if coll[1] == item
            _encode_direct(coll[2:end], [count + 1, item], pairs)
        else
            _encode_direct(coll[2:end], [1, coll[1]], vcat(pairs, Vector{}[current_pair]))
        end
    end
end

function encode_direct(coll)
    encoded = isempty(coll) ? coll : _encode_direct(coll[2:end], [1, coll[1]], [])
    map(pair -> pair[1] == 1 ? pair[2] : pair, encoded)
end

encoded = encode_direct([:a, :a, :a, :a, :b, :c, :c, :a, :a, :d, :e, :e, :e, :e])
@assert encoded[1] == [4, :a]
@assert encoded[2] == :b
@assert encoded[3] == [2, :c]
@assert encoded[4] == [2, :a]
@assert encoded[5] == :d
@assert encoded[6] == [4, :e]

@assert encode_direct([]) == []

println("Tests passed: JL-13.jl")


text:
function visualize!(vis, model, q;
		Δt = 0.1, scenario = :vertical)

    r_foot = 0.05
    r_leg = 0.5 * r_foot

	default_background!(vis)

    setobject!(vis["body"], Sphere(Point3f0(0),
        convert(Float32, 0.1)),
        MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))

    setobject!(vis["foot"], Sphere(Point3f0(0),
        convert(Float32, r_foot)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

    n_leg = 100
    for i = 1:n_leg
        setobject!(vis["leg$i"], Sphere(Point3f0(0),
            convert(Float32, r_leg)),
            MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))
    end

    p_leg = [zeros(3) for i = 1:n_leg]
    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(q)
        p_body = [q[t][1], 0.0, q[t][2]]
        p_foot = [kinematics(model, q[t])[1], 0.0, kinematics(model, q[t])[2]]

        q_tmp = Array(copy(q[t]))
        r_range = range(0, stop = q[t][4], length = n_leg)
        for i = 1:n_leg
            q_tmp[4] = r_range[i]
            p_leg[i] = [kinematics(model, q_tmp)[1], 0.0, kinematics(model, q_tmp)[2]]
        end
        q_tmp[4] = q[t][4]
        p_foot = [kinematics(model, q_tmp)[1], 0.0, kinematics(model, q_tmp)[2]]

        z_shift = [0.0; 0.0; r_foot]

        MeshCat.atframe(anim, t) do
            settransform!(vis["body"], Translation(p_body + z_shift))
            settransform!(vis["foot"], Translation(p_foot + z_shift))

            for i = 1:n_leg
                settransform!(vis["leg$i"], Translation(p_leg[i] + z_shift))
            end
        end
    end

	# if scenario == :vertical
	# 	settransform!(vis["/Cameras/default"],
	# 		compose(Translation(0.0, 0.5, -1.0),LinearMap(RotZ(-pi / 2.0))))
	# end

    MeshCat.setanimation!(vis, anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# Extra stuff

newref!(::Type{Ref{HSA.Region}}, value) = Ref{HSA.Region}(HSA.Region(value))

# Memory
memset!(ref::Ref{T}, value::T) where T<:Integer = _memset!(ref, value, T)

function _memset!(ref, value, ::Type{T}; len=1) where T
    ptr = Base.unsafe_convert(Ptr{T}, ref)
    for i in 1:len
        unsafe_store!(ptr, value, i)
    end
end

function _memset!(ref::Ref{Ptr{Nothing}}, value; len=1)
    ptr = Base.unsafe_convert(Ptr{UInt8}, ref[])
    value = UInt8(value)
    for i in 1:len
        unsafe_store!(ptr, value, i)
    end
end

# Atomic store using LLVM intrinsics
# Necessary for writing the AQL packet header to the queue
# prior to launching a kernel.
@eval atomic_store_n!(x::Ptr{UInt16}, v::UInt16) =
    Base.llvmcall($"""
    %ptr = inttoptr i$(Sys.WORD_SIZE) %0 to i16*
    store atomic i16 %1, i16* %ptr release, align 64
    ret void
    """, Cvoid, Tuple{Ptr{UInt16}, UInt16},
    Base.unsafe_convert(Ptr{UInt16}, x), v)

# Overloads for interface functions

function getinfo(agent::Agent, attribute::HSA.AgentInfo,
                 value::Union{Vector,Base.RefValue,String})
    #TODO: allocation/create Refs here
    # based on value of AgentInfo
    HSA.agent_get_info(agent, attribute, value)
end
getinfo(agent::Agent, info::HSA.AMDAgentInfo, ret::Union{String, Base.RefValue, Vector}) =
    getinfo(agent, reinterpret(HSA.AgentInfo, info), ret)

function getinfo(isa::HSA.ISA, attribute::HSA.ISAInfo,
                 value::Union{Vector,Base.RefValue,String})
    HSA.isa_get_info_alt(isa, attribute, value)
end

function getinfo(exsym::HSA.ExecutableSymbol, attribute::HSA.ExecutableSymbolInfo,
                 value::Union{Vector,Base.RefValue,String})
    HSA.executable_symbol_get_info(exsym, attribute, value)
end


text:
using Primes;
using DelimitedFiles;
using SHA;

function primeHash(value::BigInt)
  _value = unsafe_wrap(Vector{UInt8}, Ptr{UInt8}(value.d), (value.size*8,))
  res = sha256(_value)
  bn = BigInt(0)
  ccall((:__gmpz_import, :libgmp), Cvoid, (Ref{BigInt}, Csize_t, Cint, Csize_t, Cint, Csize_t, Ptr{Cvoid}), bn, 4, 1, 8, 0, 0, res)
  while !isprime(bn)
    bn-=1
  end
  return bn
end

function primeN(bn::BigInt)
  while !isprime(bn)
    bn-=1
  end
  return bn
end

let
  randstate = zeros(UInt8,32)
  global bigrand
  ccall((:__gmp_randinit_default,:libgmp),Cvoid,(Ptr{Cvoid},),randstate)
  function rand(min::BigInt, max::BigInt)
    delta = max - min
    x = BigInt()
    ccall((:__gmpz_urandomm,:libgmp),Cvoid,(Ref{BigInt},Ptr{Cvoid},Ref{BigInt}),x,randstate,delta)
    return x+min
  end

  bigrand(max::BigInt) = bigrand(BigInt(0), max)
end

function normalize(rem::BigInt, divisor::BigInt)
  return powermod(rem, -1, divisor)
end

function test()
  mantice = rand(UInt8, 2^20) .% 2
  primes = map(x -> primeN(x<<10), BigInt(1):BigInt(2^20))
  println(length(unique(primes))-2^20)
end

test()


text:
module EBayes

using Reexport

@reexport using StatsBase
@reexport using Distributions

using StructArrays

import Base:length, rand, size
import Base.Broadcast: broadcastable
import Base: iterate,
             getindex,
             firstindex,
             lastindex,
             eltype,
             zero,
             zeros

import StatsBase:predict,
                 fit,
                 dof,
                 leverage

import StatsBase:response

using Statistics
import Statistics:var

import Distributions:location

import MLJBase, MLJ # to avoid conflicts
import MLJBase:predict_mean
import MLDataPattern:kfolds,FoldsView,shuffleobs
using NLSolversBase
using Optim

include("ebayes_types.jl")
include("ebayes_samples.jl")
include("predict.jl")
include("shrinkage_locations.jl")
include("sure.jl")
include("ebcf.jl")
include("simulations.jl")

export EBayesSample,
       NormalSample,
       EBayesSamples,
       AbstractNormalSamples,
       NormalSamples,
       NormalEBayesSimulationResult,
       var, #reexport from Statistics
       predict_mean, #reexport from MLJBase
       FixedLocation,
       GrandMeanLocation,
       SURE,
       EBayesCrossFit

end # module


text:
@testset "Selection operators" begin

ss = RectSearchSpace(1, (0.0, 10.0))
fake_problem = FunctionBasedProblem(x -> 0.0, "test_problem", MinimizingFitnessScheme, ss)
fake_pop = reshape(collect(1.0:10.0), (1, 10))

@testset "SimpleSelector" begin
    @test popsize(fake_pop) == 10
    sel = BlackBoxOptim.SimpleSelector()
    for i in 1:NumTestRepetitions
        numSamples = rand(1:8)
        sampled = BlackBoxOptim.select(sel, fake_pop, numSamples)

        @test length(sampled) == numSamples

        # All sampled indices are indices into the population
        @test all(index -> in(index, 1:popsize(fake_pop)), sampled)
    end
end

@testset "RadiusLimitedSelector" begin
    sel = BlackBoxOptim.RadiusLimitedSelector(20)

    for i in 1:NumTestRepetitions
        numSamples = rand(1:sel.radius)
        sampled = BlackBoxOptim.select(sel, fake_pop, numSamples)

        @test length(sampled) == numSamples

        # All sampled indices are indices into the population
        @test all(index -> in(index, 1:popsize(fake_pop)), sampled)

        mini, maxi = minimum(sampled), maximum(sampled)
        if (maxi - mini) > max(numSamples+2, sel.radius)
            @test mini + popsize(fake_pop) <= maxi + sel.radius
        end
    end
end

@testset "Tournament" begin
    sel = BlackBoxOptim.TournamentSelector(MinimizingFitnessScheme, 3)
    fake_pop = FitPopulation(reshape(collect(1.0:10.0), (1, 10)),
                             NaN, collect(1.0:10.0))

    @testset "tournament()" begin
        for i in 1:NumTestRepetitions
            tourSize = rand(0:popsize(fake_pop))
            cand_ixs = findall(bitrand(popsize(fake_pop)))

            winner_ix = BlackBoxOptim.tournament(sel, fake_pop, cand_ixs)

            if isempty(cand_ixs)
                @test winner_ix == 0
            else
                @test in(winner_ix, 1:popsize(fake_pop))
                @test fitness(fake_pop, winner_ix) == minimum(fake_pop.fitness[cand_ixs])
            end
        end
    end

    @testset "select()" begin
        for i in 1:NumTestRepetitions
            numSamples = rand(1:3)
            sampled = BlackBoxOptim.select(sel, fake_pop, numSamples)

            @test length(sampled) == numSamples

            # All sampled indices are indices into the population
            @test all(index -> in(index, 1:popsize(fake_pop)), sampled)
        end
    end
end

end


text:
using Pkg
Pkg.add([
    "PyCall","IJulia"
])

max_stars_repo_name: kashif-se/docker-ipyjula-notebooks

max_issues_repo_name: kashif-se/docker-ipyjula-notebooks

max_forks_repo_name: kashif-se/docker-ipyjula-notebooks


text:
module HSA

import ..AMDGPU: libhsaruntime

using CEnum

include("libhsa_common.jl")
include("libhsa_enums.jl")
include("libhsa_types.jl")
include("libhsa_api.jl")

end


text:
###### conversion from optical stripe data to global optical gradient data #####
# 1. convert data to global data using inner loop data with iSARA kernel
# using temp_profile_si
# l is lengthscale, θ are rescaling parameters of iSARA kernel
# nout is number of outer temperatures to record gradients at
# T_fraction is lowest temperature as a fraction of T_max
# k is kernel used for GP inference
# only for 2d data, and backward compatability
function stripe_to_global(x::AbstractVector, y::AbstractVector{<:AbstractVector},
                          σ::Real, k, P::TemperatureProfile, T_max::Real, log10_τ::Real,
                          relevant_T, input_noise::Union{Val{true}, Val{false}} = Val(true))
    stripe_to_global(x, y, σ, k, P, (T_max, log10_τ), relevant_T, input_noise)[2:3]
end

# condition is tuple of (T_max, τ, composition, ...) conditions
function stripe_to_global(x::AbstractVector, y::AbstractVector{<:AbstractVector},
                        σ::Real, k, P::TemperatureProfile, condition::NTuple,
                        relevant_T, input_noise::Union{Val{true}, Val{false}} = Val(true))
    all(==(length(x)), length.(y)) || throw(DimensionMismatch("x and elements of y do not have same lengths: length(x) = $(length(x)) and length.(y) = $(length.(y))"))
    T_max, log10_τ = condition[1:2]
    other_conditions = condition[3:end]
    G = Gaussian(k) # define prior GP
    temperature_processes = Vector{Gaussian}([G for _ in 1:length(y)])
    @threads for j in eachindex(y)
        C = get_temperature_process(G, x, y[j], σ, P, T_max, log10_τ, input_noise)
        temperature_processes[j] = C
    end
    C = temperature_processes[1]
    Tout = relevant_T(T_max, log10_τ) # IDEA: could add composition as dimension
    conditions = tuple.(Tout, log10_τ, other_conditions...)
    gradients, uncertainty = get_global_data(temperature_processes, Tout)
    return conditions, gradients, uncertainty
end

# converts stripe of xrd measurements to integrated gradients
# x is position
# Y is integrated XRD image
function xrd_to_global(x::AbstractVector, Y::AbstractMatrix, σ, k,
                       P::TemperatureProfile, condition::NTuple, relevant_T,
                       input_noise::Union{Val{true}, Val{false}} = Val(true); rtol::Real = 1e-2)
    y = svd_coefficients(Y, rtol) # get coefficients of Y in SVD basis up to relative tolerance 1e-2
    stripe_to_global(x, y, σ, k, P, condition, relevant_T, input_noise)
end

# uses data optical data as a function of position x and the inverse
# temperature profile to compute a Gaussian process as a function of temperature
# G is prior process of position, σ is noise standard deviation
# x is vector of positions corresponding to reflectance measurements y
function get_temperature_process(G::Gaussian, x::AbstractVector,
                                 y::AbstractVector, σ::Real,
                                 P::TemperatureProfile, T_max::Real, log10_τ::Real,
                                 input_noise::Val{false})
    C = conditional(G, x, Gaussian(y, σ^2), tol = 1e-6) # condition Gaussian process in position
    inv_profile = inverse_profile(P, T_max, log10_τ)
    input_transformation(C, inv_profile) # transform the input of conditional process
end

# same as above but takes into account the uncertainty in the temperature
function get_temperature_process(G::Gaussian, x::AbstractVector,
                                 y::AbstractVector, σ_out::Real,
                                 P::TemperatureProfile, T_max::Real, log10_τ::Real,
                                 input_noise::Val{true})
    # first, get regular GP w.r.t. T
    C = get_temperature_process(G, x, y, σ_out, P, T_max, log10_τ, Val(false))
    # with input noise, adjust GP with generalized NIGP model
    T = profile(P, T_max, log10_τ).(x) # maps position of measured optical data to temperature
    d = (t->ForwardDiff.derivative(mean(C), t)).(T)
    var_in = temperature_uncertainty(P, T_max, log10_τ).(x) # input uncertainty corresponding to positions
    Σ = @. var_in * d^2 + σ_out^2
    Σ = Diagonal(Σ)
    C = conditional(G, x, Gaussian(y, Σ), tol = 1e-6) # calculate new conditional process with adjusted noise variance
    inv_profile = inverse_profile(P, T_max, log10_τ)
    input_transformation(C, inv_profile) # transform the input to temperature domain
end

# converts Gaussian processes of temperature with optical target data
# to integrated gradient information of the processes at the given outer_temperature
function get_global_data(temperature_processes::AbstractVector{<:Gaussian},
                         outer_temperature::AbstractVector{<:Real},
                         temperature_domain::NTuple{2, <:Real} = (0, 1400))
    # normalize temperature input to have a normalized scale for derivatives
    ut = unit_transform(temperature_domain...)
    iut = inv_unit_transform(temperature_domain...)
    unit_outer_temperature = ut.(outer_temperature)
    nout = length(outer_temperature)
    DT = zeros(nout, length(temperature_processes))
    UT = similar(DT) # uncertainty
    @threads for i in eachindex(temperature_processes) # for (i, G) in enumerate(temperature_processes)
        C = input_transformation(temperature_processes[i], iut) # inverse of temperature unit-scaling
        D = GaussianDistributions.derivative(C) # take the derivative w.r.t. T
        DT[:, i] = mean(D).(unit_outer_temperature) # record mean and var of derivative of each process
        UT[:, i] = var(D).(unit_outer_temperature)
    end
    # euclidean norm of temperature gradients of optical coefficients
    d, u = zeros(nout), zeros(nout)
    for i in 1:nout
        di, ui = @views DT[i, :], UT[i, :]
        d[i] = norm(di) # and first-order uncertainty propagation through "norm",
        u[i] = dot(ForwardDiff.gradient(norm, di).^2, ui)
    end
    # @. u = sqrt(u) # convert variance to std?
    return d, u # gradient values and their uncertainties
end


text:
"""
Snape's Riddle - Author JK Rowling
Danger lies before you, while safety lies behind,
Two of us will help you, whichever you would find,
One among us seven will let you move ahead,
Another will transport the drinker back instead,
Two among our number hold only nettle wine,
Three of us are killers, waiting hidden in line.
Choose, unless you wish to stay here for evermore,
To help you in your choice, we give you these clues four:
First, however slyly the poison tries to hide
You will always find some on nettle wine’s left side;
Second, different are those who stand at either end,
But if you would move onwards, neither is your friend;
Third, as you see clearly, all are different size,
Neither dwarf nor giant holds death in their insides;
Fourth, the second left and the second on the right
Are twins once you taste them, though different at first sight.
"""
julia> using AbstractLogic
Start the repl in command prompt by typing `=`.

# Codes:
# 1 is netted wine - NW
# 2 move ahead     - MA
# 3 move back      - MB
# 4 poison         - PO

# Seven Potions a,b,c,d,e,f,g each have 1 of 4 possible values
# a is left of b which is left of c etc.
a, b, c, d, e, f, g  ∈  NW, MA, MB, PO [clear]

# Two among our number hold only nettle wine, #code 1
{{i}} == 'NW' {{2}}
# Equivalent to logicalparse(["a|b|c|d|e|f|g == 1 {2}"], logicset=hpset)

# One among us seven will let you move ahead, #code 2
{{i}} == 'MA' {{1}}

#Another will transport the drinker back instead, #code 3
{{i}} == 'MB' {{1}}

# Three of us are killers #code 4
{{i}} == 'PO' {{3}}
# This last specification does not reduce the feasible set

# First, however slyly the poison tries to hide
# You will always find some on nettle wine’s left side;
# So if b is poison then a is nettle wine
# The i-1! indicates that outofbounds matches are considered "false"
{{i}} == 'NW' ==> {{i-1!}} == 'PO'
# 60 possible outcomes

# Second, different are those who stand at either end,
a != g

# But if you would move onwards, neither is your friend;
a,g != 'MA'
# 30 possible outcomes

# Third, as you see clearly, all are different size,
# Neither dwarf nor giant holds death in their insides;
# c dwarf, f are giant
c,f != 'PO'

# Fourth, the second left and the second on the right
# Are twins once you taste them, though different at first sight.
b == f

search ='MA'
# c is a match, for the third bottle

search ='MB'
# g is a match, for the last bottle


text:
# Tests for Revenue DEA Models
@testset "RevenueDEAModel" begin

    ## Test Revenue DEA Model with Cooper et al. (2007)
    # Test agains book results
    X = [3 2; 1 3; 4 6]
    Y = [3; 5; 6]
    P = [6; 6; 6]

    dearevenuecooper = dearevenue(X, Y, P, rts = :CRS)
    @test efficiency(dearevenuecooper, :Economic)   ≈ [0.9; 1; 0.6] atol = 1e-3
    @test efficiency(dearevenuecooper, :Technical)  ≈ [0.9; 1; 0.6  ] atol = 1e-3
    @test efficiency(dearevenuecooper, :Allocative) ≈ [1; 1; 1] atol = 1e-3

    ## Test Revenue DEA Model with Zofío and Prieto (2006) data.
    # Test agains results with R
    X = [5 3; 2 4; 4 2; 4 8; 7 9]
    Y = [7 4; 10 8; 8 10; 5 4; 3 6]
    P = [3 2; 3 2; 3 2; 3 2; 3 2.0]

    # Revnue CRS
    dearevenuecrs = dearevenue(X, Y, P, rts = :CRS)
    @test efficiency(dearevenuecrs) ≈ [0.4915254237;
                                    1.000;
                                    1.000;
                                    0.250;
                                    0.1735537190]
    @test efficiency(dearevenuecrs, :Technical) ≈ [0.6364;
                                1.000;
                                1.000;
                                0.250;
                                0.2609]  atol = 1e-3
    @test efficiency(dearevenuecrs, :Allocative) ≈ [0.7723970944;
                                1.000;
                                1.000;
                                1.000;
                                0.6652892562]

    # Revenue VRS
    dearevenuevrs = dearevenue(X, Y, P, rts = :VRS)
    @test efficiency(dearevenuevrs) ≈ [0.6444444444 ;
                                    1.000;
                                    1.000;
                                    0.500;
                                    0.4565217391]
    @test efficiency(dearevenuevrs, :Technical) ≈ [0.7777777778 ;
                                1.000;
                                1.000;
                                0.500 ;
                                0.6000000000]
    @test efficiency(dearevenuevrs, :Allocative) ≈ [0.8285714286;
                                1.000;
                                1.000;
                                1.000;
                                0.7608695652]

    # Check defaults
    @test efficiency(dearevenue(X, Y, P)) == efficiency(dearevenuevrs)
    @test efficiency(dearevenuevrs, :Economic) == efficiency(dearevenuevrs)

    # Print
    show(IOBuffer(), dearevenuecooper)

    # Test errors
    @test_throws ErrorException dearevenue([1; 2 ; 3], [4 ; 5], [1; 1; 1]) #  Different number of observations
    @test_throws ErrorException dearevenue([1; 2; 3], [4; 5; 6], [1; 2; 3], rts = :Error) # Invalid returns to scale
    @test_throws ErrorException dearevenue([1; 2; 3], [4; 5; 6], [1; 2; 3; 4]) # Different number of observation in prices
    @test_throws ErrorException dearevenue([1; 2; 3], [4 4; 5 5; 6 6], [4 4 4; 5 5 5; 6 6 6]) # Different number of output prices and outputs

end


text:
using BinaryProvider
import Libdl

const extension = Sys.isapple() ? "dylib" : Sys.islinux() ? "so" : Sys.iswindows() ? "dll" : ""

print_error() = error(
    "FastTransforms could not be properly installed.\n Please check that you have all dependencies installed. " *
    "Sample installation of dependencies:\n" *
    print_platform_error(platform_key_abi())
)

print_platform_error(p::Platform) = "On $(BinaryProvider.platform_name(p)), please consider opening a pull request to add support to build from source.\n"
print_platform_error(p::MacOS) = "On MacOS\n\tbrew install libomp fftw mpfr\n"
print_platform_error(p::Linux) = "On Linux\n\tsudo apt-get install libomp-dev libblas-dev libopenblas-base libfftw3-dev libmpfr-dev\n"
print_platform_error(p::Windows) = "On Windows\n\tvcpkg install openblas:x64-windows fftw3[core,threads]:x64-windows mpir:x64-windows mpfr:x64-windows\n"

ft_build_from_source = get(ENV, "FT_BUILD_FROM_SOURCE", "false")
if ft_build_from_source == "true"
    make = Sys.iswindows() ? "mingw32-make" : "make"
    flags = Sys.isapple() ? "FT_USE_APPLEBLAS=1" : Sys.iswindows() ? "FT_FFTW_WITH_COMBINED_THREADS=1" : ""
    script = """
        set -e
        set -x
        if [ -d "FastTransforms" ]; then
            cd FastTransforms
            git fetch
            git checkout master
            git pull
            cd ..
        else
            git clone https://github.com/MikaelSlevinsky/FastTransforms.git FastTransforms
        fi
        cd FastTransforms
        $make assembly
        $make lib $flags
        cd ..
        mv -f FastTransforms/libfasttransforms.$extension libfasttransforms.$extension
    """
    try
        run(`bash -c $(script)`)
    catch
        print_error()
    end
    println("FastTransforms built from source.")
else
    println("FastTransforms using precompiled binaries.")
end


text:
@defcomp co2emissions begin
    region = Index()

    e_globalCO2emissions = Variable(index=[time], unit="Mtonne/year")
    e0_baselineCO2emissions = Parameter(index=[region], unit="Mtonne/year")
    e_regionalCO2emissions = Variable(index=[time,region], unit="Mtonne/year")
    er_CO2emissionsgrowth = Parameter(index=[time,region], unit="%")

    function run_timestep(p, v, d, t)

        # eq.4 in Hope (2006) - regional CO2 emissions as % change from baseline
        for r in d.region
            v.e_regionalCO2emissions[t,r] = p.er_CO2emissionsgrowth[t,r] * p.e0_baselineCO2emissions[r] / 100
        end
        # eq. 5 in Hope (2006) - global CO2 emissions are sum of regional emissions
        v.e_globalCO2emissions[t] = sum(v.e_regionalCO2emissions[t,:])
    end
end


text:
# ---
# title: 1140. Stone Game II
# id: problem1140
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/stone-game-ii/description/>
# hidden: true
# ---
#
# Alice and Bob continue their games with piles of stones.  There are a number
# of piles  **arranged in a row** , and each pile has a positive integer number
# of stones `piles[i]`.  The objective of the game is to end with the most
# stones.
#
# Alice and Bob take turns, with Alice starting first.  Initially, `M = 1`.
#
# On each player's turn, that player can take **all the stones** in the
# **first** `X` remaining piles, where `1 <= X <= 2M`.  Then, we set `M = max(M,
# X)`.
#
# The game continues until all the stones have been taken.
#
# Assuming Alice and Bob play optimally, return the maximum number of stones
# Alice can get.
#
#
#
# **Example 1:**
#
#
#
#     Input: piles = [2,7,9,4,4]
#     Output: 10
#     Explanation:  If Alice takes one pile at the beginning, Bob takes two piles, then Alice takes 2 piles again. Alice can get 2 + 4 + 4 = 10 piles in total. If Alice takes two piles at the beginning, then Bob can take all three piles left. In this case, Alice get 2 + 7 = 9 piles in total. So we return 10 since it's larger.
#
#
# **Example 2:**
#
#
#
#     Input: piles = [1,2,3,4,5,100]
#     Output: 104
#
#
#
#
# **Constraints:**
#
#   * `1 <= piles.length <= 100`
#   * `1 <= piles[i] <= 104`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@doc raw"""
    Stiefel{n,k,𝔽} <: Manifold

The Stiefel manifold consists of all $n × k$, $n ≥ k$ unitary matrices, i.e.

````math
\operatorname{St}(n,k)\{ p ∈ 𝔽^{n × k} : p^{\mathrm{H}}p = I_k \},
````

where $𝔽 ∈ \{ℝ, ℂ\}$,
$\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian, and
$I_k ∈ ℝ^{k × k}$ denotes the $k × k$ identity matrix.

The tangent space at a point $p ∈ \mathcal M$ is given by

````math
T_p \mathcal M = \{ X ∈ 𝔽^{n × k} : p^{\mathrm{H}}X + X^{\mathrm{H}}p = 0_k\},
````

where $0_k$ is the $k × k$ zero matrix.

This manifold is modeled as an embedded manifold to the [`Euclidean`](@ref), i.e.
several functions like the [`inner`](@ref inner(::Euclidean, ::Any...)) product and the
[`zero_tangent_vector`](@ref zero_tangent_vector(::Euclidean, ::Any...)) are inherited from the embedding.

The manifold is named after
[Eduard L. Stiefel](https://en.wikipedia.org/wiki/Eduard_Stiefel) (1909–1978).

# Constructor
    Stiefel(n, k, field = ℝ)

Generate the (real-valued) Stiefel manifold of $n × k$ dimensional orthonormal matrices.
"""
struct Stiefel{n,k,𝔽} <: AbstractEmbeddedManifold{DefaultIsometricEmbeddingType} end

Stiefel(n::Int, k::Int, field::AbstractNumbers = ℝ) = Stiefel{n,k,field}()

function allocation_promotion_function(
    M::Stiefel{n,k,ℂ},
    f,
    args::Tuple,
) where {n,k}
    return complex
end

@doc raw"""
    check_manifold_point(M::Stiefel, p; kwargs...)

Check whether `p` is a valid point on the [`Stiefel`](@ref) `M`=$\operatorname{St}(n,k)$, i.e. that it has the right
[`AbstractNumbers`](@ref) type and $p^{\mathrm{H}}p$ is (approximately) the identity, where $\cdot^{\mathrm{H}}$ is the
complex conjugate transpose. The settings for approximately can be set with `kwargs...`.
"""
function check_manifold_point(M::Stiefel{n,k,𝔽}, p; kwargs...) where {n,k,𝔽}
    mpv = invoke(check_manifold_point, Tuple{supertype(typeof(M)), typeof(p)}, M, p; kwargs...)
    mpv === nothing || return mpv
    c = p' * p
    if !isapprox(c, one(c); kwargs...)
        return DomainError(
            norm(c - one(c)),
            "The point $(p) does not lie on $(M), because x'x is not the unit matrix.",
        )
    end
end

@doc raw"""
    check_tangent_vector(M::Stiefel, p, X; check_base_point = true, kwargs...)

Checks whether `X` is a valid tangent vector at `p` on the [`Stiefel`](@ref)
`M`=$\operatorname{St}(n,k)$, i.e. the [`AbstractNumbers`](@ref) fits and
it (approximately) holds that $p^{\mathrm{H}}X + X^{\mathrm{H}}p = 0$.
The optional parameter `check_base_point` indicates, whether to call [`check_manifold_point`](@ref)  for `p`.
The settings for approximately can be set with `kwargs...`.
"""
function check_tangent_vector(
    M::Stiefel{n,k,𝔽},
    p,
    X;
    check_base_point = true,
    kwargs...,
) where {n,k,𝔽}
    if check_base_point
        mpe = check_manifold_point(M, p; kwargs...)
        mpe === nothing || return mpe
    end
    mpv = invoke(
        check_tangent_vector,
        Tuple{supertype(typeof(M)), typeof(p), typeof(X)},
        M,
        p,
        X;
        check_base_point = false, # already checked above
        kwargs...
    )
    mpv === nothing || return mpv
    if !isapprox(p' * X + X' * p, zeros(k, k); kwargs...)
        return DomainError(
            norm(p' * X + X' * p),
            "The matrix $(X) is does not lie in the tangent space of $(p) on the Stiefel manifold of dimension ($(n),$(k)), since x'v + v'x is not the zero matrix.",
        )
    end
end

decorated_manifold(M::Stiefel{N,K,𝔽}) where {N,K,𝔽} = Euclidean(N, K; field = 𝔽)

embed!(::Stiefel, q, p) = (q .= p)

embed!(::Stiefel, Y, p, X) = (Y .= X)

@doc raw"""
    exp(M::Stiefel, p, X)

Compute the exponential map on the [`Stiefel`](@ref)`{n,k,𝔽}`() manifold `M`
emanating from `p` in tangent direction `X`.

````math
\exp_p X = \begin{pmatrix}
   p\\X
 \end{pmatrix}
 \operatorname{Exp}
 \left(
 \begin{pmatrix} p^{\mathrm{H}}X & - X^{\mathrm{H}}X\\
 I_n & p^{\mathrm{H}}X\end{pmatrix}
 \right)
\begin{pmatrix}  \exp( -p^{\mathrm{H}}X) \\ 0_n\end{pmatrix},
````

where $\operatorname{Exp}$ denotes matrix exponential,
$\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian, and $I_k$ and
$0_k$ are the identity matrix and the zero matrix of dimension $k × k$, respectively.
"""
exp(::Stiefel, ::Any...)

function exp!(M::Stiefel{n,k}, q, p, X) where {n,k}
    return copyto!(
        q,
        [p X] *
        exp([p'X -X' * X; one(zeros(eltype(p), k, k)) p' * X]) *
        [exp(-p'X); zeros(eltype(p), k, k)],
    )
end

@doc raw"""
    inverse_retract(M::Stiefel, p, q, ::PolarInverseRetraction)

Compute the inverse retraction based on a singular value decomposition
for two points `p`, `q` on the [`Stiefel`](@ref) manifold `M`.
This follows the folloing approach: From the Polar retraction we know that

````math
\operatorname{retr}_p^{-1}q = qs - t
````

if such a symmetric positive definite $k × k$ matrix exists. Since $qs - t$
is also a tangent vector at $p$ we obtain

````math
p^{\mathrm{H}}qs + s(p^{\mathrm{H}}q)^{\mathrm{H}} + 2I_k = 0,
````
which can either be solved by a Lyapunov approach or a continuous-time
algebraic Riccati equation as described in [^KanekoFioriTanaka2013]

This implementation follows the Lyapunov approach.

[^KanekoFioriTanaka2013]:
    > T. Kaneko, S. Fiori, T. Tanaka: "Empirical Arithmetic Averaging over the
    > Compact Stiefel Manifold", IEEE Transactions on Signal Processing, 2013,
    > doi: [10.1109/TSP.2012.2226167](https://doi.org/10.1109/TSP.2012.2226167).
"""
inverse_retract(::Stiefel, ::Any, ::Any, ::PolarInverseRetraction)

@doc raw"""
    inverse_retract(M::Stiefel, p, q, ::QRInverseRetraction)

Compute the inverse retraction based on a qr decomposition
for two points `p`, `q` on the [`Stiefel`](@ref) manifold `M` and return
the resulting tangent vector in `X`. The computation follows Algorithm 1
in [^KanekoFioriTanaka2013].

[^KanekoFioriTanaka2013]:
    > T. Kaneko, S. Fiori, T. Tanaka: "Empirical Arithmetic Averaging over the
    > Compact Stiefel Manifold", IEEE Transactions on Signal Processing, 2013,
    > doi: [10.1109/TSP.2012.2226167](https://doi.org/10.1109/TSP.2012.2226167).
"""
inverse_retract(::Stiefel, ::Any, ::Any, ::QRInverseRetraction)

function inverse_retract!(::Stiefel, X, p, q, ::PolarInverseRetraction)
    A = p' * q
    H = -2 * one(p' * p)
    B = lyap(A, H)
    return copyto!(X, q * B - p)
end
function inverse_retract!(::Stiefel{n,k}, X, p, q, ::QRInverseRetraction) where {n,k}
    A = p' * q
    R = zeros(typeof(one(eltype(p)) * one(eltype(q))), k, k)
    for i = 1:k
        b = zeros(i)
        b[i] = 1
        b[1:(end-1)] = -transpose(R[1:(i-1), 1:(i-1)]) * A[i, 1:(i-1)]
        R[1:i, i] = A[1:i, 1:i] \ b
    end
    return copyto!(X, q * R - p)
end

function isapprox(M::Stiefel, p, X, Y; kwargs...)
    return isapprox(sqrt(inner(M, p, zero_tangent_vector(M, p), X - Y)), 0; kwargs...)
end
isapprox(M::Stiefel, p, q; kwargs...) = isapprox(norm(p - q), 0; kwargs...)

@doc raw"""
    manifold_dimension(M::Stiefel)

Return the dimension of the [`Stiefel`](@ref) manifold `M`=$\operatorname{St}(n,k,𝔽)$.
The dimension is given by

````math
\begin{aligned}
\dim \mathrm{St}(n, k, ℝ) &= nk - \frac{1}{2}k(k+1)\\
\dim \mathrm{St}(n, k, ℂ) &= 2nk - k^2\\
\dim \mathrm{St}(n, k, ℍ) &= 4nk - k(2k-1)
\end{aligned}
````
"""
manifold_dimension(::Stiefel{n,k,ℝ}) where {n,k} = n * k - div(k * (k + 1), 2)
manifold_dimension(::Stiefel{n,k,ℂ}) where {n,k} = 2 * n * k - k * k
manifold_dimension(::Stiefel{n,k,ℍ}) where {n,k} = 4 * n * k - k * (2k - 1)

@doc doc"""
    project(M::Stiefel,p)

Projects `p` from the embedding onto the [`Stiefel`](@ref) `M`, i.e. compute `q`
as the polar decomposition of $p$ such that $q^{\mathrm{H}q$ is the identity,
where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed.
"""
project(::Stiefel, ::Any, ::Any)

function project!(M::Stiefel, q, p)
    s = svd(p)
    e = eigen(s.U' * s.U)
    qsinv = e.vectors * Diagonal(1 ./ sqrt.(e.values))
    q .= s.U * qsinv * e.vectors' * s.V'
    return q
end

@doc raw"""
    project(M::Stiefel, p, X)

Project `X` onto the tangent space of `p` to the [`Stiefel`](@ref) manifold `M`.
The formula reads

````math
\operatorname{proj}_{\mathcal M}(p, X) = X - p \operatorname{Sym}(p^{\mathrm{H}}X),
````

where $\operatorname{Sym}(q)$ is the symmetrization of $q$, e.g. by
$\operatorname{Sym}(q) = \frac{q^{\mathrm{H}}+q}{2}$.
"""
project(::Stiefel, ::Any...)

function project!(::Stiefel, Y, p, X)
    A = p' * X
    copyto!(Y, X - p * Hermitian((A + A') / 2))
    return Y
end

@doc raw"""
    retract(M::Stiefel, p, X, ::PolarRetraction)

Compute the SVD-based retraction [`PolarRetraction`](@ref) on the
[`Stiefel`](@ref) manifold `M`. With $USV = p + X$ the retraction reads

````math
\operatorname{retr}_p X = U\bar{V}^\mathrm{H}.
````
"""
retract(::Stiefel, ::Any, ::Any, ::PolarRetraction)

@doc raw"""
    retract(M::Stiefel, p, X, ::QRRetraction )

Compute the QR-based retraction [`QRRetraction`](@ref) on the
[`Stiefel`](@ref) manifold `M`. With $QR = p + X$ the retraction reads

````math
\operatorname{retr}_p X = QD,
````

where $D$ is a $n × k$ matrix with

````math
D = \operatorname{diag}\bigl(\operatorname{sgn}(R_{ii}+0,5)_{i=1}^k \bigr),
````

where $\operatorname{sgn}(p) = \begin{cases}
1 & \text{ for } p > 0,\\
0 & \text{ for } p = 0,\\
-1& \text{ for } p < 0.
\end{cases}$
"""
retract(::Stiefel, ::Any, ::Any, ::QRRetraction)

function retract!(::Stiefel, q, p, X, ::PolarRetraction)
    s = svd(p + X)
    return mul!(q, s.U, s.Vt)
end
function retract!(::Stiefel, q, p, X, ::QRRetraction)
    qrfac = qr(p + X)
    d = diag(qrfac.R)
    D = Diagonal(sign.(sign.(d .+ 0.5)))
    return copyto!(q, Matrix(qrfac.Q) * D)
end

@doc raw"""
    representation_size(M::Stiefel)

Returns the representation size of the [`Stiefel`](@ref) `M`=$\operatorname{St}(n,k)$,
i.e. `(n,k)`, which is the matrix dimensions.
"""
@generated representation_size(::Stiefel{n,k}) where {n,k} = (n, k)

show(io::IO, ::Stiefel{n,k,F}) where {n,k,F} = print(io, "Stiefel($(n), $(k), $(F))")