text:
module Test_sklfapps

using Random
using LinearAlgebra
using MatrixPencils
using Test

@testset "Structured Matrix Pencils Applications" begin

# test example for SLICOT subroutine AB08ND
A = [
   1.0   0.0   0.0   0.0   0.0   0.0
   0.0   1.0   0.0   0.0   0.0   0.0
   0.0   0.0   3.0   0.0   0.0   0.0
   0.0   0.0   0.0  -4.0   0.0   0.0
   0.0   0.0   0.0   0.0  -1.0   0.0
   0.0   0.0   0.0   0.0   0.0   3.0];
E = I;
B = [
   0.0  -1.0
  -1.0   0.0
   1.0  -1.0
   0.0   0.0
   0.0   1.0
  -1.0  -1.0];
C = [
   1.0   0.0   0.0   1.0   0.0   0.0
   0.0   1.0   0.0   1.0   0.0   1.0
   0.0   0.0   1.0   0.0   0.0   1.0];
D = [
   0.0   0.0
   0.0   0.0
   0.0   0.0];
M2 = [A B; C D]; N2 = [E zeros(6,2); zeros(3,8)];

# Kronecker structure
kinfo  = spkstruct(A, E, B, C, D, atol1 = 1.e-7, atol2 = 1.e-7)
@test (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [2], [2, 2], 2, 8)

kinfo  = spkstruct(A, E, view(B,:,1), C, view(D,:,1), atol1 = 1.e-7, atol2 = 1.e-7)
@test (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [1,2], [2], 2, 7)

# zeros
@time val, iz, kinfo  = spzeros(A, E, B, C, D)
@test sort(real(val)) ≈ [-1, 2, Inf,Inf] && iz == [1, 1] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [2], [2, 2], 2, 8)

@time val, iz, kinfo  = spzeros(A, E, view(B,:,1), C, view(D,:,1))
@test sort(real(val)) ≈ [-1, 2, Inf] && iz == [1] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [1, 2], [2], 2, 7)

@time val, kinfo  = speigvals(A, E, B, C, D)
@test sort(real(val)) ≈ [-1, 2, Inf, Inf, Inf, Inf] &&
(kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [2], [2, 2], 2, 8)

@time val, kinfo  = speigvals(A, E, view(B,:,1), C, view(D,:,1))
@test sort(real(val)) ≈ [-1, 2, Inf, Inf] &&
(kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [1, 2], [2], 2, 7)

@test sprank(A, E, B, C, D, fastrank = true) == 8 && sprank(A, E, B, C, D, fastrank = false) == 8

@test sprank(A, E, view(B,:,1), C, view(D,:,1), fastrank = true) == 7 && sprank(A, E, view(B,:,1), C, view(D,:,1), fastrank = false) == 7

# output decoupling zeros
kinfo  = spkstruct(A, E, missing, C, missing, atol1 = 1.e-7, atol2 = 1.e-7)
@test (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [1, 2, 2], [], 1, 6)

@time val, iz, kinfo  = spzeros(A, E, missing, C, missing)
@test val ≈ [-1] && iz == [] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [1, 2, 2], [], 1, 6)

@time val, kinfo  = speigvals(A, E, missing, C, missing, atol1 = 1.e-7, atol2 = 1.e-7)
@test val ≈ [-1] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [1, 2, 2], [], 1, 6)

@test sprank(A, E, missing, C, missing, fastrank = true) == 6 &&
      sprank(A, E, missing, C, missing, fastrank = false) == 6

# input decoupling zeros
kinfo  = spkstruct(A, E, B, missing, missing, atol1 = 1.e-7, atol2 = 1.e-7)
@test (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([2,3], [], [], 1, 6)

@time val, iz, kinfo  = spzeros(A, E, B, missing, missing)
@test val ≈ [-4] && iz == [] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([2,3], [], [], 1, 6)

@time val, kinfo  = speigvals(A, E, B, missing, missing)
@test val ≈ [-4] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([2,3], [], [], 1, 6)

@test sprank(A, E, B, missing, missing, fastrank = true) == 6 &&
      sprank(A, E, B, missing, missing, fastrank = false) == 6

@test sprank(A, E, missing, missing, missing, fastrank = true) == 6 &&
      sprank(A, E, missing, missing, missing, fastrank = false) == 6

@test sprank(missing, missing, missing, missing, D, fastrank = true) == 0 &&
      sprank(missing, missing, missing, missing, D, fastrank = false) == 0

# test example for SLICOT subroutine AB08ND
A = [
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     1     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     1     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0
    0     0     0     0     0     0     0     0     1];
E = [
    0     0     0     0     0     0     0     0     0
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     0     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     0     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0];
B = [
    -1     0     0
     0     0     0
     0     0     0
     0    -1     0
     0     0     0
     0     0     0
     0     0    -1
     0     0     0
     0     0     0];
C = [
    0     1     1     0     3     4     0     0     2
    0     1     0     0     4     0     0     2     0
    0     0     1     0    -1     4     0    -2     2];
D = [
    1     2    -2
    0    -1    -2
    0     0     0];
M2 = [A B; C D]
N2 = [E zeros(9,3); zeros(3,12)]

# poles
kinfo  = spkstruct(A, E, missing, missing, missing, atol1 = 1.e-7, atol2 = 1.e-7)
@test (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [], [3, 3, 3], 0, 9)

@time val, iz, kinfo  = spzeros(A, E, missing, missing, missing, atol1 = 1.e-7, atol2 = 1.e-7)
@test val ≈ [Inf, Inf, Inf, Inf, Inf, Inf] && iz == [2, 2, 2] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [], [3, 3, 3], 0, 9)

