text:
module NonstiffWrapper

using JavaCall

const AbstractFieldIntegrator = @jimport org.hipparchus.ode.AbstractFieldIntegrator
const AbstractFieldODEStateInterpolator = @jimport org.hipparchus.ode.sampling.AbstractFieldODEStateInterpolator
const AbstractIntegrator = @jimport org.hipparchus.ode.AbstractIntegrator
const AbstractODEStateInterpolator = @jimport org.hipparchus.ode.sampling.AbstractODEStateInterpolator
const AdamsBashforthFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.AdamsBashforthFieldIntegrator
const AdamsBashforthIntegrator = @jimport org.hipparchus.ode.nonstiff.AdamsBashforthIntegrator
const AdamsFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.AdamsFieldIntegrator
const AdamsFieldStateInterpolator = @jimport org.hipparchus.ode.nonstiff.AdamsFieldStateInterpolator
const AdamsIntegrator = @jimport org.hipparchus.ode.nonstiff.AdamsIntegrator
const AdamsMoultonFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.AdamsMoultonFieldIntegrator
const AdamsMoultonIntegrator = @jimport org.hipparchus.ode.nonstiff.AdamsMoultonIntegrator
const AdamsNordsieckFieldTransformer = @jimport org.hipparchus.ode.nonstiff.AdamsNordsieckFieldTransformer
const AdamsNordsieckTransformer = @jimport org.hipparchus.ode.nonstiff.AdamsNordsieckTransformer
const AdamsStateInterpolator = @jimport org.hipparchus.ode.nonstiff.AdamsStateInterpolator
const AdaptiveStepsizeFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator
const AdaptiveStepsizeIntegrator = @jimport org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator
const Array2DRowFieldMatrix = @jimport org.hipparchus.linear.Array2DRowFieldMatrix
const Array2DRowRealMatrix = @jimport org.hipparchus.linear.Array2DRowRealMatrix
const BracketedRealFieldUnivariateSolver = @jimport org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver
const BracketedUnivariateSolver = @jimport org.hipparchus.analysis.solvers.BracketedUnivariateSolver
const ButcherArrayProvider = @jimport org.hipparchus.ode.nonstiff.ButcherArrayProvider
const Class = @jimport java.lang.Class
const ClassicalRungeKuttaFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.ClassicalRungeKuttaFieldIntegrator
const ClassicalRungeKuttaIntegrator = @jimport org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator
const Collection = @jimport java.util.Collection
const DormandPrince54FieldIntegrator = @jimport org.hipparchus.ode.nonstiff.DormandPrince54FieldIntegrator
const DormandPrince54Integrator = @jimport org.hipparchus.ode.nonstiff.DormandPrince54Integrator
const DormandPrince853FieldIntegrator = @jimport org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator
const DormandPrince853Integrator = @jimport org.hipparchus.ode.nonstiff.DormandPrince853Integrator
const EmbeddedRungeKuttaFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaFieldIntegrator
const EmbeddedRungeKuttaIntegrator = @jimport org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator
const EquationsMapper = @jimport org.hipparchus.ode.EquationsMapper
const EulerFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.EulerFieldIntegrator
const EulerIntegrator = @jimport org.hipparchus.ode.nonstiff.EulerIntegrator
const ExpandableODE = @jimport org.hipparchus.ode.ExpandableODE
const Field = @jimport org.hipparchus.Field
const FieldButcherArrayProvider = @jimport org.hipparchus.ode.nonstiff.FieldButcherArrayProvider
const FieldEquationsMapper = @jimport org.hipparchus.ode.FieldEquationsMapper
const FieldExpandableODE = @jimport org.hipparchus.ode.FieldExpandableODE
const FieldODEEventHandler = @jimport org.hipparchus.ode.events.FieldODEEventHandler
const FieldODEIntegrator = @jimport org.hipparchus.ode.FieldODEIntegrator
const FieldODEState = @jimport org.hipparchus.ode.FieldODEState
const FieldODEStateAndDerivative = @jimport org.hipparchus.ode.FieldODEStateAndDerivative
const FieldODEStepHandler = @jimport org.hipparchus.ode.sampling.FieldODEStepHandler
const FieldOrdinaryDifferentialEquation = @jimport org.hipparchus.ode.FieldOrdinaryDifferentialEquation
const GillFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.GillFieldIntegrator
const GillIntegrator = @jimport org.hipparchus.ode.nonstiff.GillIntegrator
const GraggBulirschStoerIntegrator = @jimport org.hipparchus.ode.nonstiff.GraggBulirschStoerIntegrator
const GraggBulirschStoerStateInterpolator = @jimport org.hipparchus.ode.nonstiff.GraggBulirschStoerStateInterpolator
const HighamHall54FieldIntegrator = @jimport org.hipparchus.ode.nonstiff.HighamHall54FieldIntegrator
const HighamHall54Integrator = @jimport org.hipparchus.ode.nonstiff.HighamHall54Integrator
const JString = @jimport java.lang.String
const LutherFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.LutherFieldIntegrator
const LutherIntegrator = @jimport org.hipparchus.ode.nonstiff.LutherIntegrator
const MidpointFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.MidpointFieldIntegrator
const MidpointIntegrator = @jimport org.hipparchus.ode.nonstiff.MidpointIntegrator
const MultistepFieldIntegrator = @jimport org.hipparchus.ode.MultistepFieldIntegrator
const MultistepIntegrator = @jimport org.hipparchus.ode.MultistepIntegrator
const ODEEventHandler = @jimport org.hipparchus.ode.events.ODEEventHandler
const ODEIntegrator = @jimport org.hipparchus.ode.ODEIntegrator
const ODEState = @jimport org.hipparchus.ode.ODEState
const ODEStateAndDerivative = @jimport org.hipparchus.ode.ODEStateAndDerivative
const ODEStepHandler = @jimport org.hipparchus.ode.sampling.ODEStepHandler
const Object = @jimport java.lang.Object
const OrdinaryDifferentialEquation = @jimport org.hipparchus.ode.OrdinaryDifferentialEquation
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const RungeKuttaFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.RungeKuttaFieldIntegrator
const RungeKuttaIntegrator = @jimport org.hipparchus.ode.nonstiff.RungeKuttaIntegrator
const ThreeEighthesFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.ThreeEighthesFieldIntegrator
const ThreeEighthesIntegrator = @jimport org.hipparchus.ode.nonstiff.ThreeEighthesIntegrator

