text:
### more documentation needed; see base.jl in PowerSystems.jl, JJS 11/19/19
"""
System

A water system (municipal network).
"""
struct System <: WaterSystemType
    data::IS.SystemData
    #runchecks::Bool
    internal::InfrastructureSystemsInternal

    function System(data, internal; kwargs...)
        #runchecks = get(kwargs, :runchecks, true)
        #sys = new(data, runchecks, internal)
        sys = new(data, internal)
        return sys
    end
end

"""Construct an empty System. Useful for building a System while parsing raw data."""
function System(; kwargs...)
    return System(_create_system_data_from_kwargs(; kwargs...); kwargs...)
end

"""Construct a `System` from `InfrastructureSystems.SystemData`"""
function System(data; kwargs...)
    internal = get(kwargs, :internal, IS.InfrastructureSystemsInternal())
    return System(data, internal; kwargs...)
end

"""System constructor when components are constructed externally."""
function System(
    junctions::Vector{Junction},
    arcs::Vector{Arc},
    reservoirs::Vector{Reservoir},
    demands::Vector{<:WaterDemand}, # practically, there will always be a demand to have flow
    pipes::Vector{<:Pipe}, # theortically, a network _could_ exist without pipes, but not
                           # practically
    tanks::Union{Nothing, Vector{<:Tank}},
    pumps::Union{Nothing, Vector{Pump}},
    valves::Union{Nothing, Vector{Valve}}
    # curves? -- included with pump object
    # patterns? -- to be included as a forecasts with pump and demand objects? TBD!
    ;
    kwargs...,
)

    data = _create_system_data_from_kwargs(; kwargs...)
    sys = System(data; kwargs...)

    arrays = [junctions, arcs, reservoirs, demands, pipes]
    if !isnothing(tanks) && !isempty(tanks)
        push!(arrays, tanks)
    end
    if !isnothing(pumps) && !isempty(pumps)
        push!(arrays, pumps)
    end
    if !isnothing(valves) && !isempty(valves)
        push!(arrays, valves)
    end

    error_detected = false
    for component in Iterators.flatten(arrays)
        try
            add_component!(sys, component)
        catch e
            if isa(e, IS.InvalidRange)
                error_detected = true
            else
                rethrow()
            end
        end
    end

    runchecks = get(kwargs, :runchecks, true)

    if error_detected
        throw(IS.InvalidRange("Invalid value(s) detected"))
    end

    # if runchecks
    #     check!(sys)
    # end

    return sys
end

"""System constructor without nothing-able arguments."""
function System(
    junctions::Vector{Junction},
    arcs::Vector{Arc},
    reservoirs::Vector{Reservoir},
    demands::Vector{<:WaterDemand},
    pipes::Vector{<:Pipe}
    ;
    kwargs...,
)
    return System(
        junctions,
        arcs,
        reservoirs,
        demands,
        pipes,
        nothing, # tanks
        nothing, # pumps
        nothing # valves
        ;
        kwargs...,
    )
end

# # I'm not sure this is necessary; right now, it would conflict with existing empty System
# # creator System() defined above, JJS 1/29/20
# """System constructor with keyword arguments."""
# function System(
#     ;
#     junctions,
#     reservoirs,
#     pipes,
#     pumps,
#     demands,
#     tanks,
#     kwargs...,
# )
#     return System(
#         junctions,
#         reservoirs,
#         pipes,
#         pumps,
#         demands,
#         tanks
#         ;
#         kwargs...,
#     )
# end

function _create_system_data_from_kwargs(; kwargs...)
    validation_descriptor_file = nothing
    runchecks = get(kwargs, :runchecks, true)
    if runchecks
        validation_descriptor_file = get(kwargs, :configpath,
                                         WATER_SYSTEM_STRUCT_DESCRIPTOR_FILE)
    end

    return IS.SystemData(; validation_descriptor_file=validation_descriptor_file)
end

"""
    add_component!(sys::System, component::T; kwargs...) where T <: Component

Add a component to the system.

Throws ArgumentError if the component's name is already stored for its concrete type.

Throws InvalidRange if any of the component's field values are outside of defined valid
range.

# Examples
```julia
sys = System()

# Add a single component.
add_component!(sys, junction)

# Add many at once.
junctions = [junction1, junction2, junction3]
tanks = [tank1, tank2, tank3]
foreach(x -> add_component!(sys, x), Iterators.flatten((junctions, tanks)))
```
"""
function add_component!(sys::System, component::T; kwargs...) where T <: Component
    # check_component_addition(sys, component) # not implemented yet, JJS 1/29/20

    if Junction in fieldtypes(T) && T != Arc
        check_junction(sys, get_junction(component), component)
    end

    # if sys.runchecks && !validate_struct(sys, component)
    #     throw(InvalidValue("Invalid value for $(component)"))
    # end

    IS.add_component!(sys.data, component; kwargs...)
end

"""
Throws ArgumentError if the bus is not stored in the system.
"""
function check_junction(sys::System, junction::Junction, component::Component)
    name = get_name(junction)
    if isnothing(get_component(Junction, sys, name))
        throw(ArgumentError("$component has junction $name that is not stored in the system"))
    end
end

### below was copied from PowerSystems.jl and has not yet been checked or modified, JJS,
### 12/5/19

"""
    to_json(sys::System, filename::AbstractString)

Serializes a system to a JSON string.
"""
function to_json(sys::System, filename::AbstractString)
    IS.prepare_for_serialization!(sys.data, filename)
    return IS.to_json(sys, filename)
end

"""
    to_json(io::IO, sys::System)

Serializes a system an IO stream in JSON.
"""
function to_json(io::IO, sys::System)
    return IS.to_json(io, sys)
end

"""Constructs a System from a JSON file."""
function System(filename::String)
    sys = IS.from_json(System, filename)
    check!(sys)
    return sys
end

"""
    add_forecasts!(
                   sys::System,
                   metadata_file::AbstractString,
                   label_mapping::Dict{Tuple{String, String}, String};
                   resolution=nothing,
                  )

Adds forecasts from a metadata file or metadata descriptors.

# Arguments
- `sys::System`: system
- `metadata_file::AbstractString`: metadata file for timeseries
  that includes an array of IS.TimeseriesFileMetadata instances or a vector.
- `label_mapping::Dict{Tuple{String, String}, String}`: maps customized component field names to
  WaterSystem field names
- `resolution::DateTime.Period=nothing`: skip forecast that don't match this resolution.
"""
function add_forecasts!(
                        sys::System,
                        metadata_file::AbstractString,
                        label_mapping::Dict{Tuple{String, String}, String};
                        resolution=nothing,
                       )
    return IS.add_forecasts!(Component, sys.data, metadata_file, label_mapping;
                             resolution=resolution)
end

"""
    add_forecasts!(
                   sys::System,
                   timeseries_metadata::Vector{IS.TimeseriesFileMetadata},
                   resolution=nothing,
                  )

Adds forecasts from a metadata file or metadata descriptors.

# Arguments
- `sys::System`: system
- `timeseries_metadata::Vector{IS.TimeseriesFileMetadata}`: metadata for timeseries
- `resolution::DateTime.Period=nothing`: skip forecast that don't match this resolution.
"""
function add_forecasts!(
                        sys::System,
                        timeseries_metadata::Vector{IS.TimeseriesFileMetadata},
                        resolution=nothing
                       )
    return IS.add_forecasts!(Component, sys.data, timeseries_metadata;
                             resolution=resolution)
end

function IS.add_forecast!(
                          ::Type{<:Component},
                          data::IS.SystemData,
                          forecast_cache::IS.ForecastCache,
                          metadata::IS.TimeseriesFileMetadata;
                          resolution=nothing,
                         )
    IS.set_component!(metadata, data, WaterSystems)
    component = metadata.component
    if isnothing(component)
        return
    end

    forecast, ts_data = IS.make_forecast!(forecast_cache, metadata; resolution=resolution)
    if isnothing(forecast)
        return
    end

    if component isa LoadZones
        uuids = Set([IS.get_uuid(x) for x in get_buses(component)])
        for component_ in (load for load in IS.get_components(ElectricLoad, data)
                          if get_bus(load) |> IS.get_uuid in uuids)
            if forecast isa IS.DeterministicInternal
                forecast_ = IS.DeterministicInternal(IS.get_label(forecast), ts_data)
            else
                # TODO: others
                error("forecast type is not supported yet: $(typeof(forecast))")
            end
            IS.add_forecast!(data, component_, forecast, ts_data)
        end
    else
        IS.add_forecast!(data, component, forecast, ts_data)
    end
end

"""
    iterate_components(sys::System)

Iterates over all components.

# Examples
```julia
for component in iterate_components(sys)
    @show component
end
```

See also: [`get_components`](@ref)
"""
function iterate_components(sys::System)
    return IS.iterate_components(sys.data)
end

"""
    clear_components!(sys::System)

Remove all components from the system.
"""
function clear_components!(sys::System)
    return IS.clear_components!(sys.data)
end

"""
    remove_components!(::Type{T}, sys::System) where T <: Component

Remove all components of type T from the system.

Throws ArgumentError if the type is not stored.
"""
function remove_components!(::Type{T}, sys::System) where T <: Component
    return IS.remove_components!(T, sys.data)
end

"""
    remove_component!(sys::System, component::T) where T <: Component

Remove a component from the system by its value.

Throws ArgumentError if the component is not stored.
"""
function remove_component!(sys::System, component)
    return IS.remove_component!(sys.data, component)
end

"""
    remove_component!(
                      ::Type{T},
                      sys::System,
                      name::AbstractString,
                      ) where T <: Component

Remove a component from the system by its name.

Throws ArgumentError if the component is not stored.
"""
function remove_component!(
                           ::Type{T},
                           sys::System,
                           name::AbstractString,
                          ) where T <: Component
    return IS.remove_component!(T, sys.data, name)
end

"""
    get_component(
                  ::Type{T},
                  sys::System,
                  name::AbstractString
                 )::Union{T, Nothing} where {T <: Component}

Get the component of concrete type T with name. Returns nothing if no component matches.

See [`get_components_by_name`](@ref) if the concrete type is unknown.

Throws ArgumentError if T is not a concrete type.
"""
function get_component(::Type{T}, sys::System, name::AbstractString) where T <: Component
    return IS.get_component(T, sys.data, name)
end

"""
    get_components(
                   ::Type{T},
                   sys::System,
                  )::FlattenIteratorWrapper{T} where {T <: Component}

Returns an iterator of components. T can be concrete or abstract.
Call collect on the result if an array is desired.

# Examples
```julia
iter = WaterSystems.get_components(ThermalStandard, sys)
iter = WaterSystems.get_components(Generator, sys)
generators = WaterSystems.get_components(Generator, sys) |> collect
generators = collect(WaterSystems.get_components(Generator, sys))
```

See also: [`iterate_components`](@ref)
"""
function get_components(
                        ::Type{T},
                        sys::System,
                       )::FlattenIteratorWrapper{T} where {T <: Component}
    return IS.get_components(T, sys.data)
end

"""
    get_components_by_name(
                           ::Type{T},
                           sys::System,
                           name::AbstractString
                          )::Vector{T} where {T <: Component}

Get the components of abstract type T with name. Note that WaterSystems enforces unique
names on each concrete type but not across concrete types.

See [`get_component`](@ref) if the concrete type is known.

Throws ArgumentError if T is not an abstract type.
"""
function get_components_by_name(
                                ::Type{T},
                                sys::System,
                                name::AbstractString
                               )::Vector{T} where {T <: Component}
    return IS.get_components_by_name(T, sys.data, name)
end

"""
    add_forecast!(sys::System, component::Component, forecast::Forecast)

Adds forecast to the system.

Throws ArgumentError if the component is not stored in the system.

"""
function add_forecast!(sys::System, component::Component, forecast::Forecast)
    return IS.add_forecast!(sys.data, component, forecast)
end

"""
    add_forecast!(
                  sys::System,
                  filename::AbstractString,
                  component::Component,
                  label::AbstractString,
                  scaling_factor::Union{String, Float64}=1.0,
                 )

Add a forecast to a system from a CSV file.

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of
scaling_factor.
"""
function add_forecast!(
                       sys::System,
                       filename::AbstractString,
                       component::Component,
                       label::AbstractString,
                       scaling_factor::Union{String, Float64}=1.0,
                      )
    return IS.add_forecast!(sys.data, filename, component, label, scaling_factor)
end

"""
    add_forecast!(
                  sys::System,
                  ta::TimeSeries.TimeArray,
                  component,
                  label,
                  scaling_factor::Union{String, Float64}=1.0,
                 )

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

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of
scaling_factor.
"""
function add_forecast!(
                       sys::System,
                       ta::TimeSeries.TimeArray,
                       component,
                       label,
                       scaling_factor::Union{String, Float64}=1.0,
                      )
    return IS.add_forecast!(sys.data, ta, component, label, scaling_factor)
end

