text:
#=
This is where OpCode definitions and some methods for OpCodes reside.

Originally, opcodes and source/destination identifiers were to be Symbols,
since Symbols in Julia are interned and thus offer themselves to quick
comparison. However, I realized that "Symbols" are not actually "bits" types,
rather they are pointers to the actual symbol data that lives on the heap.

Thus, immutable structs containing Symbols contain a reference and thus seem
to be heap allocated. This is true even for small structs. I want to avoid this
because if a bunch of simulated cores are allocating a lot of memory every
cycle, this will lead to longer GC times.

DISCLAIMER: I did not benchmark an entire system with the Symbol based
implementation because that was not functional at the time. Rather, I switched
over to the current implementation before the whole system was functional.

Also, this may be fixed in Julia 0.7/1.0

--------------------------------------------------------------------------------

The current implementation still uses some of that symbol based ideas, but uses
Enums to encode opcodes and source/destination identifiers.

These enum types are automatically generated from come base constants, so adding
op codes should be easy and automatically propogate to the created Enums.

One downside is that referencing Enums now required prefixing everything with
AsapSim if used outside of this module.

That's generally not a problem because the macro based assembler can take care
of those non-scoped references.

=#

################################################################################
# Instruction Definitions
################################################################################

# This is a collection of the instructions for Asap4.

# Instructions with three operands:
#   inst dst src1 src2 [options]

@enum OpType::Int8 NOP_TYPE ALU_TYPE MAC_TYPE BRANCH_TYPE OTHER_TYPE PSUEDO_TYPE

@with_kw_noshow struct InstructionDefinition
    opcode :: Symbol
    # True if this instruction has these fields.
    dest :: Bool = true
    src1 :: Bool = true
    src2 :: Bool = true

    # Classification of operand type
    optype :: OpType
    signed :: Bool = false
    saturate :: Bool = false
end

# Convenience constructor so we don't have to write "opcode" and "optype"
# all the time.
function InstructionDefinition(opcode::Symbol, optype::OpType; kwargs...)
    return InstructionDefinition(;
        opcode = opcode,
        optype = optype,
        kwargs...
    )
end

const InstD = InstructionDefinition

# Define all of the Instructions
const AsapInstructionList = (
    # --- Instructions using all 3 operands --- #

    # Unsigned Arithmetic
    InstD(:ADDU,   ALU_TYPE),
    InstD(:ADDSU,  ALU_TYPE, saturate = true),
    InstD(:ADDCU,  ALU_TYPE),
    InstD(:ADDCSU, ALU_TYPE, saturate = true),
    InstD(:SUBU,   ALU_TYPE),
    InstD(:SUBSU,  ALU_TYPE, saturate = true),
    InstD(:SUBCU,  ALU_TYPE),
    InstD(:SUBCSU, ALU_TYPE, saturate = true),

    # Signed Arithmetic
    InstD(:ADD,    ALU_TYPE, signed = true),
    InstD(:ADDS,   ALU_TYPE, signed = true, saturate = true),
    InstD(:ADDC,   ALU_TYPE, signed = true),
    InstD(:ADDCS,  ALU_TYPE, signed = true, saturate = true),
    InstD(:SUB,    ALU_TYPE, signed = true),
    InstD(:SUBS,   ALU_TYPE, signed = true, saturate = true),
    InstD(:SUBC,   ALU_TYPE, signed = true),
    InstD(:SUBCS,  ALU_TYPE, signed = true, saturate = true),

    # Logic Operatins
    InstD(:OR,  ALU_TYPE),
    InstD(:AND, ALU_TYPE),
    InstD(:XOR, ALU_TYPE),

    # 3 operand conditional moves
    InstD(:MOVC, ALU_TYPE),
    InstD(:MOVZ, ALU_TYPE),
    InstD(:MOVCX0, ALU_TYPE),
    InstD(:MOVCX1, ALU_TYPE),

    # Shifts
    InstD(:SHL,  ALU_TYPE),
    InstD(:SHR,  ALU_TYPE),
    InstD(:SRA,  ALU_TYPE),
    InstD(:SHLC, ALU_TYPE),
    InstD(:SHRC, ALU_TYPE),
    InstD(:SRAC, ALU_TYPE),

    # Multiplication
    InstD(:MULTL,  MAC_TYPE, signed = true),
    InstD(:MULTH,  MAC_TYPE, signed = true),
    InstD(:MULTLU, MAC_TYPE),
    InstD(:MULTHU, MAC_TYPE),

    # MACC
    InstD(:MACL,   MAC_TYPE, signed = true),
    InstD(:MACH,   MAC_TYPE, signed = true),
    InstD(:MACCL,  MAC_TYPE, signed = true),
    InstD(:MACCH,  MAC_TYPE, signed = true),
    InstD(:MACLU,  MAC_TYPE),
    InstD(:MACHU,  MAC_TYPE),
    InstD(:MACCLU, MAC_TYPE),
    InstD(:MACCHU, MAC_TYPE),

    # --- Operatins without a SRC2 ---

    # Unconditional moves
    InstD(:MOVE,    ALU_TYPE, src2 = false),
    InstD(:MOVI,    ALU_TYPE, src2 = false),
    # Unary Reductions
    InstD(:ANDWORD, ALU_TYPE, src2 = false),
    InstD(:ORWORD,  ALU_TYPE, src2 = false),
    InstD(:XORWORD, ALU_TYPE, src2 = false),
    # Misc Logic
    InstD(:BTRV,    ALU_TYPE, src2 = false),
    InstD(:LSD,     ALU_TYPE, src2 = false),
    InstD(:LSDU,    ALU_TYPE, src2 = false),
    # Accumulate
    InstD(:ACCSH,   MAC_TYPE, src2 = false, signed = true),
    InstD(:ACCSHU,  MAC_TYPE, src2 = false),

    # --- Instructions without a SRC2 or DEST
    InstD(:RPT,   OTHER_TYPE, src2 = false, dest = false),
    InstD(:BR,    BRANCH_TYPE, src2 = false, dest = false),
    InstD(:BRL,   BRANCH_TYPE, src2 = false, dest = false),
    InstD(:STALL, OTHER_TYPE, src2 = false, dest = false),

    # --- Instructions with no arguments ---
    InstD(:NOP, NOP_TYPE, src1 = false, src2 = false, dest = false),

    # --- Pseudo instructions ---
    # Appear in input assembly but not in output program. Need to declare it
    # as an instruction though so it gets processed during source
    # transformations
    InstD(:END_RPT, PSUEDO_TYPE, src1 = false, src2 = false, dest = false),
)

# Make an Enum containing all the OpCodes.
#
# This macro will gather all the opcodes defined above and create an Enum with
# names corresponding to each opcode.
#
# After the invocation of this macro, key words like ADD, STALL, BR will be this
# enum type.
macro MakeOpcodeEnum()
    # Get all of the opcodes defined above, quote them into a vector of
    # Expr.
    enums = [:($(i.opcode)) for i in AsapInstructionList]

    # Splat the vector into a enum macro invokation. When @MakeOpcodeEnum is
    # called, it will return the below expression, which will in-turn get
    # expanded to form the actual Enum.
    #
    # Yeah, this is a bit strange, but it avoids having to repeat all of the
    # Opcodes.
    return esc(:(@enum AsapOpcode $(enums...)))
end
@MakeOpcodeEnum

# Create a dictionary mapping opcodes to their instruction definition
const Instruction_Dict = Dict(i.opcode => i for i in AsapInstructionList)

# Create ways for both the Symbol version and Enum version to access the
# Instruction_Dict
isop(x::Symbol) = haskey(Instruction_Dict, x)
isop(x::AsapOpcode) = haskey(Instruction_Dict, Symbol(x))
isop(x) = false

getdef(x::Symbol) = Instruction_Dict[x]
getdef(x::AsapOpcode) = Instruction_Dict[Symbol(x)]
# no generic fallback.

################################################################################
# Sources / Destinations
################################################################################
#=
All instructions need to encode where the expect to get their operands from
and where they expect to write their results to.

Make a few constant global structs defining all of the valid sources and
destinations for Asap4. These will also be turned into an Enum like we did
above with the Opcodes.
=#

# All the symbols for expected sources in the Asap4 architecture.
const SourceSymbols = (
    :DMEM,      # data memory (index needed)
    :IMMEDIATE, # immediate
    :IBUF,      # input fifo (index needed)
    :IBUF_NEXT, # input fifo, no increment (index needed)
    :BYPASS,    # bypass resigter (index needed)
    :POINTER,   # dereference hardware pointer (index needed)
    :AG,        # dereference address generator, no increment (index needed)
    :AG_PI,     # dereference address generator, increment (index needed)
    :POINTER_BYPASS, # dereference bypass register 1.
    :ACC,       # read lower 16 bits from accumulator
    :RET,       # read return address

    # --- Unimplemented --- #
    # :external_mem
    # :external_mem_next
    # :packet_router
    # :packet_router_next
)

#-------------------------------------------------------------------------------
# Destination encoding

const DestinationSymbols = (
    :DMEM,      # data memory (index needed)
    # TODO: For the dcmem, will probably be helpful to put in aliases for
    # common operations to help with the readability of assembly.
    :DCMEM,     # dynamic configuration memory (index needed)
    :POINTER,   # set hardware pointers (index needed)
    :AG,        # dereference address generator, no increment (index needed)
    :AG_PI,     # dereference address generator, increment (index needed)
    :POINTER_BYPASS, # dereference bypass register 1
    :NULL,       # No write back
    :OUTPUT,     # Write to OBUF, directions starting at 0 are:
        # east, north, west, south, right, up, left, down. (index needed)
    :OBUF,       # Write to all output broadcast directions.
    :OBUF_MASK,
    :SET_POINTER,
    :AG_START,
    :AG_STRIDE,
    :CXMASK,

    # Return from a branch.
    :BACK,

    # --- Unimplemented ---#
    # :external_mem
    # :packet_router
    # :dynamic_net
    # :done_flag
    # :start_pulse_test
)

const SourceDestSymbols = union(SourceSymbols, DestinationSymbols)

# Do the same trick as we did with Opcodes to create an ENUM for source and
# Destinations.
macro MakeSrcDestEnum()
    vals = [:($i) for i in SourceDestSymbols]
    return esc(:(@enum SrcDest $(vals...)))
end
@MakeSrcDestEnum

# Use this "Loc" type to store instruction source/dest information. It will
# contain a Symbol to identify the name of the source/dest, and an "index" field
# that will be either a true index for accesses to places like dmem or output,
# Or store other information such as immediate values.
struct Loc
    sym::SrcDest
    ind::Int32
end

Loc() = Loc(NULL, -1)
Loc(x::Int) = Loc(IMMEDIATE, x)
Loc(x::SrcDest) = Loc(x, -1)
Loc(x::Symbol) = SourceAliases[x]
Loc(x::Loc) = x

Base.:(==)(a::Loc, b::Loc) = (a.sym == b.sym) && (a.ind == b.ind)

Base.convert(::Type{Loc}, x::Union{Int,Symbol,SrcDest}) = Loc(x)

sym(x::Loc) = x.sym
ind(x::Loc) = x.ind

Base.show(io::IO, x::Loc) = print(io, "$(sym(x))[$(ind(x))]")

# Commonly used symbols for sources/destinations can be put here.
const SourceAliases = Dict(
    # Inputs
    :ibuf0 => Loc(IBUF, 0),
    :ibuf1 => Loc(IBUF, 1),
    :ibuf0_next => Loc(IBUF_NEXT, 0),
    :ibuf1_next => Loc(IBUF_NEXT, 1),
    # Address Generators
    :ag0pi => Loc(AG_PI, 0),
    :ag1pi => Loc(AG_PI, 1),
    :ag2pi => Loc(AG_PI, 2),
    :ag0   => Loc(AG, 0),
    :ag1   => Loc(AG, 1),
    :ag2   => Loc(AG, 2),
)

# Mark symbols as destination outputs so the pipeline avoids stalling these
# signals as they write.
const DestinationOutputs = (
    OUTPUT,
    OBUF,
)

function oneofin(a, b)
    for i in a
        i in b && return true
    end
    return false
end

# Enums for instruction type and conditional execution
@enum CXFlag::Int8 NO_CX CX_SET CX_TRUE CX_FALSE

