text:
include("logging.jl")
include("steps.jl")

# Empty progress function
simulate(scenario::Scenario) = simulate(nothing, scenario)

# Simulation from file name
simulate(file_name::String, context = @__MODULE__) = simulate(nothing, setup(Scenario(), file_name, context))

# Simulation from file name with progress function
simulate(progress_fcn::Function, file_name::String, context = @__MODULE__) = simulate(progress_fcn, setup(Scenario(), file_name, context))

"""
    simulate

This is the primary simulation function that begins with a set up scenario and
propagates until the end of time.

"""
function simulate(progress_fcn::Union{Function,Nothing}, scenario::Scenario, needs_init::Bool = true)

    # Create the log.
    if !isempty(scenario.sim.log_file)
        println("Creating log ", scenario.sim.log_file)
        log = HDF5Logger.Log(scenario.sim.log_file)
    else
        log = nothing
    end

    E = Dict{String,Dict{String,Tuple}}() # "These are the effects that I generate."
    U = Dict{String,Dict{String,Any}}()   # "Anyone can write to these inputs. I will consume them."
    Y = Dict{String,Dict{String,Any}}()   # "I will write to these outputs. Anyone can consume them."

    # Seed the global random number generator, which we'll only use to draw any
    # necessary seeds for other RNGs. When this is used as part of a Monte-Carlo
    # test, this global seed and all seeds below will be overridden. However,
    # we're not there yet; that's Milestone 3.
    Random.seed!(scenario.sim.seed)

    # Seed the components' RNGs.
    seed(scenario.environment.rand)
    for vehicle in scenario.vehicles
        seed(vehicle.body.rand)
        for component in vehicle.components; seed(component.rand); end;
        for computer in vehicle.computers; seed(computer.board.rand); end;
    end

    # Get the complete set of sample times.
    t, ϵ = calculate_time_steps(scenario)
    nt = length(t) # Total number of steps to take

    # Record the start time.
    start_time = time_ns()

    # If anything goes wrong, we'll want to tidy up here before rethrowing the
    # error.
    k = 1
    try

        ###########
        # Startup #
        ###########

        if scenario.environment.startup != nothing
            scenario.environment.startup(t[1], scenario.environment.constants, scenario.environment.state)
        end

        for vehicle in scenario.vehicles

            if vehicle.body.startup != nothing
                vehicle.body.startup(t[1], vehicle.body.constants, vehicle.body.state)
            end
            for component in vehicle.components
                if component.startup != nothing
                    component.startup(t[1], component.constants, component.state)
                end
            end
            for computer in vehicle.computers
                if computer.board.startup != nothing
                    computer.board.startup(t[1], computer.board.constants, computer.board.state)
                end
                for software in computer.software
                    if software.startup != nothing
                        software.startup(t[1], software.constants, software.state)
                    end
                end
            end

        end

        ########
        # Init #
        ########

        # We'll initialize the state for each model and get the inputs and
        # outputs from the appropriate components.

        # TODO: Much of the below looks redundant. Make little functions for this stuff.

        # Set up input and output buses.
        for vehicle in scenario.vehicles
            U[vehicle.name] = Dict{String,Any}()
            Y[vehicle.name] = Dict{String,Any}()
            U[vehicle.name][vehicle.body.name]  = vehicle.body.inputs
            Y[vehicle.name][vehicle.body.name] = vehicle.body.outputs
            for component in vehicle.components
                U[vehicle.name][component.name]  = component.inputs
                Y[vehicle.name][component.name] = component.outputs
            end
            for computer in vehicle.computers
                U[vehicle.name][computer.name] = Dict{String,Any}()
                Y[vehicle.name][computer.name] = Dict{String,Any}()
                U[vehicle.name][computer.name][computer.board.name] = computer.board.inputs
                Y[vehicle.name][computer.name][computer.board.name] = computer.board.outputs
                for software in computer.software
                    U[vehicle.name][computer.name][software.name] = software.inputs
                    Y[vehicle.name][computer.name][software.name] = software.outputs
                end
            end
        end

        # If the states and outputs are totally initialized already, we don't
        # need to run init.
        if needs_init

            # Initialize the environment.
            if scenario.environment.init != nothing
                scenario.environment.state =
                    scenario.environment.init(t[1],
                                              scenario.environment.constants,
                                              scenario.environment.state,
                                              draw(scenario.environment.rand, :init))
            end

            # Initialize the vehicles.
            for vehicle in scenario.vehicles

                E[vehicle.name] = Dict{String,Tuple}()

                # Initialize each vehicle body before doing components.
                if vehicle.body.init != nothing
                    vehicle.body.state =
                        vehicle.body.init(t[1],
                                          vehicle.body.constants,
                                          vehicle.body.state,
                                          draw(vehicle.body.rand, :init))
                end

                # Get the body's effects.
                if vehicle.body.effects != nothing
                    E[vehicle.name][vehicle.body.name] =
                        vehicle.body.effects(t[1],
                                             vehicle.body.constants,
                                             vehicle.body.state,
                                             draw(vehicle.body.rand, :effects),
                                             vehicle.body.implicit)
                else
                    #E[vehicle.name][vehicle.body.name] = () # TODO: Do we want to always have an entry, even if it were to produce no effects? It's a hash table, so it's not clear that this is problematic in terms of speed.
                end

                # Get the effects of the environment on the body.
                if scenario.environment.effects != nothing
                    E[vehicle.name][scenario.environment.name] =
                        scenario.environment.effects(t[1],
                                                     scenario.environment.constants,
                                                     scenario.environment.state,
                                                     draw(scenario.environment.rand, :effects),
                                                     E[vehicle.name],
                                                     E)
                else
                    #E[vehicle.name][scenario.environment.name] = ()
                end

                # TODO: Both Bodies and Environments could produce outputs.

            end

            # Now initialize each vehicle's subsystems.
            for vehicle in scenario.vehicles

                for component in vehicle.components

                    if component.init != nothing
                        component.state, Y[vehicle.name][component.name] = component.init(
                            t[1],
                            component.constants,
                            component.state,
                            draw(component.rand, :init),
                            U[vehicle.name][component.name],
                            E[vehicle.name],
                            E)
                    end

                    if component.effects != nothing
                        E[vehicle.name][component.name] = component.effects(
                            t[1],
                            component.constants,
                            component.state,
                            draw(component.rand, :effects),
                            component.implicit,
                            U[vehicle.name][component.name],
                            E[vehicle.name],
                            E) # TODO: Limit this to the body and environment effects.
                    end

                    # This counts as a "tick" of the system if this would have
                    # been an appropriate time for a step anyway.
                    if t[k] >= component.timing.t_next - ϵ
                        component.timing.count = 1
                        component.timing.t_next = component.timing.dt + component.timing.t_start
                    end

                end # for each component

                for computer in vehicle.computers

                    if computer.board.init != nothing
                        computer.board.state,
                        Y[vehicle.name][computer.name][computer.board.name] =
                            computer.board.init(
                                t[1],
                                computer.board.constants,
                                computer.board.state,
                                draw(computer.board.rand, :init),
                                U[vehicle.name][computer.name][computer.board.name],
                                E[vehicle.name],
                                E)
                    end

                    if computer.board.effects != nothing
                        E[vehicle.name][computer.name][computer.board.name] =
                            computer.board.effects(
                                t[1],
                                computer.board.constants,
                                computer.board.state,
                                draw(computer.board.rand, :effects),
                                computer.board.implicit,
                                U[vehicle.name][computer.name][computer.board.name],
                                E[vehicle.name],
                                E) # TODO: Limit this to the body and environment effects.
                    end

                    # This counts as a "tick" of the system if this would have
                    # been an appropriate time for a step anyway.
                    if t[k] >= computer.board.timing.t_next - ϵ
                        computer.board.timing.count = 1
                        computer.board.timing.t_next = computer.board.timing.dt + computer.board.timing.t_start
                    end

                    for software in computer.software

                        if software.init != nothing
                            software.state,
                            Y[vehicle.name][computer.name][software.name],
                            U[vehicle.name] =
                                software.init(
                                    t[1],
                                    software.constants,
                                    software.state,
                                    U[vehicle.name][computer.name][software.name],
                                    Y[vehicle.name],
                                    U[vehicle.name])
                        end

                        # This counts as a "tick" of the system if this would have
                        # been an appropriate time for a step anyway.
                        if t[k] >= software.timing.t_next - ϵ
                            software.timing.count = 1
                            software.timing.t_next = software.timing.dt + software.timing.t_start
                        end

                    end

                end # for each computer

            end # for each vehicle

        end # if needs_init

        # We now have the state of, inputs to, and outputs from all of pieces of
        # the simulation for the first sample. We also have all of the effects
        # resulting from the current state.

        ############
        # Log Init #
        ############

        # TODO: So that a sim can be recreated from a log, shouldn't we log the
        # sim parameters too? How do we log t_next/t_start?

        # Preallocate space for the streams we'll need in the logger.
        if log != nothing

            # TODO: Log environment.

            for vehicle in scenario.vehicles

                ## Set up the streams for the vehicle's truth.
                slug = "/" * vehicle.name * "/" * vehicle.body.name * "/state/"
                add!(log, slug, t[1], vehicle.body.state, nt)

                # Set up the streams for the components.
                for component in vehicle.components

                    # State
                    if component.state != nothing
                        slug = "/" * vehicle.name * "/" * component.name * "/state/"
                        if component.derivatives != nothing
                            num_samples = nt
                        else
                            num_samples = Int64(floor((scenario.sim.t_end - component.timing.t_start) / component.timing.dt)) + 1
                        end
                        add!(log, slug, t[1], component.state, num_samples)
                    end

                    # TODO: Inputs

                    # Outputs
                    if component.outputs != nothing
                        slug = "/" * vehicle.name * "/" * component.name * "/outputs/"
                        num_samples = Int64(floor((scenario.sim.t_end - component.timing.t_start) / component.timing.dt)) + 1
                        add!(log, slug, t[1], Y[vehicle.name][component.name], num_samples)
                    end

                end

                for computer in vehicle.computers

                    # TODO: State

                    # TODO: Inputs

                    # Outputs
                    if computer.board.outputs != nothing
                        slug = "/" * vehicle.name * "/" * computer.name * "/" * computer.board.name * "/outputs/"
                        num_samples = Int64(floor((scenario.sim.t_end - component.board.timing.t_start) / component.board.timing.dt)) + 1
                        add!(log, slug, t[1], Y[vehicle.name][computer.name][computer.board.name], num_samples)
                    end

                    for software in computer.software

                        # TODO: State

                        # TODO: Inputs

                        # Outputs
                        if software.outputs != nothing
                            slug = "/" * vehicle.name * "/" * computer.name * "/" * software.name * "/outputs/"
                            num_samples = Int64(floor((scenario.sim.t_end - software.timing.t_start) / software.timing.dt)) + 1
                            add!(log, slug, t[1], Y[vehicle.name][computer.name][software.name], num_samples)
                        end

                    end

                end

            end

        end

        ############
        # Sim Loop #
        ############

        # Loop over each index of the time array.
        for outer k = 1:nt

            # Update the progress bar.
            if progress_fcn != nothing
                if progress_fcn(k, nt) == false
                    break;
                end
            end

            # Put all of the continuous-time states, draws, and implicit variables
            # into vectors.
            X = []
            D = []
            V = []
            # TODO: Add environments; they should be able to have continuous states.
            for vehicle in scenario.vehicles
                push!(X, vehicle.body.state)
                push!(D,  draw(vehicle.body.rand, :derivatives))
                push!(V, vehicle.body.implicit)
            end
            for vehicle in scenario.vehicles
                for component in vehicle.components
                    push!(X, component.state)
                    push!(D, draw(component.rand, :derivatives))
                    push!(V, component.implicit)
                end
                for computer in vehicle.computers
                    push!(X, computer.board.state)
                    push!(D, draw(computer.board.rand, :derivatives))
                    push!(V, computer.board.implicit)
                end
            end

            ##############
            # Dimensions #
            ##############

            # On the first sample, we discover some dimensions.
            if k == 1

                # We'll count up the number of states and the number of elements in the state
                # derivative, as well as the mapping between them.
                nx  = 0
                nxd = 0
                x_to_xd = Vector{Int64}()

                # Make a dummy call to the derivatives.
                Xd = ode(t, X, V, D, U, scenario)

                # For each element of the state and its derivative...
                for c = 1:length(X)
                    nxc  = length(stackem(X[c]))  # Determine how many elements it has
                    nxdc = length(stackem(Xd[c])) # And the number of elements in the derivative.
                    if nxdc == 0 # If this state is discrete only, add a bunch of zero indices.
                        append!(x_to_xd, zeros(Int64, nxc))
                    elseif nxdc == nxc # Otherwise, map from state element to derivative elements.
                        append!(x_to_xd, nxd .+ (1:nxdc))
                    else # And if the numbers of those things don't match, it's a problem.
                        error("The number of derivatives don't match the number of states for the ", X[c], " type.")
                    end
                    nx  += nxc  # Count 'em up.
                    nxd += nxdc
                end

                # Store the results.
                scenario.dims = ScenarioDimensions(nt, nx, nxd, length(stackem(V)), x_to_xd)

            end
# display(X)
            #############
            # Propagate #
            #############

            # Propagate from k-1 to k. We skip the first sample and just log it.
            if k > 1

                #####################
                # Continuous Update #
                #####################

                # Use RK4 with a constraint solver for semi-explicity, index 1 DAE support.
                dt = t[k] - t[k-1] # Time step (s)
                Xd1, V = minor(t[k-1],         X,               V, D, U, scenario)
                Xd2, V = minor(t[k-1] + 0.5dt, X + 0.5dt * Xd1, V, D, U, scenario)
                Xd3, V = minor(t[k-1] + 0.5dt, X + 0.5dt * Xd2, V, D, U, scenario)
                Xd4, V = minor(t[k-1] +    dt, X +    dt * Xd3, V, D, U, scenario)
                X = X + dt/6. * Xd1 + dt/3. * Xd2 + dt/3. * Xd3 + dt/6. * Xd4

                # Put back all of the states.
                # TODO: Environments?
                # TODO: Why bother storing these this way?
                c = 0
                for vehicle in scenario.vehicles
                    c += 1
                    vehicle.body.state    = X[c]
                    vehicle.body.implicit = V[c]
                end
                for vehicle in scenario.vehicles
                    for component in vehicle.components
                        c += 1
                        component.state    = X[c]
                        component.implicit = V[c]
                    end
                    for computer in vehicle.computers
                        c += 1
                        computer.board.state    = X[c]
                        computer.board.implicit = V[c]
                    end
                end

                # Get the effects on the current state. Note that this may
                # actually be different from the effects we'll need at the
                # beginning of the next propagation, because the discrete
                # updates below can change their states (and hense resulting
                # effects) instantly.
                E = get_effects(t, X, V, D, U, scenario)

                ###################
                # Discrete Update #
                ###################

                # Update components and computer boards.
                for vehicle in scenario.vehicles

                    # TODO: Might bodies have discrete updates?

                    # Create an iterator over all of the effects for this vehicle for
                    # convenience and speed.
                    vehicle_effects = copy(E[vehicle.name])

                    # Update any sensor whose t_next has come.
                    for component in vehicle.components
                        if component.update != nothing && t[k] >= component.timing.t_next - ϵ

                            # Update the state and get the measurement.
                            component.state,
                            Y[vehicle.name][component.name] = component.update(
                                t[k],
                                component.constants,
                                component.state,
                                draw(component.rand, :update),
                                U[vehicle.name][component.name],
                                vehicle_effects,
                                E)

                            # Update the next hit time.
                            component.timing.count += 1
                            component.timing.t_next = component.timing.dt * component.timing.count + component.timing.t_start

                            # TODO: Log my inputs, resulting state, my outputs. Draws?
                            # Skip the state if there's a derivative function; it will be logged below.

                        end
                    end

                    # Update any sensor whose t_next has come.
                    for computer in vehicle.computers
                        if computer.board.update != nothing && t[k] >= computer.timing.t_next - ϵ

                            # Update the state and get the measurement.
                            computer.board.state,
                            Y[vehicle.name][computer.name][computer.board.name] = computer.board.update(
                                t[k],
                                computer.board.constants,
                                computer.board.state,
                                draw(computer.board.rand, :update),
                                U[vehicle.name][computer.name][computer.board.name],
                                vehicle_effects,
                                E)

                            # Update the next hit time.
                            computer.board.timing.count += 1
                            computer.board.timing.t_next = computer.board.timing.dt * computer.board.timing.count + computer.board.timing.t_start

                            # TODO: Log my inputs, resulting state, my outputs. Draws?
                            # Skip the state if there's a derivative function; it will be logged below.

                        end
                    end

                end

                ###################
                # Software Update #
                ###################

                # Update sensors, software, and actuators.
                for vehicle in scenario.vehicles

                    vehicle_outputs_km1  = deepcopy(Y[vehicle.name])
                    vehicle_outputs_temp = deepcopy(vehicle_outputs_km1)

                    # If computers were components and all components had software
                    # to run...
                    for computer in vehicle.computers # Filter on "active" components?

                        # Each computer can overwrite outputs for each software.

                        # Run each software, allowing subsequent processes to consume
                        # inputs/outputs from prior processes on this component.
                        for software in computer.software
                            if software.update != nothing && t[k] >= software.timing.t_next

                                # Any software of this computer can consume
                                # inputs/outputs from the prior components. They cannot
                                # see outputs from other components for this sample.

                                # Update the inputs bus and outputs for this software.
                                software.state,
                                vehicle_outputs_temp[computer.name][software.name], # This software's outputs
                                U[vehicle.name] = # Commands or software inputs for any other computer on this vehicle
                                    software.update(t[k],
                                                software.constants,
                                                software.state,
                                                U[vehicle.name][computer.name][software.name], # Commands or software inputs from this vehicle
                                                vehicle_outputs_temp, # Includes measurements and other software outputs
                                                U[vehicle.name]) # To write to the inputs of other components

                                # Update the next sample time.
                                software.timing.count  += 1
                                software.timing.t_next = software.timing.dt * software.timing.count + software.timing.t_start

                                # Copy the updated outputs to the complete set of
                                # updated outputs that won't be used again until the
                                # next sample.
                                Y[vehicle.name][computer.name][software.name] = vehicle_outputs_temp[computer.name][software.name]

                                # TODO: Log my inputs, resulting state, my outputs.

                            end
                        end # software

                        # Now that we're done with this component, return any of the
                        # modified outputs to their previous state.
                        vehicle_outputs_temp[computer.name] = deepcopy(vehicle_outputs_km1[computer.name])

                    end # computer

                    # The outputs bus has now been completely updated for this vehicle.

                end # for each vehicle

            end # k > 1

            ###########
            # Logging #
            ###########

            # Log the continuous states.
            if log != nothing
                for vehicle in scenario.vehicles
                    if vehicle.body.derivatives != nothing
                        slug = "/" * vehicle.name * "/" * vehicle.body.name * "/state/"
                        log!(log, slug, t[k], vehicle.body.state)
                    end
                    for component in vehicle.components
                        if component.derivatives != nothing
                            slug = "/" * vehicle.name * "/" * component.name * "/state/"
                            log!(log, slug, t[k], component.state)
                        end
                    end
                    for computer in vehicle.computers
                        if computer.board.derivatives != nothing
                            slug = "/" * vehicle.name * "/" * computer.name * "/" * computer.board.name * "/state/"
                            log!(log, slug, t[k], computer.board.state)
                        end
                    end
                end
            end

        end # sim loop

    catch err

        # TODO: Tidy up logs (get rid of unused preallocated space).

        # The sim stops upon receiving ctrl+c (ctrl+shift+c in Juno). That's not
        # really an error. All other errors really are errors though.
        if isa(err, InterruptException)
            println("Stopping sim due to user interrupt.")
            if k > 1; k -= 1; end # The last step completed was k-1.
        else
            rethrow(err)
        end

    # No matter what happens, we want to shut down when the sim has stopped
    # (either completed or errored out) so that we can tidy up resources.
    finally

        ############
        # Shutdown #
        ############

        # Shut down each vehicle first.
        for vehicle in scenario.vehicles

            if vehicle.body.shutdown != nothing
                vehicle.body.shutdown(t[k], vehicle.body.constants, vehicle.body.state)
            end
            for component in vehicle.components
                if component.shutdown != nothing
                    component.shutdown(t[k], component.constants, component.state)
                end
            end
            for computer in vehicle.computers
                if computer.board.shutdown != nothing
                    computer.board.shutdown(t[k], computer.board.constants, computer.board.state)
                end
                for software in computer.software
                    if software.shutdown != nothing
                        software.shutdown(t[k], software.constants, software.state)
                    end
                end
            end

        end

        # Close the log.
        if log != nothing
            close!(log)
        end

    end

    # Report how we did.
    run_time  = (time_ns() - start_time) * 1e-9
    println("Simulated ", t[k], " seconds in ", run_time, " seconds (", t[k]/run_time, "× real time).")

    return scenario
end

# Figure out every time step that we're going to need.
function calculate_time_steps(scenario)

    tf = scenario.sim.t_end
    dt = scenario.sim.dt

    # Timing slop. Discrete processes will occur when they're within this amount of
    # time from when they're supposed to happen.
    ϵ = eps(10 * scenario.sim.t_end)

    # Get all of the discrete time steps and start times.
    dts      = Vector{Float64}()
    t_starts = Vector{Float64}()
    for vehicle in scenario.vehicles
        # TODO: Bodies might have discrete steps.
        for component in Iterators.filter(c -> c.timing.dt != 0., vehicle.components)
            push!(dts, component.timing.dt)
            push!(t_starts, component.timing.t_start)
            component.timing.t_next = component.timing.t_start
        end
        for computer in vehicle.computers
            if computer.board.timing.dt != 0.
                push!(dts, computer.board.timing.dt)
                push!(t_starts, computer.board.timing.t_start)
                computer.board.timing.t_next = computer.board.timing.t_start
            end
            for software in Iterators.filter(s -> s.timing.dt != 0., computer.software)
                push!(dts, software.timing.dt)
                push!(t_starts, software.timing.t_start)
                software.timing.t_next = software.timing.t_start
            end
        end
    end

    # Make the time history.
    if !isempty(dts)

        t      = 0.
        ts     = Vector{Float64}() # array for time steps
        ndst   = copy(t_starts) # next discrete sample times
        counts = zeros(Float64, length(dts)) # how many times each discrete process has triggered

        # We can't need more than the sum of each process triggering on its own.
        # Create a counter so that, when it completes, we bail with an error.
        count = Int64(tf / dt + sum(tf ./ dts))

        # While we aren't already at the end of time...
        while t < tf

            # See if there's been a problem that's creating an infinite loop.
            count -= 1
            if count == 0
                error("Maximum number of steps exhausted; something went wrong when calculating time steps.")
            end

            # Find the next discrete sample time.
            tn = minimum(ndst)
            tn = min(tn, tf) # Don't step past the end of time.

            # See how many steps we'll need to take between now and then. Create
            # evenly-spaced steps between the points.
            Δt      = tn - t
            n_steps = ceil((Δt-ϵ)/dt)
            Δt_fit  = Δt / n_steps

            # Add on the evenly-spaced steps, up to the next discrete step.
            for k = 1:n_steps-1
                t += Δt_fit
                push!(ts, t) # Pushing each is better than vcating (no copies).
            end

            # Now add on the discrete step. (Doing this allows the step to occur
            # right at count * dt and aNothings the accumulation of roundoff error over
            # time.
            t = tn
            push!(ts, t)

            # Update the times for the next steps of each discrete process.
            for k = 1:length(dts)
                if t > ndst[k] - ϵ
                    counts[k] += 1
                    ndst[k] = counts[k] * dts[k] + t_starts[k]
                end
            end

        end

    else

        # Otherwise, we have no discrete times. Step at the maximum rate.
        t = collect(0.:scenario.sim.dt:scenario.sim.t_end)
        if t[end] != scenario.sim.t_end # Stop at exactly the end time.
            push!(t, senario.sim.t_end)
        end

    end

    return (ts, ϵ, counts)

end


text:
using NeRCA
using ROy
using Sockets

calib = NeRCA.read_calibration("scripts/NeRCA_00000043_20022019.detx")
calib = NeRCA.read_calibration("/home/tgal/data/detx/NeRCA_-00000001_20171212.detx")

client = NeRCA.CHClient(ip"127.0.0.1", 5553, ["IO_EVT"])

try
    for message in client
        event = NeRCA.read_io(IOBuffer(message.data), NeRCA.DAQEvent)
        hits = event.snapshot_hits
        chits = NeRCA.calibrate(hits, calib)
        thits = event.triggered_hits
        cthits = NeRCA.calibrate(thits, calib)
        println(length(cthits))
    end
catch InterruptException
    println("ok")
end

ext: jl


text:
module ImageBinarization

using Base.Iterators: repeated
using MappedArrays
using LinearAlgebra
using Polynomials
using Statistics

using ImageContrastAdjustment
using HistogramThresholding

using ImageCore
using ImageCore: GenericGrayImage
using ColorVectorSpace

# TODO: port BinarizationAPI to ImagesAPI
include("BinarizationAPI/BinarizationAPI.jl")
import .BinarizationAPI: AbstractImageBinarizationAlgorithm,
                         binarize, binarize!

include("integral_image.jl")
include("util.jl")
include("compat.jl")

# Concrete binarization algorithms

include("algorithms/adaptive_threshold.jl") # AdaptiveThreshold
include("algorithms/balanced.jl") # Balanced
include("algorithms/entropy.jl") # Entropy
include("algorithms/intermodes.jl") # Intermodes
include("algorithms/minimum.jl") # MinimumIntermodes
include("algorithms/minimum_error.jl") # MinimumError
include("algorithms/moments.jl") # Moments
include("algorithms/niblack.jl") # Niblack
include("algorithms/otsu.jl") # Otsu
include("algorithms/polysegment.jl") # Polysegment
include("algorithms/sauvola.jl") # Sauvola
include("algorithms/unimodal.jl") # UnimodalRosin
include("algorithms/yen.jl") # Yen

include("deprecations.jl")

export
    # generic API
    binarize, binarize!,

    # Algorithms
    AdaptiveThreshold, recommend_size,
    Balanced,
    Entropy,
    Intermodes,
    MinimumError,
    MinimumIntermodes,
    Moments,
    Niblack,
    Otsu,
    Polysegment,
    Sauvola,
    UnimodalRosin,
    Yen

end # module ImageBinarization


text:
# * higher order? rdiff? ctaylor? audi?
# * moulla2
# * test cases: fluid, beam
# * Hamiltonian function -> use forwarddiff for gradient and Hessian?
# * automatic coupling phs?
# * better way to simulate?
# * automatic subdivide in elements?
__precompile__()

"""
	PortHamiltonian

	This module has several tools for discretization of infinite-dimensional systems
	and manipulation of finite-dimensional port-Hamiltonian systems.

"""
module PortHamiltonian

using ForwardDiff
using LinearAlgebra

###########
# imports #
###########
# additional functionality is added to the following functions:
import LinearAlgebra: eigen
#import SparseArrays: blkdiag
import Base: show

###########
# exports #
###########
export set_constraint!,
	   removeport!,
	   frequencies,
	   constraint_elimination,
	   coupled_gyrator,
	   discrete_phs,
	   discrete_phs_closed,
	   discrete_phs2,
	   discrete_phs2_distports,
	   finelem,
	   damp,
	   blkdiag

"""
    blkdiag(x :: Array, y :: Array)

    creates a block-diagonal matrix from two arrays x and y
"""
function blkdiag(x :: Array, y :: Array)
	res = zeros(size(x,1)+size(y,1),size(x,2)+size(y,2));
	res[1:size(x,1), 1:size(x,2)] = x;
	res[1+size(x,1):end, 1+size(x,2):end] = y;
	return res;
end

include("types.jl")
"""
	blkdiag(x :: Phs, y :: Phs)

	concatenate two phs objects
"""
function blkdiag(x :: Phs, y :: Phs)
	warn("only linear phs with constraint can be concatenated yet")
	phnew = Phs(blkdiag(x.J,y.J),blkdiag(x.B,y.B),blkdiag(x.D,y.D), blkdiag(x.Q,y.Q))
	phnew.TransfMatrix = blkdiag(x.TransfMatrix,y.TransfMatrix)
	phnew.R = blkdiag(x.R,y.R)
	Nx = size(x.TransfMatrix,1)
	for key in keys(y.StatesNames)
		y.StatesNames[key] += Nx
	end
	phnew.StatesNames = merge(x.StatesNames, y.StatesNames)
	phnew.InputsNames = [x.InputsNames; y.InputsNames]
	return phnew
end

"""
	removeport!(p :: Phs, portnumber :: Integer)

	remover input/output port of Phs p
"""
function removeport!(p :: Phs, portnumber :: Integer)
	keep = 1:size(p.B,2) .!= portnumber
	p.B = p.B[:,keep]
	p.D = p.D[keep,keep]
	return p
end
include("lgwt.jl")
include("lglnodes.jl")
include("derivative.jl")
include("massmatrix.jl")
include("phdiscretization.jl")
include("phdiscretization_closed.jl")
include("eigenphs.jl")
include("constraint_elimination.jl")
include("spectral_element_phs.jl")
include("coupled_phs.jl")
include("weakdiscretization.jl")

function leg_pol(x :: Array,xi :: Array,j :: Integer)
	leg_pol(x[1],xi,j)
end
"""
	leg_pol(x :: Number,xi :: Array,j :: Integer)

	evaluate the "j" Lagrange polynomial defined using
	the collocation points "xi" at point "x"
"""
function leg_pol(x :: Number,xi :: Array,j :: Integer)

	P = 1.;
        for k = 1:length(xi)
	        if (k != j)
		        P = P*(x[1]-xi[k])/(xi[j]-xi[k])
				#println((x[1]-xi[k])/(xi[j]-xi[k]))
            end
        end
    return P
end

end

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
cursor_deleteline(buf::IO=terminal.out_stream)               = print(buf, "\x1b[2K")

cursor_hide(buf::IO=terminal.out_stream)                     = print(buf, "\x1b[?25l")
cursor_show(buf::IO=terminal.out_stream)                     = print(buf, "\x1b[?25h")

cursor_save_position(buf::IO=terminal.out_stream)            = print(buf, "\x1b[s")
cursor_restore_position(buf::IO=terminal.out_stream)         = print(buf, "\x1b[u")

cursor_move_abs(buf::IO, c::Vector{Int}=[0,0])  = print(buf, "\x1b[$(c[2]);$(c[1])H")
cursor_move_abs(c::Vector{Int}) = cursor_move_abs(terminal.out_stream, c)

# ToDo: Remove x offset after newline
function cursor_move_rel(buf::IO, c=[0,0])
    x = c[1] >= 0 ? "\x1b[$(abs(c[1]))A" : "\x1b[$(abs(c[1]))B"
    y = c[2] >= 0 ? "\x1b[$(abs(c[2]))C" : "\x1b[$(abs(c[2]))D"
    print(buf, x,y)
end
cursor_move_rel(c::Vector{Int}) = cursor_move_rel(terminal.out_stream, c)

function clear_screen()
    buf = IOBuffer()
    cursor_move_abs(buf, [0,0])
    for y in 1:50
        cursor_move_abs(buf, [0,y])
        cursor_deleteline(buf)
    end
    print(String(take!(buf)))
end

"""
    put(pos::Vector, s::String)
Put text `s` on screen at coordinates `pos`.
Does not change cursor position.
"""
function put(buf::IO, pos::Vector, s::String)
    cursor_save_position(buf)
    cursor_move_abs(buf, pos)
    print(buf, s)
    cursor_restore_position(buf)
end
put(pos::Vector, s::String) = put(terminal.out_stream, pos::Vector, s::String)

function put(buf::IO, pos::Vector, color::Crayon, s::String)
    print(buf, color)
    put(buf, pos, s)
end
put(buf::IO, pos::Vector, c::Symbol, s::String) = put(buf, pos, Crayon(foreground=c), s)
put(pos::Vector, c::Union{Symbol, Crayon}, s::String) = put(terminal.out_stream, pos, c, s)

function terminal_screen(data::Array{String}, origin::Vector{Int}=[1,1];
                colors::Array{Symbol}=fill(:white, size(data)))
    buf = IOBuffer()
    x1, y1 = origin
    dy, dx = size(data)
    row_heigths, column_widths = _row_column_sizes(data)
    for x in 1:dx, y in 1:dy
        d = data[y, x]
        s = String.(split(d, "\n"))
        for i in 1:row_heigths[y]
            line = i > length(s) ? "" : s[i]
            put(buf, [sum(column_widths[1:x-1])+3x+x1,sum(row_heigths[1:y-1])+y+y1+i], colors[y,x], line)
        end
    end
    print(String(take!(buf)))
end

function _row_column_sizes(data::Array{String})
    size_string(s::String) = length(split(s, "\n")), maximum(length.(split(s, "\n")))
    size_data = size_string.(data)
    row_heigths = maximum(getindex.(size_data, 1), 2)
    column_widths = maximum(getindex.(size_data,2), 1)
    row_heigths, column_widths
end


text:
const AttributeType = Union{Symbol, String}
const DbIndexType = Union{String, Integer}
const DbIndexArrayType = Union{Array{String}, Array{Int}}


text:
@testset "hamiltonian.jl" begin

    ## Old implementation
    using SparseArrays
    function xxzmodel_old(N, J, Δ)
        Sx = complex(σx/2); Sy = σy./2; Sz = complex(σz/2);
        res = spzeros(complex(eltype(J)), 2^N, 2^N)
        for i in 1:N
            for j in i+1:N
                res += J[i, j] * (correlator(Sx, i, j, N) + correlator(Sy, i, j, N) + Δ*correlator(Sz, i, j, N))
            end
        end
        res
    end

    @testset "fieldterm" begin
        @test fieldterm(2, σx) ≈ σx ⊗ 𝟙 + 𝟙 ⊗ σx
        @test fieldterm(5, 𝟙) ≈ 5*identity_op(5)
    end

    @testset "hopping operator" begin
        @test hopping(2,1,2) ≈ [0 0 0 0;
                                0 0 1 0;
                                0 1 0 0;
                                0 0 0 0]
        @test hopping(5,4,1) == hopping(5,1,4)
    end

    @testset "xxzmodel" begin
        @test xxzmodel(2, 1, -1)*4  ≈ [-1  0  0  0;
                                        0  1  2  0;
                                        0  2  1  0;
                                        0  0  0 -1]
        @test xxzmodel([0 1 1; 1 0 1; 1 1 0], -1) ≈ xxzmodel(3, 1, -1)

        # comparison to old implementation
        for N in 2:14
            J = Symmetric(rand(N,N))
            @test xxzmodel(J,-0.7) ≈ xxzmodel_old(N, J, -0.7)
        end
    end

    @testset "z_field" begin
        @test diag(z_field(3)) == [3,1,1,-1,1,-1,-1,-3]
        @test diag(z_field([1,2,3])) == [6,0,2,-4,4,-2,0,-6]
        @test abs(sum(z_field(rand(10)))) < 1e-13 # can be suprisingly large with ~1e-14
    end