export AdamsBashforthFieldIntegrator
export AdamsBashforthIntegrator
export AdamsFieldIntegrator
export AdamsFieldStateInterpolator
export AdamsIntegrator
export AdamsMoultonFieldIntegrator
export AdamsMoultonIntegrator
export AdamsNordsieckFieldTransformer
export AdamsNordsieckTransformer
export AdamsStateInterpolator
export AdaptiveStepsizeFieldIntegrator
export AdaptiveStepsizeIntegrator
export ButcherArrayProvider
export ClassicalRungeKuttaFieldIntegrator
export ClassicalRungeKuttaFieldStateInterpolator
export ClassicalRungeKuttaIntegrator
export ClassicalRungeKuttaStateInterpolator
export DormandPrince54FieldIntegrator
export DormandPrince54FieldStateInterpolator
export DormandPrince54Integrator
export DormandPrince54StateInterpolator
export DormandPrince853FieldIntegrator
export DormandPrince853FieldStateInterpolator
export DormandPrince853Integrator
export DormandPrince853StateInterpolator
export EmbeddedRungeKuttaFieldIntegrator
export EmbeddedRungeKuttaIntegrator
export EulerFieldIntegrator
export EulerFieldStateInterpolator
export EulerIntegrator
export EulerStateInterpolator
export FieldButcherArrayProvider
export GillFieldIntegrator
export GillFieldStateInterpolator
export GillIntegrator
export GillStateInterpolator
export GraggBulirschStoerIntegrator
export GraggBulirschStoerStateInterpolator
export HighamHall54FieldIntegrator
export HighamHall54FieldStateInterpolator
export HighamHall54Integrator
export HighamHall54StateInterpolator
export LutherFieldIntegrator
export LutherFieldStateInterpolator
export LutherIntegrator
export LutherStateInterpolator
export MidpointFieldIntegrator
export MidpointFieldStateInterpolator
export MidpointIntegrator
export MidpointStateInterpolator
export RungeKuttaFieldIntegrator
export RungeKuttaFieldStateInterpolator
export RungeKuttaIntegrator
export RungeKuttaStateInterpolator
export ThreeEighthesFieldIntegrator
export ThreeEighthesFieldStateInterpolator
export ThreeEighthesIntegrator
export ThreeEighthesStateInterpolator
export add_event_handler
export add_step_handler
export clear_event_handlers
export clear_step_handlers
export compute_derivatives
export equals
export estimate_error
export get_a
export get_b
export get_c
export get_class
export get_current_signed_stepsize
export get_current_state
export get_evaluations
export get_event_handlers
export get_field
export get_global_current_state
export get_global_previous_state
export get_instance
export get_interpolated_state
export get_max_evaluations
export get_max_growth
export get_max_step
export get_min_reduction
export get_min_step
export get_n_steps
export get_name
export get_order
export get_previous_state
export get_safety
export get_starter_integrator
export get_step_handlers
export get_step_start
export hash_code
export initialize_high_order_derivatives
export initialize_step
export integrate
export is_current_state_interpolated
export is_forward
export is_previous_state_interpolated
export notify
export notify_all
export restrict_step
export set_control_factors
export set_initial_step_size
export set_interpolation_control
export set_max_evaluations
export set_max_growth
export set_min_reduction
export set_order_control
export set_safety
export set_stability_check
export set_starter_integrator
export set_step_size_control
export single_step
export taylor
export to_string
export update_high_order_derivatives_phase1
export update_high_order_derivatives_phase2
export wait