"""
    add_forecast!(
                  sys::System,
                  df::DataFrames.DataFrame,
                  component,
                  label,
                  scaling_factor::Union{String, Float64}=1.0,
                 )

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

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of
scaling_factor.
"""
function add_forecast!(
                       sys::System,
                       df::DataFrames.DataFrame,
                       component,
                       label,
                       scaling_factor::Union{String, Float64}=1.0,
                      )
    return IS.add_forecast!(sys.data, df, component, label, scaling_factor)
end

"""
    make_forecasts(sys::System, metadata_file::AbstractString; resolution=nothing)

Return a vector of forecasts from a metadata file.

# Arguments
- `data::SystemData`: system
- `metadata_file::AbstractString`: path to metadata file
- `resolution::{Nothing, Dates.Period}`: skip any forecasts that don't match this resolution

See [`InfrastructureSystems.TimeseriesFileMetadata`](@ref) for description of what the file
should contain.
"""
function make_forecasts(sys::System, metadata_file::AbstractString; resolution=nothing)
    return IS.make_forecasts(sys.data, metadata_file, WaterSystems; resolution=resolution)
end

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

Return a vector of forecasts from a vector of TimeseriesFileMetadata values.

# Arguments
- `data::SystemData`: system
- `timeseries_metadata::Vector{TimeseriesFileMetadata}`: metadata values
- `resolution::{Nothing, Dates.Period}`: skip any forecasts that don't match this resolution
"""
function make_forecasts(sys::System, metadata::Vector{IS.TimeseriesFileMetadata};
                        resolution=nothing)
    return IS.make_forecasts(sys.data, metadata, WaterSystems; resolution=resolution)
end

"""
    generate_initial_times(sys::System, interval::Dates.Period, horizon::Int)

Generates all possible initial times for the stored forecasts. This should be used when
contiguous forecasts have been stored in chunks, such as a one-year forecast broken up into
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.
"""
function generate_initial_times(sys::System, interval::Dates.Period, horizon::Int)
    return generate_initial_times(sys.data, interval, horizon)
end

"""
    get_forecast(
                 ::Type{T},
                 component::Component,
                 initial_time::Dates.DateTime,
                 label::AbstractString,
                ) where T <: Forecast

Return a forecast for the entire time series range stored for these parameters.
"""
function get_forecast(
                      ::Type{T},
                      component::Component,
                      initial_time::Dates.DateTime,
                      label::AbstractString,
                      ) where T <: Forecast
    return IS.get_forecast(T, component, initial_time, label)
end

"""
    get_forecast(
                 ::Type{T},
                 component::Component,
                 initial_time::Dates.DateTime,
                 label::AbstractString,
                 horizon::Int,
                ) where T <: Forecast

Return a forecast for a subset of the time series range stored for these parameters.
"""
function get_forecast(
                      ::Type{T},
                      component::InfrastructureSystemsType,
                      initial_time::Dates.DateTime,
                      label::AbstractString,
                      horizon::Int,
                     ) where T <: Forecast
    return IS.get_forecast(T, component, initial_time, label, horizon)
end

function get_forecast_initial_times(
                                    ::Type{T},
                                    component::Component,
                                   ) where T <: Forecast
    return IS.get_forecast_initial_times(T, component)
end

function get_forecast_initial_times(
                                    ::Type{T},
                                    component::Component,
                                    label::AbstractString
                                   ) where T <: Forecast
    return IS.get_forecast_initial_times(T, component, label)
end

function get_forecast_labels(
                             ::Type{T},
                             component::Component,
                             initial_time::Dates.DateTime,
                            ) where T <: Forecast
    return IS.get_forecast_labels(T, component, initial_time)
end

"""
    get_forecast_values(component::Component, forecast::Forecast)

Return a TimeSeries.TimeArray where the forecast data has been multiplied by the forecasted
component field.
"""
function get_forecast_values(component::Component, forecast::Forecast)
    return IS.get_forecast_values(component, forecast)
end

"""
    get_forecast_initial_times(sys::System)::Vector{Dates.DateTime}

Return sorted forecast initial times.

"""
function get_forecast_initial_times(sys::System)::Vector{Dates.DateTime}
    return IS.get_forecast_initial_times(sys.data)
end

"""
    get_forecast_keys(component::InfrastructureSystemsType)

Return an iterable of NamedTuple keys for forecasts stored for this component.
"""
function get_forecast_keys(component::Component)
    return IS.get_forecast_keys(component)
end

"""
    get_forecasts_horizon(sys::System)

Return the horizon for all forecasts.
"""
function get_forecasts_horizon(sys::System)
    return IS.get_forecasts_horizon(sys.data)
end

"""
    get_forecasts_initial_time(sys::System)

Return the earliest initial_time for a forecast.
"""
function get_forecasts_initial_time(sys::System)
    return IS.get_forecasts_initial_time(sys.data)
end

"""
    get_forecasts_interval(sys::System)

Return the interval for all forecasts.
"""
function get_forecasts_interval(sys::System)
    return IS.get_forecasts_interval(sys.data)
end

"""
    get_forecasts_resolution(sys::System)

Return the resolution for all forecasts.
"""
function get_forecasts_resolution(sys::System)
    return IS.get_forecasts_resolution(sys.data)
end

"""
    iterate_forecasts(sys::System)

Iterates over all forecasts in order of initial time.

# Examples
```julia
for forecast in iterate_forecasts(sys)
    @show forecast
end
```
"""
function iterate_forecasts(sys::System)
    return IS.iterate_forecasts(sys.data)
end

"""
    clear_forecasts!(sys::System)

Remove all forecasts from the system.
"""
function clear_forecasts!(sys::System)
    return IS.clear_forecasts!(sys.data)
end

"""
    check_forecast_consistency(sys::System)

Throws DataFormatError if forecasts have inconsistent parameters.
"""
function check_forecast_consistency(sys::System)
    IS.check_forecast_consistency(sys.data)
end

"""
    validate_forecast_consistency(sys::System)

Return true if all forecasts have consistent parameters.
"""
function validate_forecast_consistency(sys::System)
    return IS.validate_forecast_consistency(sys.data)
end

"""
    remove_forecast!(
                     ::Type{T},
                     sys::System,
                     component::Component,
                     initial_time::Dates.DateTime,
                     label::String,
                    )

Remove the time series data for a component.
"""
function remove_forecast!(
                          ::Type{T},
                          sys::System,
                          component::Component,
                          initial_time::Dates.DateTime,
                          label::String,
                         ) where T <: Forecast
    return IS.remove_forecast!(T, sys.data, component, initial_time, label)
end

"""
    validate_struct(sys::System, value::WaterSystemType)

Validates an instance of a WaterSystemType against System data.
Returns true if the instance is valid.