end


text:
geneLen  = Dict{String, Int64}()

f = open( "hg19_gencode_basic.gtf", "r" )
n = countlines( f )  # 计算数据行数
seekstart( f )       # 返回文件起始位置

for i = 1:n
    s = split( readline(f), "\t")
    a_start = parse(Int64, s[4])
    a_end = parse(Int64, s[5])

    gene_Info = s[9]
    geneID_pre = split(gene_Info, ";")[1]
    geneID = replace(geneID_pre, "gene_id " => "")
    geneID_a = strip(geneID, '"')

    if !haskey(geneLen, geneID_a)
        geneLen[geneID_a] = a_end - a_start
    else
        geneLen[geneID_a] += a_end - a_start
    end

end

for (k,v) in geneLen

    println("$k length is $v")
end


text:
################################################################################
# setup
################################################################################

using DataBench, FastGroupBy, ShortStrings, CategoricalArrays, SortingLab,
    SortingAlgorithms, BenchmarkTools, DataFrames

using Statistics: mean

using Random
Random.seed!(1);
const N = 10_000_000; const K = 100

@time df = createSynDataFrame(N, K); #31 #40

# using JuliaDB, IterableTables
# @time dft = table(df)

# if isfile("df$(N÷1_000_000)m.feather")
#     @time df = Feather.read("df$(N÷1_000_000)m.feather")
# else
#     @time df = createSynDataFrame(N, K); #31 #40
#     @time Feather.write("df$(N÷1_000_000)m.feather", df) # save for future use
# end
# ZJ: It's actually slower to load from feather

################################################################################
# DT[, sum(v1), keyby=id1]
# short string & string
# Status: SLOWER but if converted to Categorical then it's faster
################################################################################
# test String15
# @time df[!, :id1_ss] = ShortString7.(df[!, :id1]);
# @time a = fastby(sum, df[!, :id1_ss], df[!, :v1]);
# @time a = fastby(sum, df[!,:id1_ss], df[!,:v1]);
#
# @time a = fastby(sum, df[!,:id1], df[!,:v1]);
# @time a = fastby(sum, df[!,:id1], df[!,:v1]);

################################################################################
# DT[, sum(v1), keyby=id1]
# test categorical
# Status: FASTER 9x
################################################################################
@time df[!,:id1_cate] = categorical(df[!,:id1]); #7
@time df[!,:id1_cate] = compress(df[!,:id1_cate]); # 0.5
# @time sumby(df[!,:id1_cate], df[!,:v1]);
@time fastby(sum, df[!,:id1_cate], df[!,:v1]);
@time fastby(sum, df[!,:id1_cate], df[!,:v1]);

# @time fgroupreduce(+, df[!,:id1_cate], df[!,:v1]);
# @time fgroupreduce(+, df[!,:id1_cate], df[!,:v1]);

################################################################################
# DT[, sum(v1), keyby="id1,id2"]
# test categorical
# Status: FASTER 9x
################################################################################
@time df[!,:id2_cate] = categorical(df[!,:id2]); #7
@time df[!,:id2_cate] = compress(df[!,:id2_cate]); # 0.5
# @time fgroupreduce(+, (df[!,:id1_cate], df[!,:id2_cate]), df[!,:v1]);
# @time fgroupreduce(+, (df[!,:id1_cate], df[!,:id2_cate]), df[!,:v1]);

################################################################################
# DT[, list(sum(v1),mean(v3)), keyby=id3]
# test categorical
# Status: SLOWER
# TODO: make a fast fsortandperm for integers
# TODO: make this faster than data.table
# TODO: look for mean and make it into a reduce
# TO BEAT: 0.6s for 10m; 14.11s for 100m;
#
# NOTES:
# The key issue here is that we need to sort two vectors in order to perform the
# groupby operation. R's seems to have a really fast way to doing that
# Potential make the group into radix sort instead of counting sort
#
# What about making an RLE vector type?
#
# I can radixsort the 10m refs in 0.08 seconds so still 0.52 seconds to distribute
#
# I can sortperm the 10m refs in 0.28 seconds so still
################################################################################
@time df[!,:id3_cate] = categorical(df[!,:id3]) |> compress;
# import FastGroupBy: BaseRadixSortSafeTypes, fastby

# @benchmark fastby((sum, mean), df[!,:id3_cate], (df[!,:v1], df[!,:v3]))
# BenchmarkTools.Trial:
#   memory estimate:  435.87 MiB
#   allocs estimate:  1000315
#   --------------
#   minimum time:     1.084 s (0.00% GC)
#   median time:      1.130 s (0.00% GC)  mean time:        1.133 s (0.00% GC)
#   maximum time:     1.170 s (0.00% GC)
#   --------------
#   samples:          5
#   evals/sample:     1

# BenchmarkTools.Trial:
#   memory estimate:  430.21 MiB
#   allocs estimate:  977648
#   --------------
#   minimum time:     801.169 ms (0.00% GC)
#   median time:      912.131 ms (0.00% GC)
#   mean time:        945.560 ms (0.00% GC)
#   maximum time:     1.200 s (0.00% GC)
#   --------------
#   samples:          6
#   evals/sample:     1

# alternatives
Base.zero(t::Type{T}) where T <: Tuple = ((zero(p) for p in t.parameters)...)

function fastby2(fns, byvec, valvec)
    xx = zip(valvec...) |> collect
    a = fgroupreduce((x,y)->(x[1]+y[1], x[2]+y[2], x[3] + 1), df[!,:id3_cate], xx, (0,0.0,0))
end

function fastby3(fns, byvec, valvec)
    (fgroupreduce(+, byvec, v) for (fn, v) in zip(fns, valvec)) |> collect
end

@time fastby3((sum, mean), df[!,:id3_cate], (df[!,:v1], df[!,:v3]));
# BenchmarkTools.Trial:
#   memory estimate:  885.54 MiB
#   allocs estimate:  10000828
#   --------------
#   minimum time:     25.588 s (48.27% GC)
#   median time:      26.755 s (50.79% GC)
#   mean time:        26.430 s (50.20% GC)
#   maximum time:     27.005 s (51.47% GC)
#   --------------
#   samples:          5
#   evals/sample:     1
@time fastby2((sum, mean), df[!,:id3_cate], (df[!,:v1], df[!,:v3]));
# BenchmarkTools.Trial:
#   memory estimate:  1.95 GiB  allocs estimate:  5000410  --------------
#   minimum time:     14.963 s (37.27% GC)
#   median time:      18.198 s (48.77% GC)
#   mean time:        18.539 s (48.99% GC)
#   maximum time:     21.212 s (53.51% GC)
#   --------------
#   samples:          5
#   evals/sample:     1

@benchmark fastby3((sum, mean), df[!,:id3_cate], (df[!,:v1], df[!,:v3])) samples=5 evals=1 seconds=120
# BenchmarkTools.Trial:
#   memory estimate:  100.84 MiB
#   allocs estimate:  1000671
#   --------------
#   minimum time:     640.946 ms (0.00% GC)
#   median time:      661.054 ms (0.00% GC)
#   mean time:        691.247 ms (0.00% GC)
#   maximum time:     762.335 ms (0.00% GC)
#   --------------
#   samples:          5
#   evals/sample:     1

@benchmark fastby2((sum, mean), df[!,:id3_cate], (df[!,:v1], df[!,:v3])) samples=5 evals=1 seconds=120
# BenchmarkTools.Trial:
#   memory estimate:  206.06 MiB
#   allocs estimate:  500331
#   --------------
#   minimum time:     439.494 ms (0.00% GC)
#   median time:      472.370 ms (0.00% GC)
#   mean time:        479.215 ms (0.00% GC)
#   maximum time:     558.852 ms (0.00% GC)
#   --------------
#   samples:          5
#   evals/sample:     1

################################################################################
# DT[, lapply(.SD, mean), keyby=id4, .SDcols=7:9]
# TO BEAT: 0.6s for 10m; 2.76s for 100m; 6.45s for 250m
# TODO: make CountingSort
# Status: MUCH SLOWER if the number of elements is large and saturating memory
################################################################################
df[!,:id4_32] = Int32.(df[!,:id4])
df[!,:id4_8] = Int8.(df[!,:id4])

@time a = fastby((mean, mean, mean), df[!,:id4], (df[!,:v1], df[!,:v2], df[!,:v3]));
@time a = fastby((mean, mean, mean), df[!,:id4], (df[!,:v1], df[!,:v2], df[!,:v3]));

@time a = fastby((mean, mean, mean), df[!,:id4_32], (df[!,:v1], df[!,:v2], df[!,:v3]));
@time a = fastby((mean, mean, mean), df[!,:id4_32], (df[!,:v1], df[!,:v2], df[!,:v3]));

@time a = fastby((mean, mean, mean), df[!,:id4_8], (df[!,:v1], df[!,:v2], df[!,:v3]));
@time a = fastby((mean, mean, mean), df[!,:id4_8], (df[!,:v1], df[!,:v2], df[!,:v3]));

################################################################################
# DT[, lapply(.SD, sum), keyby=id6, .SDcols=7:9]
# TO BEAT: 0.4s for 10m; 11.33s for 100m; 31.31 for 250m
################################################################################
df[!,:id6] = Int32.(df[!,:id6]); # this is needed not because it was `Int64` but `Float64`
@time a = fastby((sum, sum, sum), df[!,:id6], (df[!,:v1], df[!,:v2], df[!,:v3]));
@time a = fastby((sum, sum, sum), df[!,:id6], (df[!,:v1], df[!,:v2], df[!,:v3]));

# @time a = FastGroupBy.fastby4((sum, sum, sum), df[!,:id6], (df[!,:v1], df[!,:v2], df[!,:v3]));
# fns, byvec, valvec = (sum, sum, sum), df[!,:id6], (df[!,:v1], df[!,:v2], df[!,:v3])

using Base.Threads
function fastby2(fns, byvec, valvec)
    l = length(fns)
    for i=1:l
        fastby(fns[i], byvec, valvec[i])
    end
end

function fastby3(fns, byvec, valvec)
    l = length(fns)
    @threads for i=1:l
        fastby(fns[i], byvec, valvec[i])
    end
end

@time a = fastby2((sum, sum, sum), df[!,:id6], (df[!,:v1], df[!,:v2], df[!,:v3]));
@time a = fastby3((sum, sum, sum), df[!,:id6], (df[!,:v1], df[!,:v2], df[!,:v3]));
@time a = fastby.((sum, sum, sum), (df[!,:id6], df[!,:id6], df[!,:id6]), (df[!,:v1], df[!,:v2], df[!,:v3]));

@time a = fgroupreduce(+, df[!,:id6], df[!,:v1])

@benchmark fastby(sum, df[!,:id6], df[!,:v1]) samples=5 evals=1 seconds=60

@time aggregate(df[[:id4, :v1]], :id4, sum)

using BenchmarkTools
@benchmark fastby2((sum, sum, sum), df[!,:id6], (df[!,:v1], df[!,:v2], df[!,:v3]))
# BenchmarkTools.Trial:
#   memory estimate:  370.06 MiB
#   allocs estimate:  300124
#   --------------
#   minimum time:     246.012 ms (0.00% GC)
#   median time:      258.472 ms (0.00% GC)
#   mean time:        318.356 ms (18.04% GC)
#   maximum time:     1.192 s (77.09% GC)
#   --------------
#   samples:          16
#   evals/sample:     1
@benchmark fastby3((sum, sum, sum), df[!,:id6], (df[!,:v1], df[!,:v2], df[!,:v3]))
# BenchmarkTools.Trial:
#   memory estimate:  291.33 MiB
#   allocs estimate:  254586
#   --------------
#   minimum time:     120.978 ms (0.00% GC)
#   median time:      136.074 ms (0.00% GC)
#   mean time:        222.204 ms (37.22% GC)
#   maximum time:     1.235 s (88.69% GC)
#   --------------
#   samples:          25
#   evals/sample:     1

@time a = fastby(sum, df[!,:id6], df[!,:v1])
@which fastby(sum, df[!,:id6], df[!,:v1])

################################################################################
# plotting of results
################################################################################


text:
"""
BarycentricInterpolation
========================

This package implements the Barycentric formula for polynomial interpolation on
equispaced points and Chebyshev points of the first and second kind. The
formulae used are taken from the paper of Berrut and Trefethen, SIAM Review,
2004.

Additional Barycentric weights for Legendre points are taken from the paper of Wang,
Huybrechs, and Vandewalle, Mathematics of Computation, 83 (290) 2893-2914, 2014.

Other packages that may be of interest are FastGaussQuadrature and ApproxFun.

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2016-2021 and licensed
under the MIT license <https://opensource.org/licenses/MIT>
"""
module BarycentricInterpolation

using FastGaussQuadrature: gausslegendre

export Chebyshev1, Chebyshev2, Legendre, Equispaced, ArbitraryPolynomial

export weights, nodes, interpolate, interpolation_matrix, differentiation_matrix, degree

#-- Node distributions

abstract type AbstractPolynomial{N, T <: Number} end

for name in [:Chebyshev1, :Chebyshev2, :Legendre, :Equispaced]
    @eval struct $name{N, T <: Number} <: AbstractPolynomial{N, T}
        shift::T
        scale::T
        nodes::Vector{T}
        weights::Vector{T}
        function $name{N, T}(start::T, stop::T) where {N, T <: Number}
            _shift = T(stop + start)/2
            _scale = T(stop - start)/2
            (_nodes, _weights) = nodes_weights($name{N, T}, _shift, _scale)
            new{N, T}(_shift, _scale, _nodes, _weights)
        end
    end
    @eval $name{N, T}(start=-1, stop=1) where {N, T} = $name{N, T}(convert(T, start), convert(T, stop))
    @eval $name{N}(start=-1, stop=1) where N = $name{N, Float64}(convert(Float64, start), convert(Float64, stop))
    @eval (poly::$name)(y) = interpolate(poly, y)
    @eval (poly::$name)(y, x) = interpolate(poly, y, x)
end

struct ArbitraryPolynomial{N, T <: Number} <: AbstractPolynomial{N, T}
    nodes::Vector{T}
    weights::Vector{T}
    function ArbitraryPolynomial(nodes::Vector{T}) where T <: Number
        N = length(nodes)-1
        new{N, T}(nodes, weights(ArbitraryPolynomial{N, T}, nodes))
    end
end
(poly::ArbitraryPolynomial)(y) = interpolate(poly, y)
(poly::ArbitraryPolynomial)(y, x) = interpolate(poly, y, x)

"""
    degree(poly)

Return the degree of the polynomial specified.
"""
degree(poly::AbstractPolynomial{N}) where N = N

"""
    nodes_weights(poly)

Return the nodes and weights of the polynomial specified.
"""
function nodes_weights(poly::Type{<:AbstractPolynomial}, shift=0, scale=1)
    return (nodes(poly, shift, scale), weights(poly))
end

nodes_weights(poly::AbstractPolynomial) = (poly.nodes, poly.weights)

function nodes_weights(::Type{<:Legendre{N, Float64}}, shift=0, scale=1) where N
    (x, w) = gausslegendre(N+1)
    _nodes = [xᵢ * scale + shift for xᵢ in x]
    _weights = [(2*isodd(i)-1)*sqrt((1-x[i]^2)*w[i]) for i in eachindex(x)]
    return (_nodes, _weights)
end

"""
    weights(poly)

Return the Barycentric weights for the specified orthogonal polynomials.
"""
function weights end

weights(poly::Type{<: AbstractPolynomial{N}}) where N = [_weight(poly, j) for j = 0:N]

# Eq. (5.1)
@inline _weight(poly::Type{<: Equispaced{N, T}}, j::Integer) where {N, T} = T((2*xor(isodd(N), iseven(j))-1)*binomial(N, j))

# Eq. (5.3)
@inline _weight(poly::Type{<: Chebyshev1{N, T}}, j::Integer) where {N, T} = -(2*iseven(j)-1)*sinpi(T(2j + 1)/(2N + 2))

# Eq. (5.4)
@inline _weight(poly::Type{<: Chebyshev2{N, T}}, j::Integer) where {N, T} = -T((1.0 - 0.5*((j == 0) || (j == N)))*(2*iseven(j) - 1))

function weights(poly::Type{<: ArbitraryPolynomial{N, T}}, x::AbstractVector{T}) where {N, T}
    w = similar(x)
    for i = eachindex(w)
        xᵢ = x[i]
        wᵢ = one(T)
        for j = 1:i-1
            wᵢ *= xᵢ - x[j]
        end
        for j = i+1:N+1
            wᵢ *= xᵢ - x[j]
        end
        w[i] = 1/wᵢ
    end
    return w
end

weights(poly::AbstractPolynomial) = poly.weights

"""
    nodes(poly)

Return the nodes for the specified orthogonal polynomials.
"""
function nodes end

nodes(poly::Type{<: AbstractPolynomial{N, T}}) where {N, T} = nodes(poly, zero(T), one(T))

nodes(poly::Type{<: AbstractPolynomial{N}}, shift::Number, scale::Number) where N = [_node(poly, j)*scale + shift for j = 0:N]

nodes(poly::Type{<: Equispaced{N}}, shift::Number, scale::Number) where N = range(shift - scale, stop=shift + scale, length=N+1)

@inline _node(poly::Type{<: Chebyshev1{N, T}}, j::Integer) where {N, T} = -cospi(T(2j + 1)/(2N + 2))

@inline _node(poly::Type{<: Chebyshev2{N, T}}, j::Integer) where {N, T} = -cospi(T(j)/N)

nodes(poly::AbstractPolynomial) = poly.nodes

"""
    interpolation_matrix(poly::AbstractPolynomial, x)

Return the interpolation matrix from the nodes of `poly` to the point(s) `x`.
For example :

    P = Chebyshev2{5}()
    x = range(-1, stop=1, length=10)
    M = interpolation_matrix(P, x)

Now `y(x) ≈ M*y₀` given that `y(nodes(poly)) = y₀.
"""
function interpolation_matrix(poly::AbstractPolynomial{N, T}, x::Union{Number, AbstractVector}) where {N, T}
    w = weights(poly)
    x₀ = nodes(poly)
    M = Matrix{T}(undef, length(x), N+1)
    # Eq. (4.2)
    for j = eachindex(x)
        xx = convert(T, x[j])
        Msum = zero(T)
        exact = 0
        for i = Base.OneTo(N + 1)
            exact = ifelse(xx == x₀[i], i, exact)
            M[j, i] = w[i] / (xx - x₀[i])
            Msum += M[j, i]
        end
        if Msum == 0
            for i = Base.OneTo(N + 1)
                M[j, i] = zero(T)
            end
        elseif exact > 0
            for i = Base.OneTo(N + 1)
                M[j, i] = zero(T)
            end
            M[j, exact] = one(T)
        else
            for i = Base.OneTo(N + 1)
                M[j, i] /= Msum
            end
        end
    end
    return M
end

"""
    interpolate(poly::AbstractPolynomial, y₀, [x])

Return the value of `y(x)` given that `y(nodes(poly)) = y₀`. If the value of `x`
is not provided, return a function `y(x)` that evaluates the interpolant at any
`x`.
"""
function interpolate(poly::AbstractPolynomial{N, T}, y₀::AbstractVector{<: Number}, x::Number) where {N, T}
    w = weights(poly)
    x₀ = nodes(poly)
    xx = convert(T, x)
    # Eq. (4.2)
    numer = zero(T)
    denom = zero(T)
    exact = 0
    for j = Base.OneTo(N+1)
        xdiff = xx - x₀[j]
        temp = w[j] / xdiff
        numer += temp*convert(T, y₀[j])
        denom += temp
        exact = ifelse(xdiff == 0, j, exact)
    end
    if exact > 0
        return convert(T, y₀[exact])
    else
        return numer/denom
    end
end

interpolate(poly::AbstractPolynomial{N, T}, y₀::AbstractVector{<: Number}) where {N, T} = x -> interpolate(poly, y₀, x)
interpolate(poly::AbstractPolynomial{N, T}, y₀::AbstractVector{<: Number}, x::AbstractVector{<: Number}) where {N, T} = [interpolate(poly, y₀, xᵢ) for xᵢ in x]

"""
    differentiation_matrix(poly::AbstractPolynomial)

Return the differentiation matrix at the nodes of the polynomial specified.

    P = Chebyshev2{5}()
    D = differentiation_matrix(P)

Now dy/dx ≈ `D*y` at the nodes of the polynomial.
"""
function differentiation_matrix(poly::AbstractPolynomial{N, T}) where {N, T}
    # Eqs. (9.4) and (9.5)
    w = weights(poly)
    x = nodes(poly)
    D = Matrix{T}(undef, N+1, N+1)
    for i = Base.OneTo(N+1)
        Dsum = zero(T)
        for j = Base.OneTo(i-1)
            temp = (w[j] / w[i]) / (x[i] - x[j])
            D[i, j] = temp
            Dsum += temp
        end
        for j = i+1:N+1
            temp = (w[j] / w[i]) / (x[i] - x[j])
            D[i, j] = temp
            Dsum += temp
        end
        D[i, i] = -Dsum
    end
    return D
end

end  # module


text:
using ..DGMethods.NumericalFluxes:
    NumericalFluxFirstOrder, NumericalFluxGradient, NumericalFluxSecondOrder

import ..DGMethods: boundary_state!

export OceanBC

"""
    OceanBC(velocity    = Impenetrable(NoSlip())
            temperature = Insulating())

The standard boundary condition for OceanModel. The default options imply a "no flux" boundary condition.
"""
Base.@kwdef struct OceanBC{M, T}
    velocity::M = Impenetrable(NoSlip())
    temperature::T = Insulating()
end

"""
    boundary_state!(nf, ::HBModel, args...)

applies boundary conditions for the hyperbolic fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@inline function boundary_state!(nf, ocean::HBModel, args...)
    boundary_condition = ocean.problem.boundary_condition
    return ocean_boundary_state!(nf, boundary_condition, ocean, args...)
end

"""
    boundary_state!(nf, ::LinearHBModel, args...)

applies boundary conditions for the hyperbolic fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@inline function boundary_state!(nf, linear::LinearHBModel, args...)
    ocean = linear.ocean
    boundary_condition = ocean.problem.boundary_condition

    return ocean_boundary_state!(nf, boundary_condition, ocean, args...)
end

"""
    ocean_boundary_state!(nf, boundaries::Tuple, ::HBModel,
                          Q⁺, A⁺, n, Q⁻, A⁻, bctype)
applies boundary conditions for the first-order and gradient fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@generated function ocean_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    boundaries::Tuple,
    ocean,
    Q⁺,
    A⁺,
    n,
    Q⁻,
    A⁻,
    bctype,
    t,
    args...,
)
    N = fieldcount(boundaries)
    return quote
        Base.Cartesian.@nif(
            $(N + 1),
            i -> bctype == i, # conditionexpr
            i -> ocean_boundary_state!(
                nf,
                boundaries[i],
                ocean,
                Q⁺,
                A⁺,
                n,
                Q⁻,
                A⁻,
                t,
            ), # expr
            i -> error("Invalid boundary tag")
        ) # elseexpr
        return nothing
    end
end

"""
    ocean_boundary_state!(nf, boundaries::Tuple, ::HBModel,
                          Q⁺, A⁺, D⁺, n, Q⁻, A⁻, D⁻, bctype)
applies boundary conditions for the second-order fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@generated function ocean_boundary_state!(
    nf::NumericalFluxSecondOrder,
    boundaries::Tuple,
    ocean,
    Q⁺,
    D⁺,
    A⁺,
    n,
    Q⁻,
    D⁻,
    A⁻,
    bctype,
    t,
    args...,
)
    N = fieldcount(boundaries)
    return quote
        Base.Cartesian.@nif(
            $(N + 1),
            i -> bctype == i, # conditionexpr
            i -> ocean_boundary_state!(
                nf,
                boundaries[i],
                ocean,
                Q⁺,
                D⁺,
                A⁺,
                n,
                Q⁻,
                D⁻,
                A⁻,
                t,
            ), # expr
            i -> error("Invalid boundary tag")
        ) # elseexpr
        return nothing
    end
end

"""
    ocean_boundary_state!(nf, bc::OceanBC, ::HBModel)

splits boundary condition application into velocity and temperature conditions
"""
function ocean_boundary_state!(nf, bc::OceanBC, ocean, args...)
    ocean_velocity_boundary_state!(nf, bc.velocity, ocean, args...)
    ocean_temperature_boundary_state!(nf, bc.temperature, ocean, args...)
end

include("bc_velocity.jl")
include("bc_temperature.jl")

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
"""
测试温度流的refine
"""

using Plots
using MARY_fRG.Drawers
using MARY_fRG.Fermi
using MARY_fRG.FlowEquation

#阻止画图
ENV["GKSwstype"] = "100"

function run_rftf()
    model = common_square_lattice(0.20)
    Γ4 = TFGamma4(
        model, 8.0, 16, 800
    )
    Γ4.V .+= 1.0
    lval = 0.
    lstep = 0.01
    #
    for idx in 1:1:2501
        #if idx % 10 == 1
        #    blval = lval - 5
        #    blval = max(blval, 0.)
        #    blval = min(blval, 5.)
        #    Γ4b = TFGamma4_addition_ltris_mt(Γ4, blval; maxnum=1000000)
        #    println(length(Γ4.ltris), "->", length(Γ4b.ltris))
        #    #plt = draw_points([tri.center for tri in Γ4b.ltris])
        #    #@savepng plt "brlu"*string(lval)
        #end
        #blval = min(lval, 10.0)
        bubb_pp, bubb_fs, bubb_ex = all_bubble_tf_mt(Γ4, lval)
        dl = dl_tf_mt(Γ4, bubb_pp, bubb_fs, bubb_ex)
        Γ4.V .+= dl .* lstep
        if idx % 50 == 1
            plt = heatmap(Γ4.V[1, 1, 1, 1, :, :, 1])
            png(plt, "Gamma4"*string(idx))
        end
        lval += lstep
    end
end

run_rftf()


text:
"""
    DynamicsJacobian{n,nm,T}

Custom `n × (n+m)` matrix specifying a dynamics Jacobian for a forced dynamical system with
`n` states and `m` controls. The Jacobian is structured as `[∂x ∂u]` where `x` and `u` are
the state and control vectors, respectively.

The `DynamicsJacobian` `D` provides access to the partial derivatives `A = ∂x`  and `B = ∂u` via
direct access `D.A` and `D.B`, returning a view into the underlying `Matrix`, or

    RobotDynamics.get_A(D)
    RobotDynamics.get_B(D)

which return an `SMatrix`. Note that this method should be used with caution for systems with
large state and/or control dimensions.

# Constructors
    DynamicsJacobian(model::AbstractModel)
    DynamicsJacobian(n::Int, m::Int)
    DynamicsJacobian(D::StaticMatrix)

where `D` is a `StaticMatrix` of appropriate size. Since `DynamicsJacobian` implements the
`StaticMatrix` interface, is also supports all the constructors and operations inherent to
a `StaticMatrix`.
"""
struct DynamicsJacobian{S1,S2,T} <: StaticMatrix{S1,S2,T}
    data::SizedMatrix{S1,S2,T,2,Matrix{T}}
    A::SubArray{T,2,Matrix{T}, Tuple{Base.Slice{Base.OneTo{Int}}, UnitRange{Int}},true}
    B::SubArray{T,2,Matrix{T}, Tuple{Base.Slice{Base.OneTo{Int}}, UnitRange{Int}},true}
    function DynamicsJacobian(F::StaticMatrix{n,nm}) where {n,nm}
        m = nm-n
        data = SizedMatrix{n,nm}(F)
        A = view(data.data, :, 1:n)
        B = view(data.data, :, n .+ (1:m))
        new{n,nm,eltype(data)}(data, A, B)
    end
end

DynamicsJacobian(n::Int, m::Int) = DynamicsJacobian(SizedMatrix{n,n+m}(zeros(n,n+m)))
@inline DynamicsJacobian(model::AbstractModel) = DynamicsJacobian(state_dim(model), control_dim(model))

@inline DynamicsJacobian{S1,S2}(t::NTuple) where {S1,S2}= DynamicsJacobian(SMatrix{S1,S2}(t))
@inline DynamicsJacobian{S1,S2,T}(t::NTuple) where {S1,S2,T}= DynamicsJacobian(SMatrix{S1,S2,T}(t))
Base.@propagate_inbounds Base.getindex(D::DynamicsJacobian, i::Int) = D.data[i]
Base.@propagate_inbounds Base.setindex!(D::DynamicsJacobian, x, i::Int) = D.data[i] = x
Base.@propagate_inbounds Base.setindex!(D::DynamicsJacobian, x, inds::Vararg{Int64}) =
    setindex!(D.data, x, inds...)
@inline Base.Tuple(D::DynamicsJacobian) = Tuple(D.data)

"""
    get_A(D::DynamicsJacobian)

Return the partial derivative of the dynamics with respect to the state input as an `SMatrix`,
given the full dynamics Jacobian.
"""
@generated function get_static_A(D::DynamicsJacobian{n,nm}) where {n,nm}
    expr = [:(D[$i]) for i = 1:n*n]
    :(SMatrix{$n,$n}(tuple($(expr...))))
end

"""
    get_B(D::DynamicsJacobian)

Return the partial derivative of the dynamics with respect to the control input as an `SMatrix`,
given the full dynamics Jacobian.
"""
@generated function get_static_B(D::DynamicsJacobian{n,nm}) where {n,nm}
    m = nm - n
    inds = (1:n*m) .+ n*n
    expr = [:(D[$i]) for i in inds]
    :(SMatrix{$n,$m}(tuple($(expr...))))
end

get_A(D::DynamicsJacobian) = D.A
get_B(D::DynamicsJacobian) = D.B

@inline get_data(A::AbstractArray) = A
@inline get_data(A::SizedArray) = A.data
@inline get_data(A::DynamicsJacobian) = A.data.data


text:
add_verbose_scope(:NfOrd)
add_assert_scope(:NfOrd)

#set_verbose_level(:NfOrd, 1)

include("NfOrd/NfMaxOrd.jl")
include("NfOrd/Ideal.jl")
include("NfOrd/Zeta.jl")
include("NfOrd/FracIdl.jl")
include("NfOrd/Clgp.jl")
include("NfOrd/Unit.jl")
include("NfOrd/ResidueField.jl")
include("NfOrd/ResidueRing.jl")
include("NfOrd/ResidueRingMultGrp.jl")
include("NfOrd/FactorBaseBound.jl")
include("NfOrd/FacElem.jl")
include("NfOrd/LinearAlgebra.jl")
include("NfOrd/IdealLLL.jl")
include("NfOrd/Narrow.jl")
include("NfOrd/norm_eqn.jl")
include("NfOrd/RayClass.jl")
include("NfOrd/RayClassFacElem.jl")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""MetricStatus describes the last-read state of a single metric.

    IoK8sApiAutoscalingV2beta1MetricStatus(;
        external=nothing,
        object=nothing,
        pods=nothing,
        resource=nothing,
        type=nothing,
    )

    - external::IoK8sApiAutoscalingV2beta1ExternalMetricStatus : external refers to a global metric that is not associated with any Kubernetes object. It allows autoscaling based on information coming from components running outside of cluster (for example length of queue in cloud messaging service, or QPS from loadbalancer running outside of cluster).
    - object::IoK8sApiAutoscalingV2beta1ObjectMetricStatus : object refers to a metric describing a single kubernetes object (for example, hits-per-second on an Ingress object).
    - pods::IoK8sApiAutoscalingV2beta1PodsMetricStatus : pods refers to a metric describing each pod in the current scale target (for example, transactions-processed-per-second).  The values will be averaged together before being compared to the target value.
    - resource::IoK8sApiAutoscalingV2beta1ResourceMetricStatus : resource refers to a resource metric (such as those specified in requests and limits) known to Kubernetes describing each pod in the current scale target (e.g. CPU or memory). Such metrics are built in to Kubernetes, and have special scaling options on top of those available to normal per-pod metrics using the \"pods\" source.
    - type::String : type is the type of metric source.  It will be one of \"Object\", \"Pods\" or \"Resource\", each corresponds to a matching field in the object.
"""
mutable struct IoK8sApiAutoscalingV2beta1MetricStatus <: SwaggerModel
    external::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta1ExternalMetricStatus } # spec name: external
    object::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta1ObjectMetricStatus } # spec name: object
    pods::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta1PodsMetricStatus } # spec name: pods
    resource::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta1ResourceMetricStatus } # spec name: resource
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiAutoscalingV2beta1MetricStatus(;external=nothing, object=nothing, pods=nothing, resource=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta1MetricStatus, Symbol("external"), external)
        setfield!(o, Symbol("external"), external)
        validate_property(IoK8sApiAutoscalingV2beta1MetricStatus, Symbol("object"), object)
        setfield!(o, Symbol("object"), object)
        validate_property(IoK8sApiAutoscalingV2beta1MetricStatus, Symbol("pods"), pods)
        setfield!(o, Symbol("pods"), pods)
        validate_property(IoK8sApiAutoscalingV2beta1MetricStatus, Symbol("resource"), resource)
        setfield!(o, Symbol("resource"), resource)
        validate_property(IoK8sApiAutoscalingV2beta1MetricStatus, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiAutoscalingV2beta1MetricStatus

const _property_map_IoK8sApiAutoscalingV2beta1MetricStatus = Dict{Symbol,Symbol}(Symbol("external")=>Symbol("external"), Symbol("object")=>Symbol("object"), Symbol("pods")=>Symbol("pods"), Symbol("resource")=>Symbol("resource"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiAutoscalingV2beta1MetricStatus = Dict{Symbol,String}(Symbol("external")=>"IoK8sApiAutoscalingV2beta1ExternalMetricStatus", Symbol("object")=>"IoK8sApiAutoscalingV2beta1ObjectMetricStatus", Symbol("pods")=>"IoK8sApiAutoscalingV2beta1PodsMetricStatus", Symbol("resource")=>"IoK8sApiAutoscalingV2beta1ResourceMetricStatus", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta1MetricStatus }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta1MetricStatus))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta1MetricStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta1MetricStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta1MetricStatus }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta1MetricStatus[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta1MetricStatus)
    (getproperty(o, Symbol("type")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta1MetricStatus }, name::Symbol, val)
end


text:
@testset "count missing" begin
    for i in 1:10
        x = rand(10:30)
        rad = Array{Union{Float64, Missing}}(rand(20:100.0, x))
        missings = unique(rand(1:x,rand(1:x)))
        rad[missings] .= missing
        @test NighttimeLights.count_missing(rad) == length(missings)
    end
end

@testset "max_missing" begin
    for i in 1:10
        x = rand(10:30)
        rad = Array{Union{Float64, Missing}}(rand(20:100.0, x))
        max = findmax(rad)
        missings = rand(1:20)
        for j in 1:missings
            missing_index = rand(1:x)
            if missing_index == max[2]
                continue
            else
                rad[missing_index] = missing
            end
        end
        @test NighttimeLights.max_missing(rad) == max[1]
    end
end


text:
module WorkQueueUtils
using  ..BitWiseUtils,CUDA, Logging,..CUDAGpuUtils, ..ResultListUtils,..WorkQueueUtils, Logging,StaticArrays, ..IterationUtils, ..ReductionUtils, ..CUDAAtomicUtils,..MetaDataUtils, ..BitWiseUtils
export allocateWorkQueue,appendToWorkQueue,@appendToWorkQueueBasic,@appendToLocalWorkQueue,@appendToGlobalWorkQueue
"""
allocate memory for  work queues
    entries means
        1)xMeta
        2)yMeta
        3)zMeta
        4)isGold
 In order to prevent overwriting  we will create 8 separe work queues for each even or odd metax,metay and meta Z ...
    also initializes counters for each
    names are encoded like E is about even and O for odd so workQueueEEE will mean that all metaX,metaY and metaZ are even
   returns  workQueue,workQueueCounter

    """
function allocateWorkQueue(metaDataLength)
    queueSize = cld(metaDataLength*2,8)+2#*2 becouse of gold and segm pass divided by 8 becouse we have 8 work queues
    return (  CUDA.zeros(UInt16,4,(queueSize)),CUDA.zeros(UInt16,(1)) )

end

"""
atomically append the block linear index and information is it gold or other pass
also we need to be sure that we appended to the correct work queue based on the properties of the xMeta,yMeta,zMeta - so are they even, odd ...
"""
macro appendToLocalWorkQueue(metaX,metaY,metaZ,isGold )
    return esc(quote
    old =  atomicallyAddOne(workCounterLocalInShmem)
      #CUDA.@cuprint " ooo old $(old) metaX $($metaX) metaY $($metaY) metaZ $($metaZ) isGold $($isGold) workCounterLocalInShmem[1] $(workCounterLocalInShmem[1]) \n"
    @inbounds shmemblockData[old*4+1]= $metaX
    @inbounds shmemblockData[old*4+2]= $metaY
    @inbounds shmemblockData[old*4+3]= $metaZ
    @inbounds shmemblockData[old*4+4]= $isGold

    # workQueue[1,old]= $metaX
    # workQueue[2,old]= UInt16($metaY)
    # workQueue[3,old]= UInt16($metaZ)
    # workQueue[4,old]= UInt16($isGold)

    end)#qote
end#appendToWorkQueue

"""
the function above is block private now we need to load the work queaue into main work queue
    and additionally atomically add to general work queue counter
"""
macro appendToGlobalWorkQueue()
    return esc(quote
    if(workCounterLocalInShmem[1]>0)
        #old is a spot where we should start adding new entries
        sync_threads()
        @ifXY 1 1 workCounterHelper[1] =  atomicallyAddToSpot(workQueueCounter,1,workCounterLocalInShmem[1])-workCounterLocalInShmem[1]
        sync_threads()
        # @ifXY 1 1 CUDA.@cuprint "appendToGlobalWorkQueue old $(old)  workCounterLocalInShmem[1] $(workCounterLocalInShmem[1]) \n "
        #we multiply it by 4 so each thread will have potentially only one data entry to do
        @iterateLinearlyCheckAll(cld(workCounterLocalInShmem[1]*4,blockDimX()*blockDimY()),workCounterLocalInShmem[1]*4,begin
            @inbounds workQueue[(workCounterHelper[1])*4+i] =  shmemblockData[i]
        end)
    end

end)#qote
end#appendToWorkQueue

# macro appendToWorkQueueBasic(workQueaue,workQueauecounter, metaX,metaY,metaZ,isGold )
#     return esc(quote
#     old =  atomicallyAddOne($workQueauecounter)+1
#     # CUDA.@cuprint " ooo old $(old) metaX $($metaX) metaY $($metaY) metaZ $($metaZ) isGold $($isGold) workQueaue $(length($workQueaue)) \n"

#     $workQueaue[1,old]= $metaX
#     $workQueaue[2,old]= UInt16($metaY)
#     $workQueaue[3,old]= UInt16($metaZ)
#     $workQueaue[4,old]= UInt16($isGold)
#     end)#qote
# end#appendToWorkQueue

end#WorkQueueUtils

# basicList = [-1,0,1]
# function isAdjacent(liist)::Bool
#     for tupl in liist
#         for tuplInner in liist
#             xx=0;yy=0;zz=1
#             if( (tupl[1]+xx,tupl[2]+yy,tupl[3]+zz)== tuplInner   )
#                     return true
#                 end

#            xx=0;yy=0;zz=-1
#             if( (tupl[1]+xx,tupl[2]+yy,tupl[3]+zz)== tuplInner   )
#                     return true
#                 end

#            xx=0;yy=-1;zz=0
#             if( (tupl[1]+xx,tupl[2]+yy,tupl[3]+zz)== tuplInner   )
#                     return true
#                 end

#            xx=0;yy=1;zz=0
#             if( (tupl[1]+xx,tupl[2]+yy,tupl[3]+zz)== tuplInner   )
#                     return true
#                 end

#            xx=1;yy=0;zz=0
#             if( (tupl[1]+xx,tupl[2]+yy,tupl[3]+zz)== tuplInner   )
#                     return true
#                 end

#            xx=-1;yy=0;zz=0
#             if( (tupl[1]+xx,tupl[2]+yy,tupl[3]+zz)== tuplInner   )
#                     return true
#                 end

#         end
#     end
#     return false
# end

# even= [2,4,6]
# odd = [1,3,5]

# lA = []
# for x in even,y in even,z in even
#     push!(lA,(x,y,z)  )
# end
# lB = []
# for x in odd,y in even,z in even
#     push!(lB,(x,y,z)  )
# end
# lC = []
# for x in even,y in odd,z in even
#     push!(lC,(x,y,z)  )
# end
# lD = []
# for x in even,y in even,z in odd
#     push!(lD,(x,y,z)  )
# end
# lE = []
# for x in odd,y in odd,z in even
#     push!(lE,(x,y,z)  )
# end
# lF = []
# for x in even,y in odd,z in odd
#     push!(lF,(x,y,z)  )
# end
# lG = []
# for x in odd,y in odd,z in odd
#     push!(lG,(x,y,z)  )
# end
# lH = []
# for x in odd,y in even,z in odd
#     push!(lH,(x,y,z)  )
# end

# alll =[lA,lB,lC,lD,lE,lF,lG,lH]

# map(it-> isAdjacent(it),alll)
# sum(map(it-> length(it),alll))== 6*6*6

# evens

    #@appendToWorkQueueBasic(workQueue,workQueuecounter, $metaX,$metaY,$metaZ,$isGold )

        # if(iseven($metaX) && iseven($metaY) && iseven($metaZ) )
        #     @appendToWorkQueueBasic(workQueueEEE,workQueueEEEcounter, $metaX,$metaY,$metaZ,$isGold )
        # elseif(iseven($metaX) && isodd($metaY) && iseven($metaZ))
        #     @appendToWorkQueueBasic(workQueueEOE,workQueueEOEcounter, $metaX,$metaY,$metaZ,$isGold )
        # elseif(iseven($metaX) && iseven($metaY) && isodd($metaZ))
        #     @appendToWorkQueueBasic(workQueueEEO,workQueueEEOcounter, $metaX,$metaY,$metaZ,$isGold )
        # elseif(isodd($metaX) && iseven($metaY) && iseven($metaZ))
        #     @appendToWorkQueueBasic(workQueueOEE,workQueueOEEcounter, $metaX,$metaY,$metaZ,$isGold )
        # elseif(isodd($metaX) && isodd($metaY) && iseven($metaZ))
        #     @appendToWorkQueueBasic(workQueueOOE,workQueueOOEcounter, $metaX,$metaY,$metaZ,$isGold )
        # elseif(iseven($metaX) && isodd($metaY) && isodd($metaZ))
        #     @appendToWorkQueueBasic(workQueueEOO,workQueueEOOcounter, $metaX,$metaY,$metaZ,$isGold )
        # elseif(isodd($metaX) && iseven($metaY) && isodd($metaZ))
        #     @appendToWorkQueueBasic(workQueueOEO,workQueueOEOcounter, $metaX,$metaY,$metaZ,$isGold )
        # elseif(isodd($metaX) && isodd($metaY) && isodd($metaZ))
        #     @appendToWorkQueueBasic(workQueueOOO,workQueueOOOcounter, $metaX,$metaY,$metaZ,$isGold )
        # end

            # workQueueEEE,workQueueEEEcounter,workQueueEEO,workQueueEEOcounter
    # ,workQueueEOE,workQueueEOEcounter,workQueueOEE,workQueueOEEcounter
    # ,workQueueOOE,workQueueOOEcounter,workQueueEOO,workQueueEOOcounter
    # ,workQueueOEO,workQueueOEOcounter,workQueueOOO,workQueueOOOcounter
    # ,CUDA.zeros(UInt16,4,(queueSize)),CUDA.zeros(Int16,(1))
    # ,CUDA.zeros(UInt16,4,(queueSize)),CUDA.zeros(Int16,(1))
    # ,CUDA.zeros(UInt16,4,(queueSize)),CUDA.zeros(Int16,(1))
    # ,CUDA.zeros(UInt16,4,(queueSize)),CUDA.zeros(Int16,(1))
    # ,CUDA.zeros(UInt16,4,(queueSize)),CUDA.zeros(Int16,(1))
    # ,CUDA.zeros(UInt16,4,(queueSize)),CUDA.zeros(Int16,(1))
    # ,CUDA.zeros(UInt16,4,(queueSize)),CUDA.zeros(Int16,(1))


text:
function sbatch_test()
	@testset "sbatch" begin
		filePath = joinpath(test_file_dir(), "sbatch_test.sl");
		if isfile(filePath)
			rm(filePath);
		end
		cmdPath = joinpath(test_file_dir(), "test1.jl");
		outPath = joinpath(test_file_dir(), "test1.out");
		timeStr = "7-00";

		lineV = sbatch_commands(outPath, timeStr = timeStr);
		@test isa(lineV, Vector{String})
		@test length(lineV) > 5

		# write_sbatch_file(filePath, cmdPath, outPath,  timeStr = timeStr);
		# @test isfile(filePath)

		# # Show on screen
		# println("------------")
		# for line in eachline(filePath)
		# println(line)
	end
end


text:
######################### Lazy Difference Vector ###############################
# TODO: is it worth using the LazyArrays package?
# could be replaced by LazyVector(applied(-, x, y)), which is neat.
# lazy difference between two vectors, has no memory footprint
struct LazyDifference{T, U, V} <: AbstractVector{T}
    x::U
    y::V
    function LazyDifference(x, y)
        length(x) == length(y) || throw(DimensionMismatch("x and y do not have the same length: $(length(x)) and $(length(y))."))
        T = promote_type(eltype(x), eltype(y))
        new{T, typeof(x), typeof(y)}(x, y)
    end
end

difference(x::Number, y::Number) = x-y # avoid laziness for scalars
difference(x::AbstractVector, y::AbstractVector) = LazyDifference(x, y)
difference(x::Tuple, y::Tuple) = LazyDifference(x, y)
# this creates a type instability:
# difference(x, y) = length(x) == length(y) == 1 ? x[1]-y[1] : LazyDifference(x, y)
Base.size(d::LazyDifference) = (length(d.x),)
Base.getindex(d::LazyDifference, i::Integer) = d.x[i]-d.y[i]

################################################################################
LinearAlgebra.diagzero(D::Diagonal{<:Diagonal{T}}, i, j) where T = zeros(T, (size(D.diag[i], 1), size(D.diag[j], 2)))

# if we have Matrix valued kernels, this should be different
Base.eltype(k::AbstractKernel{T}) where {T} = T
Base.eltype(::MultiKernel{T}) where {T} = Matrix{T}

# TODO: simplify
const VecOfVec{T} = AbstractVector{<:AbstractVector{T}}
const VecOrVecOfVec{T} = AbstractVector{<:AbstractVector{T}}
const AbstractMatOrFac{T} = Union{AbstractMatrix{T}, Factorization{T}}
const AbstractVecOfVec{T} = AbstractVector{<:AbstractVector{T}}
const AbstractVecOfVecOrMat{T} = AbstractVector{<:AbstractVecOrMat{T}}

################################################################################
# euclidean distance
# IDEA: @inline?
function euclidean2(x, y)
    length(x) == length(y) || throw(DimensionMismatch("inputs have to have the same length: $(length(x)), $(length(y))"))
    val = zero(promote_type(eltype(x), eltype(y)))
    @inbounds @simd for i in eachindex(x)
        val += (x[i] - y[i])^2
    end
    return val
end
euclidean(x, y) = sqrt(euclidean2(x, y))

# energetic norm
enorm(A::AbstractMatOrFac, x::AbstractVector) = sqrt(enorm2(A, x))
enorm2(A::AbstractMatOrFac, x::AbstractVector) = dot(x, A, x)

# IDEA: maybe include in and output dimension of kernel in type?
# this makes it easier to type check admissability of input arguments
# and differentiate between data vector or vector of data
# or instead, write run-time input check
# change to a boundscheck?
function checklength(x::AbstractArray, y::AbstractArray)
    lx = length(x)
    ly = length(y)
    if lx != ly; throw(DimensionMismatch("length(x) ($lx) ≠ length(y) ($ly)")) end
    return length(x)
end

# do we need an MVector version of this?
function vector_of_static_vectors(A::AbstractMatrix)
    D = size(A, 1)
    [SVector{D}(col) for col in eachcol(A)]
end
function vector_of_static_vectors(A::AbstractVector{<:Number})
    [SVector{1}(a) for a in A]
end
vector_of_static_vectors(A::Vector{V}) where {V <: Union{SVector, MVector}} = A
# NOTE: element-vectors have
function vector_of_static_vectors(A::Vector{<:AbstractVector})
    # length(A) > 0 || return [SVector{0, eltype(A)}()]
    D = size(A[1], 1)
    [SVector{D}(a) for a in A]
end

# _Matrix(A::AbstractMatrix) = A
# _Matrix(A::Vector{<:Union{SVector, MVector}}) = A
# this first converts a colview back to the corresponding matrix
# function _Matrix(x::Vector{<:SubArray{<:Any, 1, <:Matrix, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}})
#     x[1].parent
# end

################################################################################
# is positive semi-definite
function ispsd end
ispsd(A::Number, tol::Real = 0.) = A ≥ -tol
# TODO: replace with own pivoted cholesky
ispsd(A::AbstractMatrix, tol::Real = 0.) = all(A->ispsd(A, tol), eigvals(A))
iscov(A::AbstractMatrix, tol::Real = 0.) = issymmetric(A) && ispsd(A, tol)
# import LinearAlgebraExtensions: iscov
iscov(A::AbstractMatrix{<:Real}, tol::Real = 0.) = issymmetric(A) && ispsd(A, tol)
iscov(k::MercerKernel, x = randn(32), tol = 1e-10) = iscov(gramian(k, x), tol)

####################### randomized stationarity tests ##########################
# NOTE this test only allows for euclidean isotropy ...
# does not matter for 1d tests
function isisotropic(k::AbstractKernel, x::AbstractVector)
    isiso = false
    if isstationary(k)
        isiso = true
        r = norm(x[1]-x[2])
        println(r)
        kxr = k(x_i/r, x_j/r)
        for x_i in x, x_j in x
            r = norm(x_i-x_j)
            val = k(x_i/r, x_j/r)
            if !(kxr ≈ val)
                isiso = false
            end
        end
    end
    return isiso
end

# tests if k is stationary on finite set of points x
# need to make this multidimensional
function isstationary(k::AbstractKernel, x::AbstractVector)
    n = length(x)
    d = length(x[1])
    is_stationary = true
    for i in 1:n, j in 1:n
        ε = eltype(x) <: AbstractArray ? randn(d) : randn()
        iseq = k(x[i], x[j]) ≈ k(x[i]+ε, x[j]+ε)
        if !iseq
            println(i ,j)
            println(x[i], x[j])
            println(k(x[i], x[j]) - k(x[i]+ε, x[j]+ε))
            is_stationary = false
            break
        end
    end
    K = typeof(k)
    if !is_stationary && isstationary(k)
        println("Covariance function " * string(K) * " is non-stationary but kernel k returns isstationary(k) = true.")
        return false
    end

    # if the kernel seems stationary but isn't a subtype of stationary kernel, suggest making it one
    if is_stationary && !isstationary(k)
        println("Covariance function " * string(K) * " seems to be stationary. Consider defining isstationary(k::$K) = true.")
    end
    return is_stationary
end

######################### perfect shuffle matrices #############################
# lazy representation of perfect shuffle matrices
# # X is n by m
# S::PerfectShuffle will perform S*vec(X) = vec(X')
struct PerfectShuffle <: AbstractMatrix{Bool}
    n::Int
    m::Int
end

Base.size(S::PerfectShuffle) = (n*m, n*m)
Base.size(S::PerfectShuffle, i::Int) = 0 < i ≤ 2 ? n*m : 1

SparseArrays.SparseMatrixCSC(S::PerfectShuffle) = perfect_shuffle(S.n, S.m)
function LinearAlgebra.mul!(b::AbstractVector, S::PerfectShuffle, a::AbstractVector, α::Number = 1, β::Number = 0)
    length(b) == length(a) == size(S, 1) || throw(DimensionMismatch())
    A, B = reshape(a, S.n, S.m), reshape(b, S.n, S.m)
    if β == 0
        @. B = α * A'
    else
        @. B = α * A' + β * B
    end
    return b
end
Base.:*(S::PerfectShuffle, a::AbstractVector) = mul!(zero(a), S, a)

function perfect_shuffle(n::Int)
    S = spzeros(Bool, n^2, n^2)
    for i in 1:n, j in 1:n
        S[j+n*(i-1), i+n*(j-1)] = 1
    end
    return S
end

# X is n by m
# S will perform S*vec(X) = vec(X')
function perfect_shuffle(n::Int, m::Int)
    S = spzeros(Bool, n*m, n*m)
    for i in 1:n, j in 1:m
        S[j+m*(i-1), i+n*(j-1)] = 1
    end
    return S
end

# exchange matrix (anti-diagonal identity matrix)
function exchange_matrix(r::Int)
    E = spzeros(Bool, r, r)
    for i in 1:r
        E[r-i+1, i] = 1
    end
    return E
end

# function anti_diagonal(x::AbstractVector)
#     i = tuple.(r:-1:1, 1:r)
# end

# computes all products of subsets of length(x)-1 elements of x and stores them in x
# useful for gradient product kernel structure
function leave_one_out_products!(x::AbstractArray)
    nz = sum(iszero, x)
    if nz == 0
        @. x = $prod(x) / x
    else
        @. x = 0 # with more than two zeros, all leave-one-out products will be zero
        if nz == 1
            i = findfirst(iszero, x)
            x[i] = prod(x[1:length(x) .!= i])
        end
    end
    return x
end


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

    CalibrationData.jl

    container for calibration data

    data has 3 dimensions:
    1. Well     :0 - :15
    2. Channel  1, 2
    3. Dye      water, dye1, dye2

    Author: Tom Price
    Date:   June 2019

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

import DataStructures.OrderedDict
import StaticArrays: SArray

struct CalibrationData{N <: NumberOfChannels, R <: Union{Int_T,Float_T}}
    array       ::SArray{S,R,3} where {S <: Tuple}
    ## inner constructors supply type parameters
    ## and enforce constraints on size of data array
    CalibrationData(
        ::Type{SingleChannel},
        data ::AbstractArray{R}) where{R <: Union{Int_T,Float_T}} =
            new{SingleChannel,R}(SArray{Tuple{size(data,1),1,2},R}(data))
    CalibrationData(
        ::Type{DualChannel},
        data ::AbstractArray{R}) where{R <: Union{Int_T,Float_T}} =
            new{DualChannel,R}(SArray{Tuple{size(data,1),2,3},R}(data))
end

## parser / outer constructor
function CalibrationData(calib ::Associative)
    local water, dye1, dye2 ## not enum
    num_wells = count_wells(calib)
    water     = calib[WATER_KEY][FLUORESCENCE_VALUE_KEY]
    R         = water |> first |> first |> typeof
	if (R == Int64 && Int_T != Int64)
		R = Int_T
	end

    dye1      = calib[CHANNEL_KEY * "_1"][FLUORESCENCE_VALUE_KEY]
    if length(dye1) > 1 && thing(dye1[2]) && haskey(calib, CHANNEL_KEY * "_2")
        dye2  = calib[CHANNEL_KEY * "_2"][FLUORESCENCE_VALUE_KEY]
        data  = try
            [water, dye1, dye2] |> gather(vcat) |> gather(hcat) |>
                morph(num_wells,2,3) |> mold(bless(R))
        catch
            warn(logger, "dual channel calibration data are in the wrong format")
            rethrow()
        end ## try
        return CalibrationData(DualChannel, data)
    end ## if
    ## else single channel
    data  = try
        [water, dye1] |> moose(first, hcat) |> morph(num_wells,1,2) |> mold(bless(R))
    catch
        warn(logger, "single channel calibration data are in the wrong format")
        rethrow()
    end ## try
    return CalibrationData(SingleChannel, data)
end ## CalibrationData()

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
struct Entity
    _comps::Vector{ComponentPoolIndex} # Reference to component's index in attached Kosmo._components_index
end
function Entity!(k::Kosmo, args...)
    v = Vector{ComponentPoolIndex}()
    for arg in args
        push!(v, add!(k, arg))
    end
    for set in k._component_pool._sets
        make_set(k, set, v)
    end
    Entity(v)
end

function make_set(k::Kosmo, set::Type{ComponentSet{T}}, v::Vector{ComponentPoolIndex}) where T
    add!(k, ComponentSet([filter(z->z.type==x, v)[1] for x in T]...))
end