include("adams_bashforth_field_integrator.jl")
include("adams_bashforth_integrator.jl")
include("adams_field_integrator.jl")
include("adams_field_state_interpolator.jl")
include("adams_integrator.jl")
include("adams_moulton_field_integrator.jl")
include("adams_moulton_integrator.jl")
include("adams_nordsieck_field_transformer.jl")
include("adams_nordsieck_transformer.jl")
include("adams_state_interpolator.jl")
include("adaptive_stepsize_field_integrator.jl")
include("adaptive_stepsize_integrator.jl")
include("butcher_array_provider.jl")
include("classical_runge_kutta_field_integrator.jl")
include("classical_runge_kutta_field_state_interpolator.jl")
include("classical_runge_kutta_integrator.jl")
include("classical_runge_kutta_state_interpolator.jl")
include("dormand_prince54_field_integrator.jl")
include("dormand_prince54_field_state_interpolator.jl")
include("dormand_prince54_integrator.jl")
include("dormand_prince54_state_interpolator.jl")
include("dormand_prince853_field_integrator.jl")
include("dormand_prince853_field_state_interpolator.jl")
include("dormand_prince853_integrator.jl")
include("dormand_prince853_state_interpolator.jl")
include("embedded_runge_kutta_field_integrator.jl")
include("embedded_runge_kutta_integrator.jl")
include("euler_field_integrator.jl")
include("euler_field_state_interpolator.jl")
include("euler_integrator.jl")
include("euler_state_interpolator.jl")
include("field_butcher_array_provider.jl")
include("gill_field_integrator.jl")
include("gill_field_state_interpolator.jl")
include("gill_integrator.jl")
include("gill_state_interpolator.jl")
include("gragg_bulirsch_stoer_integrator.jl")
include("gragg_bulirsch_stoer_state_interpolator.jl")
include("higham_hall54_field_integrator.jl")
include("higham_hall54_field_state_interpolator.jl")
include("higham_hall54_integrator.jl")
include("higham_hall54_state_interpolator.jl")
include("luther_field_integrator.jl")
include("luther_field_state_interpolator.jl")
include("luther_integrator.jl")
include("luther_state_interpolator.jl")
include("midpoint_field_integrator.jl")
include("midpoint_field_state_interpolator.jl")
include("midpoint_integrator.jl")
include("midpoint_state_interpolator.jl")
include("runge_kutta_field_integrator.jl")
include("runge_kutta_field_state_interpolator.jl")
include("runge_kutta_integrator.jl")
include("runge_kutta_state_interpolator.jl")
include("three_eighthes_field_integrator.jl")
include("three_eighthes_field_state_interpolator.jl")
include("three_eighthes_integrator.jl")
include("three_eighthes_state_interpolator.jl")

end


text:
function value_of(::Type{LibrationOrbitType}, arg0::JString)
    return jcall(LibrationOrbitType, "valueOf", LibrationOrbitType, (JString,), arg0)
end

function values(::Type{LibrationOrbitType})
    return jcall(LibrationOrbitType, "values", Vector{LibrationOrbitType}, ())
end

ext: jl