Users implementing this function for custom types should consider implementing
InfrastructureSystems.validate_struct instead if the validation logic only requires data
contained within the instance.
"""
function validate_struct(sys::System, value::WaterSystemType)::Bool
    return true
end

# # an analagous check is probably not needed for WaterSystems
# function check!(sys::System)
#     buses = get_components(Bus, sys)
#     slack_bus_check(buses)
#     buscheck(buses)
# end

function JSON2.read(io::IO, ::Type{System})
    raw = JSON2.read(io, NamedTuple)
    data = IS.deserialize(IS.SystemData, Component, raw.data)
    sys = System(data, float(raw.elevation); runchecks=raw.runchecks)
    return sys
end

function IS.deserialize_components(
                                   ::Type{Component},
                                   data::IS.SystemData,
                                   raw::NamedTuple,
                                  )
    # TODO: This adds components through IS.SystemData instead System, which is what should
    # happen. There is a catch-22 between creating System and SystemData.
    component_cache = Dict{Base.UUID, Component}()

    # Buses and Arcs are encoded as UUIDs.
    composite_components = [Bus]
    for composite_component in composite_components
        for component in IS.get_components_raw(IS.SystemData, composite_component, raw)
            comp = IS.convert_type(composite_component, component)
            IS.add_component!(data, comp)
            component_cache[IS.get_uuid(comp)] = comp
        end
    end

    # Skip Services this round because they have Devices.
    for c_type_sym in IS.get_component_types_raw(IS.SystemData, raw)
        c_type = getfield(WaterSystems, Symbol(IS.strip_module_name(string(c_type_sym))))
        (c_type in composite_components || c_type <: Service) && continue
        for component in IS.get_components_raw(IS.SystemData, c_type, raw)
            comp = IS.convert_type(c_type, component, component_cache)
            IS.add_component!(data, comp)
            component_cache[IS.get_uuid(comp)] = comp
        end
    end

    # Now get the Services.
    for c_type_sym in IS.get_component_types_raw(IS.SystemData, raw)
        c_type = getfield(WaterSystems, Symbol(IS.strip_module_name(string(c_type_sym))))
        if c_type <: Service
            for component in IS.get_components_raw(IS.SystemData, c_type, raw)
                comp = IS.convert_type(c_type, component, component_cache)
                IS.add_component!(data, comp)
            end
        end
    end
end

function JSON2.write(io::IO, component::T) where T <: Component
    return JSON2.write(io, encode_for_json(component))
end

function JSON2.write(component::T) where T <: Component
    return JSON2.write(encode_for_json(component))
end

"""
Encode composite buses as UUIDs.
"""
function encode_for_json(component::T) where T <: Component
    fields = fieldnames(T)
    vals = []

    for name in fields
        val = getfield(component, name)
        if val isa Bus
            push!(vals, IS.get_uuid(val))
        else
            push!(vals, val)
        end
    end

    return NamedTuple{fields}(vals)
end

function IS.convert_type(
                         ::Type{T},
                         data::NamedTuple,
                         component_cache::Dict,
                        ) where T <: Component
    @debug T data
    values = []
    for (fieldname, fieldtype)  in zip(fieldnames(T), fieldtypes(T))
        val = getfield(data, fieldname)
        if fieldtype <: Bus
            uuid = Base.UUID(val.value)
            bus = component_cache[uuid]
            push!(values, bus)
        elseif fieldtype <: Component
            # Recurse.
            push!(values, IS.convert_type(fieldtype, val, component_cache))
        else
            obj = IS.convert_type(fieldtype, val)
            push!(values, obj)
        end
    end

    return T(values...)
end

function IS.compare_values(x::System, y::System)::Bool
    match = true

    if !IS.compare_values(x.data, y.data)
        @debug "SystemData values do not match"
        match = false
    end

    if x.elevation != y.elevation
        @debug "elevation does not match" x.elevation y.elevation
        match = false
    end

    return match
end

function parse_types(mod)
    for name in names(mod)
        mod_type = getfield(mod, name)
        try
            !isstructtype(mod_type) && continue
        catch(e)
            continue
        end
        !isconcretetype(mod_type) && continue
        println("object $mod_type")
    end
    for name in names(mod)
        mod_type = getfield(mod, name)
        !isstructtype(mod_type) && continue
        !isconcretetype(mod_type) && continue
        for (fname, ftype) in zip(fieldnames(mod_type), fieldtypes(mod_type))
            if ftype in names(mod)
                println("$mod_type o-- $ftype")
            end
        end
    end
end


text:
for T1 in (Float32, Float64, Complex{Float32}, Complex{Float64})
    for T2 in (Float32, Float64, Complex{Float32}, Complex{Float64})
        @testset "multiplication with $SV: $T1 times $T2" for SV in (StridedView, UnsafeStridedView)
            d = 103
            A1 = rand(T1, (d,d))
            A2 = rand(T2, (d,d))
            T3 = promote_type(T1,T2)
            A3 = rand(T3, (d,d))
            A4 = rand(T3, (d,d))
            A1c = copy(A1)
            A2c = copy(A2)
            A3c = copy(A3)
            A4c = copy(A4)
            GC.@preserve A1c A2c A3c A4c begin
                B1 = SV(A1c)
                B2 = SV(A2c)
                B3 = SV(A3c)
                B4 = SV(A4c)

                for op1 in (identity, conj, transpose, adjoint)
                    @test op1(A1) == op1(B1)
                    for op2 in (identity, conj, transpose, adjoint)
                        @test op1(A1)*op2(A2) ≈ op1(B1)*op2(B2)
                        for op3 in (identity, conj, transpose, adjoint)
                            α = randn(T3)
                            β = randn(T3)
                            copyto!(B3, B4)
                            mul!(op3(B3), op1(B1), op2(B2), α, β)
                            @test B3 ≈ op3(β)*A4 + op3(α * op1(A1)*op2(A2)) # op3 is its own inverse
                        end
                    end
                end
            end
        end
    end
end

ext: jl


text:
#!/usr/bin/env julia

using ColorSchemes, Colors

using Base.Test

# create a colorscheme from image file

hokusai_test = extract(dirname(@__FILE__) * "/hokusai.jpg")

@test length(hokusai_test) == 10
@test typeof(hokusai_test) == Array{ColorTypes.RGB{Float64},1}

# extract colors and weights from image file located here in the test directory
c, w =  extract_weighted_colors(dirname(@__FILE__) * "/hokusai.jpg", 10, 10, 0.01; shrink = 4)

@test length(c) == 10
@test length(w) == 10

# load scheme from ./data
hok = loadcolorscheme("hokusai")

@test length(hok) == 32
@test typeof(hokusai_test) == Array{ColorTypes.RGB{Float64},1}

# test that schemes yield different values
@test colorscheme(hokusai_test, 0.0) != colorscheme(hokusai_test, 0.5)

# create weighted palette; there is some unpredictability here...
csw = colorscheme_weighted(c, w, 37)
@test 36 <= length(csw) <= 38

# default is 50
csw = colorscheme_weighted(c, w)
@test length(csw) == 50

# some possible tests would require image services such as ImageMagick:
# tempdir = mktempdir()
# cd(tempdir)
# savecolorscheme(csw, "test_hokusai.txt", "this is a test")
# csw1 = loadcolorscheme("test_hokusai") # no suffix required
# @test typeof(csw1) == Array{ColorTypes.RGB{Float64},1}
# rm(tempdir, recursive=true)


text:
using Documenter, HashDicts

makedocs(sitename="HashDicts documentation")

deploydocs(
    repo = "github.com/matbesancon/HashDicts.jl.git",
)


text:
VERSION > v"0.4-" && __precompile__()

module Lazy

############
# Utilities
############

using Compat

include("macros.jl")
include("dynamic.jl")
include("tail.jl")

export @listable

import Base: *, ==, +, -

macro listable(f)
  if typeof(f) == Expr && f.head == :tuple
    return Expr(:block, [:(@listable $f) for f in f.args]...)
  end
  f = esc(f)
  :($f(ls::List...) = map($f, ls...))
end

########
# Types
########

import Base: isempty, first, colon

export List, list, @lazy, prepend, tail

abstract List

type EmptyList <: List
end

type LinkedList <: List
  head
  tail::List
end

type LazyList <: List
  list::List
  realised::Bool
  f::Function

  LazyList(f) = new(EmptyList(), false, f)
end

# Empty List

first(::EmptyList) = nothing
tail(l::EmptyList) = l
isempty(::EmptyList) = true

# Lists

prepend(x, l::List) = LinkedList(x, l)
colon(x, xs::List) = prepend(x, xs)
colon(x,y,xs::List) = x:prepend(y,xs)

list() = EmptyList()
list(x, xs...) = x:list(xs...)

first(l::LinkedList) = l.head
tail(l::LinkedList) = l.tail

isempty(::LinkedList) = false

# Lazy Lists

function realise(xs::LazyList)
  xs.realised && return xs.list
  xs.realised = true
  xs.list = xs.f()
  # Unroll in a loop to avoid overflow
  while isa(xs.list, LazyList)
    xs.list = xs.list.f()
  end
  return xs.list
end

macro lazy(code)
  :(LazyList(() -> seq($(esc(code)))))
end

for f in [:first :isempty]
  @eval $f(l::LazyList) = $f(realise(l))
end

tail(l::LazyList) = @lazy tail(realise(l))

########
# Usage
########

include("liblazy.jl")
include("collections.jl")

# -----
# Eager
# -----

# Depending upon Julia and Compat versions, foreach may be defined as
# a core language generic function. If so, extend rather than define
# it.

if isdefined(:foreach) && isa(foreach, Function)
    if VERSION < v"0.5.0-dev+977"
        import Compat.foreach
    else
        import Base.foreach
    end
end

export dorun, doall, foreach

@rec dorun(xs::List) = isempty(xs) ? nothing : dorun(tail(xs))
doall(xs::List) = (dorun(xs); xs)

foreach(f, ls::List...) = map(f, ls...) |> dorun

# -------
# Interop
# -------

import Base: getindex, setindex!, start, next, done

getindex(l::List, i::Int) = i <= 1 ? first(l) : tail(l)[i-1]
getindex(l::List, r::UnitRange) = take(r.len, drop(r.start - 1, l))

setindex!(xs::LinkedList, v, i::Integer) = i <= 1 ? xs.first = v : (tail(xs)[i-1] = v)
setindex!(xs::LazyList, v, i::Integer) = i <= 1 ? realise(xs)[1] = v : (tail(xs)[i-1] = v)

# Iteration over a list holds on to the head
start(xs::List) = xs
done(::List, xs::List) = isempty(xs)
next(::List, xs::List) = first(xs), tail(xs)

###########
# Printing
###########

Base.show(io::IO, xs::List) =
  foreach(x->print(io,x), ["("] * interpose(xs, " ") * [")"])

@compat function show(io::IO, ::MIME"text/plain", xs::List)
  isempty(xs) && return println(io, "List()")

  print(io, "List:")
  for x in xs
    print(io, "\n  ")
    show(io, x)
  end
end

# Some example lists

@listable +, -, Base.exp, Base.log, Base.sin, Base.cos, Base.tan

const fibs = @lazy 0:big(1):(fibs + tail(fibs))

isprime(n) =
  @>> primes begin
    takewhile(x -> x<=sqrt(n))
    map(x -> n % x == 0)
    any; !
  end

const primes = filter(isprime, range(2))

end


text:
get_bounds(m::Function) = [],[]
get_bounds(m::Flux.Chain) = [reduce(vcat, [get_bounds(l)[1] for l in m.layers]), reduce(vcat, [get_bounds(l)[2] for l in m.layers])]
function get_bounds(m::Flux.Dense)
  lb = [[-20.0 for _ in 1:length(m.weight)]...,
        [-20.0 for _ in 1:length(m.bias)]...] |> f32
  ub = [[20.0 for _ in 1:length(m.weight)]...,
        [20.0 for _ in 1:length(m.bias)]...] |> f32
  return lb, ub
end

function get_bounds(m::DiffEqFlux.FastChain)
  lb = vcat([get_bounds(layer)[1] for layer in m.layers]...)
  ub = vcat([get_bounds(layer)[2] for layer in m.layers]...)
  return lb, ub
end

function get_bounds(m::DiffEqFlux.FastDense)
  lb = [[-2.0 for _ in 1:m.out*m.in]...,
        [-2.0 for _ in 1:m.out]...] |> f32
  ub = [[2.0 for _ in 1:m.out*m.in]...,
        [2.0 for _ in 1:m.out]...] |> f32
  return lb, ub
end

get_bounds(m::Flux.Recur) = get_bounds(m.cell)
#get_bounds(m::MyRecur) = get_bounds(m.cell)

function get_bounds(m::Mapper)
  lb = [[-20.1 for i in 1:length(m.W)]...,
        [-20.1 for i in 1:length(m.b)]...] |> f32

  ub = [[20.1 for i in 1:length(m.W)]...,
        [20.1 for i in 1:length(m.b)]...] |> f32
  return lb, ub
end

# function get_bounds(m::gNN)
#   lb = [[0 for _ in m.G]...,          #
#         [0.1 for _ in m.μ]...,
#         [1 for _ in m.σ]...,
#         [-2 for _ in m.E]...,] |> f32
#
#   ub = [[1 for _ in m.G]...,          #
#         [0.9 for _ in m.μ]...,
#         [10 for _ in m.σ]...,
#         [2 for _ in m.E]...,] |> f32
#   return lb, ub
# end
function get_bounds(m::gNN)
  lb = [[0.1 for _ in m.μ]...,
        [1 for _ in m.σ]...,
        [0.001 for _ in m.G]...,
        [-1 for _ in m.E]...,
        ] |> f32

  ub = [[0.9 for _ in m.μ]...,
        [10 for _ in m.σ]...,
        [1 for _ in m.G]...,
        [1 for _ in m.E]...,
        ] |> f32
  return lb, ub
end

# function get_bounds(m::LTCCell)
#   lower = [
#       get_bounds(m.mapin)[1]...,
#       get_bounds(m.mapout)[1]...,
#       get_bounds(m.sens_f)[1]...,
#       get_bounds(m.syn_f)[1]...,
#       [0.9 for _ in m.cm]...,               #
#       [0.01 for _ in m.Gleak]...,            #
#       [-2 for _ in m.Eleak]...,
#       [0.001 for _ in m.state0]...,
#   ] |> f32
#   upper = [
#     get_bounds(m.mapin)[2]...,
#     get_bounds(m.mapout)[2]...,
#     get_bounds(m.sens_f)[2]...,
#     get_bounds(m.syn_f)[2]...,
#     [4 for _ in m.cm]...,
#     [2.0 for _ in m.Gleak]...,
#     [2 for _ in m.Eleak]...,
#     [0.1 for _ in m.state0]...,
#   ] |> f32
#
#   lower, upper
# end

#function get_bounds(m::LTCCell)
#  lower = [
#      get_bounds(m.mapin)[1]...,
#      get_bounds(m.mapout)[1]...,
#      #reduce(vcat,[get_bounds(s)[1] for s in m.sens.synapses])...,
#      #reduce(vcat,[get_bounds(s)[1] for s in m.syns.synapses])...,
#      get_bounds(m.sens_f)[1]...,
#      get_bounds(m.syn_f)[1]...,
#      [0.001 for _ in m.cm]...,               #
#      [0.01 for _ in m.Gleak]...,            #
#      [-2 for _ in m.Eleak]...,
#      [-2 for _ in m.Esens]...,
#      [0 for _ in m.Gsens]...,
#      [-2 for _ in m.Esyn]...,
#      [0 for _ in m.Gsyn]...,
#      [-0.5 for _ in m.state0]...,
#  ] |> f32
#  upper = [
#    get_bounds(m.mapin)[2]...,
#    get_bounds(m.mapout)[2]...,
#    #reduce(vcat,[get_bounds(s)[2] for s in m.sens.synapses])...,
#    #reduce(vcat,[get_bounds(s)[2] for s in m.syns.synapses])...,
#    get_bounds(m.sens_f)[2]...,
#    get_bounds(m.syn_f)[2]...,
#    [5 for _ in m.cm]...,
#    [10.0 for _ in m.Gleak]...,
#    [1 for _ in m.Eleak]...,
#    [2 for _ in m.Esens]...,
#    [1 for _ in m.Gsens]...,
#    [2 for _ in m.Esyn]...,
#    [1 for _ in m.Gsyn]...,
#    [0.5 for _ in m.state0]...,
#  ] |> f32
#
#  lower, upper
#end

function get_bounds(m::LTCCell)
  lower = [
      get_bounds(m.mapin)[1]...,
      get_bounds(m.mapout)[1]...,
      reduce(vcat,[get_bounds(s)[1] for s in m.sens.synapses])...,
      reduce(vcat,[get_bounds(s)[1] for s in m.syns.synapses])...,
      #get_bounds(m.sens_f)[1]...,
      #get_bounds(m.syn_f)[1]...,
      [0.001 for _ in m.cm]...,               #
      [0.01 for _ in m.Gleak]...,            #
      [-2 for _ in m.Eleak]...,
      #[-2 for _ in m.Esens]...,
      #[0 for _ in m.Gsens]...,
      #[-2 for _ in m.Esyn]...,
      #[0 for _ in m.Gsyn]...,
      [-50 for _ in m.state0]...,
  ] |> f32
  upper = [
    get_bounds(m.mapin)[2]...,
    get_bounds(m.mapout)[2]...,
    reduce(vcat,[get_bounds(s)[2] for s in m.sens.synapses])...,
    reduce(vcat,[get_bounds(s)[2] for s in m.syns.synapses])...,
    #get_bounds(m.sens_f)[2]...,
    #get_bounds(m.syn_f)[2]...,
    [5 for _ in m.cm]...,
    [10.0 for _ in m.Gleak]...,
    [-0.1 for _ in m.Eleak]...,
    #[2 for _ in m.Esens]...,
    #[1 for _ in m.Gsens]...,
    #[2 for _ in m.Esyn]...,
    #[1 for _ in m.Gsyn]...,
    [20 for _ in m.state0]...,
  ] |> f32

  lower, upper
end

# function get_bounds(m::NCP)
#   lower = [
#       get_bounds(m.mapin)[1]...,
#       [0.9 for _ in m.ltc.cell.cm]...,               #
#       [0.0009 for _ in m.ltc.cell.Gleak]...,            #
#       [-2 for _ in m.ltc.cell.Eleak]...,
#       get_bounds(m.ltc.cell.sens)[1]...,
#       get_bounds(m.ltc.cell.syn)[1]...,
#       [-0.1 for _ in m.ltc.cell.state0]...,
#       get_bounds(m.mapout)[1]...,
#   ] |> f32
#   upper = [
#     get_bounds(m.mapin)[2]...,
#     [4 for _ in m.ltc.cell.cm]...,
#     [2.0 for _ in m.ltc.cell.Gleak]...,
#     [2 for _ in m.ltc.cell.Eleak]...,
#     get_bounds(m.ltc.cell.sens)[2]...,
#     get_bounds(m.ltc.cell.syn)[2]...,
#     [0.2 for _ in m.ltc.cell.state0]...,
#     get_bounds(m.mapout)[2]...,
#   ] |> f32
#
#   lower, upper
# end


text:
module AbstractPattern

export spec_gen, runterm, MK, RedyFlavoured, TypeObject
export and, or, literal, and, wildcard, decons,
       guard, effect, self
export PatternCompilationError, Target, PatternImpl, PComp
export APP, NoncachablePre, NoPre
export ChainDict, for_chaindict, child, for_chaindict_dup
export BasicPatterns
include("DataStructure.jl")
include("Target.jl")
include("PatternSignature.jl")
include("Print.jl")
include("structures/Print.jl")
include("structures/TypeTagExtraction.jl")
include("ADT.jl")
include("CaseMerge.jl")
include("UserSignature.jl")
include("Retagless.jl")
include("implementations/RedyFlavoured.jl")
include("∀/BasicPatterns.jl")

@nospecialize

function MK(m::Any)
    m.backend
end

function runterm(term, xs)
    points_of_view = Dict{Function, Int}(x => i for (i, x) in enumerate(xs))
    impls = Tuple(x(points_of_view) for x in xs)
    term(impls)
end

function spec_gen(branches :: Vector)
    cores = Branch[]
    ln = LineNumberNode(1, "<unknown>")
    for branch in branches
        branch isa LineNumberNode && begin
            ln = branch
            continue
        end
        (branch, cont) = branch :: Pair{F, Symbol} where F <: Function
        points_of_view = Dict{Function, Int}(tag_extract => 1, untagless => 2)
        impls = (tag_extract(points_of_view), untagless(points_of_view))
        type, pat = branch(impls)
        push!(cores, PatternInfo(pat::TagfulPattern, type::TypeObject) => (ln, cont))
    end
    split_cores = Branch[]
    case_split!(split_cores, cores)
    case_merge(split_cores)
end

end # module

max_stars_repo_name: KristofferC/AbstractPattern

max_issues_repo_name: KristofferC/AbstractPattern

max_forks_repo_name: KristofferC/AbstractPattern


text:
@testset "IArrays" begin

    A = rand(5,4,3)
    Aoneto = IArray(A);

    r1 = range(.1, stop = .2, length=5)
    r2 = ["a", "b", "c", "d"]
    r3 = 2:4

    ind1 = Index(r1)
    ind2 = Index(r2)
    ind3 = Index(r3)

    Aindices = IArray(A, (r1, r2, r3));
    Anamed = IArray(Aindices, (:a, :b, :c));

   # TODO ensure ind3 = 1 is out of bounds error
    @testset "axes(A, [d])" begin
        @test axes(Aindices, 1) == ind1
        @test axes(Aindices, 2) == ind2
        @test axes(Aindices, 3) == ind3

        @test axes(A, 1) == values(axes(Aoneto, 1))
        @test axes(A, 2) == values(axes(Aoneto, 2))
        @test axes(A, 3) == values(axes(Aoneto, 3))
    end

    @testset "Single Indexing" begin
        t1 = getindex(A, 1, 1, 1)
        t2 = getindex(A, 5, 4, 3)

        @test getindex(Aoneto, 1, 1, 1) == t1
        @test getindex(Aindices, r1[1], r2[1], r3[1]) == t1

        # Cartesian indexing bypasses any sort of AbstractIndex
        @test getindex(Aoneto, CartesianIndex(1,1,1)) == t1
        @test getindex(Aindices, CartesianIndex(1,1,1)) == t1

        @test getindex(Aoneto, 5, 4, 3) == t2
        @test getindex(Aindices, r1[5], r2[4], r3[3]) == t2

        # Cartesian indexing bypasses any sort of AbstractIndex
        @test getindex(Aoneto, CartesianIndex(5, 4, 3)) == t2
        @test getindex(Aindices, CartesianIndex(5, 4, 3)) == t2

        @test getindex(Aoneto, 1) == t1
        @test getindex(Aindices, 1) == t1
        @test getindex(Aoneto, 60) == t2
        @test getindex(Aindices, 60) == t2

        @test checkbounds(Bool, Aindices, 61) == false
        @test checkbounds(Bool, Anamed, 61) == false

        @test checkbounds(Bool, Aoneto, 2, 2, 2, 1) == true
        @test checkbounds(Bool, Aindices, r1[2], r2[2], r3[2], 1) == true

        @test checkbounds(Bool, Aoneto, 2, 2, 2, 2) == false
        @test checkbounds(Bool, Aindices, r1[2], r2[2], r3[2], 2) == false

        @test checkbounds(Bool, Aoneto, 1, 1)  == false
        @test checkbounds(Bool, Aindices, r1[1], r2[1])  == false
        @test checkbounds(Bool, Aindices, 1, 12) == false
        @test checkbounds(Bool, Aindices, 5, 12) == false
        @test checkbounds(Bool, Aindices, 1, 13) == false
        @test checkbounds(Bool, Aindices, 6, 12) == false

        @test checkbounds(Bool, Aoneto, 0, 1, 1) == false
        @test checkbounds(Bool, Aoneto, 1, 0, 1) == false
        @test checkbounds(Bool, Aoneto, 1, 1, 0) == false
        @test checkbounds(Bool, Aoneto, 6, 4, 3) == false
        @test checkbounds(Bool, Aoneto, 5, 5, 3) == false
        @test checkbounds(Bool, Aoneto, 5, 4, 4) == false
   end

    @testset "Vector indexing" begin
        @test getindex(Aindices, r1, r2, r3) == getindex(A, 1:5, 1:4, 1:3)
        @test getindex(Aoneto, 1:5, 1:4, 1:3) == getindex(A, 1:5, 1:4, 1:3)

        @test getindex(Aindices, 1:60) == getindex(Aoneto, 1:60) == getindex(A, 1:60)
        @test checkbounds(Bool, A, 2, 2, 2, 1:1) == true  # extra indices
    end

    @testset "NamedDimsExtra" begin

        @testset "filter_axes" begin
            @test filter_axes(allunique, A) == axes(A)
            @test filter_axes(allunique, Anamed) == axes(Anamed)
            @test length(filter_axes(allunique, A)) == ndims(A)
            @test length(filter_axes(i -> length(i) == length(axes(A,1)), A)) == 1
       end

       #=
        @testset "find_axes" begin
            @test find_axes(allunique, A) == ntuple(i->i, ndims(A))
            @test find_axes(allunique, Anamed) == ntuple(i->i, ndims(A))

            @test length(find_axes(allunique, A)) == ndims(A)
            @test find_axes(i -> i == axes(A,1), A) == (1,)
            @test find_axes(i -> i == axes(Anamed,1), Anamed) == (1,)
        end
        =#
    end
end


text:
module FIRRTL

export GroundType, Clock, UInteger, SInteger, Analog, BundleUnit, Bundle
export ispassive, ==, ≈

include("types.jl")

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module APITests

using Pkg, Test
import Pkg.Types.PkgError

include("utils.jl")

@testset "API should accept `AbstractString` arguments" begin
    temp_pkg_dir() do project_path
        with_temp_env() do
            Pkg.add(strip("  Example  "))
            Pkg.rm(strip("  Example "))
        end
    end
end

@testset "Pkg.activate" begin
    temp_pkg_dir() do project_path
        cd_tempdir() do tmp
            path = pwd()
            Pkg.activate(".")
            mkdir("Foo")
            cd(mkdir("modules")) do
                Pkg.generate("Foo")
            end
            Pkg.develop(Pkg.PackageSpec(path="modules/Foo")) # to avoid issue #542
            Pkg.activate("Foo") # activate path Foo over deps Foo
            @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
            Pkg.activate(".")
            rm("Foo"; force=true, recursive=true)
            Pkg.activate("Foo") # activate path from developed Foo
            @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
            Pkg.activate(".")
            Pkg.activate("./Foo") # activate empty directory Foo (sidestep the developed Foo)
            @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
            Pkg.activate(".")
            Pkg.activate("Bar") # activate empty directory Bar
            @test Base.active_project() == joinpath(path, "Bar", "Project.toml")
            Pkg.activate(".")
            Pkg.add("Example") # non-deved deps should not be activated
            Pkg.activate("Example")
            @test Base.active_project() == joinpath(path, "Example", "Project.toml")
            Pkg.activate(".")
            cd(mkdir("tests"))
            Pkg.activate("Foo") # activate developed Foo from another directory
            @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
            Pkg.activate() # activate home project
            @test Base.ACTIVE_PROJECT[] === nothing
        end
    end
end

@testset "Pkg.status" begin
    temp_pkg_dir() do project_path
        Pkg.add(["Example", "Random"])
        Pkg.status()
        Pkg.status("Example")
        Pkg.status(["Example", "Random"])
        Pkg.status(PackageSpec("Example"))
        Pkg.status(PackageSpec(uuid = "7876af07-990d-54b4-ab0e-23690620f79a"))
        Pkg.status(PackageSpec.(["Example", "Random"]))
        Pkg.status(; mode=PKGMODE_MANIFEST)
        Pkg.status("Example"; mode=PKGMODE_MANIFEST)
        @test_deprecated Pkg.status(PKGMODE_MANIFEST)
    end
end

@testset "Pkg.develop" begin
    temp_pkg_dir() do project_path; cd_tempdir() do tmpdir;
        entry = nothing
        # explicit relative path
        with_temp_env() do env_path
            cd(env_path) do
                foo_uuid = Pkg.generate("Foo")
                Pkg.develop(PackageSpec(;path="Foo"))
                manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
                entry = manifest[foo_uuid]
            end
            @test entry.path == "Foo"
            @test entry.name == "Foo"
            @test isdir(joinpath(env_path, entry.path))
        end
        # explicit absolute path
        with_temp_env() do env_path
            cd_tempdir() do temp_dir
                foo_uuid = Pkg.generate("Foo")
                absolute_path = abspath(joinpath(temp_dir, "Foo"))
                Pkg.develop(PackageSpec(;path=absolute_path))
                manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
                entry = manifest[foo_uuid]
                @test entry.name == "Foo"
                @test entry.path == absolute_path
                @test isdir(entry.path)
            end
        end
        # name
        with_temp_env() do env_path
            Pkg.develop("Example")
            manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
            for (uuid, entry) in manifest
                if entry.name == "Example"
                    @test entry.path == joinpath(Pkg.depots1(), "dev", "Example")
                    @test isdir(entry.path)
                end
            end
        end
        # name + local
        with_temp_env() do env_path
            Pkg.develop("Example"; shared=false)
            manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
            for (uuid, entry) in manifest
                if entry.name == "Example"
                    @test entry.path == joinpath("dev", "Example")
                    @test isdir(joinpath(env_path, entry.path))
                end
            end
        end
        # url
        with_temp_env() do env_path
            url = "https://github.com/JuliaLang/Example.jl"
            Pkg.develop(PackageSpec(;url=url))
            manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
            for (uuid, entry) in manifest
                if entry.name == "Example"
                    @test entry.path == joinpath(Pkg.depots1(), "dev", "Example")
                    @test isdir(entry.path)
                end
            end
        end
        # unregistered url
        with_temp_env() do env_path
            url = "https://github.com/00vareladavid/Unregistered.jl"
            Pkg.develop(PackageSpec(;url=url))
            manifest = Pkg.Types.read_manifest(joinpath(env_path, "Manifest.toml"))
            for (uuid, entry) in manifest
                if entry.name == "Unregistered"
                    @test uuid == UUID("dcb67f36-efa0-11e8-0cef-2fc465ed98ae")
                    @test entry.path == joinpath(Pkg.depots1(), "dev", "Unregistered")
                    @test isdir(entry.path)
                end
            end
        end
        # with rev
        with_temp_env() do env_path
            @test_throws PkgError Pkg.develop(PackageSpec(;name="Example",rev="Foobar"))
        end
    end end
end

end # module APITests


text:
using Random
using Distributions
using Plots

mutable struct Point
    label::Int64
    coordinates::Vector{Float64}
end

function DBSCAN(data, eps, min_cluster_size)
    cluster_counter = 0
    point_data = [Point(0, p) for p in eachcol(data)]

    for point in point_data
        if point.label == 0
            (neighbor_counter, neighbors) = find_neighbors(point_data, point, eps)
            if neighbor_counter < min_cluster_size
                point.label = -1
            else
                cluster_counter += 1
                point.label = cluster_counter
                for seed in neighbors
                    if seed.label < 1
                        seed.label = cluster_counter
                    end

                    (neighbor_counter, seed_neighbors) = find_neighbors(point_data, seed, eps)
                    if neighbor_counter >= min_cluster_size
                        union!(neighbors, seed_neighbors)
                    end
                end
            end
        end
    end

    return (point_data, cluster_counter)
end

function find_neighbors(data, point, eps)
    number_of_neighbors = 0
    neighbors = []

    for second_point in data
        if point != second_point
            if euclidian_distance_alt(point.coordinates, second_point.coordinates) < eps
                number_of_neighbors += 1
                append!(neighbors, [second_point])
            end
        end
    end

    return (number_of_neighbors, neighbors)
end

function euclidian_distance_alt(p1, p2)
    sum = 0

    for i in 1:length(p1)
        sum += abs2(p1[i] - p2[i])
    end

    return sqrt(sum)
end

function euclidian_distances(v1, v2)
    v1 = [v1[i, :] for i in 1:size(v1,1)]
    v2 = [v2[i, :] for i in 1:size(v2,1)]
    s = size(v1,1)
    d = zeros((length(v1[1]),length(v2[2])))

    for j in 1:length(v2[2])
        for k in 1:length(v1[1])
            summed = 0
            for i in 1:s
                summed += (v1[i][k] - v2[i][j])^2
            end
            d[k, j] = sqrt(summed)
        end
    end

    return d'
end

data = zeros(Float64, 3, 15000)
for i in 1:5000
    data[1, i] = rand(Normal(100, 10))
    data[2, i] = rand(Normal(40, 5))
    data[3, i] = rand(Normal(40, 10))
end

for i in 5000:10000
    data[1, i] = rand(Normal(20, 4))
    data[2, i] = rand(Normal(50, 6))
    data[3, i] = rand(Normal(50, 10))
end

for i in 10000:15000
    data[1, i] = rand(Normal(50, 5))
    data[2, i] = rand(Normal(90, 10))
    data[3, i] = rand(Normal(70, 8))
end

# display(data)
plotly()
scatter()
@time ( clustered, cluster_counter ) = DBSCAN(data, 4, 6)
for i in -1:cluster_counter
    cluster = [(x.coordinates[1], x.coordinates[2], x.coordinates[3]) for x in clustered if x.label == i]
    scatter!(cluster, markersize=1.5)
end
display(current())

# @time find_neighbors(data, data[:, 1], 10)


text:
import StatsModels

struct FeatureContrastsNotYetGenerated <: AbstractNonExistentFeatureContrasts
end

struct FitNotYetRunUnderlyingObject <: AbstractNonExistentUnderlyingObject
end

struct FitFailedUnderlyingObject <: AbstractNonExistentUnderlyingObject
end

"""
"""
mutable struct GLMModel <: AbstractEstimator
    name::T1 where T1 <: AbstractString
    isclassificationmodel::T2 where T2 <: Bool
    isregressionmodel::T3 where T3 <: Bool

    formula::T4 where T4 <: Any
    family::T5 where T5 <: Any
    link::T6 where T6 <: Any

    # parameters (learned from data):
    underlyingglm::T7 where T7 <: Any
end

"""
"""
mutable struct KnetModel <: AbstractEstimator
    name::T1 where T1 <: AbstractString
    isclassificationmodel::T2 where T2 <: Bool
    isregressionmodel::T3 where T3 <: Bool

    # hyperparameters (not learned from data):
    predict_function_source::T4 where T4 <: AbstractString
    loss_function_source::T5 where T5 <: AbstractString
    predict_function::T6 where T6 <: Any
    loss_function::T7 where T7 <: Any
    losshyperparameters::T8 where T8 <: AbstractDict
    optimizationalgorithm::T9 where T9 <: Symbol
    optimizerhyperparameters::T10 where T10 <: AbstractDict
    minibatchsize::T11 where T11 <: Integer
    maxepochs::T12 where T12 <: Integer
    printlosseverynepochs::T13 where T13 <: Integer

    # parameters (learned from data):
    modelweights::T14 where T14 <: AbstractArray
    modelweightoptimizers::T15 where T15 <: Any

    # learning state
    history::T16 where T16 <: Any
end

"""
"""
mutable struct SimplePipeline{S<:AbstractString, T<:AbstractVector} <: AbstractPipeline
    name::S
    objectsvector::T
end

struct PGFPlotsXPlot <: AbstractPlot
    underlying_object::T where T <: Any
end

"""
"""
struct ImmutablePackageMultiLabelPredictionTransformer <: AbstractEstimator
    label_names::T1 where T1 <: AbstractVector{Symbol}
end

"""
"""
struct ImmutablePackageSingleLabelPredictionTransformer <: AbstractEstimator
    single_label_name::T1 where T1 <: Symbol
end

"""
"""
struct ImmutablePackageSingleLabelPredictProbaTransformer <:
        AbstractEstimator
    single_label_name::T1 where T1 <: Symbol
end

"""
"""
struct ImmutablePredictionsSingleLabelInt2StringTransformer <:
        AbstractEstimator
    index::T1 where T1 <: Integer
    levels::T2 where T2 <: AbstractVector
end

"""
"""
struct ImmutablePredictProbaSingleLabelInt2StringTransformer <:
        AbstractEstimator
    index::T1 where T1 <: Integer
    levels::T2 where T2 <: AbstractVector
end

"""
"""
struct DataFrameFeatureContrasts <: AbstractFeatureContrasts
    columns::T1 where T1 <: AbstractVector{Symbol}
    num_df_columns::T2 where T2 <: Integer
    schema_without_intercept::T3 where T3 <: StatsModels.Schema
    formula_without_intercept::T4 where T4 <: StatsModels.AbstractTerm
    num_array_columns_without_intercept::T5 where T5 <: Integer
    schema_with_intercept::T6 where T6 <: StatsModels.Schema
    formula_with_intercept::T7 where T7 <: StatsModels.AbstractTerm
    num_array_columns_with_intercept::T8 where T8 <: Integer
end

"""
"""
mutable struct MutableDataFrame2DecisionTreeTransformer <: AbstractEstimator
    feature_names::T1 where T1 <: AbstractVector
    single_label_name::T2 where T2 <: Symbol
    levels::T3 where T3 <: AbstractVector
    dffeaturecontrasts::T4 where T4 <: AbstractFeatureContrasts
end

"""
"""
struct ImmutableDataFrame2GLMSingleLabelBinaryClassTransformer <:
        AbstractEstimator
    label::T1 where T1 <: Symbol
    positive_class::T2 where T2 <: AbstractString
end

"""
"""
mutable struct MutableDataFrame2ClassificationKnetTransformer <:
        AbstractEstimator
    feature_names::T1 where T1 <: AbstractVector
    label_names::T2 where T2 <: AbstractVector{Symbol}
    label_levels::T3 where T3 <: AbstractDict
    index::T4 where T4 <: Integer
    transposefeatures::T5 where T5 <: Bool
    transposelabels::T6 where T6 <: Bool
    dffeaturecontrasts::T7 where T7 <: AbstractFeatureContrasts
end

"""
"""
mutable struct MutableDataFrame2RegressionKnetTransformer <:
        AbstractEstimator
    feature_names::T1 where T1 <: AbstractVector
    label_names::T2 where T2 <: AbstractVector{Symbol}
    transposefeatures::T3 where T3 <: Bool
    transposelabels::T4 where T4 <: Bool
    dffeaturecontrasts::T5 where T5 <: AbstractFeatureContrasts
    function MutableDataFrame2RegressionKnetTransformer(
            feature_names::AbstractVector,
            label_names::AbstractVector{Symbol};
            transposefeatures::Bool = true,
            transposelabels::Bool = false,
            )
        result = new(
            feature_names,
            label_names,
            transposefeatures,
            transposelabels,
            )
        return result
    end
end

"""
"""
struct ImmutableFeatureArrayTransposerTransformer <: AbstractEstimator
end

"""
"""
mutable struct LIBSVMModel <: AbstractEstimator
    name::T1 where T1 <: AbstractString
    isclassificationmodel::T2 where T2 <: Bool
    isregressionmodel::T3 where T3 <: Bool

    levels::T4 where T4 <: AbstractVector

    # hyperparameters (not learned from data):
    hyperparameters::T5 where T5 <: AbstractDict

    # parameters (learned from data):
    underlyingsvm::T6 where T6 <: Any
end

"""
"""
mutable struct DecisionTreeModel <:
        AbstractEstimator
    name::T1 where T1 <: AbstractString
    isclassificationmodel::T2 where T2 <: Bool
    isregressionmodel::T3 where T3 <: Bool

    single_label_name::T4 where T4 <: Symbol
    levels::T5 where T5 <: AbstractVector

    # hyperparameters (not learned from data):
    hyperparameters::T6 where T6 <: AbstractDict

    # parameters (learned from data):
    underlyingrandomforest::T7 where T7 <: Any
end

"""
"""
struct CrossValidation{T}
    leavein::Vector{CrossValidation{T}}
    leaveout::Vector{Vector{T}}
end


text:
isdefined(Base, :__precompile__) && __precompile__(false)

module SMC
    using Dates, Distributed, Distributions
    using FileIO, HDF5, JLD2, LinearAlgebra, Random
    using ModelConstructors

    using Roots: fzero, ConvergenceFailed
    using StatsBase: sample, Weights

    import Base.<, Base.isempty, Base.min, Base.max
    import Calculus, ModelConstructors

    export
        mutation, resample, smc,
        mvnormal_mixture_draw, initial_draw!, Cloud, get_cloud,
        isempty, join_cloud, split_cloud, add_parameters_to_cloud

    const VERBOSITY   = Dict(:none => 0, :low => 1, :high => 2)
    const DATE_FORMAT = "yymmdd"

    include("particle.jl")
    include("initialization.jl")
    include("helpers.jl")
    include("util.jl")
    include("mutation.jl")
    include("resample.jl")
    include("smc_main.jl")
end


text:
using Plots, QuasiGeostrophy, Test
include(pwd() * "/test/test_utils.jl")
Ω = S¹(0,2π)
Nx = 2^3;
plot_flag = false
plot_flag2 = false
fourier_grid = FourierGrid(Nx, Ω)
fieldnames = ("ϕ_numerical", "ϕ_exact")
ϕ_numerical, ϕ_exact = create_fields(names = fieldnames, grid = fourier_grid)
∂x = create_operators(fourier_grid)
# initialize fields with nontrivial data
x = fourier_grid.grid[1]
ϕ_exact(sin.(x));
# Forward Euler
Δt_list = []
error_list = []
for scale in [1, 10^1, 10^2, 10^3, 10^4]
    ϕ_numerical(sin.(x))
    local Δt = 0.1 / scale * 2π
    push!(Δt_list, Δt)
    t = randn(1) .* 0
    for i in 1:10*scale
        local f = ∂x(ϕ_numerical)
        nϕ = ϕ_numerical + Δt * f
        ϕ_numerical.data .= nϕ.data
        if i%(scale) == 0
            if plot_flag
                display(plot(ϕ_numerical))
                sleep(0.01)
            end
        end
        t[1] += Δt
    end
    error = norm(ϕ_numerical - ϕ_exact)
    push!(error_list, error)
end
#  Heuns Method
Δt_list_2 = []
error_list_2 = []
for scale in [1, 10^1, 10^2, 10^3, 10^4]
    ϕ_numerical(sin.(x))
    local Δt = 0.1 / scale * 2π
    push!(Δt_list_2, Δt)
    t = randn(1) .* 0
    for i in 1:10*scale
        local f = ∂x(ϕ_numerical)
        nϕ = ϕ_numerical + Δt * f
        nϕ = ϕ_numerical + Δt/2 * ( f +  ∂x(nϕ) )
        ϕ_numerical.data .= nϕ.data
        if i%(scale) == 0
            if plot_flag
                display(plot(ϕ_numerical))
                sleep(0.01)
            end
        end
        t[1] += Δt
    end
    error = norm(ϕ_numerical - ϕ_exact)
    push!(error_list_2, error)
end
#
prefactor = error_list ./ Δt_list # first order convergence
prefactor2 = error_list_2 ./ Δt_list_2 .^2 # second order convergence
if plot_flag2
    xc = collect(0:0.1:4.5)
    plot(xc, -1 .* xc .+ log10(prefactor[end]), label = "1st order theoretical", xlims = (0, 4.5), linestyle = :dash, linewidth = 2, color = :orange)
    scatter!(-log10.(Δt_list), log10.(error_list), label = "1st order numerical", color = :orange)
    p1 = plot!(xlabel = "-log10(Δt)", ylabel = "log10(error)")
    plot!(xc, -2 .* xc .+ log10(prefactor2[end]), label = "2nd order theoretical", linestyle = :dash, linewidth = 2, color = :green)
    scatter!(-log10.(Δt_list_2), log10.(error_list_2), label = "2nd order numerical", color = :green)
    p2 = plot!(xlabel = "-log10(Δt)", ylabel = "log10(error)")
    display(p2)
end

@testset "TimeStepping Convergence Test" begin
    # forward euler prefactor for this test is about 40
    @test prod(prefactor[2:end] .< 40)
    # Heun's method prefactor for this test is about 12
    @test prod(prefactor2[2:end] .< 12)
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
using Docile
module JuDaisy

  export Daisy, verbose!, image!, parameters!, compute!, normalize!, getindex, descriptor_size

  push!(LOAD_PATH, dirname(Base.source_path()))
  const daisyC = "build/libdaisyC.so"

  """Maps the dummy C daisy struct"""
  immutable Daisy_struct
  end

  """Interface type to the descriptor"""
  type Daisy
    ptr::Ptr{Daisy_struct}
    descriptor_size::Cint
    function Daisy()
      ptr = ccall( (:new_daisy, daisyC), Ptr{Daisy_struct}, () )
      desc = new(ptr)
      delete_daisy!(desc::Daisy) = ccall( (:delete_daisy, daisyC), Void, (Ptr{Daisy_struct},), desc.ptr)
      finalizer(desc, delete_daisy!)
      desc
    end
  end

  """Set verbose level v=1,2,3,4"""
  function verbose!(desc::Daisy, v)
    ccall( (:daisy_verbose, daisyC), Void, (Ptr{Daisy_struct}, Csize_t), desc.ptr, v)
  end

  """Set image. im can be of types that can be casted to float.
  Expects values between 0. and 1."""
  function image!(desc::Daisy, im::Array{Float32, 2})
    ccall( (:daisy_set_float_image, daisyC), Void, (Ptr{Daisy_struct}, Ptr{Float32}, Cint, Cint), desc.ptr, im, size(im, 2), size(im, 1))
  end

  """Set parameters and initialize for Mode 1"""
  function parameters!(desc::Daisy, rad::Float64 = 15., radq::Int = 3, thq::Int = 8, histq::Int = 8)
    ccall( (:daisy_set_parameters, daisyC), Void, (Ptr{Daisy_struct}, Cdouble, Cint, Cint, Cint), desc.ptr, rad, radq, thq, histq)
    desc.descriptor_size = descriptor_size(desc)
    ccall( (:daisy_initialize_single_descriptor_mode, daisyC), Void, (Ptr{Daisy_struct},), desc.ptr)
  end

  """Compute all descriptors"""
  function compute!(desc::Daisy)
    ccall( (:daisy_compute_descriptors, daisyC), Void, (Ptr{Daisy_struct},), desc.ptr)
  end

  function normalize!(desc::Daisy)
    ccall( (:daisy_normalize_descriptors, daisyC), Void, (Ptr{Daisy_struct},), desc.ptr)
  end

  function descriptor_size(desc::Daisy)
    ccall( (:daisy_descriptor_size, daisyC), Cint, (Ptr{Daisy_struct},), desc.ptr)
  end

  import Base.getindex
  """To meet Julia convention x and y start at 1
  no bounds checking !"""
  function getindex(desc::Daisy, y::Int, x::Int)
    thor::Ptr{Cfloat} = ccall( (:daisy_get_descriptor, daisyC), Ptr{Cfloat}, (Ptr{Daisy_struct}, Cint, Cint), desc.ptr, y-1, x-1)
    pointer_to_array(thor, (desc.descriptor_size,), false)
  end

end


text:
"""
    Adam(::Vector{XVariable}; lr=1e-3, b1=0.9, b2=0.999, eps=1e-8, L1decay=0.001, L2decay=0.01)