function getoptions(args)
    # Construct an empty container for kwargs for the options.
    kwargs = Pair{Symbol,Any}[]
    # Check for nops. Iterate in the order of number of NOPS for style points.
    nops = zero(Int8)
    for (i, nop) in enumerate((:nop1, :nop2, :nop3))
        if nop in args
            nops = Int8(i)
        end
    end
    # Save nop argument if not zero
    if !iszero(nops)
        push!(kwargs, (:nops => nops))
    end
    # Check fo jump
    if :j in args
        push!(kwargs, (:jump => true))
    end

    # Check for conditional execution
    csx = oneofin((:cxs0, :cxs1, :cxt0, :cxt1, :cxf0, :cxf1), args)
    if csx
        # Set flag type
        if oneofin((:cxs0, :cxs1), args)
            push!(kwargs, (:cxflag => CX_SET))
        elseif oneofin((:cxt0, :cxt1), args)
            push!(kwargs, (:cxflag => CX_TRUE))
        else
            push!(kwargs, (:cxflag => CX_FALSE))
        end

        # Figure out the index of the conditional execution flag.
        if oneofin((:cxs0, :cxt0, :cxf0), args)
            cx_index = UInt8(0)
        else
            cx_index = UInt8(1)
        end
        push!(kwargs, (:cxindex => cx_index))
    end
    # Check double write
    if :dw in args
        push!(kwargs, (:sw => false))
    end

    return kwargs
end

# NOTE: There seems to be a problem with @with_kw where the type of the returned
# object cannot be inferred called with no arguments.
struct AsapInstruction
    # The opcode of the instruction. Set the defaultes for the "with_kw"
    # constructor to be a "nop" if constructed with no arguments.
    op :: AsapOpcode

    # Source and destination values.
    dest :: Loc
    src1 :: Loc
    src2 :: Loc

    # ------------------------ #
    # Options for instructions #
    # ------------------------ #

    # Number of no-ops to put after this instruction.
    nops :: Int8
    # Set the "jump" flag for branches
    jump :: Bool
    # Mark branch instructions as a return instruction.
    isreturn :: Bool
    # conditional execution
    cxflag :: CXFlag
    # Index of the conditional register to use.
    cxindex :: Int8
    # Options for writeback to dmem. If "true", single-write will happen
    # Default to this as it should be the common case and to bring it into
    # alignment with the c++ simulator
    sw :: Bool

    # Metadata for faster decoding during Stage 4 arithmetic or for performing
    # stall detection
    signed      :: Bool
    saturate    :: Bool
    # For fast differentiation between ALU and MAC instructions
    optype          :: OpType
    dest_is_output  :: Bool
end

AsapInstruction() = AsapInstruction(
    NOP,        # op
    Loc(),      # dest
    Loc(),      # src1
    Loc(),      # src2
    zero(Int8), # nops
    false,      # jump
    false,      # isreturn
    NO_CX,      # cxflag
    zero(Int8), # cxindex
    true,       # sw
    false,      # signed
    false,      # saturate
    NOP_TYPE,   # optype
    false,      # dest_is_output
)

# Alias "NOP" to an empty constructor, which should provide a NOP by default.
const InstNOP = AsapInstruction

function AsapInstructionKeyword(;
        op = NOP,
        dest = Loc(),
        src1 = Loc(),
        src2 = Loc(),
        nops = zero(Int8),
        jump = false,
        isreturn = false,
        cxflag = NO_CX,
        cxindex = zero(Int8),
        sw = true,
        # TODO: Move this data into this constructor.
        signed = false,
        saturate = false,
        optype = NOP_TYPE,
        #dest_is_output = false
    )

    if dest isa SrcDest
        dest_is_output = dest in DestinationOutputs
    elseif dest isa Loc
        dest_is_output = sym(dest) in DestinationOutputs
    else
        dest_is_output = false
    end

    return AsapInstruction(
        op,
        dest, src1, src2,
        nops, jump, isreturn,
        cxflag, cxindex,
        sw, signed, saturate,
        optype,
        dest_is_output,
    )
end

abstract type Mutator{T} end

struct SrcDestCollection <: Mutator{AsapInstruction}
    dest        :: Loc
    src1        :: Loc
    src2        :: Loc
end

struct InstSrc1 <: Mutator{AsapInstruction}
    src1 :: Loc
end
struct InstSrc2 <: Mutator{AsapInstruction}
    src2 :: Loc
end
struct InstDest <: Mutator{AsapInstruction}
    dest :: Loc
end

# Here's a method to change the src's and dest of an instruction quickly
# (using only the positional constructor rather than one generated by @with_kw)
#
# This is helpful for Stage 2 where pointers are dereferenced and turned into
# memory locations.
#
# Make this a @generated function so it keeps working even if the other fields
# of the AsapInstruction change.
@generated function set(a::T, b::U) where U <: Mutator{T} where T
    # Get the fieldnames for the two arguments.
    a_fields = fieldnames(T)
    b_fields = fieldnames(U)

    # Build up a list of arguments for a positional constructor of AsapInstruction
    # by iterating ove inst_fields. If a given field is in the SrcDestCollection,
    # use that instead.
    args = map(a_fields) do f
        i = findfirst(isequal(f), b_fields)
        # If "i" is in the operands fields, return that field.
        return i == nothing ? :(a.$f) : :(b.$f)
    end

    return :(T($(args...)))
end

@generated function set!(a::T, b::U) where U <: Mutator{T} where T
    # Get the fields of "b" that are also fields of "a"
    a_fields = fieldnames(T)
    b_fields = fieldnames(U)

    mutated_fields = [f for f in b_fields if f in a_fields]

    args = [:(a.$f = b.$f;) for f in mutated_fields]
    return :($(args...,))
end

# Methods on instructions - Since some values may be stored in some
# less-than-obvious places, use these methods to get these fields if needed.

# Get branch-target - This is encoded in the "dest_val" field of the branch
# instructions.
function branch_target(i::AsapInstruction)
    # For debugging, ensure that this is only called on branch instructions.
    # If not, something ahs gone wrong.
    @assert i.op == BR || i.op == BRL
    return ind(i.dest)
end

# For repeat start, use src2_index. src1_index may be used if the RPT instruction
# uses a dmem or other source
repeat_start(i::AsapInstruction) = ind(i.src2)

# Similarlym, use the dest_index for the repeat end.
repeat_end(i::AsapInstruction) = ind(i.dest)

set_branch_target(t,x) = InstDest(Loc(t,x))
set_repeat_start(x) = InstSrc2(Loc(x))
set_repeat_end(x) = InstDest(Loc(x))


text:
import Pkg; Pkg.activate(@__DIR__);
Pkg.instantiate()

using RigidBodyDynamics
using StaticArrays
using RobotDynamics
using Altro
using TrajectoryOptimization
using Parameters
using Rotations
using LinearAlgebra
using ForwardDiff

# export
#     Arthur,
#     rk4

# export
#     MPC_Params,
#     ArthurProblem,
#     ArthurHorizonProblem
    # mpc_update,
    # update_guess_trajectory!,
    # update_trajectory_fill!,
    # update_constraints!

include("../src/Arthur.jl")
include("../src/MPCUtil.jl")

ext: jl


text:
module MLDataUtils

using StatsBase
using LearnBase
using MLLabelUtils
using MLDataPattern
using DataFrames

using LearnBase: ObsDimension
import LearnBase: nobs, getobs, getobs!, datasubset, default_obsdim

using Statistics
using DelimitedFiles

eval(Expr(:toplevel, Expr(:export, setdiff(names(MLLabelUtils), [:MLLabelUtils])...)))
eval(Expr(:toplevel, Expr(:export, setdiff(names(MLDataPattern), [:MLDataPattern])...)))

export

    noisy_function,
    noisy_sin,
    noisy_poly,
    noisy_spiral,

    center!,
    rescale!,

    FeatureNormalizer,

    expand_poly,

    load_iris,
    load_line,
    load_sin,
    load_poly,
    load_spiral,

    fit,
    predict,
    predict!

include("feature_scaling.jl")
include("basis_expansion.jl")
include("noisy_function.jl")
include("datasets.jl")
include("datapattern.jl")

end


text:
module Sexagesimal

export 干支
include("干支.jl")

end # module


text:
using Test

using DataDeps: env_bool
@testset "Env Bool" begin
    withenv("A"=>"True") do
        @test env_bool("A")
    end
    withenv("B"=>"No") do
        @test !env_bool("B")
    end

    withenv("C"=>"") do
        @test !env_bool("C")
    end

    withenv("D"=>"0") do
        @test !env_bool("D")
    end

    withenv("E"=>"1") do
        @test env_bool("E")
    end
end


text:
function stringify(ssa::Array)
    ssa = [line for line in ssa if !isnothing(line)]
    return join(spacedjoin.(ssa), "\n")
end

"""
    inlinessarefs(ssa::Array)

Copypaste ssa refs into place and delete used ssa refs.
"""
function inlinessarefs(ssa::Array)
    usedrefs = []
    usedrefs!(usedrefs, ssa)

    # ssa = addparens.(ssa)
    for (i, line) in enumerate(ssa)
        ssa[i] = replacerefs(ssa, line)
    end

    for i in usedrefs
        ssa[i] = [nothing]
    end
    return ssa
end

usedrefs!(refs::Array, item) = isa(item, SSAValue) ? push!(refs, item.id) : nothing
usedrefs!(refs::Array, item::Array) = usedrefs!.((refs,), item)

replacerefs(ssa::Array, item) = isa(item, SSAValue) ? ssa[item.id] : item
replacerefs(ssa::Array, item::Array) = replacerefs.((ssa,), item)

addparens(item) = item
function addparens(items::Array)
    if isnothing(items) || length(items) == 1
        return items
    else
        return ["("; addparens.(items); ")"]
    end
end

spacedjoin(item) = item
spacedjoin(item::Array) = join(spacedjoin.(item), " ")


text:
using Statistics
using Random
using Distributions
using StaticArrays
using JSON
using POMDPs
using MCTS
using Logging
using POMDPs
using JLD2
using SCoBA
using SCoBA.SCoBADomains
using SCoBA.SCoBASolver

global_logger(SimpleLogger(stderr, Logging.Error))

const rng = MersenneTwister(5678)
const mcts_params = (disc=0.02, trials=100, explore=0.1, depth=20)
const qlearning_params = (disc=0.05, max_boxes=3)

function main()

    # Get command line input from script
    parsed_args = parse_conveyor_commandline()

    # Set up bookkeeping datastructures for stats
    trials = parsed_args["trials"]
    results_dict = Dict()
    results_dict["trials"] = trials
    lost_boxes = Int64[]
    deposited_boxes = Int64[]

    # Print the filename for tracking from the high-level bash script
    println(parsed_args["out_file_name"])

    # Iterate over different randomized trials
    for i = 1:trials

        println("Trial $(i)")
        grasp_probs = [parse(Float64, p) for p in split(parsed_args["grasp_probs"], ",")]
        server, cbelt = setup_conveyor_scenario(grasp_probs=grasp_probs,
                                            belt_speed=parsed_args["belt_speed"],
                                            new_box_prob=parsed_args["new_box_prob"])

        # Check which method was specified and benchmark accordingly
        if parsed_args["baseline"] == "mcts"

            # Create full MMDP version of problem
            mdp = create_convbelt_mdp(cbelt, server, mcts_params.disc, convert(Int64, parsed_args["timesteps"]))

            # Heuristic EDD Rollout Policy for MCTS Solver
            eddpol = EDDPolicy(mdp)
            mcts = MCTSSolver(n_iterations=mcts_params.trials, depth=mcts_params.depth,
                              rng=rng, exploration_constant=mcts_params.explore,
                              reuse_tree=false,
                              estimate_value = RolloutEstimator(eddpol))

            mcts_policy = POMDPs.solve(mcts, mdp)

            # Simulate for timesteps and collect stats
            for t = 1:parsed_args["timesteps"]

                # Look at current belt and update time window information online
                update_time_windows!(mdp.cbelt, mdp.server)

                # Current (s,a) for MCTS
                curr_state = get_current_belt_mdp_state(mdp, rng)
                curr_action = POMDPs.action(mcts_policy, curr_state)

                # Execute action on problem
                update_server_with_mdp_action!(mdp, curr_action)

                # Simulate belt and arms for one step
                propagate_belt!(mdp.cbelt, mdp.server)
                propagate_arms!(mdp.server, mdp.cbelt, rng)
                mirror_gen_step!(mdp.cbelt, rng)
            end # t = 1:timesteps
            server = mdp.server
            cbelt = mdp.cbelt

        elseif parsed_args["baseline"] == "qlearning"

            # Load up QLearning policy trained offline
            # See scripts/qlearning_conveyor_belt.jl for how it was trained
            @load parsed_args["qpol_name"] policy

            # Set up single arm MDP for Independent QLearning (IQL) to simulate on
            mdp = create_convbelt_singlearm_mdp(cbelt, server,
                                        qlearning_params.disc, qlearning_params.max_boxes)

            for t = 1:parsed_args["timesteps"]

                update_time_windows!(mdp.cbelt, mdp.server)
                vec_action = Int64[]

                # IQL; hence get individual actions separately and then concatenate
                for ag_idx = 1:length(mdp.server.agent_ordering)
                    state = get_singlearm_mdp_state(mdp, ag_idx)
                    a = POMDPs.action(policy, state)
                    push!(vec_action, a)
                end

                # Execute action on problem
                update_server_with_mdp_action!(mdp, vec_action)

                propagate_belt!(mdp.cbelt, mdp.server)
                propagate_arms!(mdp.server, mdp.cbelt, rng)
                mirror_gen_step!(mdp.cbelt, rng)

            end # t = 1:timesteps
            server = mdp.server
            cbelt = mdp.cbelt

        else

            # One of the non-MDP methods
            if parsed_args["baseline"] == "edd"
                baseline_fn = earliest_due_date!
            elseif parsed_args["baseline"] == "hungarian"
                baseline_fn = expected_hungarian!
            elseif parsed_args["baseline"] == "scoba"
                baseline_fn = scoba_conveyor!
            else
                throw(ErrorException("Incorrect baseline argument $(parsed_args["baseline"])!"))
            end

            for t = 1:parsed_args["timesteps"]

                update_time_windows!(cbelt, server)

                # Apply the method, whatever it is
                baseline_fn(server, cbelt)

                propagate_belt!(cbelt, server)
                propagate_arms!(server, cbelt, rng)

                mirror_gen_step!(cbelt, rng)
            end # t = 1:timesteps

        end

        # Update terminal counts if boxes grabbed but not deposited
        count_terminal_grabbed_boxes!(cbelt, server)

        push!(lost_boxes, cbelt.lost_boxes)
        push!(deposited_boxes, cbelt.deposited_boxes)
    end # for i=1:trials

    results_dict["lost"] = lost_boxes
    results_dict["deposited"] = deposited_boxes

    # Print to output file
    open(parsed_args["out_file_name"], "w") do f
        JSON.print(f, results_dict, 2)
    end