@time val, kinfo  = speigvals(A, E, missing, missing, missing, atol1 = 1.e-7, atol2 = 1.e-7)
@test val ≈ [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [], [3, 3, 3], 0, 9)

@test sprank(A, E, missing, missing, missing, fastrank = true) == 9 &&
      sprank(A, E, missing, missing, missing, fastrank = false) == 9

@test sprank(missing, missing, missing, missing, missing, fastrank = true) == 0 &&
      sprank(missing, missing, missing, missing, missing, fastrank = false) == 0

# output decoupling zeros

kinfo  = spkstruct(A, E, missing, C, missing, atol1 = 1.e-7, atol2 = 1.e-7)
@test (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [0, 1, 1], [1, 3, 3], 0, 9)

@time val, iz, kinfo  = spzeros(A, E, missing, C, missing)
@test val ≈ [Inf, Inf, Inf, Inf] && iz == [2, 2] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [0, 1, 1], [1, 3, 3], 0, 9)

@time val, kinfo  = speigvals(A, E, missing, C, missing, atol1 = 1.e-7, atol2 = 1.e-7)
@test val ≈ [Inf, Inf, Inf, Inf, Inf, Inf, Inf] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [0, 1, 1], [1, 3, 3], 0, 9)

@test sprank(A, E, missing, C, missing, fastrank = true) == 9 &&
      sprank(A, E, missing, C, missing, fastrank = false) == 9