Implements Adam algorithm. Refer to `Adam: A Method for Stochastic Optimization`.
"""
mutable struct Adam <: Optimizer
    xparams::Vector{XVariable}
    w1::Vector
    w2::Vector
    lr::AbstractFloat
    b1::AbstractFloat
    b2::AbstractFloat
    ϵ::AbstractFloat
    t::UInt
    b1t::AbstractFloat
    b2t::AbstractFloat
    L1decay::AbstractFloat
    L2decay::AbstractFloat
    name::String
    function Adam(xparams::Vector{XVariable}; lr=1e-3, b1=0.9, b2=0.999, eps=1e-8, L1decay=0.001, L2decay=0.01)
        num = length(xparams)
        w1  = Vector(undef,num)
        w2  = Vector(undef,num)
        for i = 1:num
            c , θ = xparams[i]
            w1[i] = Zeros(typeof(ᵛ(θ)), θ.shape)
            w2[i] = Zeros(typeof(ᵛ(θ)), θ.shape)
        end
        new(xparams,w1,w2,lr, b1, b2, eps, 0, 1.0, 1.0, L1decay, L2decay, "Adam")
    end
end

function Base.show(io::IO, O::Adam)
    print("Adam(lr=$(O.lr), β₁=$(O.b1), β₂=$(O.b2), ϵ=$(O.ϵ), L1decay=$(O.L1decay), L2decay=$(O.L2decay))");
end

function update!(O::Adam; clipfn::Function=LPInfNormClip, clipvalue=10.0)
    w₁ = O.w1
    w₂ = O.w2
    lr = O.lr
    b₁ = O.b1
    b₂ = O.b2
    ϵ  = O.ϵ
    λ₁ = O.L1decay
    λ₂ = O.L2decay
    O.t   += 1
    O.b1t *= b₁
    O.b2t *= b₂
    b₁ᵗ = O.b1t
    b₂ᵗ = O.b2t

    for i = 1:length(O.xparams)
        c , θ = O.xparams[i]
        μ = - sqrt(1-b₂ᵗ) / (1-b₁ᵗ) * lr
        ∇ = clipfn(setNanInfZero(δ(θ)), clipvalue)
        𝒗 = ᵛ(θ)
        @. w₁[i] = b₁ * w₁[i] + (1-b₁) * ∇
        @. w₂[i] = b₂ * w₂[i] + (1-b₂) * ∇ * ∇

        if c == 'w'
            if λ₁==0 && λ₂==0
                @. 𝒗 += μ * w₁[i] / sqrt(w₂[i] + ϵ)
            elseif λ₁==0 && λ₂!=0
                @. 𝒗 += μ * (w₁[i] / sqrt(w₂[i] + ϵ) + λ₂ * 𝒗)
            elseif λ₁!=0 && λ₂==0
                @. 𝒗 += μ * (w₁[i] / sqrt(w₂[i] + ϵ) + λ₁ * sign(𝒗))
            else  # λ₁!=0 && λ₂!=0
                @. 𝒗 += μ * (w₁[i] / sqrt(w₂[i] + ϵ) + λ₁ * sign(𝒗) + λ₂ * 𝒗)
            end
        else
            @. 𝒗 += μ * w₁[i] / sqrt(w₂[i] + ϵ)
        end
    end
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible ideal GLM-MHD equations

equations = IdealGlmMhdEquations3D(5/3)

initial_condition = initial_condition_convergence_test

volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
solver = DGSEM(polydeg=3, surface_flux=(flux_hll, flux_nonconservative_powell),
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = (-1.0, -1.0, -1.0)
coordinates_max = ( 1.0,  1.0,  1.0)

# Create P4estMesh with 2 x 2 x 2 trees
trees_per_dimension = (2, 2, 2)
mesh = P4estMesh(trees_per_dimension,
                 polydeg=3, initial_refinement_level=2,
                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
                 periodicity=true)

# OBS! Workaround to add a refinement patch after mesh is constructed
# Refine bottom left quadrant of each tree to level 4
function refine_fn(p8est, which_tree, quadrant)
  if quadrant.x == 0 && quadrant.y == 0 && quadrant.z == 0 && quadrant.level < 4
    # return true (refine)
    return Cint(1)
  else
    # return false (don't refine)
    return Cint(0)
  end
end

# Refine recursively until each bottom left quadrant of a tree has level 4
# The mesh will be rebalanced before the simulation starts
refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p8est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p8est_quadrant_t}))
Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

cfl = 1.0
stepsize_callback = StepsizeCallback(cfl=cfl)

glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback,
                        glm_speed_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
# Computes graph Stein discrepancy based on the Riemannian-Langevin operator.
# Our parametrization is of the form
#
# dXt = b(Xt) dt + sigma(Xt) dBt
#
# where Xt is the diffusion process and Bt is d-dimensional Brownian Motion.
# Since this is a Riemannian-Langevin diffusion, b and sigma are linked via
# the relationship
#
# a(Xt) = (1/2) sigma(Xt) sigma(Xt)^T,
# grad a(Xt) = <grad, a^T> (e.g. the grad operator applied to each column of M^T),
# b(Xt) = a(Xt) grad log p(Xt) + grad a(Xt).
#
# This ensures that the diffusion has p as its invariant distribution.
#
# In this case, the infinitesmal operator has the form
#
# T g = 2 <g, b> + 2 <a, grad g>.
#
# Args:
#  points - n x d array of sample points
#  weights - n x 1 array of real-valued weights associated with sample points
#    (default: equal weights)
#  gradlogdensity - the gradlogdensity of the target distribution
#  volatility_covariance - the matrix-valued function a = (1/2) sigma sigma^T
#  grad_volatility_covariance - the vector-valued function grad a = <grad, a^t>
#  solver - optimization program solver supported by JuMP
function riemannian_langevin_graph_discrepancy(; points=[],
                                               weights=fill(1/size(points,1), size(points,1)),
                                               gradlogdensity=nothing,
                                               volatility_covariance=nothing,
                                               grad_volatility_covariance=nothing,
                                               solver=nothing,
                                               kwargs...)
    # make sure solver is defined
    isa(solver, AbstractMathProgSolver) ||
        error("Must specify solver of type AbstractMathProgSolver")
    ## Extract inputs
    sample = SteinDiscrete(points, weights)
    points = sample.support
    weights = sample.weights
    n = length(weights)
    d = size(points,2)
    # make sure the volatility matrices are non-nothing
    isa(volatility_covariance, Function) ||
        error("For riemannian langevin methods, you must supply the volatility_covariance")

    isa(grad_volatility_covariance, Function) ||
        error("For riemannian langevin methods, you must supply the grad_volatility_covariance")

    # Now prepare the coefficients for g and grad g
    gobjcoefficients = Array{Float64}(undef, n, d)
    gradgobjcoefficients = Array{Float64}(undef, n, d, d)
    for i in 1:n
        xi = points[i,:]
        wi = weights[i]
        gradlogpxi = gradlogdensity(points[i,:])
        a = volatility_covariance(xi)
        grad_a = grad_volatility_covariance(xi)
        # Objective coefficients for each g(x_i)
        gobjcoefficients[i,:] = 2 * wi * (a * gradlogpxi + grad_a)
        # Objective coefficients for each grad g(x_i)
        gradgobjcoefficients[i,:,:] .= 2 * wi * a
    end
    # solve the optimization problem
    result = affine_graph_discrepancy(sample,
                                      gobjcoefficients,
                                      gradgobjcoefficients,
                                      solver,
                                      "riemannian-langevin";
                                      kwargs...)
    result
end


text:
using Photometry.Aperture: edges,
                           bbox,
                           oblique_coefficients

@testset "Apertures" begin

    ap_ellipse = EllipticalAperture(0, 0, 20, 10, 0)
    @test bbox(ap_ellipse) == (-20, 20, -10, 10)

    @test EllipticalAperture([0, 0], 20, 10, 0) == ap_ellipse

    ap_ellipse = EllipticalAperture(0, 0, 2, 1, 45)
    @test bbox(ap_ellipse) == (-2, 2, -2, 2)
end

@testset "Elliptical Aperture" begin

    e = EllipticalAperture(0, 0, 20, 10, 0)
    @test sprint(show, e) == "EllipticalAperture(0, 0, a=20, b=10, θ=0°)"

    e = EllipticalAperture(0, 0, 10, 10, 0)
    @test mask(e, method = :center) == mask(e, method = (:subpixel, 1))

    # test modding of angle
    @test EllipticalAperture(0, 0, 3, 4, 380).theta == 20

    @test_throws ErrorException EllipticalAperture(0, 0, -2, 4, 0)
    @test_throws ErrorException EllipticalAperture(0, 0, 2, -3, 0)

end

@testset "oblique_coefficients" begin
    @test all(oblique_coefficients(2, 2, 0) .≈ (0.25, 0.25, 0.0))
    @test all(oblique_coefficients(2, 2, 90) .≈ (0.25, 0.25, 0.0))
    @test all(oblique_coefficients(2, 1, 30) .≈ (7 / 16, 13 / 16, -6sqrt(3) / 16))
end

@testset "Elliptical Annulus" begin

    e0 = EllipticalAnnulus(0, 0, 8, 16, 4, 45)
    @test sprint(show, e0) == "EllipticalAnnulus(0.0, 0.0, a_in=8.0, a_out=16.0, b_in=2.0, b_out=4.0, θ=45.0°)"

    e = EllipticalAnnulus(0, 0, 8, 16, 4, 45)
    @test mask(e, method = :center) == mask(e, method = (:subpixel, 1))

    @test EllipticalAnnulus(0, 0, 2, 4, 3, 380).theta == 20

    @test EllipticalAnnulus([0, 0], 8, 16, 4, 45) == e

    @test_throws ErrorException EllipticalAnnulus(0, 0, -2, 4, 0, 2)
    @test_throws ErrorException EllipticalAnnulus(0, 0, 2, -3, 0, 2)
    @test_throws ErrorException EllipticalAnnulus(0, 0, 5, 4, 0, 0.5)

end

@testset "Elliptical Annulus bounding box" begin

    e = EllipticalAnnulus(0, 0, 8, 16, 4, 0)
    @test all(bbox(e) == (-16, 16, -4, 4))

end


text:
"""
linear regression function
  Compute coeff. estimates, s.e's, equation σ, Rsquared