text:
function readInitPop( mpSim::MPsim, xf::XF, filePath::String )

    if !XLSX.hassheet( xf, "Snapshot" )
        @warn "Configuration file is missing a 'Snapshoy' sheet, cannot load the initial population of the organisation."
        return
    end  # if !XLSX.hassheet( xf, "Snapshot" )

    ws = xf["Snapshot"]

    if ws["B3"] == "NO"
        return
    end  # if ws["B3"] == "NO"

    # Initial population file name.
    initPopName = ws["B4"] isa Missing ? "" : string( ws["B4"] )
    initPopName = string( initPopName,
        endswith( initPopName, ".xlsx" ) ? "" : ".xlsx" )
    initPopName = normpath( joinpath( filePath, "..", initPopName ) )

    if !ispath( initPopName )
        @warn string( "Initial population file '", initPopName,
            "' does not exist. Cannot upload initial population." )
        return
    end  # if !ispath( initPopName )

    initPopSheet = ws["B5"] isa Missing ? "" : string( ws["B5"] )
    isSheetOkay = true

    XLSX.openxlsx( initPopName ) do sxf
        isSheetOkay = XLSX.hassheet( sxf, initPopSheet )
    end  # XLSX.openxlsx( initPopName ) do sxf

    if !isSheetOkay
        @warn string( "Initial population file does not have sheet '",
            initPopSheet, "'. Cannot upload initial population." )
        return
    end  # if !isSheetOkay

    # Read special columns.
    specialCols = [ws["B6"], ws["B11"], ws["B7"], ws["B9"], ws["B12"]]

    for ii in 1:5
        if !( ( specialCols[ii] isa String ) || ( specialCols[ii] isa Int ) )
            specialCols[ii] = 0
        end  # if !( ... )
    end  # for ii in 1:5

    # Read attribute columns.
    nAttributes = ws["B16"]
    attributeCols = []

    if nAttributes > 0
        attributeCols = ws[CR( 19, 1, 18 + nAttributes, 1 )][:]
    end  # if nAttributes > 0

    attributeCols = map( attr -> ( attr isa Int ) || ( attr isa String ) ?
        attr : 0, attributeCols )
    uploadSnapshot( mpSim, initPopName, initPopSheet, attributeCols,
        Tuple( specialCols ); generateSimData=ws["B13"] == "YES" )

end  # readInitPop( mpSim, xf, filePath )


text:
"""
    find_scored_idx(collection, what::String; check::Function=_->true, score::Function=_->0)

A helper function that runs through a `collection` and select an item that
satisfies the predicate `check` with best score computed by `score`. If no
suitable item is found, an error is thrown referring to the item type as
`what`.
"""
function find_scored_idx(
    collection,
    what::String;
    check::Function = _ -> true,
    score::Function = _ -> 0,
)
    found = 0
    found_score = -Inf
    for (idx, obj) in enumerate(collection)
        if check(obj)
            s = score(obj)
            if s > found_score
                found = idx
                found_score = s
            end
        end
    end
    found > 0 || error("No suitable $what found")
    return found
end

"""
    find_physical_device(instance; args...)

Find a physical device within the instance that satisfies some given
properties; arguments are forwarded to [`find_scored_idx`](@ref).

# Example

    physical_device = find_physical_device(instance, check=props->contains(props.device_name, "MyFavVendor"))
"""
function find_physical_device(instance; args...)
    devs = unwrap(enumerate_physical_devices(instance))
    return devs[find_scored_idx(
        get_physical_device_properties.(devs),
        "physical device";
        args...,
    )]
end

"""
    find_queue_family_idx(physical_device; args...)

Find a queue family index (numbered from 0) for a physical device that
satisfies some given properties; arguments are forwarded to
[`find_scored_idx`](@ref).

# Example

    # Find a queue family with most queues that supports compute and transfer operations.
    find_queue_family_idx(physical_device,
        check=qfp->hasbits(qfp.queue_flags, QUEUE_TRANSFER_BIT | QUEUE_COMPUTE_BIT),
        score=qfp->qfp.queue_count,
    )
"""
find_queue_family_idx(physical_device; args...) =
    find_scored_idx(
        get_physical_device_queue_family_properties(physical_device),
        "queue family";
        args...,
    ) - 1

"""
    find_memory_type_idx(physical_device; check::Function=(_,_)->true, score::Function=(_,_)->0)

Find a memory type index (numbered from 0) for a physical device that satisfies
some given properties. The arguments `check` and `score` behave as with
[`find_scored_idx`](@ref), but take 2 arguments with types `MemoryType` and
`MemoryHeap`.

# Example

    # Find a host visible and device local memory type with largest heap
    find_memory_type_idx(physical_device,
        check=(mt, heap) -> hasbits(mt.property_flags, MEMORY_PROPERTY_HOST_VISIBLE_BIT|MEMORY_PROPERTY_DEVICE_LOCAL_BIT),
        score=(mt, heap) -> heap.size,
    )

"""
function find_memory_type_idx(
    physical_device;
    check::Function = (_, _) -> true,
    score::Function = (_, _) -> 0,
)
    props = get_physical_device_memory_properties(physical_device)
    return find_scored_idx(
        props.memory_types[1:props.memory_type_count],
        "memory type";
        check = mt -> check(mt, props.memory_heaps[mt.heap_index+1]),
        score = mt -> score(mt, props.memory_heaps[mt.heap_index+1]),
    ) - 1
end


text:
# Now, we need an abstract type representing a numerical discretization scheme
# for the ODE
"""
Represents a map that given (t_n,u_n) returns (t_n+1,u_n+1) and cache for the
corresponding `ODEOperator` and `NonlinearOperator`
"""
abstract type ODESolver <: GridapType end

function solve_step!(
  uF::Union{AbstractVector,Tuple{Vararg{AbstractVector}}},
  solver::ODESolver,
  op::ODEOperator,
  u0::Union{AbstractVector,Tuple{Vararg{AbstractVector}}},
  t0::Real,
  cache) # -> (uF,tF,cache)
  @abstractmethod
end

# Default API

function solve_step!(
  uF::Union{AbstractVector,Tuple{Vararg{AbstractVector}}},
  solver::ODESolver,
  op::ODEOperator,
  u0::Union{AbstractVector,Tuple{Vararg{AbstractVector}}},
  t0::Real) # -> (uF,tF,cache)
  solve_step!(uF,solver,op,u0,t0,nothing)
end

function solve(
  solver::ODESolver,
  op::ODEOperator,
  u0::T,
  t0::Real,
  tf::Real) where {T}
  GenericODESolution{T}(solver,op,u0,t0,tf)
end

# testers

function test_ode_solver(solver::ODESolver,op::ODEOperator,u0,t0,tf)
  solution = solve(solver,op,u0,t0,tf)
  test_ode_solution(solution)
end

# Specialization

include("ForwardEuler.jl")

include("ThetaMethod.jl")

include("AffineThetaMethod.jl")

include("RungeKutta.jl")

include("Newmark.jl")

include("AffineNewmark.jl")

include("ConstantNewmark.jl")

include("ConstantMatrixNewmark.jl")


text:
export Problem, print

"""
    Problem()

Define an empty mathematical polynomial optimization problem, characterized by
a polynomial `objective`, a dictionnary of `constraints` and a dictionnary of
used variables `variables`.

### Attributes
- `objective::Polynomial` : polynomial criterion.
- `constraints::SortedDict{String, Constraint}` : dictionnary of constraint name to
constraint.
- `variables::SortedDict{String, Type}` : dictionnary of variable name to type.

### Exemple
```julia
julia > x, y, z = ...
julia > pb = Problem()
julia > set_objective(pb, abs2(x+y+z))
julia > add_constraint!(pb, abs2(x) << 1)
julia > add_constraint!(pb, y << 0.5+1im)
julia > add_constraint!(pb, 0-1im << z << 1+0im)
```
"""
mutable struct Problem
  objective::Polynomial
  constraints::SortedDict{String, Constraint}
  variables::SortedDict{String, Type}
end

#############################
## Constructors
#############################
Problem() = Problem(Polynomial(), SortedDict{String, Constraint}(), SortedDict{String, Type}())

#############################
## Equality, hash
#############################
==(pb1::Problem, pb2::Problem) = (pb1.objective == pb2.objective) &&
                                 (pb1.constraints == pb2.constraints) &&
                                 (pb1.variables == pb2.variables)
!=(pb1::Problem, pb2::Problem) = !(pb1==pb2)

#############################
## Print
#############################
function Base.print(io::IO, pb::Problem)
  print(io, "▶ variables: ")
  for (varName, typ) in pb.variables
    print(io, Variable(varName, typ), " ")
  end
  println(io, "\n▶ objective: ", pb.objective)
  println(io, "▶ constraints: ")
  for (cstrName, cstr) in pb.constraints
    @printf(io, " → %10s: ", cstrName)
    println(io, cstr)
  end
end


text:
module SymbolServer

conn = stdout
(outRead, outWrite) = redirect_stdout()

module LoadingBay
end
using Serialization, Pkg, SHA
using Base: UUID
@static if VERSION < v"1.1"
    const PackageEntry = Vector{Dict{String,Any}}
else
    using Pkg.Types: PackageEntry
end
include("symbols.jl")
include("utils.jl")

server = Server(abspath(joinpath(@__DIR__, "..", "store")), Pkg.Types.Context(), Dict{Any,Any}())
if Sys.isunix()
    global const nullfile = "/dev/null"
elseif Sys.iswindows()
    global const nullfile = "nul"
else
    error("Platform not supported")
end

function write_cache(uuid, pkg)
    open(joinpath(server.storedir, "$uuid.jstore"), "w") do io
        serialize(io, pkg)
    end
end

while true
    message, payload = deserialize(stdin)
    if message == :get_context
        serialize(conn, (:success, server.context))
    elseif message == :cache_package
        for uuid in payload
            cache_package(server.context, UUID(uuid), server.depot)
        end
        out = Tuple{String,String,Bool,Bool, Bool}[] # list of saved caches
        for  (uuid, pkg) in server.depot
            overwrote = isfile(joinpath(server.storedir, "$(string(uuid)).jstore"))
            write_cache(uuid, pkg)

            isloaded = can_access(LoadingBay, Symbol(packagename(server.context, uuid))) isa Module
            issaved = isfile(joinpath(server.storedir, "$(string(uuid)).jstore"))

            push!(out, (string(uuid), packagename(server.context, uuid), isloaded, isloaded, overwrote))
        end
        serialize(conn, (:success, out))
    elseif message == :change_env
        Pkg.API.activate(payload)
        server.context = Pkg.Types.Context()
        serialize(conn, (:success, nothing))
    elseif message == :debugmessage
        out = string(eval(Meta.parse(payload)))
        serialize(conn, (:success, out))
    elseif message == :close
        break
    else
        serialize(conn, (:failure, nothing))
    end
end
end


text:
"""
    Body{T} <: Node{T}

    A rigid body object

    id: unique identifying number
    name: unique identifying name
    mass: inertial property (kilograms)
    inertia: inertia matrix (kilograms meter^2)
    state: State; representation of the system's: position, linear velocity, orientation, and angular velocity
    shape: Shape; geometry information about the Body
"""
mutable struct Body{T} <: Node{T}
    id::Int64
    name::Symbol
    mass::T
    inertia::SMatrix{3,3,T,9}
    state::State{T}
    shape::Shape{T}

    function Body(mass::T, inertia::AbstractMatrix;
        name::Symbol=:origin,
        shape::Shape=EmptyShape()) where T
        new{T}(getGlobalID(), name, mass, inertia, State{T}(), shape)
    end
end

Base.length(::Body) = 6

# warn that body has poor inertial properties
function check_body(body::Body)
    if norm(body.mass) == 0 || norm(body.inertia) == 0
        @info "Bad inertial properties detected"
    end
    #TODO check condition number and potentially adapt
end


text:
###
# Backtest using Rossman data
#
# For increasing K,
#  Look at N consecutive days of data
#  Compute various methods
#  Record output by looking at performance on next numTestDays
###
using Distributions, DelimitedFiles
include("../src/JS_SAA_main.jl")

#supp_full, ps_full are d x K matrices with true info per problem
#binned_data is just the binned_data, may have "NA"
#dates is the associated dates
#adds new subproblems in order appear in files
#s is the service level, N is the average amount of data per problem
function rollingTest(K_grid, supp_full, ps_full, binned_data_full, dates, outPath, N_grid, s;
					onlySAA = false, numTestDays=10, alpha_max = 180, alpha_len = 120,
					allCV = false)
	Kmax = maximum(K_grid)
	@assert Kmax <= size(supp_full, 2) "K_grid exceeds available subproblems"
	@assert size(supp_full) == size(ps_full) "supp_full and ps_full have incompatible dimensions"
	d = size(supp_full, 1)
	numDataPoints = size(binned_data_full, 1)

	#For safety, trim inputs to size Kmax
	supp_full = view(supp_full, 1:d, 1:Kmax)
	ps_full = view(ps_full, 1:d, 1:Kmax)
	binned_data = view(binned_data_full, 1:numDataPoints, 1:Kmax)
	dates = view(dates, 1:numDataPoints)

	p0 = ones(d) / d
	alpha_grid = range(0., stop=alpha_max, length=alpha_len)
	Gamma_grid = range(0, stop = (1-s), length=51)
	theta2_grid = range(1e-6, stop=3, length=20)
	mu_grid = range(1e-6, stop=1, length=20)

	#set up output file
	f = open("$(outPath).csv", "w")
	writedlm(f, ["StartDate" "K" "d" "N" "Method" "TruePerf" "time" "alpha"], ',')

	#generate all Kmax subproblems upfront and store in memory
	cs_full = JS.getNewsVendorCosts(supp_full, s, Kmax)
	xs_full = JS.genSSAAtrainers(supp_full, s)
	xsKS_full = JS.genKSTrainers(supp_full, s, Kmax, :crossVal)

	lam_full = ones(Kmax)

	for N in N_grid
		for ix_start = 1:N:numDataPoints
			#skip if partial window for training or test
			if ix_start + N + numTestDays - 1 > numDataPoints
				continue
			end

			#compute Training Data
			#within the window
			#Throw away NA's
			mhats_full = zeros(Int, d, Kmax)
			for k = 1:Kmax
				for i = 1:N
					if binned_data[ix_start + i - 1, k] == "NA"
						continue
					end
					mhats_full[binned_data[ix_start + i - 1, k], k] += 1
				end
			end

			Nhats_full = vec(sum(mhats_full, dims=1))
			@assert length(Nhats_full) == Kmax

			#build the out of sample training set
			mhats_out_full = zeros(Int, d, Kmax)
			for k = 1:Kmax
				for i = 1:numTestDays
					if binned_data[ix_start + N + i - 1, k] == "NA"
						continue
					else
						mhats_out_full[binned_data[ix_start + N + i - 1, k], k] += 1
					end
				end
			end

			for K in K_grid
				println("($(N), $(ix_start), $(K))")
				#Take views on everything for simplicity
				lams = view(lam_full, 1:K)
				supp = view(supp_full, 1:d, 1:K)
				Nhats = view(Nhats_full, 1:K)
				mhats = view(mhats_full, 1:d, 1:K)
				cs = view(cs_full, 1:d, 1:K)
				xs = view(xs_full, 1:K)
				xsKS = view(xsKS_full, 1:K)

				mhats_out = view(mhats_out_full, 1:d, 1:K)

				N_out = sum(mhats_out, dims=1)
				ps = mhats_out ./ N_out
				ps[:, vec(N_out) .== 0] .= 0.

				#for data-driven shrinkage anchor
				#problems with no data do not contribute to the anchor
				phat_avg = JS.get_GM_anchor(mhats)
				if !isapprox(sum(phat_avg), 1) || minimum(phat_avg) < 0
					throw("Something went terriby wrong.  Grand avg phat is not a probability distribution")
				end

				#Compute the full-info value once for reference
				t =
				  @elapsed full_info = JS.zstar(xs, cs, ps, lams)
				writedlm(f, [dates[ix_start] K d N "FullInfo" full_info t 0.], ',')

				#SAA
				t =
				  @elapsed perf_SAA = JS.zbar(xs, cs, mhats, ps, lams, (p0, 0.))
				writedlm(f, [dates[ix_start] K d N "SAA" perf_SAA t 0.0], ',')

				if onlySAA
					continue
				end

				#KS Robust value
				t =
				  @elapsed perf_KS = JS.zbar(xsKS, cs, mhats, ps, lams, (Gamma_grid, 5))
				writedlm(f, [dates[ix_start] K d N "KS" perf_KS t 0.0], ',')

				#Gen the Oracle cost with 1/d anchor
				t =
				  @elapsed alphaOR, min_indx, or_alpha_curve = JS.oracle_alpha(xs, cs, mhats, ps, lams, p0, alpha_grid)
				writedlm(f, [dates[ix_start] K d N "Oracle" or_alpha_curve[min_indx] t alphaOR], ',')

				#Gen the LOO cost with 1/d anchor
				t =
				  @elapsed alphaLOO, min_indx, looUnsc_curve = JS.loo_alpha(xs, cs, mhats, p0, alpha_grid)
				writedlm(f, [dates[ix_start] K d N "LOO_unif" or_alpha_curve[min_indx] t alphaLOO], ',')

				#Gen the CV5 cost with 1/d anchor
				t =
				  @elapsed alphaCV, min_indx, CVUnsc_curve = JS.cv_alpha(xs, cs, mhats, p0, alpha_grid, 5)
				writedlm(f, [dates[ix_start] K d N "CV5_unif" or_alpha_curve[min_indx] t alphaCV], ',')

				if allCV
					#Gen the CV2 cost with 1/d anchor
					t =
					  @elapsed alphaCV, min_indx, CVUnsc_curve = JS.cv_alpha(xs, cs, mhats, p0, alpha_grid, 2)
					writedlm(f, [dates[ix_start] K d N "CV2_unif" or_alpha_curve[min_indx] t alphaCV], ',')

					#Gen the CV10 cost with 1/d anchor
					t =
					  @elapsed alphaCV, min_indx, CVUnsc_curve = JS.cv_alpha(xs, cs, mhats, p0, alpha_grid, 10)
					writedlm(f, [dates[ix_start] K d N "CV10_unif" or_alpha_curve[min_indx] t alphaCV], ',')
				end

				##MSE version of alpha
				t =
				  @elapsed alphaMSE, min_indx = JS.mse_estimates(mhats, supp, p0, alpha_grid)
				writedlm(f, [dates[ix_start] K d N "MSE" or_alpha_curve[min_indx] t alphaMSE], ',')

				#Gen the Oracle cost with GM Anchor
				t =
				  @elapsed alphaOR_GM, min_indx, or_alpha_curve_GM = JS.oracle_alpha(xs, cs, mhats, ps, lams, phat_avg, alpha_grid)
				writedlm(f, [dates[ix_start] K d N "OraclePhat" or_alpha_curve_GM[min_indx] t alphaOR_GM], ',')

				#Gen the LOO cost with the GM Anchor
				t =
				  @elapsed alphaLOO, min_indx, looUnsc_curve = JS.loo_alpha(xs, cs, mhats, phat_avg, alpha_grid)
				writedlm(f, [dates[ix_start] K d N "LOO_avg" or_alpha_curve_GM[min_indx] t alphaLOO], ',')
				tempGM = looUnsc_curve[min_indx]

				# ## Optimizing the LOO Anchor
				# t =
				#   @elapsed optp0, alphaLOO, loo_val = JS.loo_anchor(xs, cs, mhats, init_sqrt_alpha = sqrt(alphaLOO), numClusters = (K >= 20 ? 20 : -1), info=true )
				# perf = JS.zbar(xs, cs, mhats, ps, lams, (optp0, alphaLOO))
				# writedlm(f, [dates[ix_start] K d N "OptAnchor" perf t alphaLOO], ',')
				# println("% Improve on Anchor:\t", 1-loo_val/looUnsc_curve[min_indx])

				#Gen the CV5 cost with the GM Anchor
				t =
				  @elapsed alphaCV, min_indx, CVUnsc_curve = JS.cv_alpha(xs, cs, mhats, phat_avg, alpha_grid, 5)
				writedlm(f, [dates[ix_start] K d N "CV5_avg" or_alpha_curve_GM[min_indx] t alphaCV], ',')

				if allCV
					#Gen the CV2 cost with the GM Anchor
					t =
					  @elapsed alphaCV, min_indx, CVUnsc_curve = JS.cv_alpha(xs, cs, mhats, phat_avg, alpha_grid, 2)
					writedlm(f, [dates[ix_start] K d N "CV2_avg" or_alpha_curve_GM[min_indx] t alphaCV], ',')

					#Gen the CV10 cost with the GM Anchor
					t =
					  @elapsed alphaCV, min_indx, CVUnsc_curve = JS.cv_alpha(xs, cs, mhats, phat_avg, alpha_grid, 10)
					writedlm(f, [dates[ix_start] K d N "CV10_avg" or_alpha_curve_GM[min_indx] t alphaCV], ',')
				end

				##MSE version of alpha with GM
				t =
				  @elapsed alphaMSe, min_indx = JS.mse_estimates(mhats, supp, phat_avg, alpha_grid)
				writedlm(f, [dates[ix_start] K d N "MSE_GM" or_alpha_curve_GM[min_indx] t alphaMSE], ',')

				## Optimizing a beta distribution anchor
				t =
				  @elapsed alphaLOO, p0LOO, loo_val = JS.loo_betaAnchor(xs, cs, mhats, alpha_grid, theta2_grid, mu_grid)
				true_perf = JS.zbar(xs, cs, mhats, ps, lams, (p0LOO, alphaLOO))
				writedlm(f, [dates[ix_start] K d N "BetaOptLOO" true_perf t alphaLOO], ',')
				println()
				println("----Beta %Improv LOO\t", 1- loo_val/tempGM, "------")
				println()

				t =
				  @elapsed alphaOR, p0OR, or_val = JS.oracle_betaAnchor(xs, cs, mhats, ps, lams, alpha_grid, theta2_grid, mu_grid)

				writedlm(f, [dates[ix_start] K d N "BetaOptOR" or_val t alphaOR], ',')

			end  #end K Loop
			flush(f)
		end #end window loop
	end #end N loop
	close(f)
	"$(outPath).csv"
end #end function