text:
const DEFAULT_SIMULATION_RELTOL=1e-3
const DEFAULT_SIMULATION_ABSTOL=1e-6
const DEFAULT_ALG = AutoTsit5(Rosenbrock23())

const EMPTY_PROBLEM = ODEProblem(() -> nothing, [0.0], (0.,1.))

### simulate scenario

"""
    sim(scenario::Scenario;
      parameters_upd::Vector{P}=Pair{Symbol, Float64}[],
      alg=DEFAULT_ALG,
      reltol=DEFAULT_SIMULATION_RELTOL,
      abstol=DEFAULT_SIMULATION_ABSTOL,
      kwargs...)

Simulate single `Scenario`. Returns [`SimResult`](@ref) type.

Example: `Scenario(model; tspan = (0., 200.), saveat = [0.0, 150., 250.]) |> sim`

Arguments:

- `scenario` : simulation scenario of type [`Scenario`](@ref)
- `parameters_upd` : constants, which overwrite both `Model` and `Scenario` constants. Default is empty vector.
- `alg` : ODE solver. See SciML docs for details. Default is AutoTsit5(Rosenbrock23())
- `reltol` : relative tolerance. Default is 1e-3
- `abstol` : relative tolerance. Default is 1e-6
- `kwargs...` : other solver related arguments supported by DiffEqBase.solve. See SciML docs for details
"""
function sim(
  scenario::Scenario;
  parameters_upd::Vector{P}=Pair{Symbol, Float64}[],
  alg=DEFAULT_ALG,
  reltol=DEFAULT_SIMULATION_RELTOL,
  abstol=DEFAULT_SIMULATION_ABSTOL,
  kwargs... # other solver arguments
) where P<:Pair

  prob = length(parameters_upd) > 0 ? update_init_values(scenario.prob, scenario.init_func, NamedTuple(parameters_upd)) : scenario.prob
  sol = solve(prob, alg; reltol = reltol, abstol = abstol,
    save_start = false, save_end = false, save_everystep = false, kwargs...)
  params_names = Symbol[first(x) for x in parameters_upd]
  return build_results(sol, scenario, params_names)
end

function build_results(sol::SciMLBase.AbstractODESolution, params_names::Vector{Symbol})
  params = if length(params_names) > 0
    @LArray sol.prob.p.constants[params_names] Tuple(params_names)
  else
    @LArray Float64[] ()
  end
  sv = sol.prob.kwargs[:callback].discrete_callbacks[1].affect!.saved_values
  return Simulation(sv, params, sol.retcode)
end

build_results(sol::SciMLBase.AbstractODESolution, scenario, params_names) = SimResult(build_results(sol, params_names), scenario)

### simulate scenario pairs

"""
    sim(scenario_pairs::Vector{P};
      parameters_upd::Vector{Pair{Symbol, Float64}}=Pair{Symbol, Float64}[],
      alg=DEFAULT_ALG,
      reltol=DEFAULT_SIMULATION_RELTOL,
      abstol=DEFAULT_SIMULATION_ABSTOL,
      parallel_type=EnsembleSerial(),
      kwargs...) where P<:Pair

Simulate multiple scenarios. Returns `Vector{Pair}`.

Example: `sim([:x => scn1, :y=>scn2, :z=>scn3])`

Arguments:

- `scenario_pairs` : vector of pairs containing names and scenarios of type [`Scenario`](@ref)
- `parameters_upd` : constants, which overwrite both `Model` and `Scenario` constants. Default is empty vector.
- `alg` : ODE solver. See SciML docs for details. Default is AutoTsit5(Rosenbrock23())
- `reltol` : relative tolerance. Default is 1e-3
- `abstol` : relative tolerance. Default is 1e-6
- `parallel_type` : type of multiple simulations parallelism. Default is no parallelism. See SciML docs for details
- `kwargs...` : other solver related arguments supported by DiffEqBase.solve. See SciML docs for details
"""
function sim(
  scenario_pairs::Vector{P};
  parameters_upd::Vector{Pair{Symbol, Float64}}=Pair{Symbol, Float64}[],
  alg = DEFAULT_ALG,
  reltol = DEFAULT_SIMULATION_RELTOL,
  abstol = DEFAULT_SIMULATION_ABSTOL,
  parallel_type=EnsembleSerial(),
  kwargs... # other arguments for OrdinaryDiffEq.solve()
) where P<:Pair

  isempty(scenario_pairs) && return SimResult[] # BRAKE

  progress_on = (parallel_type == EnsembleSerial()) # tmp fix
  p = Progress(length(scenario_pairs), dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, enabled=progress_on)

  function prob_func(prob,i,repeat)
    next!(p)
    prob_i = last(scenario_pairs[i]).prob
    init_func_i = last(scenario_pairs[i]).init_func
    length(parameters_upd) > 0 ? update_init_values(prob_i, init_func_i, NamedTuple(parameters_upd)) : prob_i
  end

  params_names = Symbol[first(x) for x in parameters_upd]

  function _output(sol,i)
    build_results(sol,last(scenario_pairs[i]), params_names),false
  end

  _reduction(u,batch,I) = (append!(u,batch),false)

  prob = EnsembleProblem(EMPTY_PROBLEM;
    prob_func = prob_func,
    output_func = _output,
    reduction = _reduction
  )

  solution = solve(prob, alg, parallel_type;
    trajectories = length(scenario_pairs),
    reltol = reltol,
    abstol = abstol,
    save_start = false,
    save_end = false,
    save_everystep = false,
    kwargs...
    )
  return [Pair{Symbol,SimResult}(first(cp), u) for (cp,u) in zip(scenario_pairs, solution.u)]
end

### simulate scenario array, XXX: do we need it?

"""
    sim(scenarios::AbstractVector{C}; kwargs...) where {C<:AbstractScenario}

Simulate multiple scenarios. Returns `Vector{Pair}`.

Example: `sim([scn1, scn2, scn3])`

Arguments:

- `scenarios` : `Vector` containing names and scenarios of type [`Scenario`](@ref)
- kwargs : other kwargs supported by `sim(scenario_pairs::Vector{Pair})`
"""
function sim(
  scenarios::AbstractVector{C};
  kwargs... # other arguments to sim(::Vector{Pair})
) where {C<:AbstractScenario}
  scenario_pairs = [Symbol("_$i") => scn for (i, scn) in pairs(scenarios)]
  return sim(scenario_pairs; kwargs...)
end

### simulate Platform

"""
    sim(platform::Platform;
      scenarios::Union{AbstractVector{Symbol}, Nothing} = nothing,
      kwargs...) where {C<:AbstractScenario}

Simulate scenarios included in platform. Returns `Vector{Pair}`.

Example: `sim(platform)`

Arguments:

- `platform` : platform of [`Platform`](@ref) type
- `scenarios` : `Vector` containing names of scenarios included in platform. Default value `nothing` stands for all scenarios in the platform
- `kwargs...` : other kwargs supported by `sim(scenario_pairs::Vector{Pair})`
"""
function sim(
  platform::Platform;
  scenarios::Union{AbstractVector{Symbol}, Nothing} = nothing,
  kwargs... # other arguments to sim(::Vector{Pair})
)
  @assert length(platform.scenarios) != 0 "Platform should contain at least one scenario."

  if isnothing(scenarios)
    scenario_pairs = [platform.scenarios...]
  else
    scenario_pairs = Pair{Symbol,AbstractScenario}[]
    for scn_name in scenarios
      @assert haskey(platform.scenarios, scn_name) "No scenario :$scn_name found in the platform."
      push!(scenario_pairs, scn_name=>platform.scenarios[scn_name])
    end
  end

  return sim(scenario_pairs; kwargs...)
end

function default_output(sol,i)
  sim = sol.prob.kwargs[:callback].discrete_callbacks[1].affect!.saved_values
  sim,false
end

default_reduction(u,batch,I) = (append!(u,batch),false)


text:
## Field data macros

export show_scalarlist, show_vectorlist, show_tensorlist

_negate(a::NTuple{M,Int64}) where {M} = map(x -> -x,a)
_cellshift(a::NTuple{M,Int64}) where {M} = map(x -> 0.5*(1-abs(x)),a)

# for Julia < 1.1 (since later versions have this built in)
_fieldtypes(T::Type) = ntuple(i -> fieldtype(T, i), fieldcount(T))

function _numberofcomponents(T::Type)
    n = 0
    for t in _fieldtypes(T)
        if t <: GridData
            n+=1
        end
    end
    return n
end

"""
    @griddata(wrapper,nctypes)

Create a basic set of constructors and functions for grid data wrapper type
`wrapper`. The argument `nctypes` is the number of cell types used
to parameterize this grid data type.
"""
macro griddata(wrapper, nctypes)

  ctype = Symbol[]
  for i in 1:eval(nctypes)
    push!(ctype,Symbol("C",string(i)))
  end

  return esc(quote

    export $wrapper

    celltype(::$wrapper{C}) where {C<:CellType} = C
    celltype(::Type{<:$wrapper{$(ctype...),NX,NY}}) where {$(ctype...),NX,NY} = $(ctype...)

    gridsize(::Type{<:$wrapper{$(ctype...),NX,NY}}) where {$(ctype...),NX,NY} = NX, NY

    griddatatype(::$wrapper) = $wrapper
    griddatatype(::Type{T}) where {T<:$wrapper{$(ctype...),NX,NY,T,DT} where {$(ctype...),NX,NY,T,DT}} = $wrapper

    # This allows easy construction from existing GridData on the same grid.
    $wrapper(C, ::GridData{NX,NY,T};dtype=T) where {NX, NY,T <: Number} = $wrapper(C, (NX, NY),dtype=dtype )

    $wrapper(C, nx::Int, ny::Int;dtype=Float64) = $wrapper(C,(nx,ny),dtype=dtype)
    (::Type{$wrapper{$(ctype...),NX,NY,T,DT}})() where {$(ctype...),NX,NY,T,DT} =
                  $wrapper($(ctype[1]), (NX, NY),dtype=T)

    # This constructor might be problematic? Introduced only because we have not
    # yet updated the regularization routines for the new GridData parameterization
    (::Type{$wrapper{$(ctype...),NX,NY,T}})() where {$(ctype...),NX,NY,T} =
                  $wrapper($(ctype[1]), (NX, NY),dtype=T)

    Base.similar(::$wrapper{$(ctype...),NX,NY,T,DT};element_type=T) where {$(ctype...),NX,NY,T,DT} =
                  $wrapper($(ctype[1]), (NX, NY),dtype=element_type)

    Base.parent(A::$wrapper) = A.data
    Base.size(A::$wrapper) = size(A.data)
    Base.parentindices(A::$wrapper) = parentindices(A.data)

    @propagate_inbounds Base.getindex(A::$wrapper,i::Int) = getindex(A.data,i)
    @propagate_inbounds Base.setindex!(A::$wrapper, v, i::Int) = setindex!(A.data,convert(eltype(A.data),v),i)

  end)

end

"""
    @generate_scalarlist(list)

Process a list of scalar grid types and generate an expanded form of the list,
complete with cell shifts. Used for constructing the regularization functions
and coordinate functions.
"""
macro generate_scalarlist(list)
    return quote
        newlist = []
        for (i,l) in enumerate($(esc(list)))
            wrapper, primaldn, dualdn = l

            # negative index shifts become positive values in this list
            pdn = _negate(primaldn)
            ddn = _negate(dualdn)

            # grid cell shifts: dn = 0 implies shift by half-cell, dn = 1 implies no shift
            pshift = _cellshift(pdn)
            dshift = _cellshift(ddn)
            push!(newlist,(Symbol(wrapper),:Primal,pdn...,pshift...))
            push!(newlist,(Symbol(wrapper),:Dual,ddn...,dshift...))
        end
        newlist
    end
end

show_scalarlist() = @generate_scalarlist SCALARLIST

"""
    @generate_vectorlist(list)

Process a list of vector grid types and generate an expanded form of the list,
complete with cell shifts. Used for constructing the regularization functions
and coordinate functions.
"""
macro generate_collectionlist(list)
    return quote
        newlist = []
        for (i,l) in enumerate($(esc(list)))
            wrapper, celltypes = l
            gtype = eval(wrapper){eval.(celltypes)...}

            row = (l...,)
            for ft in unique(_fieldtypes(gtype))
                if ft <: GridData
                    dn = _negate(indexshift(ft))
                    row = (row...,dn...)
                end
            end
            for ft in unique(_fieldtypes(gtype))
                if ft <: GridData
                    dn = _negate(indexshift(ft))
                    cshift = _cellshift(dn)
                    row = (row...,cshift...)
                end
            end
            push!(newlist,row)
        end
        newlist
    end
end

show_vectorlist() = @generate_collectionlist VECTORLIST

show_tensorlist() = @generate_collectionlist TENSORLIST