end # main

# Run the function and fire away!
main()

# Example call to script
# You can also run julia1.4 scripts/benchmark_conveyor_belt.jl --help to get details of the command line params
# julia1.4 scripts/benchmark_conveyor_belt.jl -g 0.6,0.6,0.6 -n 10 -t 500 scoba cbelt_grasp_pt6_scoba.json


text:
function pvcongruence(W)
	err, i, todelete, vclasses = 10^-6, 1, [], []
	verts = convert(Lar.Points, W')
	kdtree = NearestNeighbors.KDTree(verts)
	newverts = zeros(Int, size(verts,2))
	for vi in 1:size(verts,2)
		if !(vi in todelete)
			nearvs = NearestNeighbors.inrange(kdtree, verts[:,vi], err)
			push!(vclasses,nearvs)
			newverts[nearvs] .= i
			nearvs = setdiff(nearvs, vi)
			todelete = union(todelete, nearvs)
			i += 1
		end
	end
	V = zeros(3,length(vclasses))
	for (k,class) in enumerate(vclasses)
			V[:,k] = sum(W[class,:],dims=1)/length(class)
	end
	return V, vclasses
end

V, vclasses = vcongruence(W)

GL.VIEW([
	GL.GLAxis(GL.Point3d(-1,-1,-1),GL.Point3d(1,1,1))
  #GL.GLPol( V,[[1,2,3,4,5,6,7,8]], GL.COLORS[1])
	GL.GLPol( V,FV, GL.COLORS[1])
]);

function EVcongruence2(Delta_0,vclasses)
	Y = SparseVector[]
	for class in vclasses
		y = sparsevec(sum(Delta_0[:,class],dims=2))
		push!(Y,y)
	end
	copEV = hcat(Y...)

	function econgruence(copEV)
		eclasses = sort([(findnz(copEV[k,:])[1],k) for k=1:copEV.m])
		class = [eclasses[1]]
		classes = []
		for h=2:length(eclasses)
			if eclasses[h][1]!=class[end][1]
				push!(classes, class)
				class = []
			end
			push!(class, eclasses[h])
		end
		return classes
	end

	eclasses = [[idx for (edge,idx) in class] for class in econgruence(copEV)]
	lastclass = setdiff(1:Delta_0.m, cat(eclasses))
	push!(eclasses,lastclass)

	newedges = [eclass[1] for eclass in eclasses]
	newEV = hcat([copEV[k,:] for k=1:copEV.m if k in newedges]...)

	return newEV,eclasses
end


text:
"""
    Flat()

The *flat distribution* is the improper distribution of real numbers that has the improper
probability density function

```math
f(x) = 1.
```
"""
struct Flat <: ContinuousUnivariateDistribution end

Base.minimum(::Flat) = -Inf
Base.maximum(::Flat) = Inf

Base.rand(rng::Random.AbstractRNG, d::Flat) = rand(rng)
Distributions.logpdf(::Flat, x::Real) = zero(x)

# TODO: only implement `logpdf(d, ::Real)` if support for Distributions < 0.24 is dropped
Distributions.pdf(d::Flat, x::Real) = exp(logpdf(d, x))

# For vec support
Distributions.logpdf(::Flat, x::AbstractVector{<:Real}) = zero(x)
Distributions.loglikelihood(::Flat, x::AbstractVector{<:Real}) = zero(eltype(x))

"""
    FlatPos(l::Real)

The *positive flat distribution* with real-valued parameter `l` is the improper distribution
of real numbers that has the improper probability density function

```math
f(x) = \\begin{cases}
0 & \\text{if } x \\leq l, \\\\
1 & \\text{otherwise}.
\\end{cases}
```
"""
struct FlatPos{T<:Real} <: ContinuousUnivariateDistribution
    l::T
end

Base.minimum(d::FlatPos) = d.l
Base.maximum(d::FlatPos) = Inf

Base.rand(rng::Random.AbstractRNG, d::FlatPos) = rand(rng) + d.l
function Distributions.logpdf(d::FlatPos, x::Real)
    z = float(zero(x))
    return x <= d.l ? oftype(z, -Inf) : z
end

# TODO: only implement `logpdf(d, ::Real)` if support for Distributions < 0.24 is dropped
Distributions.pdf(d::FlatPos, x::Real) = exp(logpdf(d, x))

# For vec support
function Distributions.loglikelihood(d::FlatPos, x::AbstractVector{<:Real})
    lower = d.l
    T = float(eltype(x))
    return any(xi <= lower for xi in x) ? T(-Inf) : zero(T)
end

"""
    BinomialLogit(n, logitp)

The *Binomial distribution* with logit parameterization characterizes the number of
successes in a sequence of independent trials.

It has two parameters: `n`, the number of trials, and `logitp`, the logit of the probability
of success in an individual trial, with the distribution

```math
P(X = k) = {n \\choose k}{(\\text{logistic}(logitp))}^k (1 - \\text{logistic}(logitp))^{n-k}, \\quad \\text{ for } k = 0,1,2, \\ldots, n.
```

See also: [`Binomial`](@ref)
"""
struct BinomialLogit{T<:Real,S<:Real} <: DiscreteUnivariateDistribution
    n::Int
    logitp::T
    logconstant::S

    function BinomialLogit{T}(n::Int, logitp::T) where T
        n >= 0 || error("parameter `n` has to be non-negative")
        logconstant = - (log1p(n) + n * StatsFuns.log1pexp(logitp))
        return new{T,typeof(logconstant)}(n, logitp, logconstant)
    end
end

BinomialLogit(n::Int, logitp::Real) = BinomialLogit{typeof(logitp)}(n, logitp)

Base.minimum(::BinomialLogit) = 0
Base.maximum(d::BinomialLogit) = d.n

# TODO: only implement `logpdf(d, k::Real)` if support for Distributions < 0.24 is dropped
Distributions.pdf(d::BinomialLogit, k::Real) = exp(logpdf(d, k))
Distributions.logpdf(d::BinomialLogit, k::Real) = _logpdf(d, k)
Distributions.logpdf(d::BinomialLogit, k::Integer) = _logpdf(d, k)

function _logpdf(d::BinomialLogit, k::Real)
    n, logitp, logconstant = d.n, d.logitp, d.logconstant
    _insupport = insupport(d, k)
    _k = _insupport ? round(Int, k) : 0
    result = logconstant + _k * logitp - SpecialFunctions.logbeta(n - _k + 1, _k + 1)
    return _insupport ? result : oftype(result, -Inf)
end

function Base.rand(rng::Random.AbstractRNG, d::BinomialLogit)
    return rand(rng, Binomial(d.n, logistic(d.logitp)))
end
Distributions.sampler(d::BinomialLogit) = sampler(Binomial(d.n, logistic(d.logitp)))

"""
    BernoulliLogit(logitp::Real)

Create a univariate logit-parameterised Bernoulli distribution.
"""
BernoulliLogit(logitp::Real) = BinomialLogit(1, logitp)

"""
    OrderedLogistic(η, c::AbstractVector)

The *ordered logistic distribution* with real-valued parameter `η` and cutpoints `c` has the
probability mass function

```math
P(X = k) = \\begin{cases}
    1 - \\text{logistic}(\\eta - c_1) & \\text{if } k = 1, \\\\
    \\text{logistic}(\\eta - c_{k-1}) - \\text{logistic}(\\eta - c_k) & \\text{if } 1 < k < K, \\\\
    \\text{logistic}(\\eta - c_{K-1}) & \\text{if } k = K,
\\end{cases}
```
where `K = length(c) + 1`.
"""
struct OrderedLogistic{T1, T2<:AbstractVector} <: DiscreteUnivariateDistribution
    η::T1
    cutpoints::T2

    function OrderedLogistic{T1,T2}(η::T1, cutpoints::T2) where {T1,T2}
        issorted(cutpoints) || error("cutpoints are not sorted")
        return new{typeof(η), typeof(cutpoints)}(η, cutpoints)
    end
end

function OrderedLogistic(η, cutpoints::AbstractVector)
    return OrderedLogistic{typeof(η),typeof(cutpoints)}(η, cutpoints)
end

Base.minimum(d::OrderedLogistic) = 0
Base.maximum(d::OrderedLogistic) = length(d.cutpoints) + 1

# TODO: only implement `logpdf(d, k::Real)` if support for Distributions < 0.24 is dropped
Distributions.pdf(d::OrderedLogistic, k::Real) = exp(logpdf(d, k))
Distributions.logpdf(d::OrderedLogistic, k::Real) = _logpdf(d, k)
Distributions.logpdf(d::OrderedLogistic, k::Integer) = _logpdf(d, k)

function _logpdf(d::OrderedLogistic, k::Real)
    η, cutpoints = d.η, d.cutpoints
    K = length(cutpoints) + 1

    _insupport = insupport(d, k)
    _k = _insupport ? round(Int, k) : 1
    logp = unsafe_logpdf_ordered_logistic(η, cutpoints, K, _k)

    return _insupport ? logp : oftype(logp, -Inf)
end

function Base.rand(rng::Random.AbstractRNG, d::OrderedLogistic)
    η, cutpoints = d.η, d.cutpoints
    K = length(cutpoints) + 1
    # evaluate probability mass function
    ps = map(1:K) do k
        exp(unsafe_logpdf_ordered_logistic(η, cutpoints, K, k))
    end
    k = rand(rng, Categorical(ps))
    return k
end
function Distributions.sampler(d::OrderedLogistic)
    η, cutpoints = d.η, d.cutpoints
    K = length(cutpoints) + 1
    # evaluate probability mass function
    ps = map(1:K) do k
        exp(unsafe_logpdf_ordered_logistic(η, cutpoints, K, k))
    end
    return sampler(Categorical(ps))
end

# unsafe version without bounds checking
function unsafe_logpdf_ordered_logistic(η, cutpoints, K, k::Int)
    @inbounds begin
        logp = if k == 1
            -StatsFuns.log1pexp(η - cutpoints[k])
        elseif k < K
            tmp = StatsFuns.log1pexp(cutpoints[k-1] - η)
            -tmp + StatsFuns.log1mexp(tmp - StatsFuns.log1pexp(cutpoints[k] - η))
        else
            -StatsFuns.log1pexp(cutpoints[k-1] - η)
        end
    end
    return logp
end

"""
    LogPoisson(logλ)

The *Poisson distribution* with logarithmic parameterization of the rate parameter
descibes the number of independent events occurring within a unit time interval, given the
average rate of occurrence ``exp(logλ)``.

The distribution has the probability mass function

```math
P(X = k) = \\frac{e^{k \\cdot logλ}{k!} e^{-e^{logλ}}, \\quad \\text{ for } k = 0,1,2,\\ldots.
```

See also: [`Poisson`](@ref)
"""
struct LogPoisson{T<:Real,S} <: DiscreteUnivariateDistribution
    logλ::T
    λ::S

    function LogPoisson{T}(logλ::T) where T
        λ = exp(logλ)
        return new{T,typeof(λ)}(logλ, λ)
    end
end

LogPoisson(logλ::Real) = LogPoisson{typeof(logλ)}(logλ)

Base.minimum(d::LogPoisson) = 0
Base.maximum(d::LogPoisson) = Inf

function _logpdf(d::LogPoisson, k::Real)
    _insupport = insupport(d, k)
    _k = _insupport ? round(Int, k) : 0
    logp = _k * d.logλ - d.λ - SpecialFunctions.loggamma(_k + 1)

    return _insupport ? logp : oftype(logp, -Inf)
end

# TODO: only implement `logpdf(d, ::Real)` if support for Distributions < 0.24 is dropped
Distributions.pdf(d::LogPoisson, k::Real) = exp(logpdf(d, k))
Distributions.logpdf(d::LogPoisson, k::Integer) = _logpdf(d, k)
Distributions.logpdf(d::LogPoisson, k::Real) = _logpdf(d, k)

Base.rand(rng::Random.AbstractRNG, d::LogPoisson) = rand(rng, Poisson(d.λ))
Distributions.sampler(d::LogPoisson) = sampler(Poisson(d.λ))

Bijectors.logpdf_with_trans(d::NoDist{<:Univariate}, ::Real, ::Bool) = 0
Bijectors.logpdf_with_trans(d::NoDist{<:Multivariate}, ::AbstractVector{<:Real}, ::Bool) = 0
function Bijectors.logpdf_with_trans(d::NoDist{<:Multivariate}, x::AbstractMatrix{<:Real}, ::Bool)
    return zeros(Int, size(x, 2))
end
Bijectors.logpdf_with_trans(d::NoDist{<:Matrixvariate}, ::AbstractMatrix{<:Real}, ::Bool) = 0


text:
# print arrays in pretty way

function printarr(io::IO, a::AbstractArray)
    Base.with_output_limit(()->Base.showarray(io, a, header=false, repr=false))
end

printvec(io::IO, a::AbstractVector) = printarr(io, a')

printarrln(io::IO, a::AbstractArray) = (printarr(io, a); println(io))
printvecln(io::IO, a::AbstractVector) = (printvec(io, a); println(io))

# centralize

centralize(x::AbstractVector, m::AbstractVector) = (isempty(m) ? x : x - m)
centralize(x::AbstractMatrix, m::AbstractVector) = (isempty(m) ? x : x .- m)

decentralize(x::AbstractVector, m::AbstractVector) = (isempty(m) ? x : x + m)
decentralize(x::AbstractMatrix, m::AbstractVector) = (isempty(m) ? x : x .+ m)

# get a full mean vector

fullmean(d::Int, mv::Vector{T}) where T = (isempty(mv) ? zeros(T, d) : mv)

preprocess_mean(X::AbstractMatrix{T}, m) where T<:Real =
    (m === nothing ? vec(mean(X, dims=2)) : m == 0 ? T[] :  m)

# choose the first k values and columns
#
# S must have fields: values & vectors

function extract_kv(fac::Factorization{T}, ord::AbstractVector{Int}, k::Int) where T
    si = ord[1:k]
    vals = fac.values[si]::Vector{T}
    vecs = fac.vectors[:, si]::Matrix{T}
    return (vals, vecs)
end

# symmmetrize a matrix

function symmetrize!(A::Matrix)
    n = size(A, 1)
    @assert size(A, 2) == n
    for j = 1:n
        for i = 1:j-1
            @inbounds A[i,j] = A[j,i]
        end
        for i = j+1:n
            @inbounds A[i,j] = middle(A[i,j], A[j,i])
        end
    end
    return A
end

# percolumn dot

function coldot(X::AbstractMatrix{T}, Y::AbstractMatrix{T}) where T<:Real
    m = size(X, 1)
    n = size(X, 2)
    @assert size(Y) == (m, n)
    R = zeros(T, n)
    for j = 1:n
        R[j] = dot(view(X,:,j), view(Y,:,j))
    end
    return R
end

# qnormalize!

function qnormalize!(X, C)
    # normalize each column of X (say x), such that x'Cx = 1
    m = size(X, 1)
    n = size(X, 2)
    CX = C * X
    for j = 1:n
        x = view(X,:,j)
        cx = view(CX,:,j)
        rmul!(x, inv(sqrt(dot(x, cx))))
    end
    return X
end

# add_diag!

function add_diag!(A::AbstractMatrix, v::Real)
    # add v to diagonal of A
    m = size(A, 1)
    n = size(A, 2)
    @assert m == n
    if v != zero(v)
        for i = 1:n
            @inbounds A[i,i] += v
        end
    end
    return A
end

# regularize a symmetric matrix
function regularize_symmat!(A::AbstractMatrix{T}, lambda::Real) where T<:Real
    if lambda > 0
        emax = eigmax(Symmetric(A))
        add_diag!(A, emax * lambda)
    end
    return A
end

"""
    calcscattermat([covestimator::CovarianceEstimator], Z::DenseMatrix)

Calculate the scatter matrix of centered data `Z` based on a covariance
matrix calculated using covariance estimator `covestimator` (by default,
`SimpleCovariance()`).
"""
function calcscattermat(covestimator::CovarianceEstimator, Z::DenseMatrix{T}) where T<:Real
    return cov(covestimator, Z; dims=2, mean=zeros(T, size(Z, 1)))*size(Z, 2)
end

function calcscattermat(Z::DenseMatrix)
    return calcscattermat(SimpleCovariance(), Z)
end

# calculate pairwise kernel
function pairwise!(K::AbstractVecOrMat{T}, kernel::Function,
                   X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) where T<:AbstractFloat
    n = size(X, 2)
    m = size(Y, 2)
    for j = 1:m
        aj = view(Y, :, j)
        for i in j:n
            @inbounds K[i, j] = kernel(view(X, :, i), aj)[]
        end
        j <= n && for i in 1:(j - 1)
            @inbounds K[i, j] = K[j, i]   # leveraging the symmetry
        end
    end
    K
end

pairwise!(K::AbstractVecOrMat{T}, kernel::Function, X::AbstractVecOrMat{T}) where T<:AbstractFloat =
    pairwise!(K, kernel, X, X)

function pairwise(kernel::Function, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) where T<:AbstractFloat
    n = size(X, 2)
    m = size(Y, 2)
    K = similar(X, n, m)
    pairwise!(K, kernel, X, Y)
end

pairwise(kernel::Function, X::AbstractVecOrMat{T}) where T<:AbstractFloat =
    pairwise(kernel, X, X)


text:
module Bibliography

# BibInternal
import BibInternal
import BibInternal: AbstractEntry, Entry

# BibParser
import BibParser, BibParser.BibTeX

# Others
import DataStructures
import DataStructures.OrderedSet

export export_bibtex, import_bibtex
export export_web, bibtex_to_web
export select
export sort_bibliography!

include("select.jl")
include("sort_bibliography.jl")
include("bibtex.jl")
include("csl.jl")
include("staticweb.jl")

"""
    export_bibtex(target, bibliography)
Export a bibliography to BibTeX format.
"""
function export_bibtex(target, bibliography)
    data = export_bibtex(bibliography)
    if target != ""
        f = open(target, "w")
        write(f, data)
        close(f)
    end
    return data
end

"""
    bibtex_to_web(source::String)
Convert a BibTeX file to a web compatible format, specifically for the [StaticWebPages.jl](https://github.com/Humans-of-Julia/StaticWebPages.jl) package.
"""
bibtex_to_web(source) = export_web(import_bibtex(source))

end # module


text:
using HDF5
using Dates
#using Images
#https://portal.hdfgroup.org/display/HDF5/Introduction+to+HDF5
#https://imaris.oxinst.com/support/imaris-file-format
# https://github.com/imaris/ImarisWriter/blob/master/doc/Imaris5FileFormat.pdf
# https://docs.julialang.org/en/v1/stdlib/Dates/#Dates.format-Tuple{TimeType,AbstractString}
# https://cpb-us-w2.wpmucdn.com/sites.wustl.edu/dist/f/1861/files/2019/02/02-imaris-image-properties-and-edit-2jj1avj.pdf
# https://github.com/tlambert03/imarispy/blob/master/imarispy/imaris.py

"""
Version Commit
0.1		hf, first verion
0.2     hf, update API to follow upstream HDF5.jl v0.15.4
TODO: genrate thumb and fix timestamps
      add function to append image
"""

"Save 4D UInt16 Array(x,y,z,t) to .ims file"
function  save2ims(array, fname::String="myfile.ims";
              subsamp=((1, 1, 1)),
              chunks=((128, 128, 5)),
              compression=2,
              thumbsize=256,
              dx=0.108, dz=0.5, dt_min=10)

    df = DateFormat("Y-mm-dd HH:MM:SS.sss")
    RecordingDate = DateTime(2020,5,7,14,19,22,0)

    nx, ny, nz, nt = size(array)
    nc = 1
    nr = 1     #nr = length(subsamp)

    # GROUPS of Imaris file
    GROUPS = [
        "DataSetInfo",
        "Thumbnail",
        "DataSetTimes",
        "DataSetInfo/Imaris",
        "DataSetInfo/Image",
        "DataSetInfo/TimeInfo"
    ]

    # Attributes of Imaris file
o   ATTRS = [
        ("/", "ImarisDataSet", "ImarisDataSet"),
        ("/", "ImarisVersion", "5.5.0"),
        ("/", "DataSetInfoDirectoryName", "DataSetInfo"),
        ("/", "ThumbnailDirectoryName", "Thumbnail"),
        ("/", "DataSetDirectoryName", "DataSet"),
        ("DataSetInfo/Imaris", "Version", "9.5"),
        ("DataSetInfo/Imaris", "ThumbnailMode", "thumbnailMIP"),
        ("DataSetInfo/Imaris", "ThumbnailSize", thumbsize),
        ("DataSetInfo/Image", "X", nx),
        ("DataSetInfo/Image", "Y", ny),
        ("DataSetInfo/Image", "Z", nz),
        ("DataSetInfo/Image", "NumberOfChannels", nc),
        ("DataSetInfo/Image", "Noc", nc),
        ("DataSetInfo/Image", "Unit", "um"),
        ("DataSetInfo/Image", "Description", "Hi"),
        ("DataSetInfo/Image", "MicroscopeModality", "Inverted Microscope"),
        ("DataSetInfo/Image", "RecordingDate", Dates.format(RecordingDate, df)),
        ("DataSetInfo/Image", "Name", "qblab"),
        ("DataSetInfo/Image", "ExtMin0", "0"),
        ("DataSetInfo/Image", "ExtMin1", "0"),
        ("DataSetInfo/Image", "ExtMin2", "0"),
        ("DataSetInfo/Image", "ExtMax0", nx * dx),
        ("DataSetInfo/Image", "ExtMax1", ny * dx),
        ("DataSetInfo/Image", "ExtMax2", nz * dz),
        ("DataSetInfo/Image", "LensPower", "40x"),
        ("DataSetInfo/TimeInfo", "DatasetTimePoints", nt),
        ("DataSetInfo/TimeInfo", "FileTimePoints", nt)
    ];

    COLORS = ("1 1 1", "1 0 1", "1 0 0", "0 0 1")
    for c in 0:nc-1
        grp = "DataSetInfo/Channel $c"
        push!(GROUPS, grp)
        push!(ATTRS, (grp, "ColorOpacity", 1))
        push!(ATTRS, (grp, "ColorMode", "BaseColor"))
        push!(ATTRS, (grp, "Color", COLORS[3]))
        push!(ATTRS, (grp, "GammaCorrection", 1))
        push!(ATTRS, (grp, "ColorRange", "0 255"))
        push!(ATTRS, (grp, "Name", "Channel $c"))
    end

    for t in 0:nt-1
        strr = Dates.format(RecordingDate + Dates.Minute(dt_min*t), df)
        #strr = "2020-05-07 {:02d}:{:02d}:{:02d}.000".format(h, m, s)
        push!(ATTRS, ("DataSetInfo/TimeInfo", "TimePoint$t", strr))
        #TODO: Imaris don't use these timestamps as time axis
    end

    h5open(fname, "w") do file
        for grp in GROUPS
            create_group(file, grp)
        end

        for info in ATTRS
            grp, key = info[1:2]
            value = "$(info[3])"
            h5a_write_S1(file[grp], key, value)
        end

        file["Thumbnail/Data"] = UInt8.(maximum(array[:,:, 1:20,1], dims=3).>>8)
        #file["Thumbnail/Data"] = UInt8.(maximum(rawview(array[:,:, 1:20,1]), dims=3).>>8)

        # add data
        for t in 0:nt-1
            for c in 0:nc-1
                data = array[:, :, :, t+1]
                for r in 0:nr-1
                    edges, count = build_histogram(data, 256)
                    grp = create_group(file, "/DataSet/ResolutionLevel $r/TimePoint $t/Channel $c/")
                    grp["Histogram"] = UInt64.(count[1:end])
                    h5a_write_S1(grp, "HistogramMin", "1000")
                    h5a_write_S1(grp, "HistogramMax", "$(edges[end])")
                    grp["Data", chunk= chunks , compress=compression] = data
                    #grp["Data", "chunk", (256,256,4)] = data
                    h5a_write_S1(grp, "ImageSizeX", "$(size(data)[1])")
                    h5a_write_S1(grp, "ImageSizeY", "$(size(data)[2])")
                    h5a_write_S1(grp, "ImageSizeZ", "$(size(data)[3])")
                end
            end
        end
    end
	fname
end

"Append image to exits imaris5 file"
#=
function ims_append(fname, img)
    # find orignal array size
    h5open(fname, "r+") do file
        # add data
        for t in 0:nt-1
            for c in 0:nc-1
                data = array[:, :, :, t+1]
                for r in 0:nr-1
                    edges, count = build_histogram(data, 256)
                    grp = create_group(file, "/DataSet/ResolutionLevel $r/TimePoint $t/Channel $c/")
                    grp["Histogram"] = UInt64.(count[1:end])
                    h5a_write_S1(grp, "HistogramMin", "1000")
                    h5a_write_S1(grp, "HistogramMax", "$(edges[end])")
                    grp["Data", chunk= chunks , compress=compression] = data
                    #grp["Data", "chunk", (256,256,4)] = data
                    h5a_write_S1(grp, "ImageSizeX", "$(size(data)[1])")
                    h5a_write_S1(grp, "ImageSizeY", "$(size(data)[2])")
                    h5a_write_S1(grp, "ImageSizeZ", "$(size(data)[3])")
                end
            end
        end
    end

end
=#

"Encode string with ASCII S1"
function h5a_write_S1(parent, key, value)
    # https://support.hdfgroup.org/HDF5/doc1.6/UG/11_Datatypes.html
    # https://support.hdfgroup.org/HDF5/doc/RM/PredefDTypes.html
    dspace = HDF5.dataspace((length(value),))
    try
        #attr = create_attribute(parent, key, datatype(HDF5.H5T_CSET_ASCII), dspace)
        #NOTE: datatype(obj) -> return the Datatype of obj
        #      HDF5.Datatype(id) -> return Datatype responding to id
        attr = create_attribute(parent, key, HDF5.Datatype(HDF5.H5T_C_S1), dspace)
        #write(attr, "$value")
        write_attribute(attr, HDF5.Datatype(HDF5.H5T_C_S1), "$value")
    finally
        close(dspace)
    end
    nothing
end

function str2int(str)
    local tmp =""
    for i in str
        tmp=tmp*i
    end
    parse(Int, tmp)
end

"Load imaris file to multi dimension image"
function loadims(imsname::String)
    local attr_t_len = h5readattr(imsname, "DataSetInfo/TimeInfo")["DatasetTimePoints"];
    local attr_x_len = h5readattr(imsname, "DataSetInfo/Image")["X"];
    local attr_y_len = h5readattr(imsname, "DataSetInfo/Image")["Y"];
    local attr_z_len = h5readattr(imsname, "DataSetInfo/Image")["Z"];
    local t_len  = str2int(attr_t_len)
    local x_depth = str2int(attr_x_len)
    local y_depth = str2int(attr_y_len)
    local z_depth = str2int(attr_z_len)
    #local img = zeros(Gray{N0f16}, x_depth, y_depth, z_depth, t_len);
    #local img = zeros(UInt16, x_depth, y_depth, z_depth, t_len);
    local img = Array{UInt16}(undef, x_depth, y_depth, z_depth, t_len);
    h5open(imsname, "r") do imsfile
        for i in 1:t_len
            data = read(imsfile,
                        "DataSet/ResolutionLevel 0/TimePoint $(i-1)/Channel 0/Data")
            #img[:, :, :, i] = reinterpret(N0f16, data)
            img[:, :, :, i] = data
        end
    end
    return img
end

function loadims2(imsname::String)
    local attr_t_len = h5readattr(imsname, "DataSetInfo/TimeInfo")["DatasetTimePoints"];
    local t_len  = str2int(attr_t_len)
    local z_depth = 20
    local img = zeros(Gray{N0f16}, 512, 512, t_len*z_depth);
    h5open(imsname) do imsfile
        for i in 0:t_len-1
			data = imsfile["DataSet/ResolutionLevel 0/TimePoint $i/Channel 0/Data"]
            img[:, :, i*20+1:20*(i+1)] = reinterpret(N0f16,data[1:512,1:512,1:20]);
        end
    end
    img
end


text:
using NiLang

@i function loss1(out!, x, y)
    y ⊕ x
    out! -= exp(x)
    out! += exp(y)
end

@i function loss2(out!, x, y)
    y ⊕ x
    out! -= exp(x)
    x -= log(-out!)
    out! += exp(y)
end

function train(lf)
    loss = Float64[]
    y = 1.6
    for i=1:100
        out!, x = 0.0, 0.3
        @instr lf(out!, x, y)
        push!(loss, out!)
        out! = 1.0
        @instr (~lf)(out!, x, y)
    end
    loss
end

res = train(loss2)
using DelimitedFiles
writedlm("train2.dat", res)


text:
  function unique2(vec::AbstractVector)
      n = length(vec)
      id = unique(vec)
      nUnique = length(id)
      ia = zeros(Int64, nUnique)
      ic = zeros(Int64, n)
      for i = 1:nUnique
          @inbounds ia[i] = findfirst(isequal(id[i]), vec)
      end
      for i = 1:n
          @inbounds ic[i] = findfirst(isequal(vec[i]), id)
      end
      return id, ia, ic
  end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using MultidimensionalSpectroscopy
using Test

@testset "MultidimensionalSpectroscopy.jl" begin
    # Write your tests here.
end

max_stars_repo_name: lbmnky/CMDS

max_issues_repo_name: lbmnky/CMDS

max_forks_repo_name: lbmnky/CMDS


text:
using Pkg; Pkg.activate(@__DIR__)
using Comrade
using Distributions
using ComradeGalactic
using ComradeAHMC
using GalacticBBO
using GalacticOptimJL
using Plots
using StatsBase

# load eht-imaging we use this to load eht data
load_ehtim()
# To download the data visit https://doi.org/10.25739/g85n-f134
obs = ehtim.obsdata.load_uvfits(joinpath(@__DIR__, "SR1_M87_2017_096_lo_hops_netcal_StokesI.uvfits"))
obs.add_scans()
# kill 0-baselines since we don't care about
# large scale flux and make scan-average data
obs = obs.flag_uvdist(uv_min=0.1e9).avg_coherent(0.0, scan_avg=true)
# extract log closure amplitudes and closure phases
dlcamp = extract_lcamp(obs)
dcphase = extract_cphase(obs)
# form the likelihood
lklhd = RadioLikelihood(dlcamp, dcphase)
# build the model here we fit a ring with a azimuthal
#brightness variation and a Gaussian
function model(θ)
    (;rad, wid, a, b, f, sig, asy, pa, x, y) = θ
    ring = f*smoothed(stretched(MRing((a,), (b,)), μas2rad(rad), μas2rad(rad)), μas2rad(wid))
    g = (1-f)*shifted(rotated(stretched(Gaussian(), μas2rad(sig)*asy, μas2rad(sig)), pa), μas2rad(x), μas2rad(y))
    return ring + g
end
# define the priors
prior = (
          rad = Uniform(10.0, 30.0),
          wid = Uniform(1.0, 10.0),
          a = Uniform(-0.5, 0.5), b = Uniform(-0.5, 0.5),
          f = Uniform(0.0, 1.0),
          sig = Uniform((1.0), (60.0)),
          asy = Uniform(0.0, 0.9),
          pa = Uniform(0.0, 1π),
          x = Uniform(-(80.0), (80.0)),
          y = Uniform(-(80.0), (80.0))
        )
# Now form the posterior
post = Posterior(lklhd, prior, model)
# We will use HMC to sample the posterior.
# First we will find a reasonable starting location using GalacticOptim
# For optimization we need to specify what transform to use. Here we will transform to
# the unit hypercube
tpost = asflat(post)
ndim = dimension(tpost)
f = OptimizationFunction(tpost, GalacticOptim.AutoForwardDiff())
prob = OptimizationProblem(f, rand(ndim), nothing, lb=fill(-5.0, ndim), ub = fill(5.0, ndim))
sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited(); maxiters=50_000)

# Now let's get the Laplace approximation since it is cheap!
prob = OptimizationProblem(f, sol.u, nothing)
ldist = laplace(prob, LBFGS(); show_trace=true)

# transform the solution back to regular space
xopt = transform(tpost, sol)
# Let's see if the best fit looks reasonable by plotting normalized residuals for the
# log-closure amplitudes
residual(model(xopt), dlcamp)
# we can also plot the best fit model or maximum likelihood estimate (MLE)
plot(model(xopt), xlims=(-80.0,80.0),ylims=(-80.0,80.0), colorbar=nothing, title="MLE M87")

# Comrade is all about uncertainty quantification so now let's find the posterior!
# To do this we will use the `AdvancedHMC` package or rather its interface to Comrade.
metric = DiagEuclideanMetric(ndim)
chain, stats = sample(post, AHMC(;metric), 2000; nadapts=8000, init_params=xopt)
# Now let's find the mean image
images = intensitymap.(model.(sample(chain, 200)), μas2rad(160.0), μas2rad(160.0), 256, 256)
plot(mean(images), xlims=(-80.0, 80.0), ylims=(-80.0,80.0), colorbar=nothing, title="Mean M87")


text:
module Proj4

using PROJ_jll
using CEnum
using StaticArrays
using CoordinateTransformations

export Projection, # proj_types.jl
       transform, transform!,  # proj_functions.jl
       is_latlong, is_geocent, compare_datums, spheroid_params,
       xy2lonlat, xy2lonlat!, lonlat2xy, lonlat2xy!

# geodesic support
export geod_direct, geod_inverse, geod_destination, geod_distance

# re-export CoordinateTransformations methods we implement
export compose, ∘

# these are part of the deprecated PROJ API, and will be removed soon
include("projection_codes.jl") # ESRI and EPSG projection strings
include("proj_capi.jl") # low-level C-facing functions (corresponding to src/proj_api.h)
include("proj_geodesic.jl") # low-level C-facing functions (corresponding to src/geodesic.h)

# these are part of the new PROJ API
include("proj_common.jl")
include("coord.jl")
include("proj_c.jl")
include("error.jl")

function _version()
    m = match(r"(\d+).(\d+).(\d+),.+", _get_release())
    VersionNumber(parse(Int, m[1]), parse(Int, m[2]), parse(Int, m[3]))
end

"Parsed version number for the underlying version of libproj"
const version = _version()

# Detect underlying libproj support for geodesic calculations
const has_geodesic_support = true

include("proj_types.jl") # type definitions for proj objects
include("proj_functions.jl") # user-facing proj functions

"Get a global error string in human readable form"
error_message() = _strerrno()

"""
Load a null-terminated list of strings

It takes a `PROJ_STRING_LIST`, which is a `Ptr{Cstring}`, and returns a `Vector{String}`.
"""
function unsafe_loadstringlist(ptr::Ptr{Cstring})
    strings = Vector{String}()
    (ptr == C_NULL) && return strings
    i = 1
    cstring = unsafe_load(ptr, i)
    while cstring != C_NULL
        push!(strings, unsafe_string(cstring))
        i += 1
        cstring = unsafe_load(ptr, i)
    end
    strings
end

const PROJ_LIB = Ref{String}()

"Module initialization function"
function __init__()
    # register custom error handler
    funcptr = @cfunction(log_func, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Cstring))
    proj_log_func(C_NULL, funcptr)

    # point to the location of the provided shared resources
    PROJ_LIB[] = joinpath(PROJ_jll.artifact_dir, "share", "proj")
    proj_context_set_search_paths(1, [PROJ_LIB[]])