usage:

    b, seb, s, R2 = linreg(y,x)
    y = dependent variable vector
    x = matrix of independent variables (no intercept in x)
"""
function blinreg(y,x)
  # add intercept
  n = length(y)
  X = [ones(n) x]
  b = (X'*X) \ X'*y
  resids = y - X*b
  RSS = sum(resids.^2)
  s2 = RSS/n

  covb = s2.*inv(X'X)
  seb = sqrt.(diag(covb))
  k = length(b)
  ### Correct formulas for AIC and BIC
  BIC = n*log(s2) + k*log(n)
  AIC = 2*k + n*log(s2)
  odds = zeros(k)
  pvals = zeros(k)
  that = abs.(b)./seb
  for i in 1:k
    odds[i] = (1.0 + (that[i]^2)/(n-k))^((n-k+1)/2.0)
    pvals[i] = 1.0 - cdf(TDist(n-k), that[i])
  end
  println("  coeffs = ", round.(b, digits=3))
  println("   s.e's = ", round.(seb, digits=3))
  println("    odds = ", round.(odds, digits=3))
  println("p-values = ", round.(pvals, digits=4))
  println("s^2 (eqn. variance) = ",round(s2, digits=6))

  # compute R^2
  tss = sum((y .- mean(y)).^2)

  R2 = 1 - s2*n/tss
  println("Rsquared = ",round(R2, digits=3))
  println("AIC = ", round(AIC, digits = 2), "  BIC = ", round(BIC, digits = 2),)
  s = sqrt(s2)
  return b, seb, odds, pvals, s, R2, RSS, AIC, BIC
end

## example of use:
#x = randn(20)
#y = 1 + 1.*x .+ randn(20)
#b, seb, odds, pvals, s, R2, RSS, AIC, BIC = blinreg(y,x)


text:
@testset "Helper" begin
    @testset "read_varint" begin
        tests = [([0x01], 1),
                 ([0xfd, 0xd0, 0x24], 9424),
                 ([0xfe, 0x30, 0x33, 0xff, 0xb3], -1275120848),
                 ([0xff, 0x70, 0x9a, 0xeb, 0xb4, 0xbb, 0x7f, 0x00, 0x00], 140444170951280)]
        for t in tests
            @test Bitcoin.read_varint(IOBuffer(t[1])) == t[2]
        end
    end
    @testset "encode_varint" begin
        want = [0xfd, 0xfe, 0x00]
        @test Bitcoin.encode_varint(254) == want
        want = [0xfe, 0xff, 0xff, 0xff, 0xff]
        @test Bitcoin.encode_varint(0x100000000-1) == want
        want = [0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff]
        @test Bitcoin.encode_varint(0x10000000000000000-1) == want
        @test_throws UndefVarError encode_varint(0x10000000000000000)
    end
end


text:
"""
Definition of `DynamicalSystem` as well as all integrators
used in the ecosystem of DynamicalSystems.jl.