# input decoupling zeros
kinfo  = spkstruct(A, E, B, missing, missing, atol1 = 1.e-7, atol2 = 1.e-7)
@test (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([2, 2, 2], Int64[], [1, 1, 1], 0, 9)

@time val, iz, kinfo  = spzeros(A, E, B, missing, missing)
@test val ≈ Float64[] && iz == [] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([2, 2, 2], Int64[], [1, 1, 1], 0, 9)

@time val, kinfo  = speigvals(A, E, B, missing, missing, atol1 = 1.e-7, atol2 = 1.e-7)
@test val ≈ [Inf, Inf, Inf] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([2, 2, 2], Int64[], [1, 1, 1], 0, 9)

@test sprank(A, E, B, missing, missing, fastrank = true) == 9 &&
      sprank(A, E, missing, C, missing, fastrank = false) == 9

# zeros
kinfo  = spkstruct(A, E, B, C, D, atol1 = 1.e-7, atol2 = 1.e-7)
@test (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([2], [1], [1, 1, 1, 1, 3], 1, 11)

@time val, iz, kinfo  = spzeros(A, E, B, C, D)
@test val ≈ [1, Inf,Inf] && iz == [2] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([2], [1], [1, 1, 1, 1, 3], 1, 11)

@time val, kinfo  = speigvals(A, E, B, C, D, atol1 = 1.e-7, atol2 = 1.e-7)
@test val ≈ [1, Inf, Inf, Inf, Inf, Inf, Inf, Inf] &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([2], [1], [1, 1, 1, 1, 3], 1, 11)

@test sprank(A, E, B, C, D, fastrank = true) == 11 &&
      sprank(A, E, B, C, D, fastrank = false) == 11

@test sprank(missing, missing, missing, missing, D, fastrank = true) == 2 &&
      sprank(missing, missing, missing, missing, D, fastrank = false) == 2

end  # applications testset

end  # module


text:
module cywwtaip

using Base.Threads
using Distributed

using JavaCall

####################################################################################################

JAR_FILE = joinpath(@__DIR__,"../","cywwtaip.jar")

####################################################################################################

function init()
  if JavaCall.isloaded() return end
  println("--------------------------------------------------------------------------------")
  @info "JavaCall Init..."
  JavaCall.init(["-Xmx512M", "-Djava.class.path=.;$(JAR_FILE)", "-verbose:jni", "-verbose:gc"])
  println("--------------------------------------------------------------------------------")
end

function destroy()
  println("--------------------------------------------------------------------------------")
  @info "Destroy JavaCall"
  #JavaCall.destroy()
  @warn "Disabled because of this problem: If destroyed, JavaCall cannot be initalized again in same process!"
  println("--------------------------------------------------------------------------------")
end

####################################################################################################

@info "JavaCall Imports..."
#jlm = @jimport "java.lang.Math"
JNetworkClient = @jimport lenz.htw.cywwtaip.net.NetworkClient
JGraphNode = @jimport lenz.htw.cywwtaip.world.GraphNode

####################################################################################################

@info "JavaCall Classes..."

mutable struct NetworkClient
  handle::JNetworkClient
  NetworkClient(handle::JNetworkClient) = new(handle)
end
export NetworkClient

const Position = Array{Float32, 1}
const Direction = Array{Float32, 1}

mutable struct Bot
  id::Int32
  player::Int32
  position::Position
  direction::Direction
  speed::Float32

  Bot(id::Integer, player::Integer, position::Position, direction::Direction, speed::Float32) = new(id, player, position, direction, speed)

  function Bot(client::NetworkClient, id::Integer, player::Integer)
    this = new(id, player)
    this.position = getBotPosition(client, player, id)
    this.direction = getBotDirection(client, id)
    this.speed = getBotSpeed(client, id)
    this
  end
end
export Bot

mutable struct GraphNode
  id::UInt32
  handle::Union{Nothing, JGraphNode}
  position::Array{jfloat, 1}
  parent::Union{Nothing, GraphNode}
  neighbors::Array{Union{Nothing, GraphNode, JGraphNode}, 1}
  blocked::Bool
  owner::Int32

  function GraphNode(id::Integer=0; handle::Union{Nothing, JGraphNode}=nothing)
    this = new(id, handle)
    this.parent = nothing
    if handle != nothing parse(this, handle) end
    this
  end
end
export GraphNode

# TEMPLATE ARRAY
TEMPLATE_LIST = Array{GraphNode,1}(undef, 40962)
@threads for i=1:length(TEMPLATE_LIST) TEMPLATE_LIST[i] = GraphNode(i) end

const DEFAULT_NODE_FUNCTION=(x)->x

getWorldTemplate() = Array{Union{Nothing,GraphNode},1}(undef,length(TEMPLATE_LIST))
export getWorldTemplate

function reset()
  #global TEMPLATE_LIST
  #@sync @distributed for i=1:length(TEMPLATE_LIST) TEMPLATE_LIST[i] = GraphNode() end
end

parse(this::JGraphNode) = GraphNode(0; handle=this)

function parse(this::GraphNode, handle::JGraphNode)
  this.handle = handle
  this.position = getPosition(handle) #[jfield(handle, "x", jfloat), jfield(handle, "y", jfloat), jfield(handle, "z", jfloat)]
  this.blocked = Bool(jfield(handle, "blocked", jboolean))
  this.owner = jfield(handle, "owner", jint)
  this.neighbors = jfield(handle, "neighbors", Array{JGraphNode, 1})
  this
end

function parse(handles::Array{JGraphNode, 1}; f::Function=DEFAULT_NODE_FUNCTION)
  nodes = TEMPLATE_LIST #deepcopy(TEMPLATE_LIST)
  @time @sync @distributed for i=1:length(handles) f(parse(nodes[i], handles[i])) end #node in nodes, handle in handles
  nodes
end

export parse

@info "JavaCall Hooks..."

NetworkClient(hostname::String, teamName::String, winMsg::String, timeout::Int=30) = NetworkClient(JNetworkClient((JString,JString,JString,jint), hostname, teamName, winMsg, timeout))

isAlive(this::NetworkClient) = Bool(jcall(this.handle, "isAlive", jboolean, ()))
getMyPlayerNumber(this::NetworkClient) = jcall(this.handle, "getMyPlayerNumber", jint, ())+1

getScore(this::NetworkClient, player::Integer) = jcall(this.handle, "getScore", jint, (jint,), Int32(player-1))
getScore(this) = [getScore(this,1),getScore(this,2),getScore(this,3)]

getBotPosition(this::NetworkClient, player::Integer, bot::Integer) = jcall(this.handle, "getBotPosition", Array{jfloat, 1}, (jint,jint,), Int32(player-1), Int32(bot-1)) #array
getBotDirection(this::NetworkClient, bot::Integer) = jcall(this.handle, "getBotDirection", Array{jfloat, 1}, (jint,), Int32(bot-1))
getBotSpeed(this::NetworkClient, bot::Integer) = jcall(this.handle, "getBotSpeed", jfloat, (jint,), Int32(bot-1))

getBotPositions(this::NetworkClient, player::Integer) = [getBotPosition(this, player, 1),getBotPosition(this, player, 2),getBotPosition(this, player, 3)]
getBotDirections(this::NetworkClient) = [getBotDirection(this, 1),getBotDirection(this, 2),getBotDirection(this, 3)]
getBotSpeeds(this::NetworkClient) = [getBotSpeed(this, 1),getBotSpeed(this, 2),getBotSpeed(this, 3)]

getBots(this::NetworkClient, player::Integer) = [Bot(this, 1, player),Bot(this, 2, player),Bot(this, 3, player)]
getBots(this::NetworkClient) = [getBots(this, 1),getBots(this, 2),getBots(this, 3)]

function updateBot(this::NetworkClient, bot::Bot)
  bot.position = getBotPosition(this, bot.player, bot.id)
  bot.direction = getBotDirection(this, bot.id)
  bot.speed = getBotSpeed(this, bot.id)
end

updateBots(this::NetworkClient, bots::Array{Bot,1}) = for bot in bots updateBot(this,bot) end
updateBots(this::NetworkClient, playerbots::Array{Array{Bot,1},1}) = for bots in playerbots updateBots(this,bots) end

changeMoveDirection(this::NetworkClient, bot::Integer, angle::AbstractFloat) = jcall(this.handle, "changeMoveDirection", Nothing, (jint,jfloat,), Int32(bot-1), Float32(angle))

function getGraph(this::NetworkClient)
  handles = jcall(this.handle, "getGraph", Array{JGraphNode, 1}, ())
  parse(handles)
end

function update(this::GraphNode, handle::JGraphNode)
  #this.handle = handle # not used, can be ignored
  #this.position # should not change (is static anyway)
  this.blocked = Bool(jfield(handle, "blocked", jboolean)) # can change
  this.owner = jfield(handle, "owner", jint) # can change
  #this.neighbors # should not change, we replaced it previously in createWorld()
  this
end

function updateGraph(this::NetworkClient, nodes::Array{GraphNode,1})
  handles = jcall(this.handle, "getGraph", Array{JGraphNode, 1}, ())
  @threads for i=1:length(handles) update(nodes[i], handles[i]) end
end

hashCode(this::NetworkClient) = jcall(this.handle, "hashCode", jint, ())
equals(this::NetworkClient, obj::JObject) = Bool(jcall(this.handle, "equals", jboolean, (JObject,), obj))
toString(this::NetworkClient) = jcall(this.handle, "toString", JString, ())

hashCode(this::GraphNode) = jcall(this.handle, "hashCode", jint, ())
equals(this::GraphNode, obj::JObject) = Bool(jcall(this.handle, "equals", jboolean, (JObject,), obj))
toString(this::GraphNode) = jcall(this.handle, "toString", JString, ())
toInfoString(this::GraphNode) = "Node["*string(this.id)*"] owner="* string(this.owner)*", blocked="*string(this.blocked)*", pos=["*string(join(this.position,", "))*"]"

getPosition(this::JGraphNode) = [jfield(this, "x", jfloat), jfield(this, "y", jfloat), jfield(this, "z", jfloat)]

####################################################################################################

export Position
export Direction

export hashCode
export equals
export toString
export toInfoString

export isAlive
export getMyPlayerNumber
export getScore
export getBotPosition
export getBotDirection
export getBotSpeed
export getBotPositions
export getBotDirections
export getBotSpeeds
export getBots
export changeMoveDirection
export getGraph
export updateGraph
export updateBot
export updateBots
export getPosition

####################################################################################################

Base.show(io::IO, this::GraphNode) = print(io, string("$(toString(this))"))
Base.show(io::IO, this::NetworkClient) = print(io, string("$(toString(this))"))

####################################################################################################

end # cywwtaip

ext: jl

max_stars_repo_name: Gilga/AI-Cywwtaip

max_issues_repo_name: Gilga/AI-Cywwtaip

max_forks_repo_name: Gilga/AI-Cywwtaip


text:
using OrbitalTrajectories

using Test
using DifferentialEquations
using BenchmarkTools

@testset "convert_to_frame" begin
    u0 = [0.8574053516112442, 0., 0., 0., 0.47, 0.]
    system_nbp = EphemerisNBP(:earth, :moon)
    prob = State(system_nbp, SynodicFrame(), u0, (0., 3600.0*24*30))

    # Test that converting to and from inertial frames works
    inert_prob = convert_to_frame(prob, InertialFrame())
    synod_prob = convert_to_frame(inert_prob, SynodicFrame())
    @test synod_prob.prob.u0 ≈ u0

    # Test that converting to/from non-normalised frames also works
    synod_nonorm_prob = convert_to_frame(prob, SynodicFrame(false))
    synod_norm_prob = convert_to_frame(synod_nonorm_prob, SynodicFrame(true))
    @test synod_norm_prob.prob.u0 ≈ prob.prob.u0
    inert_prob2 = convert_to_frame(synod_nonorm_prob, InertialFrame())
    synod_nonorm_prob2 = convert_to_frame(inert_prob2, SynodicFrame(false))
    @test synod_nonorm_prob2.prob.u0 ≈ synod_nonorm_prob.prob.u0

    # Roughly in decreasing order of allocations (30μs~40μs, allocations 8 to 38)
    # @btime convert_to_frame(prob, InertialFrame())  # norm -> inert
    # @btime convert_to_frame(inert_prob, SynodicFrame())  # inert -> norm
    # @btime convert_to_frame(prob, SynodicFrame(false))  # norm -> unorm
    # @btime convert_to_frame(synod_nonorm_prob, SynodicFrame(true))  # unorm -> norm
    # @btime convert_to_frame(synod_nonorm_prob, InertialFrame())  # unorm -> inert
    # @btime convert_to_frame(inert_prob2, SynodicFrame(false))  # inert -> unorm

    # Solve and test that converting the solution to synodic also matches
    sol = solve(prob)
    sol2 = convert_to_frame(sol, SynodicFrame())
    @test sol2.sol.u[begin] ≈ u0
end

@testset "EphemerisNBP performance with/without CSE" begin
    u0 = [0.8574053516112442, 0., 0., 0., 0.47, 0.]
    bodies = (:earth, :moon)

    # NOTE: This test only works for EphemerisNBP (other models remain equally performant)
    model = EphemerisNBP
    sys = model(bodies...)

    tspan = (0., 3600.0*24*30)
    prob = State(sys, SynodicFrame(), u0, tspan)

    # Test default performance of prob (with CSE)
    time_new = @benchmark solve($prob)

    # Test the old performance of prob (turn off Common Subexpression Elimination)
    sys_old = model(bodies...; wrap_code=(nothing, nothing))
    prob_old = State(sys_old, SynodicFrame(), u0, tspan)
    time_old = @benchmark solve($prob_old)

    @test 10 * median(time_new.times) < median(time_old.times)  # At least a 10x improvement
end


text:
# VectorOfDecisions #
# ========================== #
struct VectorDecisionConstraintBridge{T, S <: MOI.AbstractVectorSet} <: MOIB.Constraint.AbstractBridge
    constraint::CI{MOI.VectorOfVariables, S}
    decisions::VectorOfDecisions
end

function MOIB.Constraint.bridge_constraint(::Type{VectorDecisionConstraintBridge{T,S}},
                                           model,
                                           f::VectorOfDecisions,
                                           set::S) where {T, S <: MOI.AbstractVectorSet}
    # Perform the bridge mapping manually
    mapped_variables = map(f.decisions) do decision
        bridged = MOIB.bridged_function(
            model,
            AffineDecisionFunction{T}(SingleDecision(decision)))
        bridged.decision_part.terms[1].variable_index
    end
    # Add the bridged constraint
    constraint = MOI.add_constraint(model,
                                    MOI.VectorOfVariables(mapped_variables),
                                    set)
    # Save the constraint indices and the decision to allow modifications
    return VectorDecisionConstraintBridge{T,S}(constraint, f)
end

function MOI.supports_constraint(::Type{<:VectorDecisionConstraintBridge{T}},
                                 ::Type{VectorOfDecisions},
                                 ::Type{<:MOI.AbstractVectorSet}) where T
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:VectorDecisionConstraintBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:VectorDecisionConstraintBridge{T,S}}) where {T,S}
    return [(MOI.VectorOfVariables, S), (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T})]
end
function MOIB.Constraint.concrete_bridge_type(::Type{<:VectorDecisionConstraintBridge{T}},
                              ::Type{VectorOfDecisions},
                              S::Type{<:MOI.AbstractVectorSet}) where T
    return VectorDecisionConstraintBridge{T,S}
end

MOI.get(b::VectorDecisionConstraintBridge{T,S}, ::MOI.NumberOfConstraints{MOI.VectorOfVariables, S}) where {T,S} = 1
MOI.get(b::VectorDecisionConstraintBridge{T}, ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S}) where {T,S} = [b.constraint]

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::VectorDecisionConstraintBridge)
    return bridge.decisions
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 bridge::VectorDecisionConstraintBridge)
    return MOI.get(model, MOI.ConstraintSet(), bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal,
                 bridge::VectorDecisionConstraintBridge)
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::VectorDecisionConstraintBridge)
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, ::DecisionIndex,
                 bridge::VectorDecisionConstraintBridge{T}) where T
    return bridge.constraint
end

function MOI.delete(model::MOI.ModelLike, bridge::VectorDecisionConstraintBridge)
    MOI.delete(model, bridge.constraint)
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::VectorDecisionConstraintBridge{T,S}, change::S) where {T,S}
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint, change)
    return nothing
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Luna_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Luna")
JLLWrappers.@generate_main_file("Luna", UUID("5c9374ce-fdac-595e-b98d-9d97a97a2302"))
end  # module Luna_jll


text:
module TestSerialization

using Test
using MLJBase
using Serialization
using MLJTuning
using ..Models

function test_args(mach)
    # Check source nodes are empty if any
    for arg in mach.args
        if arg isa Source
            @test arg == source()
        end
    end
end

function test_data(mach)
    @test !isdefined(mach, :old_rows)
    @test !isdefined(mach, :data)
    @test !isdefined(mach, :resampled_data)
    @test !isdefined(mach, :cache)
end

function generic_tests(mach₁, mach₂)
    test_args(mach₂)
    test_data(mach₂)
    @test mach₂.state == -1
    for field in (:frozen, :model, :old_model, :old_upstream_state, :fit_okay)
        @test getfield(mach₁, field) == getfield(mach₂, field)
    end
end

@testset "Test TunedModel" begin
    filename = "tuned_model.jls"
    X, y = make_regression(100)
    base_model = DecisionTreeRegressor()
    tuned_model = TunedModel(
        model=base_model,
        tuning=Grid(),
        range=[range(base_model, :min_samples_split, values=[2,3,4])],
    )
    mach = machine(tuned_model, X, y)
    fit!(mach, rows=1:50, verbosity=0)
    smach = MLJBase.serializable(mach)
    @test smach.fitresult isa Machine
    @test smach.report == mach.report
    generic_tests(mach, smach)

    Serialization.serialize(filename, smach)
    smach = Serialization.deserialize(filename)
    MLJBase.restore!(smach)

    @test MLJBase.predict(smach, X) == MLJBase.predict(mach, X)
    @test fitted_params(smach) isa NamedTuple
    @test report(smach) == report(mach)

    rm(filename)

    # End to end
    MLJBase.save(filename, mach)
    smach = machine(filename)
    @test predict(smach, X) == predict(mach, X)

    rm(filename)

end

end

true

max_stars_repo_name: alan-turing-institute/MLJTuning

max_issues_repo_name: alan-turing-institute/MLJTuning

max_forks_repo_name: alan-turing-institute/MLJTuning


text:
abstract type AbstractExecution end
# Execution API #
# ------------------------------------------------------------
initialize_subproblems!(ph::AbstractProgressiveHedgingSolver, subsolver::QPSolver, penaltyterm::PenaltyTerm) = initialize_subproblems!(ph, subsolver, penaltyterm, ph.execution)
iterate!(ph::AbstractProgressiveHedgingSolver) = iterate!(ph, ph.execution)
start_workers!(ph::AbstractProgressiveHedgingSolver) = start_workers!(ph, ph.execution)
close_workers!(ph::AbstractProgressiveHedgingSolver) = close_workers!(ph, ph.execution)
resolve_subproblems!(ph::AbstractProgressiveHedgingSolver) = resolve_subproblems!(ph, ph.execution)
update_iterate!(ph::AbstractProgressiveHedgingSolver) = update_iterate!(ph, ph.execution)
update_subproblems!(ph::AbstractProgressiveHedgingSolver) = update_subproblems!(ph, ph.execution)
update_dual_gap!(ph::AbstractProgressiveHedgingSolver) = update_dual_gap!(ph, ph.execution)
calculate_objective_value(ph::AbstractProgressiveHedgingSolver) = calculate_objective_value(ph, ph.execution)
fill_first_stage!(ph::AbstractProgressiveHedgingSolver, stochasticprogram::StochasticProgram, nrows::Integer, ncols::Integer) = fill_first_stage!(ph, stochasticprogram, nrows, ncols, ph.execution)
fill_submodels!(ph::AbstractProgressiveHedgingSolver, scenarioproblems::AbstractScenarioProblems, nrows::Integer, ncols::Integer) = fill_submodels!(ph, scenarioproblems, nrows, ncols, ph.execution)
# ------------------------------------------------------------
include("common.jl")
include("serial.jl")
include("distributed.jl")
include("channels.jl")
include("worker.jl")
include("synchronous.jl")
include("asynchronous.jl")


text:
#################################################################################################
#                                                                                #
#                                           L2 Loss                                             #
#                                                                                #
#################################################################################################

##################################### B-Matrix methods ##########################################

"""
    function inner_loop_BTB!(qd::QuantizerData, diag::AbstractArray, Bsum, i)

Inner loop of the BTB computation, finds the values in the diagonal by taking out the non-zero
values from ∑B.
"""
function inner_loop_BTB!(qd::QuantizerData, diag::AbstractArray, Bsum, i)
    i1_nzval = (i-1)*qd.n_dims_center*qd.n_centers+1
    i2_nzval = i1_nzval+qd.n_centers-1
    i1_diag = (i-1)*qd.n_centers+1
    i2_diag = i1_diag + qd.n_centers-1
    @inbounds diag[i1_diag:i2_diag] = @inbounds @view Bsum.nzval[i1_nzval:i2_nzval]
end

"""
    function compute_BTB(qd::QuantizerData, processing::MultiThreaded)

Helper function of the codebook update step. Computes the ∑BTB matrix by summing the sparse matrices
and deriving the values on the diagonal from the non-zero values of the result to speed up the computaiton.
"""
function compute_BTB(qd::QuantizerData, processing::SingleThreaded)
    Bsum = sum(qd.I.B)
    diag = zeros(Int, qd.n_codebooks*qd.n_centers)
    for i in 1:qd.n_codebooks
        inner_loop_BTB!(qd, diag, Bsum, i)
    end
    return Diagonal(repeat(diag, inner=qd.n_dims_center))
end

function compute_BTB(qd::QuantizerData, processing::MultiThreaded)
    Bsum = sum(qd.I.B)
    diag = zeros(Int, qd.n_codebooks*qd.n_centers)
    Threads.@threads for i in 1:qd.n_codebooks
        inner_loop_BTB!(qd, diag, Bsum, i)
    end
    return Diagonal(repeat(diag, inner=qd.n_dims_center))
end

"""
    function compute_BTx(data::AbstractMatrix, qd::QuantizerData)

Computes the ∑BTx matrix from the quantizer and the input data.
"""
function compute_BTx(data::AbstractMatrix, qd::QuantizerData)
    Bsize = (qd.n_dims_center * qd.n_codebooks * qd.n_centers)
    sumBTx = zeros(Bsize)
    for i in 1:qd.n_dp
        sumBTx+= qd.I.B[i]*@view(data[:,i])
    end
    return sumBTx
end

function compute_BMatrices(data::AbstractMatrix, qd::QuantizerData, processing::Processing)
    BTB = compute_BTB(qd, processing)
    BTx = compute_BTx(data, qd)
    return BTB, BTx
end

"""
    function compute_BMatrices(data::AbstractMatrix, qd::QuantizerData, processing::GPU)

Main helper function of the codebook update step. Calls the right functions to compute the BTB and
BTx matrix from the data and configurations. Returns them as a normal array or `CuArray` depending
on the processing setting.
"""
function compute_BMatrices(data::AbstractMatrix, qd::QuantizerData, processing::GPU)
    BTB = compute_BTB(qd, MultiThreaded())
    BTx = compute_BTx(data, qd)
    BTB = CuArray(BTB)
    BTx = CuArray(BTx)
    return BTB, BTx
end

"""
function optimisation!(BTB::AbstractMatrix, BTx::AbstractArray, update_method::UpdateMethod,
    optimisation_method::OptimizationMethod, processing::Processing)

Computes the new codebook given the BTB and BTx marix and the optimisation method. Allows
for Exact updating (BTB^-1 ⋅ BTx) and Approximate (Nesterov) updating using ∇(BTB*C ⋅ C) - (2*BTx ⋅ C).
"""
function optimisation!(BTB::AbstractMatrix, BTx::AbstractArray,
                        optimisation_method::Exact, _::Any)
    return inv(BTB)BTx
end

function optimisation!(BTB::AbstractMatrix, BTx::AbstractArray,
                            optimisation_method::Nesterov, qd::QuantizerData)
    α = 1/qd.n_dp
    C = if typeof(BTB) <: CuArray CuArray(qd.C) else qd.C end
    dist=1
    i = 0
    μ = zero(C)
    C_old = zero(C)
    while dist > 1e-3
        C_old .= C
        β = 1-(3/(i+5))
        μ.*=β
        temp_C = C .- μ
        grad = gradient(x->quantization_objective(BTB, BTx, x), temp_C)[1]
        μ .+= α.*grad
        C .-= μ
        C_old .-= C
        dist = norm(C_old)
        i+=1
    end
    return C
end
"""
    function update_codebook!(data::AbstractMatrix, qd::QuantizerData, update_method::UpdateMethod,
    processing::Processing, optimisation_method::OptimizationMethod)

Main codebook update function; selects the appropriate methods to compute the ∑BTB and ∑BTx matrices
given the configurations and updates the codebook `qd.C` in place using the optimisation method selected.

"""
function update_codebook!(data::AbstractMatrix, qd::QuantizerData, update_method::BMatrix,
                                processing::Processing, optimisation_method::OptimizationMethod)
    BTB, BTx = compute_BMatrices(data, qd, processing)
    qd.C[:] .= optimisation!(BTB, BTx, optimisation_method, qd)
end

######################################## IVF methods #############################################

function update_codebook!(data::AbstractMatrix, qd::QuantizerData, update_method::InvertedIndex,
                            processing::SingleThreaded, optimisation::Exact)
    for j in 1:qd.n_codebooks
        for k in 1:qd.n_centers
            C_i1 = ((j-1) * qd.n_dims_center * qd.n_centers) + ((k-1) * qd.n_dims_center) + 1
            C_i2 = ((j-1) * qd.n_dims_center * qd.n_centers) + ((k) * qd.n_dims_center)
            idxs = qd.I.IVF[j][k]
            if length(idxs) > 0
                qd.C[C_i1:C_i2] = dropdims(mean(@inbounds(@view(data[((j-1) * qd.n_dims_center + 1):j*qd.n_dims_center, idxs])),dims=2), dims=2)
            end
        end
    end
end

function update_codebook!(data::AbstractMatrix, qd::QuantizerData, update_method::InvertedIndex,
                            processing::MultiThreaded, optimisation::Exact)
    Threads.@threads for j in 1:qd.n_codebooks
        for k in 1:qd.n_centers
            C_i1 = ((j-1) * qd.n_dims_center * qd.n_centers) + ((k-1) * qd.n_dims_center) + 1
            C_i2 = ((j-1) * qd.n_dims_center * qd.n_centers) + ((k) * qd.n_dims_center)
            idxs = qd.I.IVF[j][k]
            filter!(x->x>1,idxs)        # Multi-threading instability seems to have a very small probability of
            filter!(x->x<qd.n_dp,idxs)  # allocating random datapoint indexes
            if length(idxs) > 0
                qd.C[C_i1:C_i2] = dropdims(mean(@inbounds(@view(data[((j-1) * qd.n_dims_center + 1):j*qd.n_dims_center, idxs])),dims=2), dims=2)
            end
        end
    end
end

#################################################################################################
#                                                                                #
#                                        Anisotropic loss                                       #
#                                                                                #
#################################################################################################

compute_inner_matrix(dp::AbstractArray, η::AnisotropicWeights) = (η.h_par - η.h_orthog)*(dp*dp')+η.h_orthog*I(length(dp))

function compute_BTx(data::AbstractMatrix, qd::QuantizerData, η::AnisotropicWeights)
    sumBTx = zeros(qd.n_dims_center * qd.n_codebooks * qd.n_centers)
    for i in 1:qd.n_dp
        sumBTx+= (η.h_par*qd.I.B[i])*@view(data[:,i])
    end
    return sumBTx
end

function compute_BTB(data::AbstractMatrix, qd::QuantizerData, η::AnisotropicWeights, threading::SingleThreaded)
    Bsize = (qd.n_dims_center * qd.n_codebooks * qd.n_centers)
    sumBTB = zeros(Float32, Bsize, Bsize)
    for i in 1:qd.n_dp
        inner = compute_inner_matrix(@view(data[:,i]), η)
        left_mul = qd.I.B[i]*inner
        a,b,c = findnz(sparse(left_mul))
        a = unique(a)
        k = 0
        for i in a
            for j in a
                k+=1
                sumBTB[i,j] += c[k]
            end
        end
    end
    return sum(sumBTB,dims=3)
end

function compute_BTB(data::AbstractMatrix, qd::QuantizerData, η::AnisotropicWeights, threading::MultiThreaded)
    Bsize = (qd.n_dims_center * qd.n_codebooks * qd.n_centers)
    sumBTB = zeros(Float32, Bsize, Bsize, Threads.nthreads())
    Threads.@threads for i in 1:qd.n_dp
        inner = compute_inner_matrix(@view(data[:,i]), η)
        left_mul = qd.I.B[i]*inner
        a,b,c = findnz(sparse(left_mul))
        a = unique(a)
        k = 0
        for i in a
            for j in a
                k+=1
                sumBTB[i,j,Threads.threadid()] += c[k]
            end
        end
    end
    return dropdims(sum(sumBTB,dims=3), dims=3)
end

function compute_BMatrices(data::AbstractMatrix, qd::QuantizerData, η::AnisotropicWeights, processing::Processing)
    BTB = compute_BTB(data, qd, η, processing)
    BTx = compute_BTx(data, qd, η)
    return BTB, BTx
end

function compute_BMatrices(data::AbstractMatrix, qd::QuantizerData, η::AnisotropicWeights, processing::GPU)
    BTB = compute_BTB(data, qd, η, MultiThreaded())
    BTx = compute_BTx(data, qd, η)
    BTB = CuArray(BTB)
    BTx = CuArray(BTx)
    return BTB, BTx
end

function update_codebook!(data::AbstractMatrix, qd::QuantizerData, η::AnisotropicWeights,
                                processing::Processing, optimisation_method::OptimizationMethod)
    BTB, BTx = compute_BMatrices(data, qd, η, processing)
    try
        qd.C[:] .= optimisation!(BTB, BTx, optimisation_method, qd)
    catch
        @warn("""Singular Matrix encountered, some clusters do not contain assignments.
        This happens when the number of clusters is relatively high ccompared to the number of data points.
        \nOptimisation method changed to Approximate""")
        qd.C[:] .= optimisation!(BTB, BTx, Nesterov(), qd)
    end
end


text:
abstract type BoxMap end
"""
Transforms a `map` defined on ℝᴺ to a `BoxMap` defined on BoxSets

`map`:              map that defines the dynamical system

`domain_points`:    the spread of test points to be mapped forward in intersection algorithms.
                    (scaled to fit a box with unit radii)

`image_points`:     the spread of test points for comparison in intersection algorithms.
                    (scaled to fit a box with unit radii)

""" # TODO: see if domain is redundant in struct
struct SampledBoxMap{F,N,T,P,I} <: BoxMap
    map::F
    domain::Box{N,T}
    domain_points::P
    image_points::I
end

function Base.show(io::IO, g::SampledBoxMap)
    center, radius = g.domain.center, g.domain.radius
    n = length(g.domain_points(center, radius))
    print(io, "BoxMap with $(n) sample points")
end

function PointDiscretizedMap(map, domain, points::AbstractArray)
    domain_points = (center, radius) -> points
    image_points = (center, radius) -> center
    return SampledBoxMap(map, domain, domain_points, image_points)
end

function BoxMap(map, domain::Box{N,T}; no_of_points::Int=20*N) where {N,T}
    points = [ tuple(2.0*rand(N).-1.0 ...) for _ = 1:no_of_points ]
    return PointDiscretizedMap(map, domain, points)
end

function BoxMap(map, P::BoxPartition{N,T}; adaptive=false, no_of_points::Int=20*N) where {N,T}
    adaptive  ?  AdaptiveBoxMap(map, P.domain) : BoxMap(map, P.domain, no_of_points=no_of_points)
end

function sample_adaptive(Df, center::SVector{N,T}) where {N,T}  # how does this work?
    D = Df(center)
    _, σ, Vt = svd(D)
    n = ceil.(Int, σ)
    h = 2.0./(n.-1)
    points = Array{SVector{N,T}}(undef, ntuple(i->n[i], N))
    for i in CartesianIndices(points)
        points[i] = ntuple(k -> n[k]==1 ? 0.0 : (i[k]-1)*h[k]-1.0, N)
        points[i] = Vt'*points[i]
    end
    @debug points
    return points
end

function AdaptiveBoxMap(f, domain::Box{N,T}) where {N,T}
    Df = x -> ForwardDiff.jacobian(f, x)
    domain_points = (center, radius) -> sample_adaptive(Df, center)

    vertices = Array{SVector{N,T}}(undef, ntuple(k->2, N))
    for i in CartesianIndices(vertices)
        vertices[i] = ntuple(k -> (-1.0)^i[k], N)
    end
    # calculates the vertices of each box
    image_points = (center, radius) -> vertices
    return SampledBoxMap(f, domain, domain_points, image_points)
end

function map_boxes(g::BoxMap, source::BoxSet)
    P, keys = source.partition, collect(source.set)
    image = [ Set{eltype(keys)}() for k = 1:nthreads() ]
    @threads for key in keys
        box = key_to_box(P, key)
        c, r = box.center, box.radius
        points = g.domain_points(c, r)
        for p in points
            fp = g.map(c.+r.*p)
            hit = point_to_key(P, fp)
            if hit !== nothing
                push!(image[threadid()], hit)
            end
        end
    end
    BoxSet(P, union(image...))
end

(g::BoxMap)(source::BoxSet) = map_boxes(g, source)