end

end # module


text:
ks = collect(setdiff(intersect(keys(io), keys(bluenote)), [:supply,:use]))
x = Drop.(:value, [missing, 0.0], "==")

for k in ks
    cols = intersect(propertynames(bluenote[k]), propertynames(io[k]))
    io_temp = edit_with(copy(io[k]), x)[:, cols]
    bluenote_temp = edit_with(copy(bluenote[k]), x)[:, cols]

    compare_values([io_temp, bluenote_temp], [:slide,:bluenote])
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using Sim6502
using BenchmarkTools

macro program(vec)
    esc(quote
        cs.ram[0x0600:(0x0600+length($vec)-1)] = $vec
    end)
end

macro p(expr)
    esc(:($expr; println(cs.cpu)))
end

function makechipset()
    cs = Chipset()

    cs.ram[0x00] = 0x01
    cs.ram[0x01] = 0x02

    cs.ram[0xa000] = 0x05
    cs.ram[0xa001] = 0x06

    # @program [0xaa, 0xa5, 0x00, 0xad, 0x01, 0xa0, 0xaa, 0xaa, 0xaa]
    @program [0xa9, 0xaa, 0xa9, 0xab, 0xa9, 0xac]

    cs
end

function showstates()
    cs = makechipset()

    @p op!(cs)
    @p op!(cs)
    @p op!(cs)

    cs
end

function runbenches()
    ref = @benchmarkable begin
        cs.cpu.A = 0xaa
        cs.cpu.A = 0xab
        cs.cpu.A = 0xac
    end setup=(cs = makechipset())

    b = @benchmarkable begin
        op!(cs)
        op!(cs)
        op!(cs)
    end setup=(cs = makechipset())

    run(ref), run(b)
end

# ref, b = runbenches()


text:
using StaticThreadSchedules
using Documenter

makedocs(;
    modules=[StaticThreadSchedules],
    authors="Zack Li",
    repo="https://github.com/xzackli/StaticThreadSchedules.jl/blob/{commit}{path}#L{line}",
    sitename="StaticThreadSchedules.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://xzackli.github.io/StaticThreadSchedules.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/xzackli/StaticThreadSchedules.jl",
)


text:
#===============================================================================

    Module QpcrAnalysis.jl

===============================================================================#

#===============================================================================
    start of module definition >>
===============================================================================#

__precompile__()

# using Base

"""
Bioinformatics module for Chai's open-source Real-Time PCR instrument.
[https://www.chaibio.com](https://www.chaibio.com)
"""
module QpcrAnalysis
    const MODULE_NAME = "QpcrAnalysis"

    # using Clustering, Combinatorics, DataFrames, DataStructures, Dierckx, Ipopt, JLD, JSON,
    # JuMP, MySQL, NamedTuples, DataArrays
    ## `NLopt` needed on BBB but not on PC ("ERROR: LoadError: Declaring __precompile__(false)
    ## is not allowed in files that are being precompiled".
    ## "ERROR: Failed to precompile NLopt to /root/.julia/lib/v0.6/NLopt.ji")
    ## In addition, `HttpServer` for "juliaserver.jl"

    ## Other functions than `include` read files from `pwd()` only instead of also `LOAD_PATH`.
    ## `pwd()` shows the present working directory in module `Main`, instead of the directory
    ## where "QpcrAnalysis.jl" is located. Therefore `LOAD_FROM_DIR` needs to be defined
    ## for those functions to find files in the directory where "QpcrAnalysis.jl" is located.
    const LOAD_FROM_DIR = LOAD_PATH[find(LOAD_PATH) do path_
        isfile("$path_/$MODULE_NAME.jl")
    end][1] ## slice by boolean vector returned a one-element vector. Assumption: LOAD_PATH is global

    ## default data width in production mode:  32 bits (BBB)
    ## default data width in development mode: 64 bits
    const DEVELOPMENT_MODE = "development"
    const PRODUCTION_MODE  = "production"
    const production_env   = (get(ENV, "JULIA_ENV", nothing) == PRODUCTION_MODE)
    ##const production_env = true

    const Float_T = Float64
	const Int_T = (Sys.WORD_SIZE == 32) ? Int32 : Int64
    const NaN_T = NaN64
    const Inf_T = Inf64

#===============================================================================
    include each script, generally in the order of workflow
===============================================================================#

    ## define constants & macros:
    include("defines/keystring_constants.jl")
    include("defines/enums.jl")
    include("defines/macros.jl")

    ## define shared functions:
    include("shared_functions.jl")

    ## define structs and types:
    include("defines/Field.jl")
    include("defines/Input.jl")
    ## calibration
    include("defines/NumberOfChannels.jl")
    include("defines/RawData.jl")
    include("defines/CalibrationData.jl")
    include("defines/CalibrationParameters.jl")
    include("defines/DeconvolutionMatrices.jl")
    ## allelic discrimination
    include("defines/ClusterAnalysisResult.jl")
    include("defines/UniqCombinCenters.jl")
    include("defines/AssignGenosResult.jl")
    ## amplification models
    include("amp_models/defines/AmpModel.jl")
    include("amp_models/defines/AmpModelFit.jl")
    include("amp_models/defines/SFC_model_definitions.jl")
    include("amp_models/defines/SFCModelDef.jl")
    ## generate amplification model definitions
    include("amp_models/generate_SFC_models.jl")
    include("amp_models/MAKx.jl")
    include("amp_models/MAKERGAULx.jl")
    include("amp_models/defines/AmpModelResults.jl")
    ## amplification
    include("amp_models/defines/AmpInput.jl")
    include("amp_models/defines/AmpOutput.jl")
    ## melting curve
    include("defines/McInput.jl")
    include("defines/PeakIndices.jl")
    include("defines/Peak.jl")
    include("defines/McPeakOutput.jl")
    include("defines/McOutput.jl")
    ## standard curve
    include("defines/StandardCurveResult.jl")
    ## thermal consistency
    include("defines/TmCheck1w.jl")
    include("defines/ThermalConsistencyOutput.jl")

    ## this code is hidden from the parser on the BeagleBone
    @static if !production_env
        ## development & testing
        import Base.Test
        using FactCheck
        include("../test/test_functions.jl")
        ## data format verification
        include("../test/verify_request.jl")
        include("../test/verify_response.jl")
    end

    ## define functions:
    ## dispatch
    include("dispatch.jl")
    ## calibration
    include("calibration.jl")
    ## amplification
    include("amplification.jl")
    include("amp_analysis.jl")
    include("amp_fit_model.jl")
    ## melting curve
    include("melting_curve.jl")
    include("mc_analysis.jl")
    include("mc_peak_analysis.jl")
    include("supsmu.jl")
    ## standard curve
    include("standard_curve.jl")
    ## custom analyses
    include("custom_analyses/optical_calibration.jl")
    include("custom_analyses/optical_test_single_channel.jl")
    include("custom_analyses/optical_test_dual_channel.jl")
    include("custom_analyses/thermal_consistency.jl")
    include("custom_analyses/thermal_performance_diagnostic.jl")
    # include("custom_analyses/your_own_analysis.jl")

#===============================================================================
    create logger >>
===============================================================================#

    ## create module level logger
    ## this can be precompiled
    using Memento
    import FactCheck.clear_results
    logger = getlogger(current_module())

#===============================================================================
    runtime initialization >>
===============================================================================#

    ## This function contains stuff that needs to happen
    ## at runtime when the module is loaded
    function __init__()
        ## create log handler
        ## NB changes to the default logger must happen at runtime
        ## otherwise segfaults are liable to occur
        push!(logger,
            DefaultHandler(
                production_env ? "/var/log/julia.log" : "/tmp/julia.log",
                DefaultFormatter("[ {date} | {level} ]: {msg}")))
        ## debug-level logging in production environment (provisional setting)
        production_env  && setlevel!(logger, "debug")
        ## debug-level logging in development environment
        !production_env && setlevel!(logger, "debug")
        ## Register the module level logger at runtime
        ## so it is accessible via `get_logger(QpcrAnalysis)`
        Memento.register(logger)
        #
        ## clear Fact Checks
        !production_env && FactCheck.clear_results()
    end ## __init__()

end ## module QpcrAnalysis

#===============================================================================
    << end of module definition
===============================================================================#

# for precompiling
using HTTP, JSON, Memento
#using Clustering, Combinatorics, DataArrays, DataFrames
#using DataStructures, Dierckx, HDF5, Ipopt, JLD, JSON, JuMP
#using MathProgBase, Nullables, SpecialFunctions, FactCheck, FunctionalData, Match, Syslogs, Memento, HttpCommon, HttpParser, IniFile

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
include( "../src/argmaxProb.jl" )
using Test

  println( "running argmaxProb unit tests." );

  vec = repeat( [23.8], outer=[3 1] )[:,1];
  append!( vec, rand( 4.0:20.0, 100 ) );
  nrRuns = 100000;
  res    = zeros( Float64, 103 );
  for i in 1:1:nrRuns
    index = argmaxUniProb.argmaxProb( vec );
    res[index] += 1;
  end
  minVal = minimum( res[1:3] );
  maxVal = maximum( res[1:3] );

  @test maxVal / minVal > 0.95
  @test maxVal / minVal < 1.05

  len        = 100;
  lmatrix    = repeat( [32.0], outer=[len 3] );
  rmatrix    = rand( 4.0:20.0, len, 3 );
  vmatrix    = hcat( lmatrix, rmatrix );
  minIndices = argmaxUniProb.argmaxProb( vmatrix );

  @test all( minIndices .>= 1 )
  @test all( minIndices .<= 3 )

  minIndices = zeros( Int64, size( vmatrix, 2 ) );
  argmaxUniProb.argmaxProb!( vmatrix, minIndices );

  @test all( minIndices .>= 1 )
  @test all( minIndices .<= 3 )


text:
function unwrap!(x, period=2π)
	y = convert(eltype(x), period)
	v = first(x)
	@inbounds for k = eachindex(x)
		x[k] = v = v + rem(x[k] - v,  y, RoundNearest)
	end
end

export unwrap!


text:
using Test
using Classes

@test superclass(Class) === nothing

@test isclass(AbstractClass) == false       # not concrete
@test isclass(Int) == false                 # not <: AbstractClass

@class Foo begin
   foo::Int

   Foo() = Foo(0)

   # Although Foo is immutable, subclasses might not be,
   # so it's still useful to define this method.
   function Foo(self::AbstractFoo)
        self.foo = 0
    end
end

@test classof(AbstractFoo) == Foo
@test classof(Foo) == Foo

@test superclass(Foo) == Class
@test_throws Exception superclass(AbstractFoo)

@class mutable Bar <: Foo begin
    bar::Int

    # Mutable classes can use this pattern
    function Bar(self::Union{Nothing, AbstractBar}=nothing)
        self = (self === nothing ? new() : self)
        superclass(Bar)(self)
        Bar(self, 0)
    end
end

@class mutable Baz <: Bar begin
   baz::Int

   function Baz(self::Union{Nothing, AbstractBaz}=nothing)
        self = (self === nothing ? new() : self)
        superclass(Baz)(self)
        Baz(self, 0)
    end
end

function sum(obj::AbstractBar)
    return obj.foo + obj.bar
end

@test Set(superclasses(Baz)) == Set([Foo, Bar, Class])
@test Set(subclasses(Foo)) == Set(Any[Bar, Baz])

x = Foo(1)
y = Bar(10, 11)
z = Baz(100, 101, 102)

@test fieldnames(Foo) == (:foo,)
@test fieldnames(Bar) == (:foo, :bar)
@test fieldnames(Baz) == (:foo, :bar, :baz)

foo(x::AbstractFoo) = x.foo

@test foo(x) == 1
@test foo(y) == 10
@test foo(z) == 100

@test sum(y) == 21

get_bar(x::AbstractBar) = x.bar

@test get_bar(y) == 11
@test get_bar(z) == 101

@test_throws Exception get_bar(x)

# Mutable
set_foo!(x::AbstractFoo, value) = (x.foo = value)
set_foo!(z, 1000)
@test foo(z) == 1000

# Immutable
@test_throws Exception x.foo = 1000

# test that where clause is amended properly
zzz(obj::AbstractFoo, bar::T) where {T} = T

@test zzz(x, :x) == Symbol
@test zzz(y, 10.6) == Float64

# Test other @class structural errors
@test_throws(LoadError, eval(Meta.parse("@class X2 x y")))
@test_throws(LoadError, eval(Meta.parse("@class (:junk,)")))

# Test that classof fails if an abstract class has multiple concrete classes
@class Blink <: Baz
struct RenegadeStruct <: AbstractBlink end

@test_throws Exception classof(AbstractBlink)

# Test that parameterized type is handled properly
@class TupleHolder{NT <: NamedTuple} begin
    nt::NT
end

nt = (foo=1, bar=2)
NT = typeof(nt)
th = TupleHolder{NT}(nt)

@test typeof(th).parameters[1] == NT
@test th.nt.foo == 1
@test th.nt.bar == 2

# Test updating using instance of parent class
bar = Bar(1, 2)
baz = Baz(100, 101, 102)

upd = Baz(bar, 555)
@test upd.foo == 1 && upd.bar == 2 && upd.baz == 555

# ...and with parameterized types
@class mutable SubTupleHolder{NT <: NamedTuple} <: Baz begin
    nt::NT
end

sub = SubTupleHolder{NT}(z, nt)
@test sub.nt.foo == 1 && sub.nt.bar == 2 && sub.foo == 1000 && sub.bar == 101 && sub.baz == 102

xyz = SubTupleHolder(sub, 10, 20, 30, (foo=111, bar=222))
@test xyz.nt.foo == 111 && xyz.nt.bar == 222 && xyz.foo == 10

@class Parameterized{T1 <: Foo, T2 <: Foo} begin
    one::T1
    two::T2
end

@class ParameterizedSub{T3 <: TupleHolder, T4 <: TupleHolder} <: Parameterized begin
    x::Float64
    y::Float64
end

# non-constrained struct type parameter
@class mutable AbstractLog{T} begin
    pos::Matrix{T}
end

obj = AbstractLog{Float64}([1. 2.; 3. 4.])

@test obj.pos == [1. 2.; 3. 4.]

# TBD: add tests on these

# Generated: needs work on parameterized types (T1, T2 are not defined)
# - convert param names to gensyms to avoid collisions
#
# function (ParameterizedSub{T3, T4}(one::T1, two::T2, x::Float64, y::Float64; )::Any) where {T3 <: TupleHolder, T4 <: TupleHolder}
#     #= /Users/rjp/.julia/packages/MacroTools/4AjBS/src/utils.jl:302 =#
#     new{T3, T4}(one, two, x, y)
# end


text:
####################################### Clustering with bins ###################
function sample_cluster(tss_bins,tss,n)
    sample_tss=[]
    tinf=length(tss[1].time_stamp)
    loops=0
    while (length(sample_tss)<n && loops<10000)
        sample_i=rand(1:tinf)
        sample=tss[1][sample_i,:].time_stamp
        ok,tss_bins=check_sample(sample,tss,tss_bins)#working on this function
        if (ok)
            println(length(sample_tss))
            push!(sample_tss,sample)#store sample set
            loops=0
        end
        loops=loops+1
    end
    return sample_tss
end

function check_sample(sample,tss,tss_bins)
    ok=true;locations=[]
    for (key0,ts_bins) in tss_bins
        for (key1,ts_bin) in ts_bins
            if (length(filter(row -> row.time_stamp == sample, ts_bin["ts"]).time_stamp)>0)
                if (ts_bin["full"]==false)
                    push!(locations,(key0,key1))
                else
                    ok=false
                    @goto bad_sample
                end
            end
        end
    end
    if (length(locations)==length(tss_bins))
        for location in locations
            if (haskey(tss_bins[location[1]][location[2]], "chosen")==false);
                push!(tss_bins[location[1]][location[2]],"chosen"=>1);
            else
                tss_bins[location[1]][location[2]]["chosen"]=tss_bins[location[1]][location[2]]["chosen"]+1
            end
            if (tss_bins[location[1]][location[2]]["elements"]-tss_bins[location[1]][location[2]]["chosen"]==0)
                tss_bins[location[1]][location[2]]["full"]=true
            end
        end
    else
        ok=false
    end
    @label bad_sample
    return ok,tss_bins
end

function cluster_probs(tss_bins,tss,n)
    tinf=length(tss[1].time_stamp)
    for (k0,ts_bins) in tss_bins
        for (k1,ts_bin) in ts_bins
            elements=length(ts_bin["ts"].time_stamp)==0 ? 0 : ceil(Int64,n*length(ts_bin["ts"].time_stamp)/tinf)
            push!(ts_bin,"elements"=>elements)
            if (elements>0); push!(ts_bin,"full"=>false);else push!(ts_bin,"full"=>true);end
        end
    end
    return tss_bins
end
#tss is an array of time series arays, n is the desired sample size
function cluster_ts(tss,n)
    number_tss=length(tss)
    tss_bins=Dict{String,Any}()
    for (key,ts) in enumerate(tss)
        push!(tss_bins,string(key)=>ts_binify(ts))#sort into bins
    end
    tss_bins=cluster_probs(tss_bins,tss,n)
    return tss_bins
end

function ts_binify(ts)
    ts_bins=Dict{String,Any}()
    nrm_ts=normalize_ts(deepcopy(ts))
    nm_bins=10
    for b=1/nm_bins:1/nm_bins:1
        for p in eachrow(nrm_ts)
            if (p.price<=b && p.price>(b-1/nm_bins))
                if (haskey(ts_bins, string(b))==false);push!(ts_bins,string(b)=>Dict{String,Any}());push!(ts_bins[string(b)],"ts"=>DataFrame("time_stamp"=>[]));end
                push!(ts_bins[string(b)]["ts"],p[1:1])
            end
        end
    end
    return ts_bins
end


text:
@block PietroVertechi [statsmakie] begin
    @cell "StatsMakie" [Group, Style, scatter, plot, wireframe, statistics] begin
        using StatsMakie
        using DataFrames, RDatasets
        using StatsMakie: linear, smooth
        N = 1000
        a = rand(1:2, N) # a discrete variable
        b = rand(1:2, N) # a discrete variable
        x = randn(N) # a continuous variable
        y = @. x * a + 0.8*randn() # a continuous variable
        z = x .+ y # a continuous variable
        @substep

        scatter(x, y, markersize = 0.2)
        @substep

        scatter(Group(a), x, y, markersize = 0.2)
        @substep

        scatter(Group(a), x, y, color = [:black, :red], markersize = 0.2)
        @substep

        scatter(Group(marker = a), x, y, markersize = 0.2)
        @substep

        scatter(Group(marker = a, color = b), x, y, markersize = 0.2)
        @substep

        scatter(Group(marker = a), Style(color = z), x, y)
        @substep

        scatter(Group(color = a), x, y, Style(markersize = z ./ 10))
        @substep

        plot(linear, x, y)
        @substep

        plot(linear, Group(a), x, y)
        @substep

        scatter(Group(a), x, y, markersize = 0.2)
        plot!(linear, Group(a), x, y)
        @substep

        plot(linear, Group(linestyle = a), x, y)
        @substep

        N = 200
        x = 10 .* rand(N)
        a = rand(1:2, N)
        y = sin.(x) .+ 0.5 .* rand(N) .+ cos.(x) .* a
        @substep

        scatter(Group(a), x, y)
        plot!(smooth, Group(a), x, y)
        @substep

        plot(histogram, y)
        @substep

        plot(histogram, x, y)
        @substep

        plot(histogram(nbins = 30), x, y)
        @substep

        wireframe(histogram(nbins = 30), x, y)
        @substep

        iris = RDatasets.dataset("datasets", "iris")
        scatter(Data(iris), Group(:Species), :SepalLength, :SepalWidth)
        @substep
        # use Position.stack to signal that you want bars stacked vertically rather than superimposed
        plot(Position.stack, histogram, Data(iris), Group(:Species), :SepalLength)
        @substep

        wireframe(
            density(trim=true),
            Data(iris), Group(:Species), :SepalLength, :SepalWidth,
            transparency = true, linewidth = 0.1
        )
        @substep

        scatter(
            Data(iris),
            Group(marker = :Species, color = bycolumn),
            :SepalLength, (:PetalLength, :PetalWidth)
        )
        @substep

        barplot(["hi", "ima", "string"], rand(3))
        @substep

    end

    @cell "Analysis"  [Group, Style, scatter, plot, wireframe, statistics] begin

        ## setup

        using StatsMakie
        using DataFrames, RDatasets # for data
        using StatsMakie: smooth, linear

        mtcars = dataset("datasets", "mtcars")    # load dataset of car statistics
        iris = dataset("datasets", "iris")

        disallowmissing!.([mtcars, iris])  # convert columns from Union{T, Missing} to T
        # We can use this because the dataset has no missing values.

        @substep

        # kde

        plot(
            density,           # the type of analysis
            Data(mtcars),
            :MPG,
            Group(color = :Cyl)
        )

        @substep

        # histogram

        plot(
            histogram,         # the type of analysis
            Data(mtcars),
            :MPG,
            Group(color = :Cyl)
        )

        @substep

        # frequency analysis

        using Distributions
        d = rand(Poisson(), 1000)
        plot(frequency, d)

        @substep

        xs = 10 .* rand(100)
        ys = sin.(xs) .+ 0.5 .* rand.()
        scatter(xs, ys)
        plot!(smooth, xs, ys)

        @substep

        scatter(xs, ys)
        plot!(linear, xs, ys)
    end

    @cell "Box plot" [boxplot, statistics] begin

        using RDatasets, StatsMakie

        d = dataset("Ecdat","Fatality")
        boxplot(Data(d), :Year, :Perinc)

    end

    @cell "Violin plot" [violin, statistics] begin

        using RDatasets, StatsMakie

        d = dataset("Ecdat", "Fatality");

        violin(Data(d), :Year, :Perinc)

    end

    @cell "Violin and box plot" [boxplot, violin, statistics] begin

        using RDatasets, StatsMakie

        d = dataset("Ecdat","Fatality");

        violin(d.Year, d.Perinc, color = :gray)
        boxplot!(d.Year, d.Perinc, color = :black)

    end

    @cell "Errorbar" [errorbar] begin
        using StatsMakie
        x = [1:4;]
        y =  [1:4;]
        Δx = fill(0.25, 4)
        Δy = fill(0.25, 4)
        p = errorbar(x,y,Δx,Δy,xcolor=:green, ycolor=:red)
    end
end


text:
__precompile__()

module Correlate

export correlate!

using Printf
using Statistics

# fmt is ignored, because the arguments to printf must be a literal string.
# Leave the door open if that ever changes

function _print_with_space(output::Any, fmt::AbstractString, v::Float64)
  if v >= 0.0
    print(output, " ")
  end
  @printf output " %.4f" v
end

# Given a matrix of differences between the current correlation and the desired
# correlation, `cdiff`, and a matrix of tolerances, return the number of columns
# out of tolerance, as well as the largest difference.
function _count_not_converged(cdiff::Array{Float64,2}, tol::Array{Float64,2})
  (nrows,ncols) = size(cdiff)

  rc = 0;
  maxdiff = 0.0

  for c1 in 1:ncols
    for c2 in (c1+1):ncols
      d = cdiff[c1,c2]
      if d > tol[c1,c2]
        if d > maxdiff
          maxdiff = d
        end
        rc += 1
      end
    end
  end

  return (rc, maxdiff)
end

# Choose two distinct items in the range [1:n].
function choose_two(n::Int64)
  while true
    r1 = rand(1:n)
    r2 = rand(1:n)
    if r1 != r2
      return (r1, r2)
    end
  end
end

function choose_two_ordered(n::Int64)
  (r1,r2) = choose_two(n)
  if (r1 < r2)
    return (r1, r2)
  else
    return (r2, r1)
  end
end

# Assumes (but does not check) that the input DATA is zero mean within each column
# This is the function that does all the work

function correlate_zm!(data::Array{Float64, 2}, target_correlation::Array{Float64, 2},
                       tol::Array{Float64, 2}; maxiter = 0, report = 0, maxsec = 0,
                       verbose = false)

# write(stderr, "$(Base.corm(data, 0.0))\n")
  (nrows,ncols) = size(data)

  # Diagonal and cross terms of the correlation computation.
  xx = zeros(Float64, ncols)
  xy = zeros(Float64, ncols, ncols)

  for c in 1:ncols
#   write(stderr, "mean of column $(c) $(mean(data[:,c]))\n")
    for r in 1:nrows
      xx[c] += (data[r,c] * data[r,c])
    end
  end

  for c in 1:ncols
    xx[c] = sqrt(xx[c])
  end

  for c1 in 1:ncols
    for c2 in (c1+1):ncols
      for r in 1:nrows
        xy[c1,c2] += (data[r,c1] * data[r,c2])
      end
    end
  end

  for c1 in 1:ncols
    for c2 in (c1+1):ncols
      xy[c2,c1] = xy[c1,c2]
    end
  end

  mycor = zeros(Float64, ncols, ncols)

  for c1 in 1:ncols
    mycor[c1,c1] = 1.0
    for c2 in (c1+1):ncols
      mycor[c1,c2] = mycor[c2,c1] = xy[c1,c2] / (xx[c1] * xx[c2])
    end
  end

# write(stderr, "My correlation $(mycor)\n")
# write(stderr, "   correlation $(cor(data,mean=0))\n")

  cdiff = abs.(target_correlation - mycor)

  (nc,md) = _count_not_converged(cdiff, tol)

  if 0 == nc
    println(stderr, "Already converged")
    return true
  end

  iter = 0

  if 0 == maxiter
    maxiter = 100 * nrows
  end

  next_report = 0 == report ? (maxiter+maxiter) : report

  tzero = time()

  switches_done = 0

  switches_last_report = 0

  last_switch = 0

  column_switched = zeros(Int64, ncols)

  converged = false

  while iter < maxiter && ! converged
    iter += 1

    if iter > next_report
      (nc,md) = _count_not_converged(cdiff, tol)
      @printf(stderr, "%d iterations, %d switches (%d this period), %d not converged, max diff %f\n", iter, switches_done, (switches_done - switches_last_report), nc, md)
      if maxsec > 0 && time() - tzero > maxsec
        @printf(stderr, "Timeout after %f seconds\n", time() - tzero)
        break
      end
      next_report += report
      switches_last_report = switches_done
    end

    (row1,row2) = choose_two_ordered(nrows)

    col_to_switch = improves_correlations(data, row1, row2, target_correlation, xx, xy, cdiff, tol)

    if col_to_switch < 1
      if iter - last_switch > div(maxiter, 10)
        @printf(stderr, "No progress in %d steps\n", iter - last_switch)
        break
      end
      continue
    end

    column_switched[col_to_switch] += 1

    last_switch = iter

    converged = true      # until we find otherwise

    for c1 in 1:ncols
      for c2 in 1:ncols
        if c1 == c2 continue end
        if c1 == col_to_switch || c2 == col_to_switch
          xy[c1,c2] += - data[row1,c1]*data[row1,c2] - data[row2,c1]*data[row2,c2] +
                       data[row1,c1]*data[row2,c2] + data[row1,c2]*data[row2,c1]
          mycor[c1,c2] = xy[c1,c2]/(xx[c1] * xx[c2])
          cdiff[c1,c2] = abs(target_correlation[c1,c2] - mycor[c1,c2])
          if converged && cdiff[c1,c2] > tol[c1,c2]
            converged = false
          end
        elseif cdiff[c1,c2] > tol[c1,c2]
          converged = false
        end
      end
    end

    data[row1,col_to_switch],data[row2,col_to_switch] = data[row2,col_to_switch],data[row1,col_to_switch]
    switches_done += 1
  end

  if verbose
    if converged
      println(stderr, "Converged. ")
#     write(stderr, "cdiff $(cdiff)")
    end
    println(stderr, "Performed $iter iterations, made $switches_done switches")
    sum_errors = 0.0
    not_converged = 0
    for c1 in 1:ncols
      for c2 in 1:ncols
        _print_with_space(stderr, " %3f", mycor[c1,c2])
        if c2 > c1
          sum_errors += abs(target_correlation[c1,c2] - mycor[c1,c2])
          if abs(target_correlation[c1,c2] - mycor[c1,c2]) > tol[c1,c2]
#           @printf(stderr, "col %d vs col %d target %f mycor %f not converged\n", c1, c2, target_correlation[c1,c2], mycor[c1,c2])
            not_converged += 1
          end
        end
      end
      @printf(stderr, "     ")
      for c2 in 1:ncols
        _print_with_space(stderr, " %3f", target_correlation[c1,c2])
      end
      println(stderr)
    end
    println(stderr, "Sum of errors $sum_errors, $not_converged items not converged")
    for c in 1:ncols
      println(stderr, "Switched column $c $(column_switched[c]) times")
    end
    @printf(stderr, "Calculation took %.2f seconds\n", time()-tzero)
  end

  return converged
end

function improves_correlations(data::Array{Float64,2}, row1::Int64, row2::Int64, target_correlation::Array{Float64, 2},
                               xx::Array{Float64, 1}, xy::Array{Float64,2}, cdiff::Array{Float64,2},
                               tol::Array{Float64,2})
  (nrows,ncols) = size(data)

  current_diff = 0.0

  for c1 in 1:ncols
    for c2 in (c1+1):ncols
      d = cdiff[c1,c2]
      if d > tol[c1,c2]
        current_diff += d
      end
    end
  end
# @printf(stderr, "Rows %d and %d, cirrent diff %f\n", row1, row2, current_diff)

  lowest = current_diff
  col_to_switch = -1

  for c in 1:ncols     # what happens if we swap column C
    sum = 0.0
    for x1 = 1:ncols      # check all correlations involving C
      for x2 = (x1+1):ncols
        if x1 == c || x2 == c
          t = xy[x1,x2] - data[row1,x1]*data[row1,x2] - data[row2,x1]*data[row2,x2] +
                          data[row1,x2]*data[row2,x1] + data[row1,x1]*data[row2,x2]
          tc = t /(xx[x1] * xx[x2])
#         jj =  xy[x1,x2]/(xx[x1] * xx[x2])

#         @printf(stderr, "switching col %d, testing %d and %d, c %f current %f target %f\n", c, x1, x2, tc, jj, target_correlation[x1,x2])
          d = abs(tc - target_correlation[x1,x2])
          if d > tol[x1,x2]
            sum += d
          end
        else
          d = cdiff[x1,x2]
          if d > tol[x1,x2]
            sum += d
          end
        end
      end
    end
#   @printf(stderr, "switching column %d, sum %f\n", c, sum)
    if sum < lowest
      lowest = sum
      col_to_switch = c
    end
  end

# @printf(stderr, "Existing diff %f, new %f, column %d\n", current_diff, lowest, col_to_switch)

  return col_to_switch
end

function improves_correlations(data::Array{Float64,2}, row1::Int64, row2::Int64, target_correlation::Array{Float64, 2},
                               xx::Array{Float64, 1}, xy::Array{Float64,2}, cdiff::Array{Float64,2})
  (nrows,ncols) = size(data)

  current_diff = 0.0

  for c1 in 1:ncols
    for c2 in (c1+1):ncols
      current_diff += cdiff[c1,c2]
    end
  end
# @printf(stderr, "Rows %d and %d, cirrent diff %f\n", row1, row2, current_diff)

  lowest = current_diff
  col_to_switch = -1

  for c in 1:ncols     # what happens if we swap column C
    sum = 0.0
    for x1 = 1:(ncols-1)      # check all correlations involving C
      for x2 = (x1+1):ncols
        if x1 == c || x2 == c
          t = xy[x1,x2] - data[row1,x1]*data[row1,x2] - data[row2,x1]*data[row2,x2] +
                          data[row1,x2]*data[row2,x1] + data[row1,x1]*data[row2,x2]
          tc = t /(xx[x1] * xx[x2])
#         jj =  xy[x1,x2]/(xx[x1] * xx[x2])

#         @printf(stderr, "switching col %d, testing %d and %d, c %f current %f target %f\n", c, x1, x2, tc, jj, target_correlation[x1,x2])
          sum += abs(tc - target_correlation[x1,x2])
        else
          sum += cdiff[x1,x2]
        end
      end
    end
#   @printf(stderr, "switching column %d, sum %f\n", c, sum)
    if sum < lowest
      lowest = sum
      col_to_switch = c
    end
  end

# @printf(stderr, "Existing diff %f, new %f, column %d\n", current_diff, lowest, col_to_switch)

  return col_to_switch
end

# Basically convert the data to zero mean, call the underlying function, then shift back

"""
    correlate!(data::Array{Float64,2}, corm::Array{Float64,2}, tolm::Array{Float64,2};
               maxiter=100*nrows, report=0, maxsec=0, verbose=false)

Enforce a correlation structure on columns of independently derived data.\n
CORM: matrix of desired correlations between columns.
TOLM: tolerance associated with each pair.

# Arguments

* maxiter: number of iterations (def 100*nrows).
* report: report progress every <report> iterations.
* maxsec: give up after <maxsec> seconds of wall clock time.
* verbose: write computation summary to stderr.

# Examples

    data = rand(1000,4);
    corm = fill(0.8, (4,4))     # diagonals should be 1, but are not checked
    tolm = fill(0.01, (4,4))    # diagonals never examined
    correlate!(data, corm, tolm, verbose=true)
    cor(data)      # check result, should show somewhere around 0.8

A variant where tolerance is the same for all correlations

    correlate!(data::Array{Float64,2}, corm::Array{Float64,2}, tol_value::Float64;
               maxiter=100*nrows, report=0, maxsec=0, verbose=false)

A variant for two columns, where correlation and tolerance are scalar

    correlate!(data::Array{Float64,2}, cor_value::Float64, tol_value::Float64;
               maxiter=100*nrows, report=0, maxsec=0, verbose=false)

Run chol(corm) on your proposed correlation matrix (must be unit diagnonals) to see if
a solution is feasible. Correlate! does not check this, but will yield a "partial"
result if given an impossible task.
"""
function correlate!(data::Array{Float64, 2}, corm::Array{Float64, 2}, tolm::Array{Float64, 2};
                    maxiter = 0, report = 0, maxsec = 0, verbose = false)
  (nrows,ncols) = size(data)

  column_mean = [mean(data[:,c]) for c in 1:ncols]

  for c in 1:ncols
    cmean = column_mean[c]
    for r in 1:nrows
      data[r,c] -= cmean
    end
  end

  rc = correlate_zm!(data, corm, tolm, maxiter=maxiter, report=report, maxsec=maxsec, verbose = verbose)

  for c in 1:ncols
    cmean = column_mean[c]
    for r in 1:nrows
      data[r,c] += cmean
    end
  end

  return rc
end

# Version with two columns, so just one value for correlation and tolerance

function correlate!(data::Array{Float64, 2}, cor_value::Float64, tol_value::Float64;
                    maxiter = 0, report = 0, maxsec = 0, verbose = false)
  (nrows,ncols) = size(data)
  if 2 != ncols
    @printf(stderr, "Correlate.correlate:version with scalar correlation and tolerance operates only on N*2 data\n")
    return false
  end

  corm = fill(cor_value, (2,2))
  tolm = fill(tol_value, (2,2))

  return correlate!(data, corm, tolm, maxiter=maxiter, report=report,maxsec=maxsec, verbose = verbose)
end

function correlate!(data::Array{Float64, 2}, corm::Array{Float64,2}, tol_value::Float64;
                   maxiter = 0, report = 0, maxsec = 0, verbose = false)
  (nrows,ncols) = size(data)

  tolm = fill(tol_value, (nrows,ncols))

  return correlate!(data, corm, tolm, maxiter=maxiter, report=report, maxsec=maxsec, verbose = verbose)
end

end # module Correlate

max_stars_repo_name: IanAWatson/Correlate

max_issues_repo_name: IanAWatson/Correlate

max_forks_repo_name: IanAWatson/Correlate


text:
function onetick()
	lsetting[:time]+=lsetting[:tmove]
	lsetting[:t]+=1
end
function unload_chunks()
	plybx=(floor(ply.x)>>6)<<6
	plyby=(floor(ply.y)>>6)<<6
	nowtime=Dates.now().instant.periods.value
	for k in dim.cs
		if -128<=k.first.first-plybx<=128 && -128<=k.first.second-plyby<=128
			k.second.tm=0
		else
			tm=k.second.tm
			if tm==0
				k.second.tm=nowtime
			elseif nowtime-tm>60
				delete!(dim.cs,k.first)
			end
		end
	end
end


text:
function _iter_xv_yv(xv::Vector{<:DateDataFrame}, yv::Vector{<:DateDataFrame})
    @assert length(xv) == length(yv)
    return zip(xv, yv)
end

# \oplus
function ⊕(xv::Vector{<:DateDataFrame}, yv::Vector{<:DateDataFrame})
    return [x .+ y for (x, y) in _iter_xv_yv(xv, yv)]
end

# \ominus
function ⊖(xv::Vector{<:DateDataFrame}, yv::Vector{<:DateDataFrame})
    return [x .- y for (x, y) in _iter_xv_yv(xv, yv)]
end

# \otimes
function ⊗(xv::Vector{<:DateDataFrame}, yv::Vector{<:DateDataFrame})
    return [x .* y for (x, y) in _iter_xv_yv(xv, yv)]
end

# \oslash, \odiv is ⨸
function ⊘(xv::Vector{<:DateDataFrame}, yv::Vector{<:DateDataFrame})
    return [x ./ y for (x, y) in _iter_xv_yv(xv, yv)]
end

# \otimes
function ⊗(xv::Vector{<:DateDataFrame}, yv::Vector{<:Real})
    return [x .* y for (x, y) in zip(xv, yv)]
end

# \otimes
function ⊗(xv::Vector{<:DateDataFrame}, y::Real)
    return [x .* y for x in xv]
end

# \oslash
function ⊘(xv::Vector{<:DateDataFrame}, yv::Vector{<:Real})
    return [x ./ y for (x, y) in zip(xv, yv)]
end

# \oslash
function ⊘(xv::Vector{<:DateDataFrame}, y::Real)
    return [x ./ y for x in xv]
end


text:
export barrierGNCG

function  barrierGN(mc,pInv::InverseParam,pMis;rho = 10.0,epsilon = 0.1*(pInv.boundsHigh - pInv.boundsLow), indCredit=[],dumpResults::Function = dummy,out::Int=2,solveGN::Function=projPCG)
	low         = pInv.boundsLow .+ 1e-6; # this is not to get confused with the bounds of projGN.
	high        = pInv.boundsHigh .- 1e-6;
	#---------------------------------------------------------------------------
	#  Initialization.
	#---------------------------------------------------------------------------

	# logBarrierReg(m,mref,M) = logBarrier(m,mref,M,low,high,epsilon);
	logBarrierReg(m,mref,M)     = logBarrierSquared(m,mref,M,low,high,epsilon);
	pInv.regularizer 	 		= [pInv.regularizer;logBarrierReg];
	pInv.alpha 		 	 		= [pInv.alpha;rho];
	mref1 						= pInv.mref;
	mref2 			 			= zeros(size(mref1,1),size(mref1,2)+1);
	mref2[:,1:size(mref1,2)] 	= mref1;
	pInv.mref = mref2;

	mc,Dc,outerFlag,His = projGN(mc,pInv,pMis,indCredit=indCredit, dumpResults = dumpResults,out=out,solveGN=solveGN);

	if length(pInv.regularizer) > 2
		pInv.regularizer = pInv.regularizer[1:end-1];
		pInv.alpha = pInv.alpha[1:end-1];
	else # if we choose only one variable in the array, julia does not make it a scalar.
		pInv.regularizer = pInv.regularizer[1];
		pInv.alpha = pInv.alpha[1];
	end
	pInv.mref = mref1;

	return mc,Dc,outerFlag,His
end  # Optimization code
barrierGNCG = barrierGN;