Also contains predefined well-known systems.
"""
module DynamicalSystemsBase

using Reexport
@reexport using DelayEmbeddings

include("dynamicalsystem.jl")
include("discrete.jl")
include("continuous.jl")
include("famous_systems.jl")

include("deprecations.jl")

export Systems, reinit!
export SVector, SMatrix, @SVector, @SMatrix

end


text:
# short be
using SortingLab, ShortStrings, SortingAlgorithms, BenchmarkTools;
N = Int(1e6);
svec = [randstring(rand(1:15)) for i=1:N];
# convert to ShortString
ssvec = ShortString.(svec);
basesort = @benchmark sort($svec)
radixsort_timings = @benchmark SortingLab.radixsort($svec)
short_radixsort = @benchmark sort($ssvec, by = x->x.size_content, alg=RadixSort)

using RCall

@rput svec;
r_timings = R"""
replicate($(length(short_radixsort.times)), system.time(sort(svec, method="radix"))[3])
""";

using Plots
bar(["Base.sort","SortingLab.radixsort","ShortStrings radix sort", "R radix sort"],
    mean.([basesort.times./1e9, radixsort_timings.times./1e9, short_radixsort.times./1e9, r_timings]),
    title="String sort performance - len: 1m, variable size 15",
    label = "seconds")
savefig("readme_string_sort.png")

# short be
using SortingLab, ShortStrings, SortingAlgorithms, BenchmarkTools;
N = Int(2e7);
svec = rand([randstring(rand(1:15)) for i=1:N÷100],N)
# convert to ShortString
ssvec = ShortString.(svec);
basesort = @benchmark sort($svec) samples = 5 seconds = 120
radixsort_timings = @benchmark SortingLab.radixsort($svec) samples = 5 seconds = 120
short_radixsort = @benchmark sort($ssvec, by = x->x.size_content, alg=RadixSort) samples = 5 seconds = 120

using RCall

@rput svec;
r_timings = R"""
replicate(max(5, $(length(short_radixsort.times))), system.time(sort(svec, method="radix"))[3])
""";

using Plots
bar(["Base.sort","SortingLab.radixsort","ShortStrings radix sort", "R radix sort"],
    mean.([basesort.times./1e9, radixsort_timings.times./1e9, short_radixsort.times./1e9, r_timings]),
    title="String sort performance - len: $(N÷1_000_000)m, fixed size: 15",
    label = "seconds")
savefig("readme_string_sort_fixed_len.png")

N = 100_000_000
@time a = ShorterString.()

fsort(s) = sort(s, by = x->x.size_content, alg=RadixSort);
@time a_sorted = fsort(a);

ax = (x->x.size_content).(a)
@time sort(ax, alg=RadixSort);


text:
#=
Sparse vector
=#

using SimpleDataStructures

ssv = SimpleSparseVector(
    100,
    [1, 10, 50],
    rand(3),
)

ssv[1]
ssv[10]

ssv[2]

ssv[101]

ssv[3] = rand()

insert!(
    ssv,
    10,
    rand(),
)

deleteat!(ssv, 20)
deleteat!(ssv, 3)

ssv = SimpleSparseVector{Float64}(100)
ssv[50] = rand()


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "UtilLinuxBuilder"
version = v"2.33.0"

# Collection of sources required to build UtilLinuxBuilder
sources = [
    "https://mirrors.edge.kernel.org/pub/linux/utils/util-linux/v2.33/util-linux-2.33.tar.xz" =>
    "f261b9d73c35bfeeea04d26941ac47ee1df937bd3b0583e748217c1ea423658a",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd util-linux-2.33/
./configure --prefix=$prefix --host=$target
make
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libuuid", :libuuid),
    LibraryProduct(prefix, "libmount", :libmount),
    LibraryProduct(prefix, "libblkid", :libblkid),
    LibraryProduct(prefix, "libsmartcols", :libsmartcols),
    LibraryProduct(prefix, "libfdisk", Symbol("lobfdisk\e[D\e[D\e[D\e[D\e[D\e[i\e[4~\e[C\e[C\e[C\e[C\e[C\e[C\e[C\e[C\e[D\e[D"))
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: hessammehr/UtilLinuxBuilder

max_issues_repo_name: hessammehr/UtilLinuxBuilder

max_forks_repo_name: hessammehr/UtilLinuxBuilder


text:
include("aerosol_activation_tests.jl")
include("microphysics_tests.jl")
include("gpu_tests.jl")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Testing

@testset "read inp file" begin
    model = open(parse_abaqus, Pkg.dir("JuliaFEM")*"/geometry/3d_beam/palkki.inp")
    @test length(model["nodes"]) == 298
    @test length(model["elements"]) == 120
    @test length(model["elsets"]["Body1"]) == 120
    @test length(model["nsets"]["SUPPORT"]) == 9
    @test length(model["nsets"]["LOAD"]) == 9
    @test length(model["nsets"]["TOP"]) == 83
end

@testset "test read element section" begin
    data = """*ELEMENT, TYPE=C3D10, ELSET=BEAM
    1,       243,       240,       191,       117,       245,       242,       244,
    1,         2,       196
    2,       204,       199,       175,       130,       207,       208,       209,
    3,         4,       176
    """
    data = split(data, "\n")
    model = Dict{AbstractString, Any}()
    model["nsets"] = Dict{AbstractString, Vector{Int}}()
    model["elsets"] = Dict{AbstractString, Vector{Int}}()
    model["elements"] = Dict{Integer, Any}()
    parse_section(model, data, :ELEMENT, 1, 5, Val{:ELEMENT})
    @test length(model["elements"]) == 2
    @test model["elements"][1]["connectivity"] == [243, 240, 191, 117, 245, 242, 244, 1, 2, 196]
    @test model["elements"][2]["connectivity"] == [204, 199, 175, 130, 207, 208, 209, 3, 4, 176]
    @test model["elsets"]["BEAM"] == [1, 2]
end

@testset "parse nodes from abaqus .inp file to Mesh" begin
    fn = Pkg.dir("JuliaFEM") * "/test/testdata/cube_tet4.inp"
    mesh = abaqus_read_mesh(fn)
    info(mesh.surfaces)
    @test length(mesh.nodes) == 10
    @test length(mesh.elements) == 17
    @test haskey(mesh.elements, 1)
    @test mesh.elements[1] == [8, 10, 1, 2]
    @test mesh.element_types[1] == :Tet4
    @test haskey(mesh.node_sets, :SYM12)
    @test haskey(mesh.element_sets, :CUBE)
    @test haskey(mesh.surfaces, :LOAD)
    @test length(mesh.surfaces[:LOAD]) == 2
    @test mesh.surfaces[:LOAD][1] == (16, :S1)
    @test mesh.surface_types[:LOAD] == :ELEMENT
end


text:
module AbstractPlotting

using Observables, GeometryTypes, StaticArrays, ColorTypes, Colors, IntervalSets
using ColorBrewer
using Markdown # documentation
using Serialization # serialize events

# Text related packages
using Packing
using SignedDistanceFields
using FreeType, FreeTypeAbstraction, UnicodeFun
using LinearAlgebra, Statistics
import FileIO

using Base: RefValue
using Base.Iterators: repeated, drop
import Base: getindex, setindex!, push!, append!, parent, get, get!, delete!, haskey

include("utilities/quaternions.jl")
include("types.jl")
include("utilities/utilities.jl")
include("utilities/logging.jl")
include("utilities/texture_atlas.jl")
include("interaction/nodes.jl")

# Basic scene/plot/recipe interfaces + types
include("scenes.jl")
include("theming.jl")
include("recipes.jl")
include("interfaces.jl")
include("conversions.jl")

# camera types + functions
include("camera/projection_math.jl")
include("camera/camera.jl")
include("camera/camera2d.jl")
include("camera/camera3d.jl")

# some default recipes
include("basic_recipes/basic_recipes.jl")
# layouting of plots
include("layouting/transformation.jl")
include("layouting/data_limits.jl")
include("layouting/layouting.jl")
include("layouting/boundingbox.jl")
# more default recipes
include("basic_recipes/buffers.jl")
include("basic_recipes/axis.jl")
include("basic_recipes/legend.jl")

#
include("interaction/events.jl")
include("interaction/gui.jl")

# documentation and help functions
include("documentation.jl")
include("display.jl")

# help functions and supporting functions
export help, help_attributes, help_arguments

# Abstract/Concrete scene + plot types
export AbstractScene, SceneLike, Scene, AbstractScreen
export AbstractPlot, Combined, Atomic, Axis

# Theming, working with Plots
export Attributes, Theme, attributes, default_theme, theme, set_theme!

# Node/Signal related
export Node, node, lift, map_once, to_value, on

# utilities and macros
export @recipe, @extract, @extractvalue, @key_str, @get_attribute
export broadcast_foreach, to_vector, replace_automatic!

# conversion infrastructure
export @key_str, convert_attribute, convert_arguments
export to_color, to_colormap, to_rotation, to_font, to_align, to_textsize
export to_ndim, Reverse

# Transformations
export translated, translate!, transform!, scale!, rotate!, grid, Accum, Absolute
export boundingbox, insertplots!, center!, translation
export hbox, vbox

# camera related
export AbstractCamera, EmptyCamera, Camera, Camera2D, Camera3D, cam2d!, cam2d
export campixel!, campixel, cam3d!, update_cam!, rotate_cam!, translate_cam!
export pixelarea, plots, cameracontrols, cameracontrols!, camera, events
export to_world
# picking + interactive use cases + events
export mouseover, ispressed, onpick, pick, Events, Keyboard, Mouse
export register_callbacks
export window_area
export window_open
export mouse_buttons
export mouse_position
export mousedrag
export scroll
export keyboard_buttons
export unicode_input
export dropped_files
export hasfocus
export entered_window
export disconnect!, must_update, force_update!

# gui
export slider, button, playbutton

# Raymarching algorithms
export RaymarchAlgorithm, IsoValue, Absorption, MaximumIntensityProjection, AbsorptionRGBA, IndexedAbsorptionRGBA
export Billboard

# Reexports of
# Color/Vector types convenient for 3d/2d graphics
export Quaternion, Quaternionf0, qrotation
export RGBAf0, RGBf0, VecTypes, RealVector, FRect, FRect2D, IRect2D
export FRect3D, IRect3D, Rect3D, Transformation
export IRect, FRect, Rect, Sphere, Circle
export Vec4f0, Vec3f0, Vec2f0, Point4f0, Point3f0, Point2f0
export Vec, Vec2, Vec3, Vec4, Point, Point2, Point3, Point4
#export (..) # reexport interval

export plot!, plot

end # module


text:
module SfGravity

using PhysicalConstants.CODATA2018: c_0, g_n, G, StefanBoltzmannConstant, ħ, k_B
using Unitful, UnitfulAstro, Documenter

export vis_viva, gravity, planetary_mass, planetary_radius, escape_velocity, hill_sphere,
	orbital_velocity, orbital_period, gravitational_binding_energy, roche_limit, gravity_tug_mass,
	barycentric_distance, tidal_acceleration,
	radial_orbit_time, radial_orbit_displacement

DocMeta.setdocmeta!(SfGravity, :DocTestSetup, :(using Unitful, UnitfulAstro, ..SfGravity); recursive=true)

"""
    vis_viva(parent_mass::Unitful.Mass, semimajor_axis::Unitful.Length, current_radius::Unitful.Length)

Return the square of the relative velocities of a body with negligible mass and orbit with `semimajor_axis` and `current_radius` from a central body of `parent_mass`.

``v^2 = GM\\left( \\frac{2}{r} - \\frac{1}{a} \\right)``

where ``G`` is the gravitational constant, ``M`` is the mass of the central body, ``r`` is the current separation of the bodies and ``a`` is the semi-major axis of the orbit.
"""
function vis_viva(parent_mass::Unitful.Mass, semimajor_axis::Unitful.Length, current_radius::Unitful.Length)
	G*parent_mass*(2/current_radius - 1/semimajor_axis) |> u"m^2/s^2"
end

"""
	gravity(m::Unitful.Mass, r::Unitful.Length)

Calculate acceleration due to gravity at distance `r` from point mass `m`

# Example
```jldoctest
julia> gravity(1u"Mearth", 1u"Rearth")
9.798398133669465 m s^-2
```"""
gravity(m::Unitful.Mass, r::Unitful.Length) = G * m / r^2 |> u"m/s/s"

"""
	gravity(m1::Unitful.Mass, m2::Unitful.Mass, r::Unitful.Length)

Calculate the strength of gravity at a distance `r` between two masses, `m1` and `m2`.

# Example
```jldoctest
julia> gravity(1u"Mearth", 100u"kg", 1u"Rearth")
979.8398133669466 N
```"""
gravity(m1::Unitful.Mass, m2::Unitful.Mass, r::Unitful.Length) = G * m1 * m2 / r^2 |> u"N"

"""
	planetary_mass(g::Unitful.Acceleration, r::Unitful.Length)

Calculate the mass of a body required to produce gravitational acceleration `g` at a distance of `r`.
"""
planetary_mass(g::Unitful.Acceleration, r::Unitful.Length) = (g * r^2) / G |> u"kg"

"""
	planetary_mass(ve::Unitful.Velocity, r::Unitful.Length)

Calculate the mass of a body required to produce escape velocity `ve` at a distance of `r`.
"""
planetary_mass(ve::Unitful.Velocity, r::Unitful.Length) = r * ve^2 / 2G |> u"kg"

"""
	planetary_radius(m::Unitful.Mass, g::Unitful.Acceleration)

Radius of a spherical planet of mass `m` with surface gravitational acceleration of `g`.
"""
planetary_radius(m::Unitful.Mass, g::Unitful.Acceleration) = sqrt((G * m) / g) |> u"m"

"""
	planetary_radius(m::Unitful.Mass, g::Unitful.Acceleration)

Radius of a spherical planet of mass `m` with surface escape velocity `ve`.
"""
planetary_radius(m::Unitful.Mass, ve::Unitful.Velocity) = 2G * m / ve^2 |> u"m"

"""
	planetary_radius(ρ::Unitful.Density, g::Unitful.Acceleration)

Radius of a spherical planet of density `ρ` with surface gravitational acceleration of `g`.
"""
planetary_radius(ρ::Unitful.Density, g::Unitful.Acceleration) = 3g / (4G * π * ρ) |> u"m"

"""
	planetary_radius(ρ::Unitful.Density, ve::Unitful.Velocity)

Radius of a spherical planet of density `ρ` with surface escape velocity `ve`.
"""
planetary_radius(ρ::Unitful.Density, ve::Unitful.Velocity) = sqrt(3ve^2 / (8G * π * ρ)) |> u"m"

# kepler's third

"""
	orbital_period(m::Unitful.Mass, r::Unitful.Length)

Calculate the period of an orbit with semimajor axis `r` about a body with mass `m`.

``P = \\sqrt{\\frac{4π^2r^3}{Gm}}``
"""
orbital_period(m::Unitful.Mass, r::Unitful.Length) = sqrt((4π^2 * r^3) / (G * m)) |> u"s"

"""
	orbital_radius(m::Unitful.Mass, t::Unitful.Time)

Calculate the semimajor axis of an orbit with period `t` about a body with mass `m`.

``a = \\sqrt[3]{\\frac{t^2Gm}{4π^2}}``
"""
orbital_radius(m::Unitful.Mass, t::Unitful.Time) = cbrt((t^2 * G * m)/(4π^2)) |> u"m"

"""
    orbital_radius(m::Unitful.Mass, v::Unitful.Velocity)

Approximate the radius of a circular orbit with velocity `v` about a body with mass `m`.
"""
orbital_radius(m::Unitful.Mass, v::Unitful.Velocity) = (G * m) / v^2 |> u"m"

"""
	orbital_velocity(sma::Unitful.Length, t::Unitful.Time)

Calculate the average orbital velocity of a circular orbit of period `t`.

The orbit is assumed to be circular.
"""
orbital_velocity(sma::Unitful.Length, t::Unitful.Time) = 2π * sma / t |> u"km/s"

"""
    orbital_velocity(sma::Unitful.Length, t::Unitful.Time, e)

Approximate the average orbital velocity of a orbit with eccentricity `e` and period `t`.

``s_o = \\frac{2πa}{t} \\left( 1 - \\frac{e^2}{4} - \\frac{3e^4}{64} - \\frac{5e^6}{256} - \\frac{175e^8}{16384} \\right)``
"""
function orbital_velocity(sma::Unitful.Length, t::Unitful.Time, e)
	(2π * sma / t) * (1 - e^2 / 4 - 3e^4 / 64 - 5e^6 / 256 - 175e^8 / 16384) |> u"km/s"
end

"""
    orbital_velocity(parent_mass::Unitful.Mass, semimajor_axis::Unitful.Length, current_radius::Unitful.Length)

Approximate the orbital velocity of a body with an eccentric orbit with `semimajor_axis` and current orbital radius `current_radius` about a body of mass `parent_mass`
"""
orbital_velocity(parent_mass::Unitful.Mass, semimajor_axis::Unitful.Length, current_radius::Unitful.Length) = sqrt(vis_viva(parent_mass, sem, current_radius)) |> u"km/s"

"""
    orbital_velocity(parent_mass::Unitful.Mass, radius::Unitful.Length)

Approximate the orbital velocity of a body with a circular orbit of radius `radius` about a body with mass `parent_mass`.
"""
orbital_velocity(parent_mass::Unitful.Mass, radius::Unitful.Length) = sqrt(G * parent_mass / radius) |> u"km/s"

"""
	planetary_mass(orbital_radius::Unitful.Length, orbital_period::Unitful.Time)

Calculate the mass of a body orbited by a satellite with semimajor axis `sma` and period `orbital_period`.

``M_p = \\frac{4π^2a^3}{GT^2}``

where ``a`` is the semi-major axis and ``T`` is the orbital period.
"""
planetary_mass(sma::Unitful.Length, orbital_period::Unitful.Time) = (4π^2 * sma^3) / (G * orbital_period^2) |> u"kg"

"""
	escape_velocity(m::Unitful.Mass, r::Unitful.Length)

Calculate the velocity required to escape a body with mass `m` from distance `r`.
"""
escape_velocity(m::Unitful.Mass, r::Unitful.Length) = sqrt(2G * m / r) |> u"km/s"

"""
	hill_sphere(m_parent::Unitful.Mass, m::Unitful.Mass, sma::Unitful.Length, e = 0)

Calculate the approximate Hill sphere radius at periapse of a body of mass `m` that orbits a body of mass `m_parent` at a distance of `sma` with orbital eccentricity `e`.

``h_r = a(1-e)\\sqrt[3]{\\frac{m}{3M}}``

where ``a`` is the semi-major axis, `m` is the mass of the smaller body and ``M`` is the mass of the central body.
"""
function hill_sphere(m_parent::Unitful.Mass, m::Unitful.Mass, sma::Unitful.Length, e = 0)
	sma * (1-e) * cbrt(m / 3m_parent) |> Unitful.unit(sma)
end

"""
	gravitational_binding_energy(m::Unitful.Mass, r::Unitful.Length)

Approximate the gravitational binding energy of a body with mass `m` and radius `r`.

``E_g = \\frac{3Gm^2}{5r}``
"""
gravitational_binding_energy(m::Unitful.Mass, r::Unitful.Length) = (3G*m^2) / 5r |>u"J"

"""
    roche_limit(r_primary::Unitful.Length, ρ_primary::Unitful.Density, ρ_satellite::Unitful.Density)

Rigid body approximation of the Roche limit for a body with radius `r_primary` and density `ρ_primary`, approached by a body with density `ρ_satellite`.

``r_r = r\\sqrt[2]{\\frac{2ρ_p}{ρ_s}}``

where ``r`` and ``ρ_p`` are the radius and density of the central body, and ρ_s is the density of the orbiting body.
"""
roche_limit(r_primary::Unitful.Length, ρ_primary::Unitful.Density, ρ_satellite::Unitful.Density) = r_primary * cbrt(2ρ_primary / ρ_satellite) |> u"km"

"""
    gravity_tug_mass(force::Unitful.Force, separation::Unitful.Length, body_mass::Unitful.Mass)

Compute the mass required to apply `force` to a body of mass `body_mass` given their barycenters are at a distance of `separated`.
"""
gravity_tug_mass(force::Unitful.Force, separation::Unitful.Length, body_mass::Unitful.Mass) = (force * separation^2)/(G * body_mass) |> u"kg"

"""
    barycentric_distance(m_1::Unitful.Mass, m_2::Unitful.Mass, a::Unitful.Length)

Compute the distance from the barycenter of `m_1` to the barycenter of the `m_1`-`m_2` system, where the masses have an average separation of `a`.
"""
barycentric_distance(m1::Unitful.Mass, m2::Unitful.Mass, a::Unitful.Length) = (a * m2) / (m1 + m2) |> u"km"

"""
    tidal_acceleration(body_radius::Unitful.Length, parent_mass::Unitful.Mass, separation::Unitful.Length)

Approximate tidal acceleration felt by a body of radius `body_radius` at a distance of `separation` from a body of mass `parent_mass`.

Changes in tidal acceleration are associated with effects like heating, etc.

``a_t = \frac{2GMr}{s^3}``

where ``M`` is the mass of the central body, ``r`` is the radius of the body being affected and ``s`` is the distance between the two bodies.
"""
tidal_acceleration(body_radius::Unitful.Length, parent_mass::Unitful.Mass, separation::Unitful.Length) = 2body_radius * G * parent_mass / separation^3

function radial_orbit_time(m1::Unitful.Mass, m2::Unitful.Mass, x0::Unitful.Length, v0::Unitful.Velocity, x::Unitful.Length)
	μ = G * (m1 + m2)
	w = 1/x0 - v0^2/2μ

	return (asin(sqrt(w * x)) - sqrt(w * x * (1 - w * x))) / sqrt(2μ * w^3)
end

function radial_orbit_displacement(m1::Unitful.Mass, m2::Unitful.Mass, x0::Unitful.Length, v0::Unitful.Velocity, t::Unitful.Time)
	μ = G * (m1 + m2)
	w = 1/x0 - v0^2/2μ
	p = cbrt(9μ * t^2 / 2)

	return p - w*p^2 / 5 - 3w^2 * p^3 / 175 - 23w^3 * p^4/7875 - 1894w^4 * p^5 / 3931875 - 3293w^5 * p^6 / 21896875 - 2418092w^6 * p^7/62077640625
end

function radial_orbit_displacement2(m1::Unitful.Mass, m2::Unitful.Mass, x0::Unitful.Length, v0::Unitful.Velocity, t::Unitful.Time)
	μ = G * (m1 + m2)
	w = 1/x0 - v0^2/2μ
	p = cbrt(9μ * t^2 / 2)

	return p - w*p^2 / 5 - 3w^2 * p^3 / 175 - 23w^3 * p^4/7875 - 1894w^4 * p^5 / 3931875 - 3293w^5 * p^6 / 21896875 - 2418092w^6 * p^7/62077640625
end

function radial_orbit_displacement(m1::Real, m2::Real, x0::Real, v0::Real, t::Real)
	μ = ustrip(G) * (m1 + m2)
	w = 1/x0 - v0^2/2μ
	p = cbrt(9μ * t^2 / 2)

	return p - w*p^2 / 5 - 3w^2 * p^3 / 175 - 23w^3 * p^4/7875 - 1894w^4 * p^5 / 3931875 - 3293w^5 * p^6 / 21896875 - 2418092w^6 * p^7/62077640625
end

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics