text:
const TIME_SERIES_STORAGE_FILE = "time_series_storage.h5"

mutable struct SystemData <: InfrastructureSystemsType
    components::Components
    forecast_metadata::ForecastMetadata
    validation_descriptors::Vector
    time_series_storage::TimeSeriesStorage
    time_series_storage_file::Union{Nothing, String}  # only valid during serialization
    internal::InfrastructureSystemsInternal
end

function SystemData(; validation_descriptor_file=nothing, time_series_in_memory=false)
    if isnothing(validation_descriptor_file)
        validation_descriptors = Vector()
    else
        validation_descriptors = read_validation_descriptor(validation_descriptor_file)
    end

    ts_storage = make_time_series_storage(; in_memory=time_series_in_memory)
    components = Components(ts_storage, validation_descriptors)
    return SystemData(components, ForecastMetadata(), validation_descriptors, ts_storage,
                      nothing, InfrastructureSystemsInternal())
end

function SystemData(forecast_metadata, validation_descriptors, time_series_storage,
                    internal)
    components = Components(time_series_storage, validation_descriptors)
    return SystemData(components, forecast_metadata, validation_descriptors,
                      time_series_storage, nothing, internal)
end

"""
    SystemData(filename::AbstractString)

Deserialize SystemData from a JSON file.
"""
function SystemData(filename::AbstractString)
    return from_json(SystemData, filename)
end

"""
    add_forecasts!(
                   ::Type{T},
                   data::SystemData,
                   metadata_file::AbstractString,
                   resolution=nothing,
                  ) where T <: InfrastructureSystemsType

Adds forecasts from a metadata file or metadata descriptors.

# Arguments
- `::Type{T}`: forecasted component type; may be abstract
- `data::SystemData`: system
- `metadata_file::AbstractString`: metadata file for timeseries
  that includes an array of TimeseriesFileMetadata instances or a vector.
- `resolution::DateTime.Period=nothing`: skip forecast that don't match this resolution.
"""
function add_forecasts!(
                        ::Type{T},
                        data::SystemData,
                        metadata_file::AbstractString;
                        resolution=nothing,
                       ) where T <: InfrastructureSystemsType
    metadata = read_time_series_metadata(metadata_file)
    return add_forecasts!(T, data, metadata; resolution=resolution)
end

"""
    add_forecasts!(
                   data::SystemData,
                   timeseries_metadata::Vector{TimeseriesFileMetadata};
                   resolution=nothing,
                  )

Adds forecasts from a metadata file or metadata descriptors.

# Arguments
- `data::SystemData`: system
- `timeseries_metadata::Vector{TimeseriesFileMetadata}`: metadata for timeseries
- `resolution::DateTime.Period=nothing`: skip forecast that don't match this resolution.
"""
function add_forecasts!(
                        ::Type{T},
                        data::SystemData,
                        timeseries_metadata::Vector{TimeseriesFileMetadata};
                        resolution=nothing
                       ) where T <: InfrastructureSystemsType
    forecast_cache = ForecastCache()

    for metadata in timeseries_metadata
        add_forecast!(T, data, forecast_cache, metadata; resolution=resolution)
    end
end

"""
    add_forecast!(data::SystemData, forecast)

Add a forecast.

# Arguments
- `data::SystemData`: infrastructure
- `forecast`: Any object of subtype forecast

Throws ArgumentError if the forecast's component is not stored in the system.

"""
function add_forecast!(
                       data::SystemData,
                       component::InfrastructureSystemsType,
                       forecast::Forecast,
                      )
    ts_data = TimeSeriesData(get_data(forecast))
    forecast_internal = make_internal_forecast(forecast, ts_data)
    add_forecast!(data, component, forecast_internal, ts_data)
end

function add_forecast!(
                       data::SystemData,
                       component::InfrastructureSystemsType,
                       forecast::ForecastInternal,
                       ts_data::TimeSeriesData,
                      )
    _validate_component(data, component)
    check_add_forecast!(data.forecast_metadata, forecast)
    add_forecast!(component, forecast)
    # TODO: can this be atomic with forecast addition?
    add_time_series!(data.time_series_storage, get_uuid(component), get_label(forecast),
                     ts_data)
end

"""
    add_forecast!(
                  data::SystemData,
                  filename::AbstractString,
                  component::InfrastructureSystemsType,
                  label::AbstractString,
                  scaling_factor::Union{String, Float64}=1.0,
                 )

Add a forecast from a CSV file.

See [`TimeseriesFileMetadata`](@ref) for description of scaling_factor.
"""
function add_forecast!(
                       data::SystemData,
                       filename::AbstractString,
                       component::InfrastructureSystemsType,
                       label::AbstractString,
                       scaling_factor::Union{String, Float64}=1.0,
                      )
    component_name = get_name(component)
    ts = read_time_series(filename, component_name)
    timeseries = ts[Symbol(component_name)]
    _add_forecast!(data, component, label, timeseries, scaling_factor)
end

"""
    add_forecast!(
                  data::SystemData,
                  ta::TimeSeries.TimeArray,
                  component::InfrastructureSystemsType,
                  label::AbstractString,
                  scaling_factor::Union{String, Float64}=1.0,
                 )

Add a forecast to a system from a TimeSeries.TimeArray.

See [`TimeseriesFileMetadata`](@ref) for description of scaling_factor.
"""
function add_forecast!(
                       data::SystemData,
                       ta::TimeSeries.TimeArray,
                       component::InfrastructureSystemsType,
                       label::AbstractString,
                       scaling_factor::Union{String, Float64}=1.0,
                      )
    timeseries = ta[Symbol(get_name(component))]
    _add_forecast!(data, component, label, timeseries, scaling_factor)
end

"""
    add_forecast!(
                  data::SystemData,
                  df::DataFrames.DataFrame,
                  component::InfrastructureSystemsType,
                  label::AbstractString,
                  scaling_factor::Union{String, Float64}=1.0;
                  timestamp=:timestamp,
                 )

Add a forecast to a system from a DataFrames.DataFrame.

See [`TimeseriesFileMetadata`](@ref) for description of scaling_factor.
"""
function add_forecast!(
                       data::SystemData,
                       df::DataFrames.DataFrame,
                       component::InfrastructureSystemsType,
                       label::AbstractString,
                       scaling_factor::Union{String, Float64}=1.0;
                       timestamp=:timestamp,
                      )
    timeseries = TimeSeries.TimeArray(df; timestamp=timestamp)
    add_forecast!(data, timeseries, component, label, scaling_factor)
end

function add_forecast!(
                       ::Type{T},
                       data::SystemData,
                       forecast_cache::ForecastCache,
                       metadata::TimeseriesFileMetadata;
                       resolution=nothing,
                      ) where T <: InfrastructureSystemsType
    set_component!(metadata, data, InfrastructureSystems)
    component = metadata.component

    forecast, ts_data = make_forecast!(forecast_cache, metadata; resolution=resolution)
    if !isnothing(forecast)
        add_forecast!(data, component, forecast, ts_data)
    end
end

"""
    remove_forecast!(
                     ::Type{T},
                     data::SystemData,
                     component::InfrastructureSystemsType,
                     initial_time::Dates.DateTime,
                     label::String,
                    ) where T <: Forecast

Remove the time series data for a component.
"""
function remove_forecast!(
                          ::Type{T},
                          data::SystemData,
                          component::InfrastructureSystemsType,
                          initial_time::Dates.DateTime,
                          label::String,
                         ) where T <: Forecast
    type_ = forecast_external_to_internal(T)
    forecast = get_forecast(type_, component, initial_time, label)
    uuid = get_time_series_uuid(forecast)
    # TODO: can this be atomic?
    remove_forecast_internal!(type_, component, initial_time, label)
    remove_time_series!(data.time_series_storage, uuid, get_uuid(component), label)
end

"""
    make_forecast!(timeseries_metadata::TimeseriesFileMetadata;
                   resolution=nothing)

Return a vector of forecasts from TimeseriesFileMetadata.

# Arguments
- `timeseries_metadata::TimeseriesFileMetadata`: metadata
- `resolution::{Nothing, Dates.Period}`: skip any forecasts that don't match this resolution
"""
function make_forecast!(
                        forecast_cache::ForecastCache,
                        timeseries_metadata::TimeseriesFileMetadata;
                        resolution=nothing,
                       )
    forecast_info = add_forecast_info!(forecast_cache, timeseries_metadata)
    return _make_forecast(forecast_info, resolution)
end

function _add_forecast!(
                        data::SystemData,
                        component::InfrastructureSystemsType,
                        label::AbstractString,
                        timeseries::TimeSeries.TimeArray,
                        scaling_factor,
                       )
    timeseries = handle_scaling_factor(timeseries, scaling_factor)
    # TODO: This code path needs to accept a metdata file or parameters telling it which
    # type of forecast to create.
    ts_data = TimeSeriesData(timeseries)
    forecast = DeterministicInternal(label, ts_data)
    add_forecast!(data, component, forecast, ts_data)
end

function _make_forecasts(forecast_cache::ForecastCache, resolution)
    forecasts = Vector{Forecast}()

    for forecast_info in forecast_cache.forecasts
        forecast = _make_forecast(forecast_info)
        if !isnothing(forecast)
            push!(forecasts, forecast)
        end
    end

    return forecasts
end

function _make_forecast(forecast_info::ForecastInfo, resolution)
    len = length(forecast_info.data)
    @assert len >= 2
    timestamps = TimeSeries.timestamp(forecast_info.data)
    res = timestamps[2] - timestamps[1]
    if !isnothing(resolution) && res != resolution
        @debug "Skip forecast with resolution=$res; doesn't match user=$resolution"
        return nothing, nothing
    end

    timeseries = forecast_info.data[Symbol(get_name(forecast_info.component))]
    timeseries = handle_scaling_factor(timeseries, forecast_info.scaling_factor)
    forecast_type = get_forecast_type(forecast_info)
    ts_data = TimeSeriesData(timeseries)
    forecast = forecast_type(forecast_info.label, ts_data)
    @debug "Created $forecast"
    return forecast, ts_data
end

function add_forecast_info!(forecast_cache::ForecastCache, metadata::TimeseriesFileMetadata)
    timeseries = _add_forecast_info!(forecast_cache, metadata.data_file,
                                     metadata.component_name)
    forecast_info = ForecastInfo(metadata, timeseries)
    @debug "Added ForecastInfo" metadata
    return forecast_info
end

"""
    are_forecasts_contiguous(data::SystemData)

Return true if forecasts are stored contiguously.

Throws ArgumentError if there are no forecasts stored.
"""
function are_forecasts_contiguous(data::SystemData)
    for component in iterate_components_with_forecasts(data.components)
        if has_forecasts(component)
            return are_forecasts_contiguous(component)
        end
    end

    throw(ArgumentError("no forecasts are stored"))
end

"""
    generate_initial_times(
                           data::SystemData,
                           interval::Dates.Period,
                           horizon::Int;
                           initial_time::Union{Nothing, Dates.DateTime}=nothing,
                          )

Generates all possible initial times for the stored forecasts. This should return the same
result regardless of whether the forecasts have been stored as one contiguous array or
chunks of contiguous arrays, such as one 365-day forecast vs 365 one-day forecasts.

Throws ArgumentError if there are no forecasts stored, interval is not a multiple of the
system's forecast resolution, or if the stored forecasts have overlapping timestamps.

# Arguments
- `data::SystemData`: system
- `interval::Dates.Period`: Amount of time in between each initial time.
- `horizon::Int`: Length of each forecast array.
- `initial_time::Union{Nothing, Dates.DateTime}=nothing`: Start with this time. If nothing,
  use the first initial time.
"""
function generate_initial_times(
                                data::SystemData,
                                interval::Dates.Period,
                                horizon::Int;
                                initial_time::Union{Nothing, Dates.DateTime}=nothing,
                               )
    for component in iterate_components_with_forecasts(data.components)
        if has_forecasts(component)
            return generate_initial_times(component, interval, horizon;
                                          initial_time=initial_time)
        end
    end

    throw(ArgumentError("no forecasts are stored"))
end

"""
Checks that the component exists in data and the UUID's match.
"""
function _validate_component(
                             data::SystemData,
                             component::T,
                            ) where T <: InfrastructureSystemsType
    comp = get_component(T, data.components, get_name(component))
    if isnothing(comp)
        throw(ArgumentError("no $T with name=$(get_name(component)) is stored"))
    end

    user_uuid = get_uuid(component)
    ps_uuid = get_uuid(comp)
    if user_uuid != ps_uuid
        throw(ArgumentError(
            "comp UUID doesn't match, perhaps it was copied?; " *
            "$T name=$(get_name(component)) user=$user_uuid system=$ps_uuid"
        ))
    end
end

function get_component_types_raw(::Type{SystemData}, raw::NamedTuple)
    return get_component_types_raw(Components, raw.components)
end

function get_components_raw(
                            ::Type{SystemData},
                            ::Type{T},
                            raw::NamedTuple,
                           ) where T <: InfrastructureSystemsType
    return get_components_raw(Components, T, raw.components)
end

function compare_values(x::SystemData, y::SystemData)::Bool
    match = true
    for key in keys(x.components.data)
        if !compare_values(x.components.data[key], y.components.data[key])
            @debug "System components do not match"
            match = false
        end
    end

    if !compare_values(x.forecast_metadata, y.forecast_metadata)
        @debug "System forecasts do not match"
        match = false
    end

    return match
end

function remove_component!(::Type{T}, data::SystemData, name) where T
    return remove_component!(T, data.components, name)
end

function remove_component!(data::SystemData, component)
    remove_component!(data.components, component)
end

function remove_components!(::Type{T}, data::SystemData) where T
    return remove_components!(T, data.components)
end

function clear_forecasts!(data::SystemData)
    clear_forecasts!(data.components)
    clear_time_series!(data.time_series_storage)
end

function iterate_forecasts(data::SystemData)
    Channel() do channel
        for component in iterate_components_with_forecasts(data.components)
            for forecast in iterate_forecasts(component)
                time_series = get_time_series(data.time_series_storage,
                                              get_time_series_uuid(forecast))
                put!(channel, make_public_forecast(forecast, time_series))
            end
        end
    end
end

"""
Return the time delta between the first two stored forecasts.
if less than two are stored, return Dates.Second(0).
"""
function get_forecasts_interval(data::SystemData)
    initial_times = get_forecast_initial_times(data)
    if length(initial_times) <= 1
        return Dates.Second(0)
    end

    return initial_times[2] - initial_times[1]
end

"""
    get_forecast_counts(data::SystemData)

Return a tuple of counts of components with forecasts and total forecasts.
"""
function get_forecast_counts(data::SystemData)
    component_count = 0
    forecast_count = 0
    for component in iterate_components_with_forecasts(data.components)
        component_count += 1
        forecast_count += get_num_forecasts(component)
    end

    return component_count, forecast_count
end

"""
    prepare_for_serialization!(data::SystemData, filename::AbstractString)

Parent object should call this prior to serialization so that SystemData can store the
appropriate path information for the time series data.
"""
function prepare_for_serialization!(data::SystemData, filename::AbstractString)
    dir = dirname(filename)
    base = splitext(basename(filename))[1]
    data.time_series_storage_file = joinpath(dir, base * "_" * TIME_SERIES_STORAGE_FILE)
end

function JSON2.write(io::IO, data::SystemData)
    return JSON2.write(io, encode_for_json(data))
end

function JSON2.write(data::SystemData)
    return JSON2.write(encode_for_json(data))
end

function encode_for_json(data::SystemData)
    if isnothing(data.time_series_storage_file)
        data.time_series_storage_file = TIME_SERIES_STORAGE_FILE
    end

    json_data = Dict()
    for field in (:components, :forecast_metadata, :validation_descriptors,
                  :time_series_storage_file, :internal)
        json_data[string(field)] = getfield(data, field)
    end

    serialize(data.time_series_storage, data.time_series_storage_file)
    json_data["time_series_storage_type"] = string(typeof(data.time_series_storage))
    return json_data
end

function JSON2.read(io::IO, ::Type{SystemData})
    raw = JSON2.read(io, NamedTuple)
    sys = deserialize(SystemData, InfrastructureSystemsType, raw)
    return sys
end

function deserialize(
                     ::Type{SystemData},
                     ::Type{T},
                     raw::NamedTuple,
                    ) where T <: InfrastructureSystemsType
    forecast_metadata = convert_type(ForecastMetadata, raw.forecast_metadata)

    if strip_module_name(raw.time_series_storage_type) == "InMemoryTimeSeriesStorage"
        hdf5_storage = Hdf5TimeSeriesStorage(raw.time_series_storage_file)
        time_series_storage = InMemoryTimeSeriesStorage(hdf5_storage)
    else
        time_series_storage = from_file(Hdf5TimeSeriesStorage, raw.time_series_storage_file)
    end

    # OPT: This looks odd and is wasteful.
    # JSON2 creates NamedTuples recursively. JSON creates dicts, which is what we need.
    # Could be optimized.
    validation_descriptors = JSON.parse(JSON2.write(raw.validation_descriptors))

    internal = convert_type(InfrastructureSystemsInternal, raw.internal)
    sys = SystemData(forecast_metadata, validation_descriptors, time_series_storage,
                     internal)
    deserialize_components(T, sys, raw)
    return sys
end

"""
Deserializes components defined in InfrastructureSystems. Parent modules should override
this by changing the component type and module.
"""
function deserialize_components(
                                ::Type{InfrastructureSystemsType},
                                sys::SystemData,
                                raw::NamedTuple,
                               )
    for c_type_sym in get_component_types_raw(SystemData, raw)
        c_type = getfield(InfrastructureSystems,
                          Symbol(strip_module_name(string(c_type_sym))))
        for component in get_components_raw(SystemData, c_type, raw)
            comp = convert_type(c_type, component)
            add_component!(sys, comp)
        end
    end

    return
end

# Redirect functions to Components and Forecasts

add_component!(data::SystemData, component; kwargs...) = add_component!(
    data.components, component; kwargs...
)
iterate_components(data::SystemData) = iterate_components(data.components)

get_component(::Type{T}, data::SystemData, args...) where T = get_component(
    T, data.components, args...
)
get_components(::Type{T}, data::SystemData) where T = get_components(T, data.components)
get_components_by_name(::Type{T}, data::SystemData, args...) where T =
    get_components_by_name(T, data.components, args...)

#get_component_forecasts(::Type{T}, data::SystemData, args...) where T =
#    get_component_forecasts(T, data.forecasts, args...)
#get_forecasts(::Type{T}, data::SystemData, component, args...) where T = get_forecasts(
#    T, component, args...
#)
get_forecast_initial_times(data::SystemData) = get_forecast_initial_times(data.components)
get_forecasts_initial_time(data::SystemData) = get_forecasts_initial_time(data.components)
get_forecasts_last_initial_time(data::SystemData) = get_forecasts_last_initial_time(data.components)
get_forecasts_horizon(data::SystemData) = get_forecasts_horizon(data.forecast_metadata)
get_forecasts_resolution(data::SystemData) = get_forecasts_resolution(data.forecast_metadata)
clear_components!(data::SystemData) = clear_components!(data.components)
set_component!(metadata::TimeseriesFileMetadata, data::SystemData, mod::Module) =
    set_component!(metadata, data.components, mod)
check_forecast_consistency(data::SystemData) = check_forecast_consistency(data.components)
validate_forecast_consistency(data::SystemData) = validate_forecast_consistency(data.components)


text:
include("src/RK.jl")
include("src/parameters.jl")

# Если аргументы не переданы или переданы некорректно, то рисуем линзу Люнеберга
if (length(ARGS) > 0 && ARGS[1] == "maxwell")
  lens = "maxwell"
  include("src/maxwell.jl")
else
  lens = "luneburg"
  include("src/luneburg.jl")
end

include("src/polar_eikonal.jl")

# Создаем директорию, куда будем сохранять результаты расчетов
isdir(dir[lens]) ? true : mkdir(dir[lens])

# x_0, y_0, p1_0, p2_0
xn_0 = Vector{Float64}(undef, 4)

# Источник точечный. Из него под углом α расходятся лучи
for (i, α) in enumerate(range(α_min, stop=α_max, length=30))
  # Лучи исходят из одной точки, расходясь под углом α к оси Ох
  xn_0[1:2] = [r_0, phi_0]
  xn_0[3:4] = [-cos(α)*n(r_0), -sin(α)*n(r_0)]

  tn, xn = RK.RKp6n1(eikonal, xn_0, h, t_start, t_stop)

  file = open("./$(dir[lens])/data$i.txt", "w")
  for (t, x, y) in zip(tn, xn[:, 1], xn[:, 2])
    write(file, "$t,$x,$y\n")
  end
  close(file)

  # для отображения прогресса
  print("=")
end
println()

ext: jl


text:
using Gnuplot, RDatasets
dset = dataset("datasets", "iris")
byCat = dset.Species
categ = unique(byCat)
# aesthetics
@gp linetypes(:tab10) :-
@gp :- xlab ="Sepal Length" ylab = "Sepal Width" "set auto fix"
@gp :- title = "Iris dataset" :-
@gp :- "set offsets graph .05, graph .05, graph .05, graph .05" :-
@gp :- "set key b r font ',12' tit 'Species' box" :-
# the actual plot
for c in categ
    indc = findall(x->x == c, byCat)
    @gp :-  dset.SepalLength[indc] dset.SepalWidth[indc]  "w p tit '$(c)' pt 7 ps 1.4 "
end
save(term ="pngcairo size 600,400", output ="Categorical_iris.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
#' ---
#' title: Power Systems
#' ---

#' **Originally Contributed by**: Yury Dvorkin and Miles Lubin

#' This notebook demonstrates how to formulate basic power systems engineering models in JuMP using a 3 bus example.
#' We will consider basic "economic dispatch" and "unit commitment" models without taking into account transmission constraints.

#' ## Illustrative example
#' In the following notes for the sake of simplicity, we are going to use a three bus example mirroring the interface between Western and Eastern Texas. This example is taken from R. Baldick, "[Wind and Energy Markets: A Case Study of Texas](http://dx.doi.org/10.1109/JSYST.2011.2162798)," IEEE Systems Journal, vol. 6, pp. 27-34, 2012.

#'

#' For this example, we set the following characteristics of generators, transmission lines, wind farms and demands:

#' | Quantity         | Generator 1 | Generator 2 |
#' |:-----------------|:-----------:|------------:|
#' | $g_{min}$, MW    | 0           | 300         |
#' | $g_{max}$, MW    | 1000        | 1000        |
#' | $c^g$, \$/MWh    | 50          | 100         |
#' | $c^{g0}$, \$/MWh | 1000        | 0           |

#' | Quantity      | Line 1 | Line 2 |
#' |:--------------|:------:|-------:|
#' | $f^{max}$, MW | 100    | 1000   |
#' | x, p.u.       | 0.001  | 0.001  |

#' | Quantity        | Wind farm 1 | Wind farm 2 |
#' |:----------------|:-----------:|------------:|
#' | $w^{f}$, MW     | 150         | 50          |
#' | $c^{w}$, \$/MWh | 50          | 50          |

#' | Quantity | Bus 1 | Bus 2 | Bus 3 |
#' |:---------|:-----:|:-----:|------:|
#' | $d$, MW  | 0     | 0     | 15000 |

#' ## Economic dispatch
#' Economic dispatch (ED) is an optimization problem that minimizes the cost of supplying energy demand subject to operational constraints on power system assets. In its simplest modification, ED is an LP problem solved for an aggregated load and wind forecast and for a single infinitesimal moment. Mathematically, the ED problem can be written as follows:

#' $$
#' \min \sum_{i \in I} c^g_{i} \cdot g_{i} + c^w \cdot w,
#' $$

#' where $c_{i}$ and $g_{i}$ are the incremental cost ($\$/MWh$) and power output ($MW$) of the $i^{th}$ generator, respectively, and $c^w$ and $w$ are the incremental cost ($\$/MWh$) and wind power injection ($MW$), respectively.

#' Subject to the constraints:

#'
 Minimum ($g^{\min}$) and maximum ($g^{\max}$) limits on power outputs of generators:

#' $$
#' g^{\min}_{i} \leq g_{i} \leq g^{\max}_{i}.
#' $$

#'
Constraint on the wind power injection:

#' $$
#' 0 \leq w \leq w^f,
#' $$

#' where $w$ and $w^f$ are the wind power injection and wind power forecast, respectively.

#'
Power balance constraint:

#' $$
#' \sum_{i \in I} g_{i} + w = d^f,
#' $$

#' where $d^f$ is the demand forecast.

#' Further reading on ED models can be found in A. J. Wood, B. F. Wollenberg, and G. B. Sheblé, "Power Generation, Operation and Control", Wiley, 2013.

#' ## JuMP Implementation of Economic Dispatch

using JuMP, GLPK, LinearAlgebra, DataFrames

#+

# Define some input data about the test system
# Maximum power output of generators
g_max = [1000, 1000];
# Minimum power output of generators
g_min = [0, 300];
# Incremental cost of generators
c_g = [50, 100];
# Fixed cost of generators
c_g0 = [1000, 0]
# Incremental cost of wind generators
c_w = 50;
# Total demand
d = 1500;
# Wind forecast
w_f = 200;

#+

# In this cell we create function solve_ed, which solves the economic dispatch problem for a given set of input parameters.
function solve_ed(g_max, g_min, c_g, c_w, d, w_f)
    #Define the economic dispatch (ED) model
    ed = Model(GLPK.Optimizer)

    # Define decision variables
    @variable(ed, 0 <= g[i = 1:2] <= g_max[i]) # power output of generators
    @variable(ed, 0 <= w <= w_f) # wind power injection

    # Define the objective function
    @objective(ed, Min, dot(c_g, g) + c_w * w)

    # Define the constraint on the maximum and minimum power output of each generator
    @constraint(ed, [i = 1:2], g[i] <= g_max[i]) #maximum
    @constraint(ed, [i = 1:2], g[i] >= g_min[i]) #minimum

    # Define the constraint on the wind power injection
    @constraint(ed, w <= w_f)

    # Define the power balance constraint
    @constraint(ed, sum(g) + w == d)

    # Solve statement
    optimize!(ed)

    # return the optimal value of the objective function and its minimizers
    return value.(g), value(w), w_f - value(w), objective_value(ed)
end

# Solve the economic dispatch problem
(g_opt, w_opt, ws_opt, obj) = solve_ed(g_max, g_min, c_g, c_w, d, w_f);

println("\n")
println("Dispatch of Generators: ", g_opt, " MW")
println("Dispatch of Wind: ", w_opt, " MW")
println("Wind spillage: ", w_f - w_opt, " MW")
println("\n")
println("Total cost: ", obj, "\$")

#' ### Economic dispatch with adjustable incremental costs

#' In the following exercise we adjust the incremental cost of generator G1 and observe its impact on the total cost.

c_g_scale_df = DataFrame(Symbol("Dispatch of Generator 1(MW)") => Float64[],
               Symbol("Dispatch of Generator 2(MW)") => Float64[],
               Symbol("Dispatch of Wind(MW)") => Float64[],
               Symbol("Spillage of Wind(MW)") => Float64[],
               Symbol("Total cost(\$)") => Float64[])
for c_g1_scale = 0.5:0.1:3.0
    c_g_scale = [c_g[1] * c_g1_scale, c_g[2]] # update the incremental cost of the first generator at every iteration
    local g_opt, w_opt, ws_opt, obj = solve_ed(g_max, g_min, c_g_scale, c_w, d, w_f) # solve the ed problem with the updated incremental cost
    push!(c_g_scale_df, (g_opt[1], g_opt[2], w_opt, ws_opt, obj))
end

#+

ENV["COLUMNS"]=250 # Helps us display the complete table
c_g_scale_df

#' ## Modifying the JuMP model in place

#' Note that in the previous exercise we entirely rebuilt the optimization model at every iteration of the internal loop, which incurs an additional computational burden. This burden can be alleviated if instead of re-building the entire model, we modify a specific constraint(s) or the objective function, as it shown in the example below.

#' Compare the computing time in case of the above and below models.

function solve_ed_inplace(c_w_scale)
    start = time()
    obj_out = Float64[]
    w_out = Float64[]
    g1_out = Float64[]
    g2_out = Float64[]

    ed = Model(GLPK.Optimizer)

    # Define decision variables
    @variable(ed, 0 <= g[i = 1:2] <= g_max[i]) # power output of generators
    @variable(ed, 0 <= w <= w_f ) # wind power injection

    # Define the objective function
    @objective(ed, Min, dot(c_g, g) + c_w * w)

    # Define the constraint on the maximum and minimum power output of each generator
    @constraint(ed, [i = 1:2], g[i] <= g_max[i]) #maximum
    @constraint(ed, [i = 1:2], g[i] >= g_min[i]) #minimum

    # Define the constraint on the wind power injection
    @constraint(ed, w <= w_f)

    # Define the power balance constraint
    @constraint(ed, sum(g) + w == d)

    optimize!(ed)

    for c_g1_scale = 0.5:0.01:3.0
        @objective(ed, Min, c_g1_scale*c_g[1]*g[1] + c_g[2]*g[2] + c_w_scale*c_w*w)
        optimize!(ed)
        push!(obj_out, objective_value(ed))
        push!(w_out, value(w))
        push!(g1_out, value(g[1]))
        push!(g2_out, value(g[2]))
    end
    elapsed = time() - start
    print(string("elapsed time: ", elapsed, " seconds"))
    return obj_out, w_out, g1_out, g2_out
end

solve_ed_inplace(2.0);

#' Adjusting specific constraints and/or the objective function is faster than re-building the entire model.

#' ## A few practical limitations of the economic dispatch model

#' ### Inefficient usage of wind generators

#' The economic dispatch problem does not perform commitment decisions and, thus, assumes that all generators must be dispatched at least at their minimum power output limit. This approach is not cost efficient and may lead to absurd decisions. For example, if $ d = \sum_{i \in I} g^{\min}_{i}$, the wind power injection must be zero, i.e. all available wind generation is spilled, to meet the minimum power output constraints on generators.

#' In the following example, we adjust the total demand and observed how it affects wind spillage.

demandscale_df = DataFrame(Symbol("Dispatch of Generators(MW)") => Float64[],
               Symbol("Dispatch of Generator 2(MW)") => Float64[],
               Symbol("Dispatch of Wind(MW)") => Float64[],
               Symbol("Spillage of Wind(MW)") => Float64[],
               Symbol("Total cost(\$)") => Float64[])

for demandscale = 0.2:0.1:1.5
    local g_opt,w_opt,ws_opt,obj = solve_ed(g_max, g_min, c_g, c_w, demandscale*d, w_f)

    push!(demandscale_df, (g_opt[1], g_opt[2], w_opt, ws_opt, obj))
end

#+

demandscale_df

#' This particular drawback can be overcome by introducing binary decisions on the "on/off" status of generators. This model is called unit commitment and considered later in these notes.

#' For further reading on the interplay between wind generation and the minimum power output constraints of generators, we refer interested readers to R. Baldick, "Wind and Energy Markets: A Case Study of Texas," IEEE Systems Journal, vol. 6, pp. 27-34, 2012.

#' ### Transmission-infeasible solution

#' The ED solution is entirely market-based and disrespects limitations of the transmission network. Indeed, the flows in transmission lines would attain the following values:

#' $$f_{1-2} = 150 MW \leq f_{1-2}^{\max} = 100 MW $$

#' $$f_{2-3} = 1200 MW \leq f_{2-3}^{\max} = 1000 MW $$

#' Thus, if this ED solution was enforced in practice, the power flow limits on both lines would be violated. Therefore, in the following section we consider the optimal power flow model, which amends the ED model with network constraints.

#' The importance of the transmission-aware decisions is emphasized in E. Lannoye, D. Flynn, and M. O'Malley, "Transmission, Variable Generation, and Power System Flexibility," IEEE Transactions on Power Systems, vol. 30, pp. 57-66, 2015.

#' ## Unit Commitment model
#' The Unit Commitment (UC) model can be obtained from ED model by introducing binary variable associated with each generator. This binary variable can attain two values: if it is "1", the generator is synchronized and, thus, can be dispatched, otherwise, i.e. if the binary variable is "0", that generator is not synchronized and its power output is set to 0.

#' To obtain the mathematical formulation of the UC model, we will modify the constraints of the ED model as follows:

#' $$
#' g^{\min}_{i} \cdot u_{t,i} \leq g_{i} \leq g^{\max}_{i} \cdot u_{t,i},
#' $$

#' where $ u_{i} \in \{0;1\}. $ In this constraint, if $ u_{i} = 0$, then $g_{i}  = 0$. On the other hand, if $ u_{i} = 1$, then $g^{max}_{i} \leq g_{i}   \leq g^{min}_{i}$.

#' For further reading on the UC problem we refer interested readers to G. Morales-Espana, J. M. Latorre, and A. Ramos, "Tight and Compact MILP Formulation for the Thermal Unit Commitment Problem," IEEE Transactions on Power Systems, vol. 28, pp. 4897-4908, 2013.
#' In the following example we convert the ED model explained above to the UC model.

# In this cell we introduce binary decision u to the economic dispatch problem (function solve_ed)
function solve_uc(g_max, g_min, c_g, c_w, d, w_f)
    #Define the unit commitment (UC) model
    uc = Model(GLPK.Optimizer)

    # Define decision variables
    @variable(uc, 0 <= g[i=1:2] <= g_max[i]) # power output of generators
    @variable(uc, u[i = 1:2], Bin) # Binary status of generators
    @variable(uc, 0 <= w <= w_f ) # wind power injection

    # Define the objective function
    @objective(uc, Min, dot(c_g, g) + c_w * w)

    # Define the constraint on the maximum and minimum power output of each generator
    @constraint(uc, [i = 1:2], g[i] <= g_max[i]) #maximum
    @constraint(uc, [i = 1:2], g[i] >= g_min[i]) #minimum

    # Define the constraint on the wind power injection
    @constraint(uc, w <= w_f)

    # Define the power balance constraint
        @constraint(uc, sum(g) + w == d)

    # Solve statement
    optimize!(uc)

    status = termination_status(uc)
    if status != MOI.OPTIMAL
        return status, zeros(length(g)), 0.0, 0.0, zeros(length(u)), Inf
    end
    return status, value.(g), value(w), w_f - value(w), value.(u), objective_value(uc)
end

# Solve the economic dispatch problem
status, g_opt, w_opt, ws_opt, u_opt, obj = solve_uc(g_max, g_min, c_g, c_w, d, w_f);

println("\n")
println("Dispatch of Generators: ", g_opt[:], " MW")
println("Commitments of Generators: ", u_opt[:])
println("Dispatch of Wind: ", w_opt, " MW")
println("Wind spillage: ", w_f - w_opt, " MW")
println("\n")
println("Total cost: ", obj, "\$")

#' ### Unit Commitment as a function of demand

#' After implementing the UC model, we can now assess the interplay between the minimum power output constraints on generators and wind generation.

uc_df = DataFrame(Symbol("Commitment of Generator 1(MW)") => Float64[],
               Symbol("Commitment of Generator 2(MW)") => Float64[],
               Symbol("Dispatch of Generator 1(MW)") => Float64[],
               Symbol("Dispatch of Generator 2(MW)") => Float64[],
               Symbol("Dispatch of Wind(MW)") => Float64[],
               Symbol("Spillage of Wind(MW)") => Float64[],
               Symbol("Total cost(\$)") => Float64[])

for demandscale = 0.2:0.1:1.5
    local status, g_opt, w_opt, ws_opt, u_opt, obj = solve_uc(g_max, g_min, c_g, c_w, demandscale*d, w_f)

    if status == MOI.OPTIMAL
        push!(uc_df, (u_opt[1], u_opt[2], g_opt[1], g_opt[2], w_opt, ws_opt, obj))
    else
        println("Status: $status for demandscale = $demandscale \n")
    end
end

#+

uc_df


text:
module motion_planning

using LinearAlgebra, ForwardDiff, FiniteDiff, StaticArrays, SparseArrays
using MathOptInterface, Ipopt
using Distributions, Interpolations
using JLD2
using BenchmarkTools
using ModelingToolkit

using Colors
using CoordinateTransformations
using FileIO
using GeometryBasics
using MeshCat, MeshIO
using Rotations
# using RigidBodyDynamics, MeshCatMechanisms
Point
include("indices.jl")
include("utils.jl")

include("time.jl")
include("model.jl")
include("integration.jl")

include("problem.jl")

include("objective.jl")
include("objectives/quadratic.jl")
include("objectives/penalty.jl")

include("constraints.jl")
include("constraints/dynamics.jl")

include("moi.jl")
include("solvers/snopt.jl")
include("solvers/newton.jl")
include("solvers/levenberg_marquardt.jl")

include("lqr.jl")
include("unscented.jl")

# direct policy optimization
function include_dpo()
    include(joinpath(pwd(), "src/direct_policy_optimization/dpo.jl"))
end

# differential dynamic programming
function include_ddp()
    include(joinpath(pwd(), "src/differential_dynamic_programming/ddp.jl"))
end

end # module


text:
function m_sqrt(num)
    guess = num / 2
    tmp = 0
    while true
        tmp = 0.5 * (guess + (num / guess))
        if tmp != guess
            guess = tmp
        else
            return guess
        end
    end
end

print(m_sqrt(10), '\n')

max_stars_repo_name: zmaas/scratch

max_issues_repo_name: zmaas/scratch

max_forks_repo_name: zmaas/scratch


text:
"""
    XNLI

A collection of 5,000 test and 2,500 dev pairs for the `MultiNLI` corpus.

For details, see the [2018 paper](https://www.aclweb.org/anthology/papers/D/D18/D18-1269/) "XNLI: Evaluating Cross-lingual Sentence Representations."

```julia
XNLI.dev_tsv()
XNLI.dev_jsonl()
XNLI.test_tsv()
XNLI.test_jsonl()
```
"""
module XNLI

import ...register_data
using DataDeps: @datadep_str

xnlifile(filename) = joinpath(datadep"XNLI", "XNLI-1.0", filename)

dev_tsv()    = xnlifile("xnli.dev.tsv")
dev_jsonl()  = xnlifile("xnli.dev.jsonl")
test_tsv()   = xnlifile("xnli.test.tsv")
test_jsonl() = xnlifile("xnli.test.jsonl")

function __init__()
    register_data(
        "XNLI",
        """
        The Cross-lingual Natural Language Inference (XNLI) corpus is a crowd-sourced collection of 5,000 test and 2,500 dev pairs for the `MultiNLI` corpus.
        """,
        "http://www.nyu.edu/projects/bowman/xnli/XNLI-1.0.zip",
        "4ba1d5e1afdb7161f0f23c66dc787802ccfa8a25a3ddd3b165a35e50df346ab1")
end

end


text:
# This file includes an example file for HInterp1D interpolation.

using FractalTools
using Makie

# Generate data
xi, dx, xf = 0., 0.1, 10.
x = collect(xi : dx : xf)
y = sin.(x)
z = cos.(x)
pts = collect.(zip(x, y, z))
freevars = [0.01 * rand() * ones(2,2) for i in 1 : length(x) - 1]
interp = interpolate(pts, HInterp1D(freevars))

# Calculate interpolant
xt = collect(xi : 0.1dx : xf)
_yt = interp.(xt)
yt = getindex.(_yt, 1)

# Plot interpolation
fig, ax, plt = lines(xt, yt, color=yt, linewidth=5)
scatter!(ax, x, y, markersize=5, color=:black)
fig


text:
module LoadAllPackages

using Base: PkgId, UUID
using Pkg: TOML

function as_project_file(project)
    isfile(project) && return project
    isdir(project) || error("`$project` is not a file or directory")
    candidates = joinpath.(project, ("JuliaProject.toml", "Project.toml"))
    i = findfirst(isfile, candidates)
    if i === nothing
        error("`$project` does not have `JuliaProject.toml` or `Project.toml`")
    end
    return candidates[i]
end

"""
    LoadAllPackages.loadall(project = Base.active_project())

Load all packages in `project`.
"""
function loadall(project = Base.active_project())
    @info "Loading packages from `$project`..."
    project = as_project_file(project)
    deps = get(TOML.parsefile(project), "deps", nothing)
    deps === nothing && return
    @info "Loading $(length(deps)) packages..."
    t0 = time_ns()
    for (name, uuid) in deps
        pkg = PkgId(UUID(uuid), name)
        @debug "Loading `$pkg`..."
        Base.require(pkg)
        @debug "`$pkg` loaded."
    end
    sec = (time_ns() - t0) / 10^9
    @info "$(length(deps)) packages loaded in $sec seconds"
    return
end

_loadall(cwd, path) = loadall(joinpath(cwd, path))

"""
    LoadAllPackages.@loadall(project)

Load all packages in `project`.  If `project` is a relative path,
interpret it as the relative path from `@__DIR__` rather than `pwd()`.
"""
macro loadall(project)
    :(_loadall($Base.@__DIR__, $(esc(project))))
end

end


text:
using Markdown # hide

# ## Parallel construction and parallel search
# By Eric S. Téllez
#

ENV["GKSwstype"] = "nul" # hide

using SimilaritySearch
# Similarity search on very large datasets and high-dimensional datasets require high computational
# resources. In this example we show how to parallelize both the construction and search to be able to
# handle this kind of databases.

const dim = 16

db = MatrixDatabase(randn(Float32, dim, 10^5))
Q = MatrixDatabase(randn(Float32, dim, 30))
dist = SqL2Distance()
G = SearchGraph(; dist, db)

# The `SearchGraph` construction algorithm is incremental:
# - If the index is empty, an element is inserted just inserting it into the index
# - If the index is not empty, the element is inserted and connected to its nearest neighbors (looking into the current index)
#
# the parallel construction is made with `index!` or `append!`, for this matter these functions accept
# a `parallel_block` argument, that controls how many elements are inserted at once, i.e., looking for its nearest neighbors
# in parallel and connected also in parallel.
#
# As in the sequential version, a minimum number of elements must exists to work, and therefore, the `parallel_minimum_first_block`
# argument can also be specified. By default, it is equal to `parallel_block`.
# The `parallel_block` argument should be set to at least the number of available threads, and perhaps multiplying it by
# a small constant is also a good approach.
#
# This example didn't run in parallel due to the document generation pipeline but all demonstrations actually do it.
# Nonetheless, they ran with the parallel API in single thread mode.
index!(G; parallel_block=512)
nothing # hide

# Note that you can't call `push!`, `append!`, or `index!` from several threads, instead you must use the `parallel_block` argument
# to let the algorithm take advantage of the multiple threads.
#
# ## Searching
# Once the index is constructed, you can solve batches in parallel and also single queries.
# In contrast with append, these functions can be called in multithreading algorithms.
# However, you must pause the searching requests while perform insertions (parallel or sequential).

I, D = searchbatch(G, Q, 10; parallel=true)

Threads.@threads for i in eachindex(Q)
    res, cost = search(G, Q[i], KnnResult(10))
    println(res.id) # do something with `res`
end

# ## About searching pools
# Searching and also insertion methods (since they also perform searches) make use of several objects that are cached
# and reused to reduce memory allocations. These caches are named pools and are of type `SearchGraphPools`, their
# default values should be enough for multithreading applications, but special usage will require to review them.
# Among the special usages that require an explicit handling of pools we can think on objects indexing objects
# that are also indexes, or distance functions that can perform context switches as tasks (?).
#
#- Using BeutifulMakie style
# ## Dependencies
using Pkg # hide
Pkg.status(["SimilaritySearch"]) # hide

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
"Z-normalize the input array"
function znorm(x::AbstractVector)
    x = x .- mean(x)
    x ./= std(x, mean=0, corrected=false)
end

function znorm!(x::AbstractVector)
    x .= x .- mean(x)
    x ./= std(x, mean=0, corrected=false)
end

"""
    meanstd(x::AbstractArray)

Return mean and std of `x`. Is faster but less accurate than built-in functions.
"""
function meanstd(x::AbstractArray{T}) where T
    s = ss = zero(float(T))
    n = length(x)
    @avx for i in eachindex(x)
        s  += x[i]
        ss += x[i]^2
    end
    m   = s/n
    sig = √(ss/n - m^2)
    m, sig
end

import Base: @propagate_inbounds, @boundscheck, getindex
import LinearAlgebra: normalize

abstract type AbstractNormalizer{T,N} <: AbstractArray{T,N} end

advance!(x) = 0 # Catches everything that's not a normalizer

setup_normalizer(n::Val{Nothing}, q, y) = n, q, y

normalize(::Val{Nothing}, q) = q

normalize(T::Type, q) = normalize(Val(T), q)

abstract type AbstractZNormalizer{T,N} <: AbstractNormalizer{T,N} end

"""
    ZNormalizer{T} <: AbstractNormalizer{T, 1}

Utility object that normalizes an input array on the fly.
Works like a vector, index into window with `z[i]`, index into normalized window with `z[!,i]`.

# Arguments:
- `x::Vector{T}`: The vecetor to operate on
- `n::Int`: The lenght of a window
- `μ::T`: mean over window
- `σ::T`: std over window
- `s::T`: sum over window
- `ss::T`: sum^2 over window
- `i::Int`: the current index
- `buffer::Vector{T}`: stores calculated normalized values, do not access this, it might not be complete
- `bufi::Int`: index of last normalized value
"""
mutable struct ZNormalizer{T} <: AbstractZNormalizer{T,1}
    x::Vector{T}
    n::Int
    μ::T
    σ::T
    s::T
    ss::T
    i::Int
    buffer::Vector{T}
    bufi::Int
end

function ZNormalizer(x::AbstractVector{T}, n) where T
    @assert length(x) >= n
    s = ss = zero(T)
    @avx for i in 1:n
        s += x[i]
        ss += x[i]^2
    end
    μ = s/n
    σ = sqrt(ss/n - μ^2)
    buffer = similar(x, n)
    ZNormalizer(x, n, μ, σ, s, ss, 0, buffer, 0)
end

function normalize(::Val{ZNormalizer}, q::AbstractVector)
    q = q .- mean(q)
    q ./= std(q, corrected=false, mean=0)
end

setup_normalizer(z::Val{ZNormalizer}, q::AbstractVector, y::AbstractVector) = z, normalize(z, q), ZNormalizer(y, length(q))

@propagate_inbounds function advance!(z::ZNormalizer{T}) where T

    if z.i == 0
        return z.i = 1
    end
    @boundscheck if z.i + z.n > length(z.x)
        return z.i += 1
    end
    z.bufi = 0

    # Remove old point
    x = z.x[z.i]
    z.s -= x
    z.ss -= x^2

    # Add new point
    x = z.x[z.i+z.n]
    z.s += x
    z.ss += x^2
    z.μ = z.s/z.n
    z.σ = sqrt(z.ss/z.n - z.μ^2)
    z.i += 1
end

@inline @propagate_inbounds function getindex(z::ZNormalizer, i)
    @boundscheck 1 <= i <= z.n || throw(BoundsError(z,i))
    xi = i+z.i-1
    @boundscheck xi <= length(z.x) || throw(BoundsError(z,i))
    z.x[xi]
end

@inline @propagate_inbounds function getindex(z::ZNormalizer, ::typeof(!), i::Int, inorder = i == z.bufi + 1)
    y = (z[i]-z.μ) / z.σ
    if inorder
        z.bufi = i
        z.buffer[i] = y
    end
    y
end

@inline @propagate_inbounds function getindex(z::ZNormalizer, ::typeof(!), i::AbstractRange)
    @boundscheck (i[1] == z.i && length(i) == z.n) || throw(ArgumentError("ZNormalizers can only be indexed by ranges corresponding to their current state. Got range $i but state was $(z.i) corresponding to range $(z.i):$(z.i+z.n-1)"))
    z
end

Statistics.mean(z::AbstractZNormalizer) = z.μ
Statistics.std(z::AbstractZNormalizer) = z.σ

SlidingDistancesBase.lastlength(z::AbstractZNormalizer) = z.n
Base.length(z::AbstractZNormalizer) = z.n
Base.size(z::ZNormalizer) = (z.n,)
actuallastlength(x::AbstractZNormalizer) = lastlength(x.x)

# Multi dim ==================================================================================
mutable struct IsoZNormalizer{T} <: AbstractZNormalizer{T,2}
    x::Array{T,2}
    n::Int
    μ::Array{T,1}
    σ::Array{T,1}
    s::Array{T,1}
    ss::Array{T,1}
    i::Int
    buffer::Array{T,2}
    bufi::Int
end

function IsoZNormalizer(x::AbstractArray{T,2}, n) where T
    @assert length(x) >= n
    s  = zeros(T, size(x,1))
    ss = zeros(T, size(x,1))
    @inbounds @simd for i in 1:n
        s  .+= x[!, i]
        ss .+= x[!, i].^2
    end
    μ = s./n
    σ = sqrt.(ss./n .- μ.^2)
    buffer = similar(x, size(x,1), n)
    IsoZNormalizer(x, n, μ, σ, s, ss, 0, buffer, 0)
end

function normalize(::Val{IsoZNormalizer}, q::AbstractMatrix)
    q = q .- mean(q, dims=2) # TODO: this will cause a ton of allocations
    q ./= std(q, dims=2, corrected=false)
end

setup_normalizer(z::Val{IsoZNormalizer}, q, y) = z, normalize(z, q), IsoZNormalizer(y, lastlength(q))
setup_normalizer(z::Val{ZNormalizer}, q::AbstractMatrix, y::AbstractMatrix) = setup_normalizer(Val(IsoZNormalizer), q, y) # Only need to expose ZNormalizer to the user

@propagate_inbounds function advance!(z::IsoZNormalizer{T}) where T

    if z.i == 0
        return z.i = 1
    end
    @boundscheck if z.i + z.n > length(z.x)
        return z.i += 1
    end
    z.bufi = 0

    # Remove old point
    x = z.x[!, z.i]
    @avx z.s .-= x
    @avx z.ss .-= x.^2

    # Add new point
    x = z.x[!, z.i+z.n]
    @avx z.s .+= x
    @avx z.ss .+= x.^2
    @avx z.μ .= z.s./z.n
    @avx z.σ .= sqrt.(z.ss./z.n .- z.μ.^2)
    z.i += 1
end

@inline @propagate_inbounds function getindex(z::IsoZNormalizer, i::Union{Number, AbstractRange}, j)
    @boundscheck 1 <= j <= z.n || throw(BoundsError(z,j))
    @boundscheck 1 <= 1 <= size(z.x, 1) || throw(BoundsError(z,i))
    xj = j+z.i-1
    @boundscheck xj <= lastlength(z.x) || throw(BoundsError(z,j))
    z.x[i, xj]
end

@inline @propagate_inbounds function getindex(z::IsoZNormalizer, ::typeof(!), i, inorder = i == z.bufi + 1)
    j = inorder ? i : z.n
    xj = z.i + i - 1
    @avx for k = 1:size(z.x, 1)
        z.buffer[k, j] = (z.x[k, xj] - z.μ[k]) / z.σ[k]
    end
    if inorder
        z.bufi = i
    end
    z.buffer[!, j]
end

@inline @propagate_inbounds function getindex(z::IsoZNormalizer, ::typeof(!), i::AbstractRange)
    @boundscheck (i[1] == z.i && length(i) == z.n) || throw(ArgumentError("ZNormalizers can only be indexed by ranges corresponding to their current state. Got range $i but state was $(z.i) corresponding to range $(z.i):$(z.i+z.n-1)"))
    z
end

Base.Matrix(z::IsoZNormalizer) = z.x[:,z.i:z.i+z.n-1]

Base.length(z::IsoZNormalizer) = size(z.x,1) * z.n
Base.size(z::IsoZNormalizer) = (size(z.x,1), z.n)

for T in [ZNormalizer, IsoZNormalizer]
    @eval @inline @propagate_inbounds function normalize(::Val{$T}, z::$T)
        if z.bufi == z.n
            return z.buffer
        end
        for i = z.bufi+1:z.n
            z[!, i, true] # This populates the buffer
        end
        return z.buffer
    end
end


text:
using Test

global flog=open("log.testFacto","a+")

#@testset "Logico Tests" begin
 @testset "subst Tests" begin
  @test subst([:v1,:v2,:v3],:like,[:v1,:Akira,:v3]) == :like
  @test subst([:v1,:v2,:v3],:v1,[:v1,:Akira,:v3]) == :v1
  @test subst([:v1,:v2,:v3],:v2,[:v1,:Akira,:v3]) == :Akira
  @test subst([:v1,:v2,:v3],:v2,[:v2,:v3,:v3]) == :v3

 # this should not be happen
#@test_skipsubst([:v1,:v2,:v3],:v2,[:v2,:Qbey,:v3])==:Qbey

 end

 @testset "substitution Tests" begin
  @test substitution([:v1,:v2,:v3],[:v2],[:v1,:Akira,:v3]) == [:Akira]
  @test substitution([:v1,:v2,:v3],[:v1 :v2 :v3 :v2],[:v1,:Akira,:v3]) == [:v1 :Akira :v3 :Akira]
  @test substitution([:v1,:v2,:v3],[:v1 :v2 :Taeko :v2],[:Akira,:Akira,:v2]) == [:Akira :Akira :Taeko :Akira]
  @test substitution([:v1,:v2,:v3],[:v1 :v2 :v3],[:Akira,:v2,:v2]) == [:Akira :v2 :v2]

 ## belows are should not happen
 # right side's var should be same in left
 # @test_skip substitution([:v1, :v2, :v3], [:v1 :v2 :v3], [:v2, :Akira, :v2]) == [:Akira :Akira :Akira]
 # @test_skip substitution([:v1, :v2, :v3], [:v1 :v2 :v3], [:v2, :Akira, :v3]) == [:Akira :Akira :v3]
 # @test_skip substitution([:v1, :v2, :v3], [:v1 :v2 :v3], [:Akira, :v1, :v3]) == [:Akira :Akira :Akira]
 end

 @testset "substitutionQ Tests" begin
  @test substitutionQ([:v1, :v2, :v3], [[:v1],[:v2],[:v3]], [:A, :v2, :B]) == [[:A],[:v2],[:B]]
  @test substitutionQ([:v1, :v2, :v3], [[:v1 :v1 :v2],[:v1 :v2 :v1],[:v3 :v3 :v3]], [:v1, :A, :v3]) == [[:v1 :v1 :A],[:v1 :A :v1],[:v3 :v3 :v3]]
 end

 @testset "unifytt Tests" begin
  @test unifytt([:x,:y],:a,:b) == false

  @test unifytt([:x,:y],:x,:x) == :x
  @test unifytt([:x,:y],:x,:a) == :a
  @test unifytt([:x,:y],:a,:x) == :a
  @test unifytt([:x,:y],:x,:y) == :y
  @test unifytt([:x,:y],:a,:a) == :a
 end

 @testset "unify Tests" begin
  @test unifyqf([:x,:y],[:a :y],[:b,:b]) == false

  @test unifyqf([:x,:y],[:x :y],[:x,:y]) == [:x,:y]
  @test unifyqf([:x,:y],[:x :y],[:a,:y]) == [:a,:y]
  @test unifyqf([:x,:y],[:x :y],[:a,:b]) == [:a,:b]
  @test unifyqf([:x,:y],[:x :y],[:a,:b]) == [:a,:b]

  @test unifyqf([:x,:y],[:a :y],[:a,:b]) == [:x,:b]
  @test unifyqf([:x,:y],[:x :b :y],[:a,:b,:a]) == [:a,:a]
  @test unifyqf([:x,:y],[:x :b :y],[:a,:b,:c]) == [:a,:c]

  @test  unifyqf([:x,:y],[:a :c],[:a,:b]) == false

 end

 @testset "unifymm Tests" begin
  @test unifymm([:x,:y], [:x, :y], [:x,:b]) == [:x,:b]
  @test unifymm([:x,:y], [:a, :b], [:a,:b]) == [:a,:b]
  @test unifymm([:x,:y], [:a, :b], [:x,:y]) == [:a,:b]
  @test unifymm([:x,:y], [:a, :y], [:b,:y]) == false
  @test unifymm([:x,:y], [:x, :b], [:b,:c]) == false

 end

 @testset "unifyMM Tests" begin
  @test unifyMM([:x,:y], [[:x, :y]], [[:x,:b]]) == [[:x,:b]]
  @test unifyMM([:x,:y], [[:a, :y]], [[:x,:b],[:x,:c]]) == [[:a,:b],[:a,:c]]
  @test unifyMM([:x,:y,:z,:w], [[:a,:y,:b,:w]],[[:x,:b,:b,:d],[:x,:a,:b,:a]]) == [[:a,:b,:b,:d],[:a,:a,:b,:a]]
  @test unifyMM([:x,:y,:z,:w], [[:a,:y,:b,:w],[:d,:y,:z,:w]],[[:x,:b,:b,:d]]) == [[:a,:b,:b,:d],[:d,:b,:b,:d]]
  @test unifyMM([:x,:y,:z,:w], [[:a,:y,:z,:w],[:d,:y,:b,:w]],[[:x,:b,:b,:d]]) == [[:a,:b,:b,:d],[:d,:b,:b,:d]]
  @test unifyMM([:x,:y,:z,:w], [[:a,:y,:b,:w],[:d,:y,:e,:w]],[[:x,:b,:b,:d],[:d,:y,:e,:f]]) == [[:a,:b,:b,:d],[:d,:y,:e,:f]]

 end

#==
 @testset "makemodels Tests" begin
  mm = makemodels([:x,:y],[[:P :x],[:P :b],[:Q :x :y]],[[:P :a],[:P :b],[:Q :a :b],[:Q :b :a]]) ==
# these are irregal array..
  @test mm[1,:]= [[:a,:y] false false]
  @test mm[2,:]= [[:b,:y][:x,:y]false]
  @test mm[3,:]= [false false [:a,:b]]
  @test mm[4,:]= [false false [:b,:a]]]
 end
==#

 @testset "reduction Tests" begin
  mm=makemodels([:x,:y],[[:P :x :x]],[[:P :a :b], [:P :a :b], [:P :b :b], [:P :c :d]])
  @test reduction([:x,:y],[:x,:y],mm,1) == [[:b,:y]]

  mm=makemodels([:x,:y],[[:P :x],[:P :b],[:Q :x :y]],[[:P :a],[:P :b],[:Q :a :b],[:Q :b :a]])
  @test reduction([:x,:y],[:x :y],mm,1) == [[:a,:b],[:b,:a]]
 end

 @testset "solveq Tests" begin
  @test solveq([:x, :y],[:have :x :orange], [[:have :Hanako :ringko]]) == false
  @test solveq([:x, :y],[:have :x :ringo], [[:have :Hanako :ringo]])==[[:Hanako, :y]]
  @test solveq([:x, :y],[:have :x :y], [[:have :Hanako :ringo],[:have :Kiro :jiro]])==[[:Hanako, :ringo],[:Kiro, :jiro]]
  @test solveq([:x, :y],[:love :x :Hanako],[[:love :Taro :Hanako],[:love :Kenji :Michiko],[:love :Jiro :Hanako],[:kill :Moji :Do],[:love :Hanako :Hanako]]) == [[:Taro, :y],[:Jiro,:y],[:Hanako,:y]]
 end

 @testset "solveQ Tests" begin
  @test solveQ([:x,:y],[[:P :x :x]],[[:P :a :b],[:P :a :b],[:P :b :b],[:P :c :d]]) == [[:b,:y]]
  @test solveQ([:x,:y],[[:P :x :y]],[[:P :a :b],[:P :c :d]]) == [[:a,:b],[:c,:d]]
  @test solveQ([:x,:z],[[:P :x],[:Q :x :z],[:R :z :u]],[[:P :r],[:Q :r :o],[:R :o :u]]) == [[:r,:o]]

## ↓[R :n :*]はどこにもないが、否定できないので[[:r,:n]]が答えになる。いいのか? 否定の扱いの問題
## K's elements are all true. isn't is?
  @test solveQ([:x,:z],[[:P :x],[:Q :x :z],[:R :z :u]],[[:P :r],[:Q :r :n],[:R :o :u]]) == []
  @test solveQ([:x,:z],[[:P :x],[:Q :x :z],[:R :z :u]],[[:P :r],[:Q :r :n],[:R :n :u]]) == [[:r,:n]]

  @test solveQ([:x,:z,:w],[[:P :x],[:Q :x :z],[:R :z :w]],[[:P :r],[:P :s],[:Q :r :o],[:Q :s :u],[:R :u :v],[:R :o :u]]) == [[:r,:o,:u],[:s,:u,:v]]
  @test solveQ([:x,:y,:z], [[:P :a :x], [:Q :x :b :z], [:R :z :u]], [[:P :a :r], [:Q :r :b :o], [:R :o :u]]) == [[:r,:y,:o]]
  @test solveQ([:x,:y,:z],[[:P :a :x],[:Q :x :b :z],[:R :z :w]],[[:P :a :r],[:P :a :s],[:Q :r :b :o],[:Q :r :c :w],[:Q :s :b :h],[:R :o :u],[:R :w :v]] ) == []
  @test solveQ([:x,:y,:z],[[:P :a :x],[:Q :x :b :z],[:R :z :w]],[[:P :a :r],[:P :a :s],[:Q :r :b :o],[:Q :r :c :w],[:Q :s :b :h],[:R :o :w],[:R :w :v]] ) == [[:r,:y,:o]]

# case no contradiction but no model. what is this?
  @test solveQ([:x,:y,:z],[[:P :a :x],[:Q :x :b :z],[:R :z :w]],[[:P :a :r],[:P :a :s],[:Q :r :b :o],[:Q :r :c :w],[:Q :s :b :h],[:R :o :u],[:R :w :v]] ) == []
  @test solveQ([:x,:y],[[:P :x],[:R :y :v],[:Q :x :y]],[[:P :r],[:R :o :u],[:R :w :v],[:Q :r :o],[:Q :r :w],[:Q :s :h]] ) == [[:r,:w]]
  @test solveQ([:x,:y],[[:P :x],[:R :y :v],[:Q :x :y]],[[:P :r],[:P :s],[:R :o :u],[:R :w :v],[:Q :r :o],[:Q :r :w],[:Q :s :w]] ) == [[:r,:w],[:s,:w]]
  @test solveQ([:x,:y],[[:P :x],[:Q :x :y],[:R :y :v]],[[:P :r],[:P :s],[:R :o :u],[:R :w :v],[:Q :r :o],[:Q :r :w],[:Q :s :h]] ) == [[:r,:w]]
  @test solveQ([:x,:y],[[:P :x],[:Q :x :y],[:R :y :u]],[[:P :r],[:P :s],[:Q :r :o],[:Q :r :w],[:Q :s :h],[:R :o :u],[:R :w :v]] ) == [[:r,:o]]

 end

 kb=[[:H :know :K],[:H :know :T],[:K :has :ringo],[:T :wear :jacket],[:ringo :is :aka],[:jacket :is :aka]]
 @testset "Final answer Tests" begin
  @test solveQ([:who,:what],[[:H :love :who],[:who :has :what],[:what :is :aka]], [[:H :love :M],[:M :has :ringo],[:ringo :is :aka]]) ==[[:M, :ringo]]

  @test solveQ([:who,:act,:what],[[:H :act :who],[:who :has :what],[:what :is :aka]], [[:H :love :M],[:M :has :ringo],[:ringo :is :aka],[:H :kill :Majo],[:Majo :has :cake],[:cake :is :sweet]]) ==[[:M,:love,:ringo]]

 end

#end

@testset "END OF testFacto" begin end

close(flog)

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
using Dates
using Base.Test

include("types.jl")
include("periods.jl")
include("accessors.jl")
include("query.jl")
include("arithmetic.jl")
include("conversions.jl")
include("ranges.jl")
include("adjusters.jl")
include("io.jl")

#TODO
 #Timezones.jl

#NEED TESTS

#IDEAS
 #research JSR-310, PHP? javascript? go? C#? for API completeness
 #round(dt,period)
 #add(dt,y,m,d,h,mi,s,ms); sub(dt,y,m,d,h,mi,s,ms); many ariths at once?
 #conversions: ruby? python? javascript? etc.
 #make TimeStamp fully parameterized (Instant, Calendar)
  #have datetime field + nanosecond + timezone field?
 #d"2014-01-01"
 #dt"2014-01-01"


text:
include("../ranger.jl")

using .Ranger

const RGeo = Ranger.Geometry
const RAni = Ranger.Animation
const RNodes = Ranger.Nodes

function test_basic_create_tween()
    println("Testing basic create tween")

    node = Nodes.NodeNil()

    tween = Animation.create_tween_to(node, 1, 10.0)

    println("Passed")
end

function test_tween_pool()
    println("Testing tween pool")

    node = Nodes.NodeNil()

    tween = Animation.create_tween_to(node, 1, 10.0)
    @assert Animation.pool_size() == 0 "Expected 0 item(s) in pool"

    Animation.free!(tween)
    @assert Animation.pool_size() == 1 "Expected 1 item(s) in pool"

    tween2 = Animation.create_tween_to(node, 1, 10.0)

    Animation.free!(tween2)
    @assert Animation.pool_size() == 1 "Expected 1 item(s) in pool"

    tween3 = Animation.create_tween_to(node, 1, 10.0)
    @assert Animation.pool_size() == 0 "Expected 0 item(s) in pool"

    tween4 = Animation.create_tween_to(node, 1, 10.0)
    @assert Animation.pool_size() == 0 "Expected 0 item(s) in pool"

    Animation.free!(tween3)
    Animation.free!(tween4)
    @assert Animation.pool_size() == 2 "Expected 2 item(s) in pool"

    tween5 = Animation.create_tween_to(node, 1, 10.0)
    @assert Animation.pool_size() == 1 "Expected 1 item(s) in pool"

    tween6 = Animation.create_tween_to(node, 1, 10.0)
    @assert Animation.pool_size() == 0 "Expected 0 item(s) in pool"

    println("Passed")
end

function test_tween_manager()
    println("Testing tween manager")

    man = Animation.TweenManager(false)

    node = Nodes.NodeNil()

    seq_timeline_tween = Animation.TimeLine(Animation.TIMELINE_SEQUENCE)

    Animation.add!(man, seq_timeline_tween)

    @assert seq_timeline_tween.base.started == true "Expected a started timeline"

    Animation.update!(man, 1.0)

    println("Passed")
end

# A test fixture to testing the test_tweens
mutable struct TestAccessor <: RAni.AbstractTweenAccessor
    FIELD_POSITION::Int64

    function TestAccessor()
        o = new()
        o.FIELD_POSITION = 1
        o
    end
end

function get_values!(accessor::TestAccessor, tween::RAni.Tween, return_buffer::Array{Float64})
    if tween.value_type == accessor.FIELD_POSITION
        println("getting FIELD_POSITION")
        return_buffer[1] = tween.target.transform.position.x
        return_buffer[2] = tween.target.transform.position.y
        return 2 # specify how many values are to be changed.
    end
end

function set_values!(accessor::TestAccessor, tween::RAni.Tween, values::Array{Float64})
    if tween.value_type == accessor.FIELD_POSITION
        println("setting FIELD_POSITION")
        tween.target.transform.position.x = values[1]
        tween.target.transform.position.y = values[2]
        # Mark node dirty. etc...
    end
end

function test_tweens()
    println("Testing tweens")

    accessor = TestAccessor()

    RAni.register_accessor(RAni.CLASS_NODE, accessor)

    node = Nodes.TrivialNode()
    man = Animation.TweenManager(false)

    tween = RAni.create_tween_to(node, accessor.FIELD_POSITION, 10.0)
    RAni.set_target_values!(tween, [20.0, 0.0])
    RAni.start!(man, tween)

    println("Passed")
end

test_tweens()

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using PackageCompiler
create_sysimage([:Pluto, :PlutoUI, :Plots, :LaTeXStrings, :Unitful, :DataFrames, :LinearAlgebra];
                precompile_execution_file = "warmup.jl",
                replace_default = true,
                cpu_target = PackageCompiler.default_app_cpu_target())


text:
module AutoInteract

using Requires
using DataStructures
using AxisArrays
using Plots

@require InteractNext begin
    using InteractNext
    using Observables
    using WebIO
    using CSSUtil
end
@require Interact begin
    using Interact
    using Reactive
end
@require DataFrames begin
    using DataFrames
end
@require StaticArrays begin
    using StaticArrays
end

include("interface.jl")
include("signals.jl")
include("api.jl")
include("plots.jl")

export set_defaults, steprange, box, make_widgets, make_interface, make_plottables, spreadwidgets
export plot_all, plot_all!
export get_signals, apply_signals!

end # module


text:
module JLD2
using DataStructures, Reexport
import Base.sizeof
using MacroTools
using Printf
using Mmap
using TranscodingStreams
@reexport using FileIO: load, save

export jldopen, @load, @save, save_object, load_object, printtoc
export jldsave

const OBJECT_HEADER_SIGNATURE = htol(0x5244484f) # "OHDR"

# Currently we specify that all offsets and lengths are 8 bytes
const Length = UInt64

struct UnsupportedVersionException <: Exception
    msg::String
end
struct UnsupportedFeatureException <: Exception
    msg::String
end
struct InvalidDataException <: Exception
    msg::String
end
struct InternalError <: Exception
    msg::String
end

# In the future a more descriptive error should be returned
UnsupportedVersionException() = UnsupportedVersionException("")
UnsupportedFeatureException() = UnsupportedFeatureException("")
InvalidDataException() = InvalidDataException("")
InternalError() = InternalError("")

# Due to custom overrides we do not use Base functions directly
# but define our own to avoid type piracy
jlwrite(io, x) = Base.write(io, x)
jlread(io, x) = Base.read(io, x)
jlsizeof(x) = Base.sizeof(x)
jlunsafe_store!(p, x) = Base.unsafe_store!(p, x)
jlunsafe_load(p) = Base.unsafe_load(p)

include("file_header.jl")
include("Lookup3.jl")
include("mmapio.jl")
include("bufferedio.jl")
include("misc.jl")

"""
    RelOffset

Represents an HDF5 relative offset. This differs from a file offset (used elsewhere) in
that it is relative to the superblock base address. In practice, this means that
`FILE_HEADER_LENGTH `has been subtracted. `fileoffset` and `h5offset` convert between
`RelOffsets` and file offsets.
"""
struct RelOffset
    offset::UInt64
end
define_packed(RelOffset)
Base.:(==)(x::RelOffset, y::RelOffset) = x === y
Base.hash(x::RelOffset) = hash(x.offset)

const UNDEFINED_ADDRESS = RelOffset(0xffffffffffffffff)
const NULL_REFERENCE = RelOffset(0)

"""
    JLDWriteSession{T}

A JLDWriteSession keeps track of references to serialized objects. If `T` is a Dict,
`h5offset` maps an object ID (returned by calling `objectid`) to th `RelOffset` of the
written dataset. If it is `Union{}`, then references are not tracked, and objects
referenced multiple times are written multiple times.
"""
struct JLDWriteSession{T<:Union{Dict{UInt,RelOffset},Union{}}}
    h5offset::T
    objects::Vector{Any}

    JLDWriteSession{T}() where T = new()
    JLDWriteSession{T}(h5offset, objects) where T = new(h5offset, objects)
end
JLDWriteSession() = JLDWriteSession{Dict{UInt,RelOffset}}(Dict{UInt,RelOffset}(), Any[])

"""
    GlobalHeap

Represents an HDF5 global heap structure.
"""
mutable struct GlobalHeap
    offset::Int64
    length::Length
    free::Length
    objects::Vector{Int64}
end

"""
    H5Datatype

Supertype of all HDF5 datatypes.
"""
abstract type H5Datatype end

struct CommittedDatatype <: H5Datatype
    header_offset::RelOffset
    index::Int
end

"""
    ReadRepresentation{T,ODR}

A type encoding both the Julia type `T` and the on-disk (HDF5) representation `ODR`.
"""
struct ReadRepresentation{T,ODR} end

"""
    CustomSerialization{T,S}

On-disk representation for data that is written as if it were of Julia type `T`, but is
read as type `S`.
"""
struct CustomSerialization{T,S} end

"""
    Group{T}

JLD group object.
"""
mutable struct Group{T}
    f::T
    last_chunk_start_offset::Int64
    continuation_message_goes_here::Int64
    last_chunk_checksum_offset::Int64
    next_link_offset::Int64
    est_num_entries::Int
    est_link_name_len::Int
    unwritten_links::OrderedDict{String,RelOffset}
    unwritten_child_groups::OrderedDict{String,Group}
    written_links::OrderedDict{String,RelOffset}

    Group{T}(f; est_num_entries::Int=4, est_link_name_len::Int=8) where T =
        new(f, -1, -1, -1, -1, est_num_entries, est_link_name_len,
        OrderedDict{String,RelOffset}(), OrderedDict{String,Group}())

    Group{T}(f, last_chunk_start_offset, continuation_message_goes_here,
             last_chunk_checksum_offset, next_link_offset,
             est_num_entries, est_link_name_len,
             unwritten_links, unwritten_child_groups,
             written_links) where T =
        new(f, last_chunk_start_offset, continuation_message_goes_here,
            last_chunk_checksum_offset, next_link_offset,
            est_num_entries, est_link_name_len,
            unwritten_links, unwritten_child_groups,
            written_links)
end

"""
    JLDFile{T<:IO}

JLD file object.
"""
mutable struct JLDFile{T<:IO}
    io::T
    path::String
    writable::Bool
    written::Bool
    compress#::Union{Bool,Symbol}
    mmaparrays::Bool
    n_times_opened::Int
    datatype_locations::OrderedDict{RelOffset,CommittedDatatype}
    datatypes::Vector{H5Datatype}
    datatype_wsession::JLDWriteSession{Dict{UInt,RelOffset}}
    jlh5type::IdDict
    h5jltype::IdDict
    jloffset::Dict{RelOffset,WeakRef}
    end_of_data::Int64
    global_heaps::Dict{RelOffset,GlobalHeap}
    global_heap::GlobalHeap
    loaded_groups::Dict{RelOffset,Group}
    root_group_offset::RelOffset
    root_group::Group
    types_group::Group

    function JLDFile{T}(io::IO, path::AbstractString, writable::Bool, written::Bool,
                        compress,#::Union{Bool,Symbol},
                        mmaparrays::Bool) where T
        f = new(io, path, writable, written, compress, mmaparrays, 1,
            OrderedDict{RelOffset,CommittedDatatype}(), H5Datatype[],
            JLDWriteSession(), IdDict(), IdDict(), Dict{RelOffset,WeakRef}(),
            Int64(FILE_HEADER_LENGTH + jlsizeof(Superblock)), Dict{RelOffset,GlobalHeap}(),
            GlobalHeap(0, 0, 0, Int64[]), Dict{RelOffset,Group}(), UNDEFINED_ADDRESS)
        finalizer(jld_finalizer, f)
        f
    end
end
JLDFile(io::IO, path::AbstractString, writable::Bool, written::Bool, compress, mmaparrays::Bool) =
    JLDFile{typeof(io)}(io, path, writable, written, compress, mmaparrays)

"""
    fileoffset(f::JLDFile, x::RelOffset)

Converts an offset `x` relative to the superblock of file `f` to an absolute offset.
"""
fileoffset(f::JLDFile, x::RelOffset) = Int64(x.offset + FILE_HEADER_LENGTH)

"""
    h5offset(f::JLDFile, x::RelOffset)

Converts an absolute file offset `x` to an offset relative to the superblock of file `f`.
"""
h5offset(f::JLDFile, x::Int64) = RelOffset(x - FILE_HEADER_LENGTH)

#
# File
#

openfile(::Type{IOStream}, fname, wr, create, truncate, fallback::Nothing = nothing) =
    open(fname, read = true, write = wr, create = create,
         truncate = truncate, append = false)
openfile(::Type{MmapIO}, fname, wr, create, truncate, fallback::Nothing = nothing) =
    MmapIO(fname, wr, create, truncate)

function openfile(T::Type, fname, wr, create, truncate, fallback::Type)
    try
         openfile(T, fname, wr, create, truncate, nothing)
    catch
        @warn "Opening file with $T failed, falling back to $fallback"
        openfile(fallback, fname, wr, create, truncate, nothing)
    end
end

# default fallback behaviour : MmapIO -> IOStream -> failure
FallbackType(::Type{MmapIO}) = IOStream
FallbackType(::Type{IOStream}) = nothing

# The delimiter is excluded by default
read_bytestring(io::IOStream) = String(readuntil(io, 0x00))

const OPEN_FILES = Dict{String,WeakRef}()
const OPEN_FILES_LOCK = ReentrantLock()
function jldopen(fname::AbstractString, wr::Bool, create::Bool, truncate::Bool, iotype::T=MmapIO;
                 fallback::Union{Type, Nothing} = FallbackType(iotype),
                 compress=false,
                 mmaparrays::Bool=false
                 ) where T<:Union{Type{IOStream},Type{MmapIO}}
    mmaparrays && @warn "mmaparrays keyword is currently ignored" maxlog=1
    verify_compressor(compress)
    exists = ispath(fname)

    lock(OPEN_FILES_LOCK)
    f = try
        if exists
            rname = realpath(fname)
            # catch existing file system entities that are not regular files
            if !isfile(rname)
                throw(ArgumentError("not a regular file: $fname"))
            end
            if haskey(OPEN_FILES, rname)
                ref = OPEN_FILES[rname]
                f = ref.value
                if ref.value !== nothing
                    if truncate
                        throw(ArgumentError("attempted to truncate a file that was already open"))
                    elseif !isa(f, JLDFile{iotype})
                        throw(ArgumentError("attempted to open file with $iotype backend, but already open with a different backend"))
                    elseif f.writable != wr
                        current = wr ? "read/write" : "read-only"
                        previous = f.writable ? "read/write" : "read-only"
                        throw(ArgumentError("attempted to open file $(current), but file was already open $(previous)"))
                    elseif f.compress != compress
                        throw(ArgumentError("attempted to open file with compress=$(compress), but file was already open with compress=$(f.compress)"))
                    elseif f.mmaparrays != mmaparrays
                        throw(ArgumentError("attempted to open file with mmaparrays=$(mmaparrays), but file was already open with mmaparrays=$(f.mmaparrays)"))
                    end

                    f = f::JLDFile{iotype}
                    f.n_times_opened += 1
                    return f
                end
            end
        end

        io = openfile(iotype, fname, wr, create, truncate, fallback)
        created = !exists || truncate
        rname = realpath(fname)
        f = JLDFile(io, rname, wr, created, compress, mmaparrays)
        OPEN_FILES[rname] = WeakRef(f)
        f
    catch e
        rethrow(e)
    finally
        unlock(OPEN_FILES_LOCK)
    end
    if f.written
        f.root_group = Group{typeof(f)}(f)
        f.types_group = Group{typeof(f)}(f)
    else
        verify_file_header(f)

        seek(f.io, FILE_HEADER_LENGTH)
        superblock = jlread(f.io, Superblock)
        f.end_of_data = superblock.end_of_file_address
        f.root_group_offset = superblock.root_group_object_header_address
        f.root_group = load_group(f, superblock.root_group_object_header_address)

        if haskey(f.root_group.written_links, "_types")
            types_group_offset = f.root_group.written_links["_types"]
            f.types_group = f.loaded_groups[types_group_offset] = load_group(f, types_group_offset)
            i = 0
            for offset in values(f.types_group.written_links)
                f.datatype_locations[offset] = CommittedDatatype(offset, i += 1)
            end
            resize!(f.datatypes, length(f.datatype_locations))
        else
            f.types_group = Group{typeof(f)}(f)
        end
    end

    f
end

"""
    jldopen(fname::AbstractString, mode::AbstractString; iotype=MmapIO, compress=false)

Opens a JLD2 file at path `fname`.

`"r"`: Open for reading only, failing if no file exists
`"r+"`: Open for reading and writing, failing if no file exists
`"w"`/`"w+"`: Open for reading and writing, overwriting the file if it already exists
`"a"`/`"a+"`: Open for reading and writing, creating a new file if none exists, but
              preserving the existing file if one is present
"""
function jldopen(fname::AbstractString, mode::AbstractString="r"; iotype=MmapIO, kwargs...)
    (wr, create, truncate) = mode == "r"  ? (false, false, false) :
                             mode == "r+" ? (true, false, false) :
                             mode == "a" || mode == "a+" ? (true, true, false) :
                             mode == "w" || mode == "w+" ? (true, true, true) :
                             throw(ArgumentError("invalid open mode: $mode"))
    jldopen(fname, wr, create, truncate, iotype; kwargs...)
end

"""
    load_datatypes(f::JLDFile)

Populate f.datatypes and f.jlh5types with all of the committed datatypes from a file. We
need to do this before writing to make sure we reuse written datatypes.
"""
function load_datatypes(f::JLDFile)
    dts = f.datatypes
    cdts = f.datatype_locations
    @assert length(dts) == length(cdts)
    i = 1
    for cdt in values(cdts)
        !isassigned(dts, i) && jltype(f, cdt)
        i += 1
    end
end

"""
    prewrite(f::JLDFile)

Check that a JLD file is actually writable, and throw an error if not. Sets the `written`
flag on the file.
"""
function prewrite(f::JLDFile)
    f.n_times_opened == 0 && throw(ArgumentError("file is closed"))
    !f.writable && throw(ArgumentError("file was opened read-only"))
    !f.written && load_datatypes(f)
    f.written = true
end

Base.read(f::JLDFile, name::AbstractString) = f.root_group[name]
#Base.write(f::JLDFile, name::AbstractString, obj, wsession::JLDWriteSession=JLDWriteSession()) =
#    write(f.root_group, name, obj, wsession)

Base.getindex(f::JLDFile, name::AbstractString) = f.root_group[name]
Base.setindex!(f::JLDFile, obj, name::AbstractString) = (f.root_group[name] = obj; f)
Base.haskey(f::JLDFile, name::AbstractString) = haskey(f.root_group, name)
Base.isempty(f::JLDFile) = isempty(f.root_group)
Base.keys(f::JLDFile) = filter!(x->x != "_types", keys(f.root_group))
Base.keytype(f::JLD2.JLDFile) = String
Base.get(default::Function, f::Union{JLDFile, Group}, name::AbstractString) =
    haskey(f, name) ? f[name] : default()
Base.get(f::Union{JLDFile, Group}, name::AbstractString, default) =
    haskey(f, name) ? f[name] : default
Base.get!(f::Union{JLDFile, Group}, name::AbstractString, default) =
    get!(() -> default, f, name)
function Base.get!(default::Function, f::Union{JLDFile, Group}, name::AbstractString)
    if haskey(f, name)
        return f[name]
    else
        default_value = default()
        f[name] = default_value
        return default_value
    end
end

function Base.close(f::JLDFile)
    if f.n_times_opened != 1
        f.n_times_opened == 0 && return
        f.n_times_opened -= 1
        return
    end

    io = f.io
    if f.written
        # Save any groups we know of that have been modified
        for group in values(f.loaded_groups)
            save_group(group)
        end
        if !isempty(f.types_group) && !haskey(f.root_group, "_types")
            f.root_group["_types"] = f.types_group
        end
        res = save_group(f.root_group)
        if f.root_group_offset == UNDEFINED_ADDRESS
            f.root_group_offset = res
        end

        # Write JLD2 header
        seek(io, 0)
        jlwrite(io, FILE_HEADER)

        # Write superblock
        seek(io, FILE_HEADER_LENGTH)
        jlwrite(io, Superblock(0, FILE_HEADER_LENGTH, UNDEFINED_ADDRESS,
              f.end_of_data, f.root_group_offset))

        truncate_and_close(io, f.end_of_data)
    else
        close(io)
    end
    lock(OPEN_FILES_LOCK) do
        delete!(OPEN_FILES, f.path)
        f.n_times_opened = 0
    end
    nothing
end

"""
    jld_finalizer(f::JLDFile)

When a JLDFile is finalized, it is possible that the `MmapIO` has been munmapped, since
Julia does not guarantee finalizer order. This means that the underlying file may be closed
before we get a chance to write to it.
"""
function jld_finalizer(f::JLDFile{MmapIO})
    f.n_times_opened == 0 && return
    if f.written && !isopen(f.io.f)
        f.io.f = open(f.path, "r+")
    end
    f.n_times_opened = 1
    close(f)
end

function jld_finalizer(f::JLDFile{IOStream})
    f.n_times_opened == 0 && return
    if f.written && !isopen(f.io)
        f.io = openfile(IOStream, f.path, true, false, false)
    end
    f.n_times_opened = 1
    close(f)
end

# Display functions

# simple one-line display (without trailing line break)
function Base.show(io::IO, f::JLDFile)
    print(io, "JLDFile $(f.path)")
    if get(io, :compact, false)
        return
    else
        print(io, f.writable ? " (read/write)" : " (read-only)")
    end
end

# fancy multi-line display (unless an IOContext requests compactness)
function Base.show(io::IO, ::MIME"text/plain", f::JLDFile)
    show(io, f)
    if !get(io, :compact, false)
        print(io, "\n")
        printtoc(io, f; numlines = get(io, :jld2_numlines, 10))
    end
end

"""
    printtoc([io::IO,] f::JLDFile [; numlines])
Prints an overview of the contents of `f` to the `IO`.

Use the optional `numlines` parameter to restrict the amount of items listed.
"""
printtoc(f::JLDFile; kwargs...) = printtoc(Base.stdout, f; kwargs...)
printtoc(io::IO, f::JLDFile; numlines = typemax(Int64)) =
    show_group(io, f.root_group, numlines, " ", true)

include("superblock.jl")
include("object_headers.jl")
include("groups.jl")
include("dataspaces.jl")
include("attributes.jl")
include("datatypes.jl")
include("datasets.jl")
include("global_heaps.jl")

include("data/type_defs.jl")
include("data/specialcased_types.jl")
include("data/number_types.jl")
include("data/custom_serialization.jl")
include("data/writing_datatypes.jl")
include("data/reconstructing_datatypes.jl")

include("dataio.jl")
include("loadsave.jl")
include("stdlib.jl")
include("backwards_compatibility.jl")
include("inlineunion.jl")
include("fileio.jl")
include("compression.jl")

if Base.VERSION >= v"1.6.0"
    include("precompile.jl")
    _precompile_()
end

end # module


text:
# This file holds the definition of the functions pertaining to the History
#   type.

# The functions of the History type require the HistoryEntry type.
requiredTypes = [ "historyEntry", "history" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

# This function adds an entry to the history.
export addEntry!
function addEntry!( hist::History, timestamp::T, newState ) where T <: Real
    # Check if there's no entry with this timestamp already.
    ii = searchsortedlast( hist.history, timestamp )

    if ( ii > 0 ) && ( hist.history[ ii ] == timestamp )
        hist.history[ ii ].newVal = newState
        return
    end  # if ( ii <= length( hist.history ) ) ...

    # Add the new entry at the end.
    ii = length( hist.history )
    push!( hist.history, HistoryEntry( timestamp, newState ) )

    # Place the entry where it belongs.
    while ( ii >= 1 ) && ( hist.history[ ii ] > hist.history[ ii + 1 ] )
        hist.history[ ii ], hist.history[ ii + 1 ] =
            hist.history[ ii + 1 ], hist.history[ ii ]
        ii -= 1
    end  # while ( ii >= 1 ) ...
end  # addEntry!( hist, timestamp, newState )

# This function returns the value of the attribute at the requested timestamp.
# If the timestamp is earlier than the first time in the history, or if the
#   history is empty, nothing gets returned.
function Base.getindex( hist::History, timestamp::T ) where T <: Real
    ii = searchsortedlast( hist.history, timestamp )

    return ii == 0 ? nothing : hist.history[ ii ].newVal
end  # getindex( hist, timestamp )

function Base.setindex!( hist::History, data, timestamp::T ) where T <: Real
    addEntry!( hist, timestamp, data )
end

function Base.show( io::IO, hist::History )
    print( io, "Attribute: $(hist.attribute)" )

    for entry in hist.history
        print( io, "\n$entry" )
    end  # for entry in hist.history
end  # show( io, hist )


text:
# using StaticArrays

# Emitter Types ================================================================
abstract type Emitter end

struct Ray{T <: Real} <: Emitter
    src :: SVector{3,T}
    trj :: SVector{3,T}
end

Ray(src::Vector{Float64}, trj::Vector{Float64}) = Ray(
    SVector{3,Float64}(src),
    SVector{3,Float64}(trj),
)
Ray() = Ray(
    SVector{3,Float64}([0.0, 0.0, 0.0]),
    SVector{3,Float64}([0.0, 0.0, 1.0]),
)

struct PointSource{T <: Real} <: Emitter
    src :: SVector{3,T}
end
PointSource(src::Vector{Float64}) = PointSource(SVector{3,Float64}(src))
PointSource() = PointSource(SVector{3,Float64}([0.0, 0.0, 0.0]))

# Functions ===================================================================
function random_traj()
    cosθ = 2.0*rand() - 1.0
    sinθ = sqrt(1.0 - cosθ^2) # sinθ is always positive
    sinψ, cosψ = sincos(2π * rand())
    return [sinθ * cosψ, sinθ * sinψ, cosθ]
end

makephoton(s::Ray) = Photon(s.src, s.trj)
makephoton(s::PointSource) = Photon(s.src,  SVector{3,Float64}(random_traj()) )


text:
###
### LongSequence
###
###
### A general purpose biological sequence representation.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# About Internals
# ---------------
#
# The `data` field of a `LongSequence{A}` object contains binary representation
# of a biological character sequence. Each character is encoded with an encoder
# corresponding to the alphabet `A` and compactly packed into `data`. To extract
# a character from a sequence, you should decode this binary sequence with a
# decoder that is a pair of the encoder. The length of encoded binary bits is
# fixed, and hence a character at arbitrary position can be extracted in a
# constant time. To know the exact location of a character at a position, you
# can use the `bitindex(seq, i)` function, which returns a pair of element's index
# containing binary bits and bits' offset. As a whole, character extraction
# `seq[i]` can be written as:
#
#     j = bitindex(seq, i)
#     decode(A, (seq.data[index(j)] >> offset(j)) & mask(A))
#
#  index :        index(j) - 1       index(j)       index(j) + 1
#   data : ....|xxxxx...........|xxXxxxxxxxxxxxxx|............xxxx|....
# offset :                          |<-offset(j)-|
#  width :      |<---- 64 ---->| |<---- 64 ---->| |<---- 64 ---->|
#
#  * '.' : unused (4 bits/char)
#  * 'x' : used
#  * 'X' : used and pointed by index `i`

"""
Biological sequence data structure indexed by an alphabet type `A`.
"""
mutable struct LongSequence{A <: Alphabet} <: BioSequence{A}
    data::Vector{UInt64}  # encoded character sequence data
    len::Int

    function LongSequence{A}(data::Vector{UInt64}, len::Int) where {A <: Alphabet}
        return new(data, len)
    end
end

const LongNucleotideSequence = LongSequence{<:NucleicAcidAlphabet}
const LongDNASeq       = LongSequence{DNAAlphabet{4}}
const LongRNASeq       = LongSequence{RNAAlphabet{4}}
const LongAminoAcidSeq = LongSequence{AminoAcidAlphabet}
const LongCharSeq      = LongSequence{CharAlphabet}

###
### Required type traits and methods
###

"Gets the alphabet encoding of a given BioSequence."
BioSymbols.alphabet(::Type{LongSequence{A}}) where {A} = alphabet(A)
Alphabet(::Type{LongSequence{A}}) where {A <: Alphabet} = A()
Base.length(seq::LongSequence) = seq.len
bindata(seq::LongSequence) = seq.data
Base.eltype(::Type{LongSequence{A}}) where {A} = eltype(A)

@inline seq_data_len(s::LongSequence{A}) where A = seq_data_len(A, length(s))
@inline function seq_data_len(::Type{A}, len::Integer) where A <: Alphabet
	iszero(bits_per_symbol(A())) && return 0
    return cld(len, div(64, bits_per_symbol(A())))
end

@inline function encoded_data(seq::LongSequence)
    return seq.data
end

include("seqview.jl")
include("indexing.jl")
include("constructors.jl")
include("printing.jl")
include("copying.jl")
include("conversion.jl")
include("stringliterals.jl")
include("transformations.jl")
include("operators.jl")
include("counting.jl")


text:
using Colors, ColorSchemes
using CairoMakie, FileIO
using CSV, DataFrames
using Statistics

imgpath = joinpath("Molenaar", "Figure 4", "Plots")
resultspath = joinpath("Molenaar", "Results")

x = [res["sub_out"] for res in stats]
a = [res["cat_flux"] for res in stats]
b = [res["met_flux"] for res in stats]

fig = Figure()
ax = Axis(fig[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)

ax.xlabel = "Substrate concentration"
ax.ylabel = "Fraction of the total substarte flux"

fig[1, 1] = Legend(
    fig,
    [[scatter1, line1], [scatter2, line2]],
    ["MetEf_flux", "CatEf_flux"],
    tellheight = false,
    tellwidth = false,
    halign = :right,
    valign = :center,
    labelsize = 14,
)
fig

FileIO.save(joinpath(imgpath, "Figure_4A.pdf"), fig)


text:
struct DiffRuleOperation
    diffrule::DiffRule
    operations::OffsetVector{Operation,Vector{Operation}}
end
function DiffRuleOperation(dr::DiffRule, nparents::Int = num_parents(dr))
    ops = OffsetArray{Operation}(undef, -nparents:last(last(dr.sections)))
    DiffRuleOperation(
        dr, ops
    )
end

function DiffRuleOperation(op::Operation, newops::Vector{Operation})
    vparents = parents(op)
    nargs = length(vparents)
    dr = DERIVATIVERULES[InstructionArgs(instruction(op), nargs)]
    ops = OffsetArray{Operation}(undef, -nargs:last(last(dr.sections)))
    for i ∈ -nargs:-1
        ops[i] = newops[identifier(vparents[i + nargs + 1])]
    end
    DiffRuleOperation( dr, ops )
end
function LoopVectorization.parents(dro::DiffRuleOperation)
    nparents = num_parents(dro.diffrule)
    view(dro.operations[-nparents:-1])
end
LoopVectorization.instruction(dro::DiffRuleOperation, i::Int) = instruction(dro.diffrule, i)
LoopVectorization.operations(dro::DiffRuleOperation) = dro.operations
function returned_ind(dro::DiffRuleOperation, i::Int)
    i == 2 && return typemin(Int)
    dro.diffrule.returns[i - (i > 2)]
end
returned_inds(dro::DiffRuleOperation) = dro.diffrule.returns
section(dro::DiffRuleOperation, i::Int) = dro.diffrule.sections[i]
sections(dro::DiffRuleOperation) = dro.diffrule.sections
dependencies(dro::DiffRuleOperation) = dro.diffrule.dependencies
dependencies(dro::DiffRuleOperation, i::Int) = dro.diffrule.dependencies[i]

function Base.getindex(dro::DiffRuleOperation, i::Int)
    # i -= first(dro.operations.offsets)
    instruction(dro, i), dependencies(dro, i)
end


text:
@testset "pure state" begin
    dim = DIM

    @test CoherentState(α(2., π/4), dim=dim) == displace!(VacuumState(dim=dim), α(2., π/4))
    @test CoherentState(α(2., π/4), dim=dim, rep=StateMatrix) ==
        displace!(VacuumState(dim=dim, rep=StateMatrix), α(2., π/4))
    @test SqueezedState(ξ(2., π/4), dim=dim) == squeeze!(VacuumState(dim=dim), ξ(2., π/4))
    @test SqueezedState(ξ(2., π/4), dim=dim, rep=StateMatrix) ==
        squeeze!(VacuumState(dim=dim, rep=StateMatrix), ξ(2., π/4))
end

@testset "mixed state" begin
    dim = DIM

    n̄ = 0.5
    n = 5
    @test QSB.bose_einstein(n̄)(n) == n̄^n / (1 + n̄)^(n+1)

    @test ThermalState(n̄, dim=dim) ==
        StateMatrix(diagm(QSB.bose_einstein(n̄).(ComplexF64.(0:dim-1))), dim)
    @test SqueezedThermalState(ξ(1., π/4), n̄, dim=dim) ==
        squeeze!(ThermalState(ComplexF64, n̄, dim=dim), ξ(1., π/4))
end


text:
"""
Array type of rank 3 (taking three scalar indices) that is used for the storage
of discrete convolution operators that happen to be the tensor product of a spatial
part and a temporal part.

Using a dedicated storage type for these operators saves both on memory and computation
time for convolutions with space-time data.
"""
mutable struct SpaceTimeDyad{A,B,T} <: AbstractArray{3,T}
    spatial_factor::A
    temporal_factor::B
end

# Implementation of a minimalistic array interface
# Note that setindex! is not supported. There is no way to enforce maintaining
# the dyadic structure during element modification.
import Base: size, getindex, linearindexing
size(X::SpaceTimeDyad) = (size(X.spatial_factor)..., length(X.temporal_factor))
linearindexing(::Type{T}) where {T<:SpaceTimeDyad} = Base.LinearSlow()
function getindex(X::SpaceTimeDyad, m::Int, n::Int, k::Int)
    T = promote_type(eltype(X.spatial_factor), eltype(X.temporal_factor))
    k > length(X.spacetial_factor) && return zero(T)
    return X.spatial_factor[m,n] * X.temporal_factor[k]
end

function convolve(X::SpaceTimeDyad, x::Array{T,2})

end


text:
abstract type AbstractSearchQuery end

struct CASQuery <: AbstractSearchQuery
    ID::String
    cas::Tuple{Int32,Int16,Int16}
end

function CASQuery(str::String)
    cas = cas_parse(str)
    return CASQuery(str,cas)
end

function CASQuery(val::NTuple{3,Integer})
    n1 = val[1]
    n2 = val[2]
    n3 = val[3]
    @assert n1 >=0
    @assert n2 >=0
    @assert n3 >=0
    v1 = convert(Int32,n1)
    v2 = convert(Int16,n2)
    v3 = convert(Int16,n3)
    cas = (v1,v2,v3)
    return CASQuery("",cas)
end

function tuple_to_casstr(n1,n2,n3)

    return string(n1) * "-" * string(n2) * "-" * string(n3)
end

tuple_to_casstr(val::NTuple{3,Integer}) = tuple_to_casstr(val...)
cas(id::CASQuery) = id.cas

struct InChIQuery <: AbstractSearchQuery
    ID::String
end

struct InChIKeyQuery <: AbstractSearchQuery
    ID::String
end

struct PubChemIDQuery <: AbstractSearchQuery
    ID::String
    val::Int
end

function PubChemIDQuery(str::String)
    val = parse(Int,str)
    @assert val >0
    return PubChemIDQuery(str,val)
end

function PubChemIDQuery(val::Integer)
    @assert val >0
    return PubChemIDQuery("",val)
end

id_num(ID::PubChemIDQuery) = ID.val

struct SMILESQuery <: AbstractSearchQuery
    ID::String
end

struct FormulaQuery <: AbstractSearchQuery
    ID::String
end

struct ElementQuery <: AbstractSearchQuery
    ID::String
end

struct AnyQuery <: AbstractSearchQuery
    ID::String
end

value(id::AbstractSearchQuery)::String = strip(id.ID)

function value(id::SMILESQuery)::String
    id = strip(id.ID)
    return replace(id,"SMILES=" =>"")
end

function value(id::InChIQuery)::String
    id_raw = strip(id.ID)
    id_lower = lowercase(id_raw)
    re1 = r"^inchi=1s/"
    re2 =  r"^inchi=1/"
    t1 = occursin(re1,id_lower)
    t2 = occursin(re2,id_lower)
    if t1
        return chop(id_raw,head=9,tail=0)
    elseif t2
        return chop(id_raw,head=8,tail=0)
    else
        throw("incorrect InChI passed")
    end
end

"""
    is_cas(str)::Bool
check if a given string is a cas number and returns true or false accordingly.

"""
function is_cas(str)
    str = strip(str)
    #regex from https://gist.github.com/KhepryQuixote/00946f2f7dd5f89324d8#file-pypubchemxtractor-py-L22
    cas_regex = r"^[1-9][0-9]{1,6}\\-[0-9]{2}\\-[0-9]$"
    return occursin(cas_regex,str)
end

"""
    is_inchi(str)::Bool
check if a given string is an InChI name and returns true or false accordingly.

"""
function is_inchi(str)
    str = strip(str)
    #regex from https://chemistry.stackexchange.com/a/86892
    inchi_regex = r"^InChI\=1S?\/[^\s]+(\s|$)"
    return occursin(inchi_regex,str)
end

"""
    is_inchikey(str)::Bool
check if a given string is an InChI key and returns true or false accordingly.

"""
function is_inchikey(str)
    str = strip(str)
    #14(A-Z)-10(A-Z)-1(A-Z)
    inchikey_regex  =r"[A-Z]{14}-[A-Z]{10}-[A-Z]{1}"
    return occursin(inchikey_regex,str)
end

"""
    is_element(str)::Bool
check if a given string is an element symbol and returns true or false accordingly.

"""
function is_element(str)
    return false
end

"""
    is_pubchemid(str)::Bool
check if a given string is a PubChem ID and returns true or false accordingly.
A PubChem ID is a positive integer.
"""
function is_pubchemid(str)
    str = strip(str)
    res = tryparse(Int,str)
   if res === nothing
    return false
   elseif res <= 0
    return false
   else
    return true
   end
end

"""
    is_smiles(str)::Bool
check if a given string is an element symbol and returns true or false accordingly.

"""
function is_smiles(str)
    str = strip(str)
    str = uppercase(str)
    smiles_regex  =r"^SMILES="
    return occursin(smiles_regex,str)
end


text: include("src/test.jl")


text:
module CenteredRBMs
    import Random
    import Statistics
    import LinearAlgebra
    import ValueHistories
    import Flux
    import FillArrays
    using ValueHistories: MVHistory
    import RestrictedBoltzmannMachines as RBMs
    using RestrictedBoltzmannMachines: RBM, AbstractLayer, visible, hidden, weights
    using RestrictedBoltzmannMachines: Binary, Spin, Potts, Gaussian, ReLU, dReLU, pReLU, xReLU
    using RestrictedBoltzmannMachines:  gradnorms, update!, minibatches

    include("centered_rbm.jl")
    include("centering.jl")
    include("binary.jl")
    include("layers.jl")
    include("sampling.jl")
    include("train/pcd.jl")
    include("train/centered_gradient.jl")
    include("train/initialization.jl")
end


text:
@testset "ff f4 certify" begin
    R, (x, y, z) = PolynomialRing(GF(2^31 - 1), ["x","y","z"], ordering=:lex)

    @test Groebner.groebner([x, y], certify=true) == [y, x]
    @test Groebner.groebner([y, x], certify=true) == [y, x]

    fs = [x^2 + y, x*y]
    @test Groebner.groebner(fs, certify=true) == [y^2, x*y, x^2 + y]
end

@testset "qq f4 certify" begin

    root = Groebner.change_ordering(Groebner.rootn(3, ground=QQ), :degrevlex)
    x1, x2, x3 = gens(parent(first(root)))
    gb = Groebner.groebner(root, certify=true)
    @test gb == [x1 + x2 + x3, x2^2 + x2*x3 + x3^2, x3^3 - 1]

    root = Groebner.change_ordering(Groebner.rootn(6, ground=QQ), :deglex)
    x1, x2, x3, x4, x5, x6 = gens(parent(first(root)))
    gb = Groebner.groebner(root, certify=true)
    @test gb == [
        x1 + x2 + x3 + x4 + x5 + x6,
        x2^2 + x2*x3 + x3^2 + x2*x4 + x3*x4 + x4^2 +
        x2*x5 + x3*x5 + x4*x5 + x5^2 + x2*x6 + x3*x6 + x4*x6 +
        x5*x6 + x6^2, x3^3 + x3^2*x4 + x3*x4^2 + x4^3 +
        x3^2*x5 + x3*x4*x5 + x4^2*x5 + x3*x5^2 + x4*x5^2 +
        x5^3 + x3^2*x6 + x3*x4*x6 + x4^2*x6 + x3*x5*x6 +
        x4*x5*x6 + x5^2*x6 + x3*x6^2 + x4*x6^2 + x5*x6^2 +
        x6^3, x4^4 + x4^3*x5 + x4^2*x5^2 + x4*x5^3 + x5^4 +
        x4^3*x6 + x4^2*x5*x6 + x4*x5^2*x6 + x5^3*x6 + x4^2*x6^2 +
        x4*x5*x6^2 + x5^2*x6^2 + x4*x6^3 + x5*x6^3 + x6^4, x5^5 +
        x5^4*x6 + x5^3*x6^2 + x5^2*x6^3 + x5*x6^4 + x6^5,
        x6^6 - 1]

    ku = Groebner.change_ordering(Groebner.ku10(ground=QQ), :degrevlex)
    x1, x2, x3, x4, x5, x6, x7, x8, x9, x10 = gens(parent(first(ku)))
    gb = Groebner.groebner(ku, certify=true)

    @test gb == [
        x9 + 230965//116706*x10 - 697789//116706,
        x8 + 92386//335511*x10 + 578636//335511,
        x7 + 5173616//8563059*x10 - 30862793//8563059,
        x6 + 87951472//192015891*x10 + 488096201//192015891,
        x5 + 131980//18759*x10 - 56944//18759,
        x4 + 32995//6084*x10 - 63415//6084,
        x3 - 263960//229671*x10 + 493631//229671,
        x2 - 65990//17043*x10 + 151205//17043,
        x1 - 26396//3273*x10 + 19850//3273,
        x10^2 - 71551//26396*x10 + 45155//26396]

end