"""
    @wrapparay(wrapper,field,N)

Basic macro to develop any AbstractArray data type into a proper wrapper, with
indexing and other array-type operations.
"""
macro wraparray(wrapper, field, N)
    S = eval(wrapper)
    @assert S <: AbstractArray "Wrapped type must be a subtype of AbstractArray"

    quote
        Base.parent(A::$wrapper) = A.$field
        Base.size(A::$wrapper) = size(A.$field)
        parentindices(A::$wrapper) = parentindices(A.$field)

        if $N > 1
          function Base.show(io::IO, m::MIME"text/plain", A::$wrapper)
            println(io, "$(typeof(A)) data")
            println(io, "Printing in grid orientation (lower left is (1,1))")
            show(io,m, reverse(transpose(A.$field),dims=1))
          end
          #function Base.summary(io::IO, A::$wrapper)
          #  println(io, "$(typeof(A)) data")
          #  print(io, "Printing in grid orientation (lower left is (1,1))")
          #end
        end

        @propagate_inbounds Base.getindex(A::$wrapper, i::Int) = A.$field[i]
        @propagate_inbounds Base.setindex!(A::$wrapper, v, i::Int) = A.$field[i] = convert(eltype(A.$field), v)
        if $N > 1
          @propagate_inbounds Base.getindex(A::$wrapper, I::Vararg{Int, $N}) = A.$field[I...]
          @propagate_inbounds Base.setindex!(A::$wrapper, v, I::Vararg{Int, $N}) = A.$field[I...] = convert(eltype(A.$field), v)
        end
    end
end

ext: jl

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis
using LinearAlgebra

# Parameters
h = .1
r = 1.

# Links
origin = Origin{Float64}()
link1 = Cylinder(r, h, h*r, color = RGBA(1., 0., 0.))

# Constraints
joint1 = EqualityConstraint(Floating(origin, link1))

links = [link1]
constraints = [joint1]

mech = Mechanism(origin, links, constraints, g = 0.)

axis = [0;0;1.]
speed = 20pi #*0
setVelocity!(link1, ω = speed*axis)

function controller!(mechanism, k)
    if k==1
        setForce!(mechanism.bodies[1], F = [0;0;2.], p=[0;1.;0])
    else
        setForce!(mechanism.bodies[1], F = [0;0;0.], p=[0;0.0;0])
    end
    return
end

storage = simulate!(mech, 10., controller!, record = true)
visualize(mech, storage)


text:
const MAX_PAYLOAD_SIZE = 1460

bytes_in_flight(sender::TCPSender) = sender.bytes_in_flight

next_seqno_absolute(sender::TCPSender) = sender.next_seqno

next_seqno(sender::TCPSender) = wrap(sender.next_seqno |> UInt, sender.isn)

segments_out(sender::TCPSender) = sender.segments_out

stream_in(sender::TCPSender) = sender.stream

function fill_window!(sender::TCPSender)
    (sender.next_seqno == 0) &&
        send_seg!(sender, TCPSegment(TCPHeader(syn=true), ""))
    win = sender.win_size > 0 ? sender.win_size : 1
    remain = 0
    # (win < sender.next_seqno - sender.last_acked) && return
    while (remain = win - (sender.next_seqno - sender.last_acked)) > 0 && !sender.fin_flg
        sz = min(MAX_PAYLOAD_SIZE, remain)
        seg = TCPSegment(read!(sender.stream, sz))
        if length_in_sequence_space(seg) < win && eof(sender.stream)
            seg.header.fin = sender.fin_flg = true
        end
        (length_in_sequence_space(seg) == 0) && break
        send_seg!(sender, seg)
    end
end

function send_seg!(sender::TCPSender, seg::TCPSegment)
    header = seg.header
    header.seqno = next_seqno(sender)
    header.syn = sender.next_seqno == 0

    len = length_in_sequence_space(seg)
    sender.next_seqno += len
    sender.bytes_in_flight += len

    enqueue!(sender.segments_out, seg)
    push!(sender.outstanding_segs, seg)
    sender.timer_on = true
end

function ack_received!(sender::TCPSender, ackno::WrappingInt32, win_size::UInt16)
    abs_ackno = unwrap(ackno, sender.isn, sender.last_acked)
    sender.win_size = win_size
    (abs_ackno <= sender.last_acked || abs_ackno > sender.next_seqno) && return
    i = 0
    acc = sender.last_acked
    len = length(sender.outstanding_segs)
    while i <= len && acc < abs_ackno
        i += 1
        @inbounds acc += length_in_sequence_space(sender.outstanding_segs[i])
    end
    if acc > abs_ackno
        acc -= length_in_sequence_space(sender.outstanding_segs[i])
        i -= 1
    end
    (i == 0) && return
    sender.bytes_in_flight -= (acc - sender.last_acked)
    deleteat!(sender.outstanding_segs, 1:i)
    sender.last_acked = acc
    sender.timer_on = false
    sender.RTO = sender.initial_retransmission_timeout
    sender.ticks = 0

    fill_window!(sender)
    !isempty(sender.outstanding_segs) && (sender.timer_on = true)
end

function tick!(sender::TCPSender, ms_since_last_tick::Int)
    !sender.timer_on && return
    sender.ticks += ms_since_last_tick
    if sender.ticks >= sender.RTO && !isempty(sender.outstanding_segs)
        enqueue!(sender.segments_out, sender.outstanding_segs[1])
        if sender.last_acked == 0 || sender.win_size != 0
            sender.consecutive_retransmissions += 1
            sender.RTO <<= 1
        end
        sender.ticks = 0
    end
    # if isempty(sender.outstanding_segs)
    #     sender.timer_on = false
    #     sender.ticks = 0
    # end
end

consecutive_retransmissions(sender::TCPSender) = sender.consecutive_retransmissions

send_empty_segment!(sender::TCPSender) =
    enqueue!(sender.segments_out, TCPSegment(TCPHeader(seqno=wrap(sender.next_seqno, sender.isn)) ,""))

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
using DataDeps

register(DataDep(
  "spectf-heart",
  "http://archive.ics.uci.edu/ml/datasets/SPECTF+Heart",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/spect/SPECTF.train",
  "9c837e4fd4b7931d9be5304989e39f284f35f767518ffd4dfd42a6a519a9f30c",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:45,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
mutable struct MPSState <: QState
	sites::Vector{ITensor}
	map::QubitSiteMap
	llim::Int
	rlim::Int

	MPSState(sites::Vector{ITensor},
		     map::QubitSiteMap,
		     llim::Int,
		     rlim::Int) = new(sites, map, llim, rlim)
end #struct

function MPSState(N::Int, init::Vector{T}) where {T <:Number}
	sites = ITensor[]
	rightlink,leftlink = nothing, nothing
	for i =1:N
		if i ==1
			rightlink = Index(1, "Link, l=$(i)")
			push!(sites, ITensor(init, rightlink,Index(2, "Site, q=$(i)")))
		elseif i == N
			push!(sites, ITensor(init, leftlink, Index(2, "Site, q=$(i)")))
		else
			rightlink = Index(1, "Link, l=$(i)")
			push!(sites, ITensor(init, leftlink, rightlink, Index(2, "Site, q=$(i)")))
		end
		leftlink = rightlink
	end
	MPSState(sites,QubitSiteMap(N),0,N+1)
end
MPSState(N::Int) = MPSState(N, [1,0])

getindex(state::MPSState, n::Int) = getindex(state.sites,n)
getinds(state::MPSState, inds::Vector{Int}) = [getindex(state, i) for i ∈ inds]
setindex!(state::MPSState, T::ITensor, n::Integer) = setindex!(state.sites, T, n)
length(state::MPSState) = length(state.sites)
size(state::MPSState) = size(state.sites)
iterate(state::MPSState, itstate::Int=1) = iterate(state.sites,itstate)
copy(state::MPSState) = MPSState(copy.(state.sites), copy(state.map), state.llim, state.rlim)
numQubits(state::MPSState) = length(state)

function isapprox(state1::MPSState, state2::MPSState)
	exact1 = toExactState(state1)
	exact2 = toExactState(state2)
	return exact1 ≈ exact2
end

siteForQubit(state::MPSState, i::Int) = siteForQubit(state.map,i)
sitesForQubits(state::MPSState, inds::Vector{Int}) = sitesForQubits(state.map, inds)
qubitAtSite(state::MPSState, i::Int) = qubitAtSite(state.map, i)
qubitsAtSites(state::MPSState, inds::Vector{Int}) = qubitsAtSites(state.map, inds)

updateMap!(state::MPSState; kwargs...) = updateMap!(state.map; kwargs...)

function updateLims!(state::MPSState, leftEnd::Int, rightEnd::Int)
	state.llim >= leftEnd && (state.llim = leftEnd-1)
	state.rlim <= rightEnd && (state.rlim = rightEnd+1)
	return state
end
updateLims!(state::MPSState, touchedSite::Int) = updateLims!(state, touchedSite, touchedSite)

getQubit(state::MPSState, i::Int) = state.sites[siteForQubit(state,i)]
getQubits(state::MPSState, inds::Vector{Int}) = [getQubit(state, i) for i ∈ inds]
siteIndex(state::MPSState, i::Int) = findindex(state[i], "Site")
siteInds(state::MPSState, inds::Vector{Int}) = IndexSet([siteIndex(state, i) for i ∈ inds]...)

# TODO: Contract in a binary tree order?
toExactState(state::MPSState) = ExactState(prod(state.sites))

function normalize!(state::MPSState; kwargs...)
	center::Int = get(kwargs, :center,-1)
	if center ==-1
		center = state.llim+1
	end
	centerAtSite!(state,center)
end
# apply!(state::MPSState, gate::MeasureGate) = collapseQubits!(state, qubits(gate); reset=reset(gate))
function measure!(state::MPSState, qubits::Vector{Int}, shots::Int; kwargs...)
	binary = get(kwargs, :binary, true)
	probability = get(kwargs, :probability, false)
	if length(qubits) > 1
		localizeQubitsInOrder!(state, qubits; kwargs...)
	end
	sites = sitesForQubits(state, qubits)
	centerAtSite!(state, sites[1])
	if !probability
		if binary
			results = zeros(Int, shots, length(qubits))
			for i = 1:shots
				results[i,:] = oneShot(state, sites; kwargs...)
			end
		else
			results = zeros(Int, shots)
			for i = 1:shots
				results[i] = oneShot(state, sites; kwargs...)
			end
		end
		return results
	else
		prob = zeros(shots)
		if binary
			results = zeros(Int, shots, length(qubits))
			for i = 1:shots
				results[i,:], prob[i] = oneShot(state, sites; kwargs...)
			end
		else
			results = zeros(Int, shots)
			for i = 1:shots
				results[i], prob[i] = oneShot(state, sites; kwargs...)
			end
		end
		return results, prob
	end
end

function measure!(state::MPSState, qubit::Int, shots::Int; kwargs...)
	seed = get(kwargs, :seed, 1)
	site = siteForQubit(state, qubit)
	centerAtSite!(state, site)
	results = zeros(Int, shots)
	for i = 1: shots
		results[i] = oneShot(state, [site]; kwargs...)[1]
	end
	return results
end

function orderQubits!(state::MPSState; kwargs...)
	for q = 1:length(state)
		moveQubit!(state, q, q; kwargs...)
	end
end

function moveSite!(state::MPSState, s::Int, t::Int; kwargs...)
	while s < t
		swapSites!(state, s, s+1; kwargs...)
		s += 1
	end
	while s > t
		swapSites!(state, s-1, s; kwargs...)
		s -= 1
	end
end

moveQubit!(state::MPSState, q::Int, s::Int; kwargs...) = moveSite!(state, siteForQubit(state,q), s; kwargs...)

function localizeQubits!(state::MPSState, actpos::Vector{Int}; kwargs...)
	sortedSites = sort!(sitesForQubits(state, actpos))

	for i = 2: length(sortedSites)
		moveSite!(state, sortedSites[i], sortedSites[i-1]+1; kwargs...)
		sortedSites[i] = sortedSites[i-1]+1
	end

	return state
end

function localizeQubitsInOrder!(state::MPSState, qubits::Vector{Int}; kwargs...)
	for i = 1:length(qubits)
		moveQubit!(state, qubits[i], i; kwargs...)
	end
end

function centerAtSite!(state::MPSState, s::Int)
	while state.llim < s-1
		currIndex = state.llim+1
		curr = state[currIndex]
		if currIndex != 1
			Q,R = qr(curr, findindex(curr,"Site"), findindex(curr, "l=$(currIndex-1)"))
		else
			Q,R = qr(curr, findindex(curr,"Site"))
		end
		state[currIndex] = replacetags!(Q, "qr", "l=$(currIndex)")
		state[currIndex+1] = replacetags!(state[currIndex+1]*R , "qr", "l=$(currIndex)")
		state.llim+=1
	end

	while state.rlim > s+1
		currIndex = state.rlim-1
		curr = state[currIndex]
		if currIndex != length(state)
			Q,R = qr(curr, findindex(curr, "Site"), findindex(curr, "l=$(currIndex)"))
		else
			Q,R = qr(curr, findindex(curr, "Site"))
		end
		state[currIndex] = replacetags!(Q, "qr", "l=$(currIndex-1)")
		state[currIndex-1] = replacetags!(state[currIndex-1]*R , "qr", "l=$(currIndex-1)")
		state.rlim-=1
	end

	state.llim = s-1
	state.rlim = s+1
	state[s] /= norm(state[s]) # normalization
	return state
end

centerAtQubit!(state::MPSState, q::Int) = centerAtSite!(state, siteForQubit(state, q))

function swapSites!(state::MPSState, s1::Int, s2::Int, decomp= "svd"; kwargs...)
	s1 > s2 && ((s1, s2) = (s2, s1))
	(s2 != s1+1) && error("swapSites can only swap two neighbor sites")
	(q1, q2) = (qubitAtSite(state, s1), qubitAtSite(state, s2))
	product = state[s1] * state[s2]
	if decomp =="svd"
		if s1 > 1
			leftinds = IndexSet(findindex(product, "l=$(s1-1)"),
								findindex(product, "q=$(q2)"))
			U,S,V,u,v = svd(product, leftinds; kwargs...)
		else
			U,S,V,u,v = svd(product, findindex(product, "q=$(q2)"); kwargs...)
		end
		state[s1] = replacetags(U,"u", "l=$(s1)")
		state[s2] = replacetags(S*V,"u", "l=$(s1)")
	else #qr
		if s1 > 1
			leftinds = IndexSet(findindex(product, "l=$(s1-1)"),
								findindex(product, "q=$(q2)"))
			Q,R = qr(product, leftinds)
		else
			Q,R = qr(product, findindex(product, "q=$(q2)"))
		end
		state[s1] = replacetags(Q,"u", "l=$(s1)")
		state[s2] = replacetags(R,"u", "l=$(s1)")
	end
	updateMap!(state; s=s1, q=q2)
	updateMap!(state; s=s2, q=q1)
	state.llim >= s1 && (state.llim = s1-1)
	state.rlim <= s2 && (state.rlim = s2+1)
end

function projector(i::Int, ind::Index)
		data = zeros(2)
		data[i] = 1
		return ITensor(data, ind)
end

function oneShot(state::MPSState, sites::Vector{Int}; kwargs...)
	binary = get(kwargs, :binary, true)
	probability = get(kwargs, :probability, false)
	rng = get(kwargs, :rng, Random.GLOBAL_RNG)
	logprob = get(kwargs, :logprob, false)
	if binary
		sample = zeros(Int, length(sites))
	else
		sample = 0
	end
	logprob ? prob = 0 : prob = 1
	clamped = nothing

	for i =1:length(sites)
		ψ = state[sites[i]]
		clamped != nothing && (ψ *= clamped)
		ψ0 = ψ*projector(1, findindex(ψ, "Site"))
		prob0 = Real(scalar(ψ0 * dag(ψ0)))
		ψ1 = ψ*projector(2, findindex(ψ, "Site"))
		prob1 = Real(scalar(ψ1 * dag(ψ1)))
		total = (prob0+prob1)
		prob0 /= total
		prob1 /= total
		if rand(rng) > prob0
			if binary
				sample[i] = 1
			else
				sample += 2^(i-1)
			end
			clamped = ψ1
			logprob ? prob += log(prob1) : prob *= prob1
		else
			clamped = ψ0
			logprob ? prob += log(prob0) : prob *= prob0
		end

	end

	if probability
		return sample, prob
	else
		return sample
	end
end

function collapse!(state::MPSState, qubit::Int; kwargs...)
	reset = get(kwargs, :reset, false)
	rng = get(kwargs, :rng, Random.GLOBAL_RNG)

	site = siteForQubit(state, qubit)
	centerAtSite!(state, site)
	ψ = state[site]
	proj0 = projector(1, findindex(ψ, "Site"))
	ψ0 = ψ*proj0
	proj1 = projector(2, findindex(ψ, "Site"))
	ψ1 = ψ*proj1
	prob0 = Real(scalar(ψ0 * dag(ψ0)))
	prob1 = Real(scalar(ψ1 * dag(ψ1)))
	if rand(rng) < prob0
		result = 0
		state[site] = ψ1*proj0
	else
		result = 1
		if reset
			state[site] = ψ1*proj0
		else
			state[site] = ψ1*proj1
		end
	end
	updateLims!(state, site)
	return result
end

function collapse!(state::MPSState, qubits::Vector{Int}; kwargs...)
	results = zeros(Int, length(qubits))
	for i = 1:length(qubits)
		results[i] = collapse!(state, qubits[i]; kwargs...)
	end
	return results
end

function dag(state::MPSState)
	sites_dag = Vector{ITensor}(undef,0)
	for i =1:length(state)
		push!(sites_dag, dag(state[i]))
	end
	return MPSState(sites_dag, copy(state.map), 0, length(state)+1) # llim, rlim may be optimized
end

function ρ(state::MPSState, config::Vector{Int})
	if any(config .> 1) || any(config .<0)
		error("Invalid configuration input.")
	end
	clampedMPS = Vector{ITensor}(undef,0)
	for s = 1:length(config)
		q = qubitAtSite(state, s)
		siteindex = findindex(state[s], "Site")
		push!(clampedMPS, clamp(state[s], siteindex, config[q]+1))
	end
	return scalar(prod(clampedMPS))
end

function logprobability(state::MPSState, config::Vector{Int})
	if length(state) != length(config)
		throw(ArgumentError("state and config has different lengths"))
	end
	clamped = nothing
	logprob = 0
	centerAtSite!(state, 1)
	for i =1:length(config)
		ψ = state[i]
		clamped != nothing && (ψ *= clamped)
		ψ0 = ψ*projector(1, findindex(ψ, "Site"))
		prob0 = Real(scalar(ψ0 * dag(ψ0)))
		ψ1 = ψ*projector(2, findindex(ψ, "Site"))
		prob1 = Real(scalar(ψ1 * dag(ψ1)))
		total = (prob0+prob1)
		prob0 /= total
		prob1 /= total
		if config[qubitAtSite(state,i)]== 1
			logprob+=log(prob1)
			clamped = ψ1
		else
			logprob+=log(prob0)
			clamped = ψ0
		end
	end
	return logprob
end

function show(io::IO, state::MPSState)
	println(io, "$(length(state))-qubit MPSState")
end
function showStructure(io::IO, state::MPSState)
	printstyled(io, "MPSState\n"; bold=true, color=:red)
	for i =1: length(state)
		printstyled(io, "Site $(i):"; bold=true, color=43)
		print(io, IndexSet(state.sites[i]),"\n")
	end
	printstyled(io, "siteForQubit:"; bold=true, color=:blue)
	println(io, state.map.siteForQubit)
	printstyled(io, "qubitAtSite:"; bold=true, color=:blue)
	println(io, state.map.qubitAtSite)
	printstyled(io, "llim = $(state.llim), rlim = $(state.rlim)\n", bold=true, color=:blue)
end

function showData(io::IO,state::MPSState)
	printstyled("--- MPS state: ---\n"; bold=true, color=5)
	for i =1: length(state)
		printstyled(io, "Site $(i):\n"; bold=true, color=43)
		print(io, state[i],"\n")
	end
end
showData(state::MPSState) = showData(stdout, state)

function project!(state1::MPSState, state2::MPSState)
	if length(state1) != length(state2)
		error("Two states do not have same length")
	end
	normalize!(state1)
	normalize!(state2)
	product = contract(state1[1], dag(state2[1]), findindex(state1[1], "Site"), findindex(state2[1], "Site"))
	for location = 2: length(state1)
		product = product* state1[location]
		product = contract(product, dag(state2[location]), findindex(product, "Site"), findindex(state2[location], "Site"))
	end
	return scalar(product)
end


text:
export
    PolicyTable,
    get_utilities,
    CoordinationTable,
    make_coordination_table,
    get_coordination_sense,
    get_action_sense

#=
This file contains types that wrap qvalue and coordination tables
and related functions.
=#

#=
PolicyTable type contains all information needed to use a Qvalue table.
=#
type PolicyTable
    q_table::Matrix{Float64}
    grid::RectangleGrid
    num_states::Int64
    function PolicyTable(q_table::Matrix{Float64}, grid::RectangleGrid, num_states::Int64)
        return new(q_table, grid, num_states)
    end
end

#=
Description:
Use the policy table and belief state to calcualte the utilities of the possible joint actions

Parameters:
- policy_table: the policy table giving Q values
- belief: a sparse matrix of probabilities over states

Return Value:
- utilities: the utilities of the joint actions
    note that these utilities correspond to the joint actions
    based on their order.
=#
function get_utilities(policy_table::PolicyTable, belief_state::SparseMatrixCSC{Float64, Int64})
    # check to make sure that the indices of the belief state correspond
    # with rows of the q_table. If the last row index of the belief state
    # is in the q_table, then all belief indices must be in the table.
    if belief_state.rowval[end] > length(policy_table.q_table)
        q_table_size = size(policy_table.q_table)
        msg = "policy_table.q_table of size: $(q_table_size) at index: $(belief_state.rowval[end])"
        throw(BoundsError(msg))
    end

    # create container for the actions
    num_actions = size(policy_table.q_table, 2)
    action_utilities = zeros(num_actions)

    if belief_state[end, 1] != 1
        # iterate over nonzero rows of belief state, multiplying the belief probability
        # that we are in that state by all the values of the q_table
        # summing over these values for the entire table gives the expected value
        # of each action, which we call here the action utilities
        for ridx in belief_state.rowval
            action_utilities += reshape(belief_state[ridx] * policy_table.q_table[ridx, :], num_actions)
        end
    end
    return action_utilities
end

#=
CoordinationTable type contains all information needed to use a coordination table
=#
type CoordinationTable
    c_table::Array{Float64}
    grid::RectangleGrid
    num_states::Int64
    function CoordinationTable(c_table::Array{Float64}, grid::RectangleGrid, num_states::Int64)
        return new(c_table, grid, num_states)
    end
end

#=
Description:
When a coordination table is made from one policy it is created online in the following
section. If a coordination table is made from multiple polies, it is created offline as
presented in offline_coordination_tables.jl. This section builds a coordination table using
a q table and joint actions matrix. Does this by finding the optimal joint action for each row,
and checking if the planes are turning in the same direction or opposite directions.

Parameters:
- q_table: table of Q values shape (num_states, num_joint_actions)
- joint_actions: matrix where each row corresponds to a uav and each
    column to a possible joint action

Return Value:
- returns the coordination table
=#
function make_coordination_table(q_table::Matrix{Float64}, joint_actions::Matrix{Float64})
    # allocate coordination table
    c_table = zeros(size(q_table, 1), 1)

    # go through each row setting the value
    for ridx in 1:size(q_table, 1)
        utilities = q_table[ridx, :]

        # make it so that we do not select joint actions
        # with straight or COC as one of the options
        for aidx in 1:size(joint_actions, 2)
            if  joint_actions[1, aidx] == 0. || joint_actions[2, aidx] == 0. ||
                joint_actions[1, aidx] == COC_ACTION || joint_actions[2, aidx] == COC_ACTION
                utilities[aidx] = -Inf
            end
        end

        aidx = indmax(utilities)
        actions = joint_actions[:, aidx]

        # if the signs of the two actions match then
        # signal 1 for requiring same sense actions
        # note that we removed COC and straight above so
        # does not come up here therefore case SAME includes:
        # right + right
        # left + left
        if sign(actions[1]) == sign(actions[2])
            c_table_value = 1

        # otherwise, signal to only do opposite sense
        # therefore case DIFF includes:
        # right + left
        # left + right
        else
            c_table_value = -1
        end

        c_table[ridx] = c_table_value
    end

    return c_table
end

#=
Description:
Determines which sense should be used from the coordination table

Parameters:
- coordination_table: the coordination table to derive the sense from
- belief_state: belief state over the coordination table

Return Value:
- sense: symbol denoting sense
=#
function get_coordination_sense(coordination_table::CoordinationTable,
                        belief_state::SparseMatrixCSC{Float64, Int64})
    # decide which sense has the most weight on it in the belief state
    sense_same_weight = 0
    sense_diff_weight = 0
    sense_none_weight = 0
    for ridx in belief_state.rowval

        # index into coordination table to see
        # which sense this belief state corresponds to
        cur_sense = coordination_table.c_table[ridx, 1]

        # get the weight on this state
        belief_weight = belief_state[ridx]

        # accumulate it for the corresponding sense
        if cur_sense == -1
            sense_diff_weight += belief_weight
        elseif cur_sense == 1
            sense_same_weight += belief_weight
        else
            sense_none_weight += belief_weight
        end
    end

    # and choose whichever sense has highest weight
    if sense_same_weight > sense_diff_weight && sense_same_weight > sense_none_weight
        sense = :same_sense
    elseif sense_diff_weight > sense_none_weight
        sense = :different_sense
    else
        sense = :neither_sense
    end

    return sense
end

#=
Description:
Given two actions, return their sense

Parameters:
- action_1: first action
- action_2: second action

Return Value:
- sense: symbol denoting the sense of the actions
=#
function get_action_sense(action_1::Float64, action_2::Float64)
    # joint actions involving coc and straight are neither sense
    # because no matter what, they are always considered when
    # deciding which action to take and signal to sen
    if action_1 == COC_ACTION || action_2 == COC_ACTION ||
            action_1 == 0.0 || action_2 == 0.0
        sense = :neither_sense

    # same sense: right + right or left + left
    elseif sign(action_1) == sign(action_2)
        sense = :same_sense

    # left + right or right + left
    else
        sense = :different_sense
    end
    return sense
end


text:
module Shamir

using Polynomials

function l(x, j, k)
"""
    Create a Lagrange basis polynomial
    Reference: https://wikimedia.org/api/rest_v1/media/math/render/svg/6e2c3a2ab16a8723c0446de6a30da839198fb04b
"""
    polys = []
    for m in 1:k
        if m != j
            d = x[j] - x[m]
            r = Poly([-1 * x[m], 1]) / d
            push!(polys, r)
        end
    end

    #println(polys)

    return prod(polys)
end

function L(x, y, k)
"""
    Create a linear combination of Lagrange basis polynomials
    Reference: https://wikimedia.org/api/rest_v1/media/math/render/svg/d07f3378ff7718c345e5d3d4a57d3053190226a0
"""

    s = []

    for j in 1:k
        r = y[j] * l(x, j, k)
        push!(s, r)
    end

    #println(sum(s))

    return sum(s)
end

function construct_shares(n, production_poly)
"""
    Create shares of the secret

    Parameters:
    n: Total number of shares
    production_poly : The Polynomial created with the secret
                    to produce secret shares
"""

    share = []

    for i in 1:n
        push!(share, [i, production_poly(i)])
    end

    return share
end

function recover_secret(shares, n, k, p)
"""
    Recover the secret by finding the coefficient of x^0 i.e.
    when the x=0 in the polynomial.

    Parameters:
    shares : An array of shares of the individuals.
    n : total number of shares
    k: minimum number of shares required to unravel the secret
    p: Field number (to restrict the computation space)
"""
    if length(shares) < k
        throw("Need more parties!")
    end

    x = []
    y = []
    for i in 1:n
        push!(x, shares[i][1])
        push!(y, shares[i][2])
    end

    f = L(x, y, k)
    f = mod(f(0), p)

    return f
end

export Shamir

end


text:
function nonconv2sumofconv_v2(ps)

    V = copy(ps.Vertices[1])
    setPuntos = 1:size(V,1)
    flagConvex = polyShape.isPolyConvex(ps)[1]

    matPuntosConvHull = poly2D.verticesConvHull(V)
    setPuntosConvHull = sort(unique(matPuntosConvHull))
    setPuntosNoConvexos = setdiff(setPuntos,setPuntosConvHull)

    set_r = setPuntos
    cond = true
    ps_sub = PolyShape([],0)
    ps_r = PolyShape([],0)
    contador = 0
    V_r = copy(V)
    while cond
        if length(setPuntosNoConvexos) >= 1
            contador = contador + 1
            pnc = maximum(setPuntosNoConvexos)
            set_r = setdiff(set_r, pnc)
            set_np = copy(pnc)
            largo_r = length(set_r)
            for i = 1:largo_r
                if pnc > maximum(set_r)
                    p_adic = minimum(set_r)
                else
                    p_adic = minimum(set_r[set_r.>pnc])
                end
                set_np = union(set_np, p_adic)
                if length(set_np) >=3
                    V_np = copy(V_r[set_np,:])
                    flagConvex_np = poly2D.checkConvex(V_np)
                    flagInscrito_np = polyShape.isPolyInPoly(PolyShape([V_np],1), ps)
                    flagFactible = (flagConvex_np && flagInscrito_np)
                    if !flagFactible
                        set_np = setdiff(set_np, p_adic)
                    end
                end
                set_r = setdiff(set_r,p_adic)
            end

            ps_sub.Vertices = push!(ps_sub.Vertices, V_r[set_np,:])
            ps_sub.NumRegions = length(ps_sub.Vertices)

            ps_r = polyShape.polyDifference_v3(ps, ps_sub)
            V_r = copy(ps_r.Vertices[1])
            set_r = 1:size(V_r,1)
            matPuntosConvHull = poly2D.verticesConvHull(V_r)
            setPuntosConvHull = sort(unique(matPuntosConvHull))
            setPuntosNoConvexos = setdiff(set_r,setPuntosConvHull)

        else
            cond = false
            ps_sub.Vertices = push!(ps_sub.Vertices, V_r[set_r,:])
            ps_sub.NumRegions = length(ps_sub.Vertices)
        end

    end

    SP  = Array{SubPoly,1};
    for i = 1:ps_sub.NumRegions
        V_i = ps_sub.Vertices[i]
        set_j = []
        for j = 1:size(V_i,1)
            p_j = V_i[j,:]
            for k = 1:size(V,1)
                p_k = V[k,:]
                dist_jk = sqrt(sum((p_j .- p_k).^2))
                if dist_jk < .1
                    set_j = union(set_j,k)
                end
            end
        end
        if i==1
            SP = Array{SubPoly,1}([SubPoly(set_j,[])]);
        else
            push!(SP, SubPoly(set_j,[]));
        end
    end

    for i=1:ps_sub.NumRegions
        vecLadoComun=zeros(ps_sub.NumRegions,1);
        for j=1:ps_sub.NumRegions
            if i != j
                si = SP[i].points;
                sj = SP[j].points;
                aux=ismembern(si, sj);
                lc=[];
                if sum(aux)==2
                    vecPos=1:length(aux);
                    vecPos_=vecPos[aux];
                    if aux[1]==1 && aux[end]==1
                        lc=[vecPos_[end] vecPos_[1]];
                    else
                        lc=vecPos_;
                    end
                    vecLadoComun[j]=lc[1];
                end
            end
        end
        setfield!(SP[i], :ladoComun, vecLadoComun);
    end

    return ps_sub, SP

end


text:
# This file is a part of ParallelProcessingTools.jl, licensed under the MIT License (MIT).

_run_on_procs(f, procsel::Integer) = remotecall_fetch(f, procsel)

_run_on_procs(f, procsel::AbstractArray) =
    fetch.([remotecall(f, pid) for pid in procsel])

"""
    @onprocs procsel expr

Executes `expr` in parallel on all processes in `procsel`. Waits until
all processes are done. Returns all results as a vector (or as a single
scalar value, if `procsel` itself is a scalar).

Example:

```julia
using Distributed
addprocs(2)
workers() == @onprocs workers() myid()
```
"""
macro onprocs(procsel, expr)
    f = esc(:(()->($expr)))
    quote
        let procsel = $(esc(procsel)), f = $f
            _run_on_procs(f, procsel)
        end
    end
end
export @onprocs

"""
    @mp_async expr

Run `expr` asynchronously on a worker process.

Compatible with `@sync`.

Equivalent to `Distributed.@spawn expr` on Julia <= v1.2, equivalent to
`Distributed.@spawn :any expr` on Julia >= v1.3.
"""
macro mp_async(expr)
    # Code taken from Distributed.@spawn:
    thunk = esc(:(()->($expr)))
    var = esc(Base.sync_varname)
    quote
        local ref = Distributed.spawn_somewhere($thunk)
        if $(Expr(:isdefined, var))
            push!($var, ref)
        end
        ref
    end
end
export @mp_async


text:
# Split a list of ratings into a training and test set, with at most
# target_percentage * length(ratings) in the test set. The property we want to
# preserve is: any user in some rating in the original set of ratings is also
# in the training set and any item in some rating in the original set of ratings
# is also in the training set. We preserve this property by iterating through
# the ratings in random order, only adding an item to the test set only if we
# haven't already hit target_percentage and we've already seen both the user
# and the item in some other ratings.
function split_ratings(ratings::Array{Rating,1},
                       target_percentage=0.10)
    seen_users = Set()
    seen_items = Set()
    training_set = (Rating)[]
    test_set = (Rating)[]
    shuffled = shuffle(ratings)
    for rating in shuffled
        if in(rating.user, seen_users) && in(rating.item, seen_items) && length(test_set) < target_percentage * length(shuffled)
            push!(test_set, rating)
        else
            push!(training_set, rating)
        end
        push!(seen_users, rating.user)
        push!(seen_items, rating.item)
    end
    return training_set, test_set
end

# RMSE measures the average difference between predicted and actual ratings.
# This function measures RMSE on a test subset of the rating set that isn't used
# for training the model. As an optional final parameter, you can pass in a rank
# that's less than the model's rank to see if you should have stopped training
# features earlier.
function rmse(rating_set::RatingSet, model::RatingsModel; rank=nothing)
    total = 0.0
    if rank == nothing
        rank = size(model.S,1)
    else
        rank = min(size(model.S,1), rank)
    end
    for rating in rating_set.test_set
        predicted = 0.0
        for i=1:rank
            predicted += model.U[rating.user,i] * model.S[i] * model.V[rating.item,i]
        end
        total += (predicted - rating.value) ^ 2
    end
    sqrt(total/length(rating_set.test_set))
end

# Reduce the rank of a model. If you have a model of rank 40 but find that
# rmse(rating_set, model, rank=30) < rmse(rating_set, model, rank=40), it may
# make sense to throw away the extra 10 dimensions by running truncate_model!(model, 30).
function truncate_model!(model::RatingsModel, rank)
    current_rank = size(model.U, 2)
    if current_rank > rank
        model.U = model.U[:,1:rank]
        model.V = model.V[:,1:rank]
    end
end

# Cosine similarity between two vectors x and y. The higher the cosine
# similarity, the smaller the angle between x and y.
function cosine_similarity(x,y)
    sum, norm1, norm2 = 0.0, 0.0, 0.0
    for (a,b) in zip(x,y)
        sum += a * b
        norm1 += a * a
        norm2 += b * b
    end
    return sum / (sqrt(norm1) * sqrt(norm2))
end

# All items that the model was trained with.
function items(model::RatingsModel)
    keys(model.item_to_index)
end

# All users that the model was trained with.
function users(model::RatingsModel)
    keys(model.user_to_index)
end

# Extract the feature vector for an item.
function item_features(model::RatingsModel, item)
    model.V[model.item_to_index[item],:]
end

# Extract the feature vector for a user.
function user_features(model::RatingsModel, user)
    model.U[model.user_to_index[user],:]
end

# Returns a list of items, sorted by their values for the kth feature. This can
# help you interpret what the kth feature actually means.
function show_items_by_feature(model::RatingsModel, k)
    sort([i for i in items(model)], by=(m -> model.V[model.item_to_index[m],k]))
end

# Search for a particular item by a case-insensitive substring.
function item_search(model::RatingsModel, text)
    results = (AbstractString)[]
    r = Regex(".*?$(lowercase(text)).*")
    for item in items(model)
        if ismatch(r, lowercase(item))
            push!(results, item)
        end
    end
    return results
end

# Return the most similar items to a given item based on cosine similarity.
function similar_items(model::RatingsModel, item; max_results=10)
    features = item_features(model, item)
    sort([i for i in items(model)], by=(m -> -cosine_similarity(item_features(model, m), features)))[1:max_results]
end

# Return the most similar users to a given user based on cosine similarity.
function similar_users(model::RatingsModel, user; max_results=10)
    features = user_features(model, user)
    sort([i for i in users(model)], by=(u -> -cosine_similarity(user_features(model, u), features)))[1:max_results]
end

# Return all of the ratings in the training set for a given user as (item, rating) tuples.
function user_ratings(rating_set::RatingSet, user)
    user_id = rating_set.user_to_index[user]
    id_to_item = [id => item for (item, id) in rating_set.item_to_index]
    extract_rating = rating -> (id_to_item[rating.item], rating.value)
    is_user = rating -> rating.user == user_id
    map(extract_rating, sort(filter(is_user, rating_set.training_set), by=r -> -r.value))
end

# Return the the model's predicted rating for a given user and item.
function get_predicted_rating(model::RatingsModel, user, item)
    sum([model.U[model.user_to_index[user],r] * model.S[r] * model.V[model.item_to_index[item],r] for r=1:size(model.S,1)])
end


text:
field(k::QuoteNode, ::Type{T}) where T = :($T(kwargs[$k]))
field(k::QuoteNode, ::Type{Vector{T}}) where T = :($T.(kwargs[$k]))

function symbolize(dict::Dict{String, Any})
    new = Dict{Symbol, Any}()
    for (k, v) in dict
        new[Symbol(k)] = symbolize(v)
    end
    new
end
symbolize(t::Any) = t

macro constructors(T)
    TT = eval(T)
    args = map(f -> field(QuoteNode(f), fieldtype(TT, f)), fieldnames(TT))
    quote
        function $(esc(T))(; kwargs...)
            $(esc(T))($(args...))
        end
        $(esc(T))(d::Dict{Symbol, Any}) = $(esc(T))(; d...)
        $(esc(T))(d::Dict{String, Any}) = $(esc(T))(; symbolize(d)...)
        $(esc(T))(d::$(esc(T))) = d
        function StructTypes.keyvaluepairs(x::$(esc(T)))
            d = Dict{Symbol, Any}()
            for k = fieldnames(typeof(x))
                d[k] = getfield(x, k)
            end
            d
        end
        StructTypes.StructType(::Type{$(esc(T))}) = StructTypes.DictType()
    end
end

include("items.jl")
include("player.jl")

@constructors Recipe
@constructors Item
@constructors Effect
@constructors InventoryItem


text:
mutable struct Tensor
    inner::Ptr{primitivTensor_t}

    function Tensor(inner::Ptr{primitivTensor_t})
        @assert inner != C_NULL
        tensor = new(inner)
        finalizer(tensor, obj -> @call(DeleteTensor, (Ptr{primitivTensor_t},), obj.inner))
        tensor
    end

    function Tensor()
        tensor_ptr_ptr = Ptr{primitivTensor_t}[0]
        @call(CreateTensor, (Ptr{Ptr{primitivTensor_t}},), tensor_ptr_ptr)
        Tensor(tensor_ptr_ptr[1])
    end
end

"""
    valid(tensor)

Check whether the object is valid or not.

# Arguments
- `tensor::Tensor` tensor object
"""
function valid(tensor::Tensor)
    retval_ptr = Ref{PRIMITIV_C_BOOL}(0)
    @call(IsValidTensor, (Ptr{primitivTensor_t}, Ptr{Cuint}), tensor.inner, retval_ptr)
    retval_ptr[] == PRIMITIV_C_TRUE
end

"""
    shape(tensor)

Returns the shape of the Tensor.

# Arguments
- `tensor::Tensor` tensor object
"""
function shape(tensor::Tensor)
    shape_ptr_ptr = Ptr{primitivShape_t}[0]
    @call(GetTensorShape, (Ptr{primitivTensor_t}, Ptr{Ptr{primitivShape_t}}), tensor.inner, shape_ptr_ptr)
    Shape(shape_ptr_ptr[])
end

"""
    to_float(tensor)

Retrieves one internal value in the tensor.

# Arguments
- `tensor::Tensor` tensor object
"""
function to_float(tensor::Tensor)
    retval_ptr = Ref{Float32}(0)
    @call(EvaluateTensorAsFloat, (Ptr{primitivTensor_t}, Ptr{Cfloat}), tensor.inner, retval_ptr)
    retval_ptr[]
end

"""
    reset!(tensor)

Reset internal values using a constant.

# Arguments
- `tensor::Tensor` tensor object
"""
function reset!(tensor::Tensor, k::AbstractFloat)
    @call(ResetTensor, (Ptr{primitivTensor_t}, Cfloat), tensor.inner, k)
end


text:
"""
General Green's function.
"""

abstract type TimeDomain end
struct ImTime <: TimeDomain end
struct ReTime <: TimeDomain end
struct ImFreq <: TimeDomain end
struct ReFreq <: TimeDomain end
struct DLRFreq <: TimeDomain end
const IMTIME = ImTime()
const RETIME = ReTime()
const IMFREQ = ImFreq()
const REFREQ = ReFreq()
const DLRFREQ = DLRFreq()

abstract type InterpMethod end
struct DefaultInterp <: InterpMethod end
struct LinearInterp <: InterpMethod end
struct DLRInterp <: InterpMethod end
const DEFAULTINTERP = DefaultInterp()
const LINEARINTERP = LinearInterp()
const DLRINTERP = DLRInterp()

InterpMethod(GT::Type{<:CompositeGrids.AbstractGrid}, MT::Type{<:LinearInterp}) = CompositeGrids.Interp.LinearInterp()
InterpMethod(GT::Type{<:CompositeGrids.AbstractGrid}, MT::Type{<:DefaultInterp}) = CompositeGrids.Interp.InterpStyle(GT)
InterpMethod(GT, MT::Type{<:DLRInterp}) = DLRInterp()

"""
Green's function with two external legs that has in-built Discrete Lehmann Representation.
#Parameters:
- 'T': type of data
- 'TType': type of time domain, TType<:TimeDomain
- 'TGT': type of time grid
- 'SGT': type of space grid

#Members:
- 'name': Name of green's function
- 'color': Number of different species of Green's function (such as different spin values)
- 'dlrGrid': In-built Discrete Lehmann Representation
- 'timeGrid': Time or Frequency grid
- 'spaceType': Whether the Green's function is in coordinate space/momentum space
- 'spaceGrid': Coordinate or momentum grid
- 'instant': Instantaneous part of Green's function that is proportional to δ(τ) in τ space.
- 'dynamic': Dynamic part of Green's function
- 'instantError': Error of instantaneous part
- 'dynamicError': Error of dynamic part
"""
mutable struct Green2DLR{T<:Number,Type<:TimeDomain,TGT,SGT}
    name::Symbol
    # isFermi::Bool
    # β::Float64
    color::Int
    dlrGrid::DLRGrid

    #########    Mesh   ##############

    # timeType::DataType
    # timeSymmetry::Symbol
    timeGrid::TGT

    spaceType::Symbol
    spaceGrid::SGT

    ###########     data   ###########
    instant::Array{T,3}
    dynamic::Array{T,4}

    ####### statistical error handling #####
    instantError::Array{T,3}
    dynamicError::Array{T,4}

    """
        function Green2DLR{T}(name::Symbol, timeType::TT, β, isFermi::Bool, Euv, spaceGrid, color::Int = 1;
            timeSymmetry::Symbol = :none, rtol = 1e-8, kwargs...
        ) where {T<:Number, TT<:TimeDomain}

    Create two-leg Green's function on timeGrid, spaceGrid and color, with in-built DLR.

    #Arguements
    - 'name': Name of green's function.
    - 'timeType': type of time domain, TT<:TimeDomain
    - 'β': Inverse temperature
    - 'isFermi': Particle is fermi or boson
    - 'Euv': the UV energy scale of the spectral density
    - 'spaceGrid': k/x grid
    - 'color': Number of different species of Green's function (such as different spin values). Default: 1 color
    - 'timeSymmetry': Whether the Green's function has particle-hole symmetry, anti-particle-hole symmetry or none of them
    - 'rtol': tolerance absolute error
    #Optional Arguments
    - 'timeGrid': τ/n/ω  grid. Default: DLR grid in timeType (:τ/:n/:ω)
    - 'instant': Instantaneous part of Green's function that is proportional to δ(τ) in τ space. Default: 0 everywhere
    - 'dynamic': Dynamic part of Green's function. Default: 0 everywhere
    - 'instantError': Error of instantaneous part. Default: 0 everywhere
    - 'dynamicError': Error of dynamic part. Default: 0 everywhere
    """
    function Green2DLR{T}(name::Symbol, timeType::TT, β, isFermi::Bool, Euv, spaceGrid, color::Int = 1;
        timeSymmetry::Symbol = :none, rtol = 1e-8, kwargs...
    ) where {T<:Number, TT<:TimeDomain}
        # @assert spaceType == :k || spaceType == :x
        @assert timeSymmetry == :ph || timeSymmetry == :pha || timeSymmetry == :none

        spaceType = :k #TODO: replace it with spaceGrid.type after we have a spaceGrid package

        dlrGrid = DLRGrid(Euv, β, rtol, isFermi, timeSymmetry)

        # println(keys(kwargs))

        if :timeGrid in keys(kwargs)
            givenTimeGrid = kwargs[:timeGrid]
            if givenTimeGrid isa AbstractVector
                timeGrid = CompositeGrids.SimpleG.Arbitrary{eltype(givenTimeGrid)}(givenTimeGrid)
            elseif givenTimeGrid isa CompositeGrids.AbstractGrid
                timeGrid = givenTimeGrid
            else
                error("Input timeGrid has to be Vector or Composite grid")
            end
        else
            if TT == ImFreq
                bareTimeGrid = dlrGrid.n
            elseif TT == ImTime
                bareTimeGrid = dlrGrid.τ
            elseif TT == DLRFreq
                bareTimeGrid = dlrGrid.ω
            else
                error("$TimeType is not supported!")
            end
            timeGrid = CompositeGrids.SimpleG.Arbitrary{eltype(bareTimeGrid)}(bareTimeGrid)
        end
        #if timeGrid is nothing, then set it to be the DLR grid, which is a vector of Integer or Float64

        if TT == DLRFreq
            #@assert length(timeGrid.grid) == dlrGrid.size "The size of the DLR grid should match the DLR rank = $(dlrGrid.size)."
        elseif TT == ImFreq
            @assert eltype(timeGrid.grid) <: Int "Matsubara frequency grid is expected to be integers!"
        end

        if spaceGrid isa AbstractVector
            compSpaceGrid = CompositeGrids.SimpleG.Arbitrary{eltype(spaceGrid)}(spaceGrid)
        elseif spaceGrid isa CompositeGrids.AbstractGrid
            compSpaceGrid = spaceGrid
        else
            error("Input spaceGrid has to be Vector or Composite grid")
        end

        instant = Array{T,3}(undef, 0, 0, 0)
        dynamic = Array{T,4}(undef, 0, 0, 0, 0)
        instantError = Array{T,3}(undef, 0, 0, 0)
        dynamicError = Array{T,4}(undef, 0, 0, 0, 0)

        gnew = new{T,TT,typeof(timeGrid),typeof(compSpaceGrid)}(
            name, color, dlrGrid,
            timeGrid,
            spaceType, compSpaceGrid,
            instant, dynamic,
            instantError, dynamicError)
        return set!(gnew; kwargs...)
    end
end

function Base.getproperty(obj::Green2DLR{T,TT,TGT,SGT}, sym::Symbol) where {T,TT,TGT,SGT}
    if sym === :isFermi
        return obj.dlrGrid.isFermi
    elseif sym === :β
        return obj.dlrGrid.β
    elseif sym === :timeType
        return TT
    elseif sym === :timeSymmetry
        return obj.dlrGrid.symmetry
    else # fallback to getfield
        return getfield(obj, sym)
    end
end

function Base.size(green::Green2DLR)
    return (green.color, green.color, size(green.spaceGrid), size(green.timeGrid))
end

function set!(green::Green2DLR; kwargs...)
    dynamicSize = size(green)
    instantSize = size(green)[1:3]
    if :dynamic in keys(kwargs) && isempty(kwargs[:dynamic]) == false
        green.dynamic = reshape(kwargs[:dynamic], Tuple(dynamicSize))
    end
    if :instant in keys(kwargs) && isempty(kwargs[:instant]) == false
        green.instant = reshape(kwargs[:instant], Tuple(instantSize))
    end
    if :dynamicError in keys(kwargs) && isempty(kwargs[:dynamicError]) == false
        green.dynamicError = reshape(kwargs[:dynamicError], Tuple(dynamicSize))
    end
    if :instantError in keys(kwargs) && isempty(kwargs[:instantError]) == false
        green.instantError = reshape(kwargs[:instantError], Tuple(instantSize))
    end
    return green
end

"""
    function toTau(green::Green2DLR, targetGrid =  green.dlrGrid.τ)
Convert Green's function to τ space by Fourier transform.
If green is already in τ space then it will be interpolated to the new grid.

#Arguements
- 'green': Original Green's function
- 'targetGrid': Grid of outcome Green's function. Default: DLR τ grid
"""
function toTau(green::Green2DLR, targetGrid = green.dlrGrid.τ)

    if targetGrid isa AbstractVector
        targetGrid = CompositeGrids.SimpleG.Arbitrary{eltype(targetGrid)}(targetGrid)
    end

    # do nothing if the domain and the grid remain the same
    if green.timeType == ImTime && length(green.timeGrid.grid) ≈ length(targetGrid.grid) && green.timeGrid.grid ≈ targetGrid.grid
        return green
    end
    if isempty(green.dynamic) # if dynamic data has not yet been initialized, there is nothing to do
        return green
    end

    if (green.timeType == ImTime)
        dynamic = tau2tau(green.dlrGrid, green.dynamic, targetGrid.grid, green.timeGrid.grid; axis = 4)
    elseif (green.timeType == ImFreq)
        dynamic = matfreq2tau(green.dlrGrid, green.dynamic, targetGrid.grid, green.timeGrid.grid; axis = 4)
    elseif (green.timeType == DLRFreq)
        dynamic = dlr2tau(green.dlrGrid, green.dynamic, targetGrid.grid; axis = 4)
    end

    return Green2DLR{eltype(dynamic)}(
        green.name, IMTIME,green.β, green.isFermi, green.dlrGrid.Euv, green.spaceGrid, green.color;
        timeSymmetry = green.timeSymmetry, timeGrid = targetGrid, rtol = green.dlrGrid.rtol,
        dynamic = dynamic, instant = green.instant)
end

"""
    function toMatFreq(green::Green2DLR, targetGrid =  green.dlrGrid.n)
Convert Green's function to matfreq space by Fourier transform.
If green is already in matfreq space then it will be interpolated to the new grid.

#Arguements
- 'green': Original Green's function
- 'targetGrid': Grid of outcome Green's function. Default: DLR n grid
"""
function toMatFreq(green::Green2DLR, targetGrid = green.dlrGrid.n)

    if targetGrid isa AbstractVector
        targetGrid = CompositeGrids.SimpleG.Arbitrary{eltype(targetGrid)}(targetGrid)
    end

    # do nothing if the domain and the grid remain the same
    if green.timeType == ImFreq && length(green.timeGrid.grid) ≈ length(targetGrid.grid) && green.timeGrid.grid ≈ targetGrid.grid
        return green
    end
    if isempty(green.dynamic) # if dynamic data has not yet been initialized, there is nothing to do
        return green
    end

    if (green.timeType == ImFreq)
        dynamic = matfreq2matfreq(green.dlrGrid, green.dynamic, targetGrid.grid, green.timeGrid.grid; axis = 4)
    elseif (green.timeType == ImTime)
        dynamic = tau2matfreq(green.dlrGrid, green.dynamic, targetGrid.grid, green.timeGrid.grid; axis = 4)
    elseif (green.timeType == DLRFreq)
        dynamic = dlr2matfreq(green.dlrGrid, green.dynamic, targetGrid.grid; axis = 4)
    end

    return Green2DLR{eltype(dynamic)}(
        green.name, IMFREQ, green.β, green.isFermi, green.dlrGrid.Euv, green.spaceGrid, green.color;
        timeSymmetry = green.timeSymmetry, timeGrid = targetGrid, rtol = green.dlrGrid.rtol,
        dynamic = dynamic, instant = green.instant)

end

"""
    function toDLR(green::Green2DLR)
Convert Green's function to dlr space.

#Arguements
- 'green': Original Green's function
"""
function toDLR(green::Green2DLR)

    # do nothing if the domain and the grid remain the same
    if green.timeType == DLRFreq
        return green
    end
    if isempty(green.dynamic) # if dynamic data has not yet been initialized, there is nothing to do
        return green
    end

    if (green.timeType == ImTime)
        dynamic = tau2dlr(green.dlrGrid, green.dynamic,  green.timeGrid.grid; axis = 4)
    elseif (green.timeType == ImFreq)
        dynamic = matfreq2dlr(green.dlrGrid, green.dynamic,  green.timeGrid.grid; axis = 4)
    end

    return Green2DLR{eltype(dynamic)}(
        green.name,DLRFREQ, green.β, green.isFermi, green.dlrGrid.Euv, green.spaceGrid, green.color;
        timeSymmetry = green.timeSymmetry, timeGrid = green.dlrGrid.ω, rtol = green.dlrGrid.rtol,
        dynamic = dynamic, instant = green.instant)

end

"""
    function instant(green::Green2DLR{DT,TT,TGT,SGT}, space, color1::Int, color2::Int=color1; spaceMethod::SM = DEFAULTINTERP) where {DT,TT,TGT,SGT,SM}

Find value of Green's function's instant part at given color and k/x by interpolation.
Interpolation method is by default depending on the grid, but could also be chosen to be linear.

#Argument
- 'green': Green's function
- 'space': Target k/x point
- 'color1': Target color1
- 'color2': Target color2
- 'spaceMethod': Method of interpolation for space.
"""
function instant(green::Green2DLR{DT,TT,TGT,SGT}, space, color1::Int, color2::Int=color1; spaceMethod::SM = DEFAULTINTERP) where {DT,TT,TGT,SGT,SM}
    if isempty(green.instant)
        error("Instant Green's function can not be empty!")
    else
        IM = InterpMethod(SGT,SM)
        nei = CompositeGrids.Interp.findneighbor(IM, green.spaceGrid,space)
        instant_x = view(green.instant, color1, color2, nei.index)
        return CompositeGrids.Interp.interpsliced(nei,instant_x)
    end
end

"""
    function dynamic(green::Green2DLR{DT,TT,TGT,SGT}, time, space, color1::Int, color2::Int, timeMethod::TM , spaceMethod::SM) where {DT,TT,TGT<:CompositeGrids.AbstractGrid,SGT<:CompositeGrids.AbstractGrid,TM,SM}

Find value of Green's function's dynamic part at given color and k/x by interpolation.
Interpolation method is by default depending on the grid, but could also be chosen to be linear.

#Argument
- 'green': Green's function
- 'time': Target τ/ω_n point
- 'space': Target k/x point
- 'color1': Target color1
- 'color2': Target color2
- 'timeMethod': Method of interpolation for time
- 'spaceMethod': Method of interpolation for space
"""
function dynamic(green::Green2DLR{DT,TT,TGT,SGT}, time, space, color1::Int, color2::Int, timeMethod::TM , spaceMethod::SM) where {DT,TT,TGT<:CompositeGrids.AbstractGrid,SGT<:CompositeGrids.AbstractGrid,TM,SM}
    # for double composite
    if isempty(green.dynamic)
        error("Dynamic Green's function can not be empty!")
    else
        TIM = InterpMethod(TGT,SM)
        SIM = InterpMethod(SGT,SM)

        spaceNeighbor = CompositeGrids.Interp.findneighbor(SIM, green.spaceGrid, space)
        timeNeighbor = CompositeGrids.Interp.findneighbor(TIM, green.timeGrid, time)
        dynamic_slice = view(green.dynamic, color1, color2, spaceNeighbor.index, timeNeighbor.index)

        dynamic_slice_xint = CompositeGrids.Interp.interpsliced(spaceNeighbor,dynamic_slice, axis=1)
        result = CompositeGrids.Interp.interpsliced(timeNeighbor,dynamic_slice_xint, axis=1)
    end
    return result
end

function dynamic(
    green::Green2DLR{DT,TT,TGT,SGT}, time, space, color1::Int, color2::Int,
    timeMethod::LinearInterp = LINEARINTERP , spaceMethod::LinearInterp = LINEARINTERP
    ) where {DT,TT,TGT<:CompositeGrids.AbstractGrid,SGT<:CompositeGrids.AbstractGrid}
    # for double composite and double linear
    if isempty(green.dynamic)
        error("Dynamic Green's function can not be empty!")
    else
        dynamic_slice = view(green.dynamic, color1, color2, :,:)
        result = CompositeGrids.Interp.linear2D(dynamic_slice, green.spaceGrid,green.timeGrid,space,time)
    end
    return result
end

function dynamic(green::Green2DLR{DT,TT,TGT,SGT}, time, space, color1::Int, color2::Int, timeMethod::DLRInterp , spaceMethod::SM) where {DT,TT,TGT<:CompositeGrids.AbstractGrid,SGT<:CompositeGrids.AbstractGrid,SM}
    # for composite space and dlr time
    if isempty(green.dynamic)
        error("Dynamic Green's function can not be empty!")
    else
        SIM = InterpMethod(SGT,SM)
        spaceNeighbor = CompositeGrids.Interp.findneighbor(SIM, green.spaceGrid, space)
        dynamic_slice = view(green.dynamic, color1, color2, spaceNeighbor.index,:)

        dynamic_slice_xint = CompositeGrids.Interp.interpsliced(spaceNeighbor,dynamic_slice, axis=1)
        if green.timeType == ImFreq
            result = (matfreq2matfreq(green.dlrGrid, dynamic_slice_xint, [time,], green.timeGrid.grid))[1]
        elseif green.timeType == ImTime
            result = (tau2tau(green.dlrGrid, dynamic_slice_xint, [time,], green.timeGrid.grid))[1]
        end
    end
    return result
end


text:
"""
MDP with all state uncertainty modeled by outcome uncertainty
"""
struct OutcomeMDP{M} <: MDP{MLPhysicalState, MLAction}
    mdp::M
end

function generate_sr(m::OutcomeMDP, s::MLPhysicalState, a::MLAction, rng::AbstractRNG)
    cars = CarState[]
    for c in s.cars
        cs = CarState(c, rand(rng, m.mdp.dmodel.behaviors))
        push!(cars, cs)
    end
    mls = MLState(s, cars)
    mlsp, r = generate_sr(m.mdp, mls, a, rng)
    return MLPhysicalState(mlsp), r
end

actions(m::OutcomeMDP, s::MLPhysicalState) = actions(m.mdp, s)
discount(m::OutcomeMDP) = discount(m.mdp)

struct OutcomeSolver <: Solver
    solver::Solver
end

struct OutcomePlanner{P} <: Policy
    planner::P
end

function solve(sol::OutcomeSolver, mdp)
    omdp = OutcomeMDP(mdp)
    return OutcomePlanner(solve(sol.solver, omdp))
end

function action_info(p::OutcomePlanner, s::MLState)
    ps = MLPhysicalState(s)
    return action_info(p.planner, ps)
end
action(p::OutcomePlanner, s::MLState) = first(action_info(p, s))

Base.srand(p::OutcomePlanner, s) = srand(p.planner, s)


text:
module TimeWindowStats

using Dates

export Sum, update!, Pct, Mean

abstract TimeWindowStatistic

type Sum<:TimeWindowStatistic
    tm::Int64
    size::Int64 # size of the array including empty cells
    len::Int64 # number of values in the calculation
    tail::Int64 # index of last value
    head::Int64 # index of last value
    tms::Vector{DateTime}
    xs::Vector
    v # value
    Sum(time::Int64, size::Int64=100) = new(time, size, 0, 0, 0, fill(DateTime(1900,1,1,0,0,0),size),[])
end

function update!(s::Sum, v, tm::DateTime)

    # initalization
    reducedCnt = 0
    reducedAmt = zero(typeof(v))
    arrayFull　= false
    nextHead = s.head+1 <= s.size? s.head+1:1

    if s.tail != 0
        i = s.tail
#         msg = tm - s.tms[i] > s.tm? "reduce":"no-reduce"

        while tm - s.tms[i] > s.tm && reducedCnt < s.len
            reducedAmt += s.xs[i]
            reducedCnt += 1
            i = i < s.size?i+1:1
#             println("  reduced $(s.len) $(int64(tm - s.tms[i-1])) reducedCnt $reducedCnt i=$(i-1) $(s.tms[i-1]) $(s.xs[i-1])")
        end
#         println("     $msg i=$i reduced $reducedAmt")
        s.tail = i

    else  # only happens during initialization
        s.tail = 1
        s.xs = zeros(typeof(v), s.size)
        s.v = zero(typeof(v))
    end

    s.len -= reducedCnt
    arrayFull = s.len + 1 >= s.size
#     println("    arrayFull $arrayFull")

    if !arrayFull
        s.xs[nextHead] = v
        s.tms[nextHead] = tm
        s.head = nextHead
    else
        newSize = iceil(s.size*1.05) # increase size by 5%
        newValues = zeros(typeof(v),newSize)
        newDates  = fill(DateTime(1900,1,1,0,0,0),newSize)
        s.size = newSize
#         println("Array full! size=$(s.size) newSize=$(newSize) head=$(s.head) tail=$(s.tail) nextHead=$(nextHead)")
        if  s.head > s.tail
            newValues[1:s.len] = s.xs[s.tail:s.head]
            newDates[1:s.len] = s.tms[s.tail:s.head]
        else
#             println("Reversed $(s.tail)  $(s.head)  $(s.len)")
            newValues[1:s.len] = vcat(s.xs[s.tail:end],s.xs[1:s.head])
            newDates[1:s.len]  = vcat(s.tms[s.tail:end],s.tms[1:s.head])
#             println("Reversed done")
        end
        s.xs = newValues
        s.tms = newDates
        s.tail = 1
        s.head = s.len+1
        s.xs[s.head] = v
        s.tms[s.head] = tm
#         println("Array full! tail=$(s.tail) head=$(s.head) nextHead=$(nextHead) newSize=$newSize")
    end

    s.v += v - reducedAmt # calculate and cache latest value
    s.len += 1
    return s
end

type Mean <:TimeWindowStatistic
    tm::Int64
    n::Sum # nominator
    v::Float64 # value
    Mean(time::Int64) = new(time, Sum(time), 0)
end

function update!(s::Mean, n, tm::DateTime)
    update!(s.n, n, tm)
    s.v = s.n.v / s.n.len
    return s
end

type Pct<:TimeWindowStatistic
    tm::Int64
    n::Sum # nominator
    d::Sum # denominator
    v::Float64 # value
    Pct(time::Int64, size::Int64=100) = new(time, Sum(time,size), Sum(time,size), 0)
end

function update!(s::Pct, n, d, tm::DateTime)
    update!(s.n, n, tm)
    update!(s.d, d, tm)
    s.v = s.d ==0?: 0: s.n.v/s.d.v
    return s
end

end # module

max_stars_repo_name: Andy-P/StreamingUtil

max_issues_repo_name: Andy-P/StreamingUtil

max_forks_repo_name: Andy-P/StreamingUtil


text:
using ClimateMachine
using ClimateMachine.Atmos
using ClimateMachine.BalanceLaws
using ClimateMachine.ConfigTypes
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.GenericCallbacks
using ClimateMachine.Mesh.Geometry
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.Orientations
using ClimateMachine.SystemSolvers
using Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates
using ClimateMachine.VTK

using CLIMAParameters
using CLIMAParameters.Planet: kappa_d
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using MPI, Logging, StaticArrays, LinearAlgebra, Printf, Dates, Test

include("isentropicvortex_setup.jl")

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

const output_vtk = false

function main()
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 4
    numlevels = integration_testing ? 4 : 1

    expected_error = Dict()
    expected_error[Float64, MRIGARKERK33aSandu, 1] = 2.3357934866477940e+01
    expected_error[Float64, MRIGARKERK33aSandu, 2] = 5.3328129440361121e+00
    expected_error[Float64, MRIGARKERK33aSandu, 3] = 1.2991232057877919e-01
    expected_error[Float64, MRIGARKERK33aSandu, 4] = 6.1056067013518876e-03
    expected_error[Float64, MRIGARKERK45aSandu, 1] = 2.3207510164213900e+01
    expected_error[Float64, MRIGARKERK45aSandu, 2] = 5.2787446598866872e+00
    expected_error[Float64, MRIGARKERK45aSandu, 3] = 1.2151170640665301e-01
    expected_error[Float64, MRIGARKERK45aSandu, 4] = 2.1001271191583956e-03

    @testset "$(@__FILE__)" begin
        for FT in (Float64,), dims in 2
            for mrigark_method in (MRIGARKERK33aSandu, MRIGARKERK45aSandu)
                @info @sprintf """Configuration
                                  ArrayType  = %s
                                  mrigark_method = %s
                                  FT     = %s
                                  dims       = %d
                                  """ ArrayType "$mrigark_method" "$FT" dims

                setup = IsentropicVortexSetup{FT}()
                errors = Vector{FT}(undef, numlevels)

                for level in 1:numlevels
                    numelems =
                        ntuple(dim -> dim == 3 ? 1 : 2^(level - 1) * 5, dims)
                    errors[level] = test_run(
                        mpicomm,
                        ArrayType,
                        polynomialorder,
                        numelems,
                        setup,
                        FT,
                        mrigark_method,
                        dims,
                        level,
                    )

                    @test errors[level] ≈
                          expected_error[FT, mrigark_method, level]
                end

                rates = @. log2(
                    first(errors[1:(numlevels - 1)]) /
                    first(errors[2:numlevels]),
                )
                numlevels > 1 && @info "Convergence rates\n" * join(
                    [
                        "rate for levels $l → $(l + 1) = $(rates[l])"
                        for l in 1:(numlevels - 1)
                    ],
                    "\n",
                )
            end
        end
    end
end

function test_run(
    mpicomm,
    ArrayType,
    polynomialorder,
    numelems,
    setup,
    FT,
    mrigark_method,
    dims,
    level,
)
    brickrange = ntuple(dims) do dim
        range(
            -setup.domain_halflength;
            length = numelems[dim] + 1,
            stop = setup.domain_halflength,
        )
    end

    topology = BrickTopology(
        mpicomm,
        brickrange;
        periodicity = ntuple(_ -> true, dims),
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = polynomialorder,
    )

    problem =
        AtmosProblem(boundaryconditions = (), init_state_prognostic = setup)

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = IsentropicVortexReferenceState{FT}(setup),
        turbulence = ConstantDynamicViscosity(FT(0)),
        moisture = DryModel(),
    )
    model = AtmosModel{FT}(
        AtmosLESConfigType,
        physics;
        problem = problem,
        orientation = NoOrientation(),
        source = (),
    )
    # The linear model has the fast time scales
    fast_model = AtmosAcousticLinearModel(model)

    dg = DGModel(
        model,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )
    fast_dg = DGModel(
        fast_model,
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient();
        state_auxiliary = dg.state_auxiliary,
    )
    slow_dg = remainder_DGModel(dg, (fast_dg,))

    timeend = FT(2 * setup.domain_halflength / setup.translation_speed)
    # determine the slow time step
    elementsize = minimum(step.(brickrange))
    slow_dt =
        2 * elementsize / soundspeed_air(param_set, setup.T∞) /
        polynomialorder^2
    nsteps = ceil(Int, timeend / slow_dt)
    slow_dt = timeend / nsteps

    # arbitrary and not needed for stability, just for testing
    fast_dt = slow_dt / 3

    Q = init_ode_state(dg, FT(0), setup)

    fastsolver = LSRK144NiegemannDiehlBusch(fast_dg, Q; dt = fast_dt)

    ode_solver = mrigark_method(slow_dg, fastsolver, Q, dt = slow_dt)

    eng0 = norm(Q)
    dims == 2 && (numelems = (numelems..., 0))
    @info @sprintf """Starting refinement level %d
                      numelems  = (%d, %d, %d)
                      slow_dt   = %.16e
                      fast_dt   = %.16e
                      norm(Q₀)  = %.16e
                      """ level numelems... slow_dt fast_dt eng0

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            runtime = Dates.format(
                convert(DateTime, now() - starttime[]),
                dateformat"HH:MM:SS",
            )
            @info @sprintf """Update
                              simtime = %.16e
                              runtime = %s
                              norm(Q) = %.16e
                              """ gettime(ode_solver) runtime energy
        end
    end
    callbacks = (cbinfo,)

    if output_vtk
        # create vtk dir
        vtkdir =
            "vtk_isentropicvortex_mrigark" *
            "_poly$(polynomialorder)_dims$(dims)_$(ArrayType)_$(FT)" *
            "_$(FastMethod)_level$(level)"
        mkpath(vtkdir)

        vtkstep = 0
        # output initial step
        do_output(mpicomm, vtkdir, vtkstep, dg, Q, Q, model)

        # setup the output callback
        outputtime = timeend
        cbvtk = EveryXSimulationSteps(floor(outputtime / slow_dt)) do
            vtkstep += 1
            Qe = init_ode_state(dg, gettime(ode_solver))
            do_output(mpicomm, vtkdir, vtkstep, dg, Q, Qe, model)
        end
        callbacks = (callbacks..., cbvtk)
    end

    solve!(Q, ode_solver; timeend = timeend, callbacks = callbacks)

    # final statistics
    Qe = init_ode_state(dg, timeend)
    engf = norm(Q)
    engfe = norm(Qe)
    errf = euclidean_distance(Q, Qe)
    @info @sprintf """Finished refinement level %d
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e
    """ level engf engf / eng0 engf - eng0 errf errf / engfe
    errf
end

function do_output(
    mpicomm,
    vtkdir,
    vtkstep,
    dg,
    Q,
    Qe,
    model,
    testname = "isentropicvortex_mrigark",
)
    ## name of the file that this MPI rank will write
    filename = @sprintf(
        "%s/%s_mpirank%04d_step%04d",
        vtkdir,
        testname,
        MPI.Comm_rank(mpicomm),
        vtkstep
    )

    statenames = flattenednames(vars_state(model, Prognostic(), eltype(Q)))
    exactnames = statenames .* "_exact"

    writevtk(filename, Q, dg, statenames, Qe, exactnames)

    ## Generate the pvtu file for these vtk files
    if MPI.Comm_rank(mpicomm) == 0
        ## name of the pvtu file
        pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

        ## name of each of the ranks vtk files
        prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
            @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
        end

        writepvtu(
            pvtuprefix,
            prefixes,
            (statenames..., exactnames...),
            eltype(Q),
        )

        @info "Done writing VTK: $pvtuprefix"
    end
end

main()


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

# Wrapper around C++ pointers used by chemfiles, adding automatic memory
# management and constness checking.
mutable struct CxxPointer{T}
    __ptr::Ptr{T}
    __is_const::Bool
    function CxxPointer(ptr::Ptr{T}; is_const::Bool=true) where T
        @assert Int(ptr) != 0
        this = new{T}(ptr, is_const)
        finalizer(o -> lib.chfl_free(Ptr{Cvoid}(__const_ptr(o))), this)
        return this
    end
end

# Get a mutable C pointer from a CxxPointer, checking that we do have mutable
# access to the object.
function __ptr(ptr::CxxPointer)
    if ptr.__is_const
        throw(ChemfilesError("This object is immutable"))
    else
        return ptr.__ptr
    end
end

# Get a non mutable C pointer from a CxxPointer.
function __const_ptr(ptr::CxxPointer)
    return ptr.__ptr
end

function __get_last_error()
    message = ""
    try
        message = last_error()
    catch error
        @warn "$error while getting chemfiles error message"
    end
    return message
end

function __check(result::Integer)
    if result != 0
        throw(ChemfilesError(__get_last_error()))
    end
end

macro __check_ptr(expr)
    quote
        local ptr = $(esc(expr))
        if Int(ptr) == 0
            message = __get_last_error()
            if occursin("internal error: pointer at", message) && occursin("is already managed by shared_allocator", message)
                @warn "Got an internal error from chemfiles::shared_allocator, running GC and retrying"
                GC.gc()
                ptr = $(esc(expr))
                if Int(ptr) == 0
                    # This is a legit failure, bail out
                    throw(ChemfilesError(__get_last_error()))
                end
            else
                throw(ChemfilesError(message))
            end
        end
        ptr
    end
end

# Remove the final NULL ('\0') character in string `s` coming from C
function __strip_null(s)
    for i in 1:length(s)
        if s[i] == '\0'
            return s[1:i-1]
        end
    end
    throw(ChemfilesError("A C string is not NULL terminated"))
end

# Call the `callback` with a growing string buffer until the return string
# fits in the buffer.
function __call_with_growing_buffer(callback::Function, initial_size=64)
    function __buffer_was_big_enough(buffer)
        @assert length(buffer) > 2
        return buffer[end-2] == '\0'
    end

    size = initial_size
    buffer = repeat("\0", size)

    callback(pointer(buffer), UInt64(size))
    while !__buffer_was_big_enough(buffer)
        # Grow the buffer and retry
        size *= 2
        buffer = repeat("\0", size)
        callback(pointer(buffer), UInt64(size))
    end

    return __strip_null(buffer)
end


text:
export parameteropt

#Optimize the hyperparameters of a GP, given a data set
#For now, all parameters are in the covariance function
function parameteropt(initp::Vector, X, y)
    sig_n2 = 0.001
    function func(params)
        #println("Func call: ", params)
        llfunc(params, sig_n2, X, y)
    end
    function gradf!(params, xgrad)
        #println("Grad call: ", params)
        llgrad!(xgrad, params, sig_n2, X, y)
    end
    function nlopt_func(params, grad)
        llfunc_and_grad!(grad, params, sig_n2, X, y)
    end

    #NLOpt solving
    optprob = NLopt.Opt(:LD_MMA, length(initp))
    NLopt.min_objective!(optprob, nlopt_func)
    NLopt.lower_bounds!(optprob, 0.01)
    NLopt.upper_bounds!(optprob, 2.0)
    NLopt.ftol_abs!(optprob, 1e-3)
    NLopt.maxeval!(optprob, 1000)

    (optf,optx,ret) = NLopt.optimize(optprob, initp)
    optx
end

function parameteropt{N}(gp::GaussianProcess{N}, X, y)
    init_params = [gp.K.sig_f; gp.K.lambdas]
    parameteropt(init_params, X, y)
end
parameteropt{N}(gp::GaussianProcess{N}) = parameteropt(gp, gp.Xtrn, gp.ftrn)

#Simultaneously calculate log-likelihood and its gradient
#In-place update xgrad and return the llfunc value
function llfunc_and_grad!(xgrad, params, sig_n2, X, y)
    K = SECov(params)
    tgp = SEGP(K, sig_n2, X, y)

    Linv = inv(tgp.L)
    Kinv = Linv'*Linv
    for i=1:length(xgrad)
        dK_i = grad(tgp.K, Val{i-1}, X)
        xgrad[i] = -0.5*(y'*Kinv*dK_i*Kinv*y)[1] +
                    0.5*trace(Kinv*dK_i)
    end
    -loglikelihood(tgp)
end

#Calc log-likelihood
function llfunc(params, sig_n2, X, y)
    K = SECov(params)
    tgp = SEGP(K, sig_n2, X, y)
    -loglikelihood(tgp)
end

#Calc log-likelihood gradient
function llgrad!(xgrad, params, sig_n2, X, y)
    #println("grad params: ", params)
    Kt = SECov(params)
    tgp = SEGP(Kt, sig_n2, X, y)
    Linv = inv(tgp.L)
    Kinv = Linv'*Linv

    Np = length(params)
    for i=1:Np
        dK_i = grad(tgp.K, Val{i-1}, X)
        xgrad[i] = -0.5*(y'*Kinv*dK_i*Kinv*y)[1] +
                    0.5*trace(Kinv*dK_i)
    end
    Ngrad = norm(xgrad)
    for i=1:Np
        xgrad[i] = xgrad[i] / Ngrad
    end
    xgrad
end

#Debugging function
function checkgrad{N}(params, gp::GaussianProcess{N}, X, y)
    Np = length(params)
    dx = 1e-6
    out = zeros(Np,2)

    f0 = llfunc(params, gp.sig_n2, X, y)
    for i=1:Np
        pdiff = copy(params)
        pdiff[i] = params[i] + dx
        fpi = llfunc(pdiff, gp.sig_n2, X, y)
        pdiff[i] = params[i] - dx
        fni = llfunc(pdiff, gp.sig_n2, X, y)
        out[i,:] = [(fpi-f0) (f0-fni)]/dx
    end

    xgrad = zeros(Np)
    llgrad!(xgrad, params, gp.sig_n2, X, y)
    hcat(xgrad, out)
end


text:
using Revise
using JLD2
using Test
using DelimitedFiles

@load "data/C.jld2" C VX
CC = C
VXX = VX
@load "data/C0.jld2" C VX
C0 = C
VX0 = VX
@load "data/C1.jld2" C VX
C1 = C
VX1 = VX
@load "data/C2.jld2" C VX
C2 = C
VX2 = VX
@load "data/C3.jld2" C VX
C3 = C
VX3 = VX
CMPI = [C0 C1;C2 C3]
VXMPI = [VX0 VX1; VX2 VX3]

writedlm(stdout, CC)
println("-------------------------------------------------------------")
writedlm(stdout, CMPI)
println("-------------------------------------------------------------")
writedlm(stdout, CMPI-CC)
println("-------------------------- Relative ----------------------------------")
writedlm(stdout, abs.(CMPI-CC) ./ (1e-10 .+abs.(CC)))


text:
module util
using Distributions, RCall
import Cytof5

# Import R libraries
# TODO: Remake rcommon and cytof3 for this.
R"require(rcommon)";
R"require(cytof3)";

# Import R plotting functions
plot = R"plot";
boxplot = R"boxplot";
ari = R"cytof3::ari";
rgba = R"rcommon::rgba";
density = R"density";
lines = R"lines";
plotPost = R"rcommon::plotPost";
plotPosts = R"rcommon::plotPosts";
myImage = R"cytof3::my.image";
plotPdf = R"pdf";
# plotPng = R"png";
R"""
plotPng = function(fname, s=10, w=480, h=480, ps=12, ...) {
  png(fname, w=w*s, h=h*s, pointsize=ps*s, ...)
}
""";
plotPng = R"plotPng";
devOff = R"dev.off";
blueToRed = R"cytof3::blueToRed";
greys = R"cytof3::greys";
plot_dat = R"cytof3::plot_dat";
yZ_inspect = R"cytof3::yZ_inspect";
yZ= R"cytof3::yZ";
abline = R"abline";
addErrbar = R"rcommon::add.errbar";
hist = R"hist";
colorBtwn = R"color.btwn";
myQQ = R"my.qqplot";
# addCut = R"add.cut";
addCut = R"function(clus, s=1) abline(h=cumsum(table(clus)) + .5, lwd=3*s, col='yellow')";
addGridLines = R"function(Z, s=1) abline(v=1:NCOL(Z) + .5, h=1:NROW(Z) + .5, col='grey', lwd=1*s)";

function quantile_vm(xs, p)
  I, J = size(xs[1])
  for x in xs
    @assert size(x) == (I, J)
  end
  return [quantile([x[i, j] for x in xs], p) for i in 1:I, j in 1:J]
end

function getPosterior(sym::Symbol, monitor)
  return [ m[sym] for m in monitor]
end

nrow(X) = size(X, 1)
ncol(X) = size(X, 2)

function pairwiseAlloc(Z, W, i)
  J, K = size(Z)
  A = zeros(J, J)
  for r in 1:J
    for c in 1:J
      for k in 1:K
        if Z[r, k] == 1 && Z[c, k] == 1
          A[r, c] += W[i, k]
          A[c, r] += W[i, k]
        end
      end
    end
  end

  return A
end # pairwiseAlloc

function estimate_ZWi_index(monitor, i)
  As = [pairwiseAlloc(m[:Z], m[:W], i) for m in monitor]

  Amean = mean(As)
  mse = [ mean((A - Amean) .^ 2) for A in As]

  return argmin(mse)
end

function plotProbMiss(beta, i; xlim=[-10, 5],
                      ygrid=range(xlim[1], stop=xlim[2], length=300))

  p = [Cytof5.Model.prob_miss(yy, beta[:, i]) for yy in ygrid]
  plot(ygrid, p, main="i: $i", typ="l", xlab="y", ylab="prob of missing",
       xlim=xlim)
end

# For Posterior Predictive QQ
function y_obs_range(y)
  I = length(y)
  y_vec = vcat([vec(y[i]) for i in 1:I]...)
  y_obs = filter(yy -> !isnan(yy), y_vec)
  return quantile(y_obs, [0, 1])
end

function idx_observed_ij(y, i::Int, j::Int)
  return findall(yij -> !isnan(yij), y[i][:, j])
end

function get_mu_observed(i::Int, j::Int, y, o, o2)
  idx_observed = idx_observed_ij(y, i, j)
  gam_ij = o2[rand(1:end)][:gam][i][:, j]

  result = [begin
              k = o[:lam][i][n]
              if k == 0 # noisy class
                0.0
              else
                z = o[:Z][j, k]
                l = gam_ij[n]
                sum_delta = sum(o[:delta][z][1:l])

                z == 0 ? -sum_delta : sum_delta
              end
            end for n in idx_observed]
  return result
end

function gen_post_pred(i, j, y, out)
  y_pp = [rand.(Normal.(get_mu_observed(i, j, y, o, out[2]), sqrt(o[:sig2][i])))
          for o in out[1]]
  return hcat(y_pp...)
end

function qq_yobs_postpred(y, i::Int, j::Int, out)

  y_obs = y[i][idx_observed_ij(y, i, j), j]
  y_pp = gen_post_pred(i, j, y, out)

  return y_obs, y_pp
end

function reorder_lami(ord, lami)
  lami_new = fill(NaN, length(lami))
  K = length(ord)
  for k in 1:K
    lami_new[lami .== ord[k]] .= k
  end
  return lami_new
end

function get_common_celltypes(Wi; thresh=.9, filter_by_min_presence=false)
  ord = sortperm(Wi, rev=true)
  cumsum_Wi = cumsum(Wi[ord])

  if filter_by_min_presence
    # filter celltypes by min threshold
    return ord[Wi[ord] .> thresh]
  else
    # filter celltypes by a cumulative presence
    k = minimum(findall(cumsum_Wi .> thresh))
    return ord[1:k]
  end
end

end # util


text:
# This file will contain the SOCP TrajOpt Module
"""
    TrajOptSOCPs

A native Julia library to solve trajectory optimization problems that
contain second-order cone constraints.

Author: Daniel Neamati (Summer 2020)

Advisor: Prof. Zachary Manchester (REx Lab at Stanford University)

Funding graciously provided by Caltech SURF program and the Homer J. Stewart
Fellowship.
"""
module TrajOptSOCPs

using LinearAlgebra, SparseArrays

# Declare Abstract Types
abstract type constraint end
abstract type objectiveFunc end
abstract type constraintManager end

# Export the abstract types
export constraint, objectiveFunc, constraintManager

# Basic Test
include("extra.jl")
export fCheck2

# Include all of the files
include("rocket/rocket-setup.jl")
include("rocket/ground.jl")
include("rocket/maxThrust.jl")
include("rocket/maxThrustAngle.jl")

include("dynamics/trajectory-setup.jl")
include("objective/LQR_objective.jl")
include("constraints/constraintManager.jl")
include("constraints/constraintManagerDynamics.jl")
include("auglag/auglag-core.jl")

include("solver/AL-Primal-Main-Solver.jl")
include("solver/QP-AffineEquality-Solver.jl")

include("other_utils/parsePrimalDual.jl")
include("other_utils/trajectoryParsing.jl")
include("other_utils/montecarlo.jl")

# # Export Constraint Structs
export AL_AffineEquality, AL_AffineInequality,
        AL_simpleCone, AL_Multiple_simpleCone,
        AL_simpleAngleCone, AL_Multiple_AngleCone

# Export Top Level Structs
export augLag, solverParams, constraintManager_Dynamics,
        LQR_QP_Referenced, rocket_simple

# Export other Helper structs
export parseTrajectory, primal_dual

# Export Top Level Functions
export ALPrimalNewtonMain, solParamPrint,
        initializeTraj, initializeHoverOnlyTraj,
        rocketDynamicsFull

# Export Helper Set-Up Functions
export makeLQR_TrajReferenced, makeLQR_TrajRefTerminal,
        makeMaxThrustConstraint, makeGroundConstraint, makeMaxAngleConstraint

# Export Evaluation functions
export evalConstraints, evalGradConstraints, evalHessConstraints,
        evalAL, evalGradAL, evalHessAl, evalAffineEq

# Export Helper Utils
export parsePrimalDualVec, getParseTrajectory, splitDimensions,
        getConstraintViolationList, safeNorm

end  # module SOCP_TrajOpt

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
using MPIMeasurements
using PyPlot

include("config.jl")

scanner = MPIScanner(conf)
daq = getDAQ(scanner)

params = toDict(daq.params)

params["studyName"] = "Test"
params["studyDescription"] = "A very cool measurement"
params["scannerOperator"] = "Tester"
params["dfStrength"] = [1e-3]
params["dfWaveform"] = "SINE"
params["acqNumFrames"] = 10
params["acqNumAverages"] = 10

u, uSlowADC = measurement(daq, params, controlPhase=true)
#filename = measurement(daq, params, MDFStore, controlPhase=true)

figure(1)
clf()
plot(vec(u))
savefig("images/im1.png")

scanner = MPIScanner(conf) #what???
daq = getDAQ(scanner)
params["dfWaveform"] = "TRIANGLE"
u, uSlowADC = measurement(daq, params, controlPhase=false)
figure(1)
clf()
plot(vec(u))
savefig("images/im2.png")


text:
import FlatBuffers

include("monster_generated.jl")

# require the files generated from the schema
import .MyGame.Sample.Weapon
import .MyGame.Sample.Monster
import .MyGame.Sample.Vec3
import .MyGame.Sample.Color
import .MyGame.Sample.Equipment

sword = Weapon(;name="Sword", damage=3)
axe = Weapon(;name="Axe", damage=5)
orc = Monster(;
    name="Orc",
    pos=Vec3(1.0, 2.0, 3.0),
    hp=300,
    inventory=collect(1:10),
    color=MyGame.Sample.ColorRed,
    weapons=[sword, axe],
    equipped=axe
)

# Get the flatbuffer as a string containing the binary data
io = IOBuffer()
FlatBuffers.serialize(io, orc)
bytes = take!(io)

# Build the Monster from the raw bytes
mon = Monster(bytes)

@assert mon.mana == 150
@assert mon.hp == 300
@assert mon.name == "Orc"
@assert mon.color == MyGame.Sample.ColorRed
@assert mon.pos.x == 1.0
@assert mon.pos.y == 2.0
@assert mon.pos.z == 3.0
@assert mon.inventory == collect(1:10)

@assert mon.equipped_type == FlatBuffers.typeorder(Equipment, Weapon)
@assert mon.equipped.name == "Axe"
@assert mon.equipped.damage == 5

@assert mon.weapons[1].name == "Sword"
@assert mon.weapons[1].damage == 3

@assert mon.weapons[2].name == "Axe"
@assert mon.weapons[2].damage == 5

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia