text:
# # Usage Guide

# In this example, we will present the basics of using Ripserer. We start by loading some
# packages.

using Distances
using Plots
using Ripserer
using Random # hide
Random.seed!(1337) # hide
gr() # hide
nothing # hide

# ## Using Ripserer With Point Cloud Data

# Let's start with generating some points, randomly sampled from a noisy circle.

function noisy_circle(n; r=1, noise=0.1)
    points = NTuple{2,Float64}[]
    for _ in 1:n
        θ = 2π * rand()
        push!(points, (r * sin(θ) + noise * rand(), r * cos(θ) + noise * rand()))
    end
    return points
end

circ_100 = noisy_circle(100)
scatter(circ_100; aspect_ratio=1, legend=false, title="Noisy Circle")

# !!! tip "Point-like data types"
#     Ripserer can interpret various kinds of data as point clouds. The limitation is that
#     the data set should be an `AbstractVector` with elements with the following
#     properties:
#     * all elements are collections numbers;
#     * all elements have the same length.
#     Examples of element types that work are `Tuple`s,
#     [`SVector`](https://github.com/JuliaArrays/StaticArrays.jl)s, and
#     [`Point`](https://github.com/JuliaGeometry/GeometryBasics.jl)s.

# To compute the Vietoris-Rips persistent homology of this data set, run the
# following.

ripserer(circ_100)

# You can use the `dim_max` argument to set the maximum dimension persistent homology is
# computed in.

result_rips = ripserer(circ_100; dim_max=3)

# The result can be plotted as a persistence diagram or as a barcode.

plot(result_rips)
barcode(result_rips)
plot(plot(result_rips), barcode(result_rips)) # hide

# We can also plot a single diagram or a subset of all diagrams in the same manner. Keep in
# mind that the result is just a vector of [`PersistenceDiagram`](@ref)s. The
# zero-dimensional diagram is found at index 1.

plot(result_rips[2])
barcode(result_rips[2:end]; linewidth=2)
plot(plot(result_rips[2]), barcode(result_rips[2:end]; linewidth=2)) # hide

# Plotting can be further customized using the standard attributes from
# [Plots.jl](http://docs.juliaplots.org/latest/).

plot(result_rips; markeralpha=1, markershape=:star, color=[:red, :blue, :green, :purple])

# ## Changing Filtrations

# By default, calling [`ripserer`](@ref) will compute persistent homology with the
# [`Rips`](@ref) filtration. To use a different filtration, we have two options.

# The first option is to pass the filtration constructor as the first argument. Any keyword
# arguments the filtration accepts can be passed to [`ripserer`](@ref) and it will be
# forwarded to the constructor.

ripserer(EdgeCollapsedRips, circ_100; threshold=1, dim_max=3, metric=Euclidean())

# The second option is to initialize the filtration object first and use that as an argument
# to [`ripserer`](@ref). This can be useful in cases where constructing the filtration takes
# a long time.

collapsed_rips = EdgeCollapsedRips(circ_100; threshold=1, metric=Euclidean())
ripserer(collapsed_rips; dim_max=3)

# ## Distance Matrix Inputs

# In the previous example, we got our result by passing a collection of points to
# [`ripserer`](@ref). Under the hood, [`Rips`](@ref) and [`EdgeCollapsedRips`](@ref)
# actually work with distance matrices. Let's define a distance matrix of the shortest
# paths on a [regular icosahedron](https://en.wikipedia.org/wiki/Regular_icosahedron) graph.

# ```@raw html
# <img src="https://upload.wikimedia.org/wikipedia/commons/8/83/Icosahedron_graph.svg" height="200" width="200">
# ```

icosahedron = [
    0 1 2 2 1 2 1 1 2 2 1 3
    1 0 3 2 1 1 2 1 2 1 2 2
    2 3 0 1 2 2 1 2 1 2 1 1
    2 2 1 0 3 2 1 1 2 1 2 1
    1 1 2 3 0 1 2 2 1 2 1 2
    2 1 2 2 1 0 3 2 1 1 2 1
    1 2 1 1 2 3 0 1 2 2 1 2
    1 1 2 1 2 2 1 0 3 1 2 2
    2 2 1 2 1 1 2 3 0 2 1 1
    2 1 2 1 2 1 2 1 2 0 3 1
    1 2 1 2 1 2 1 2 1 3 0 2
    3 2 1 1 2 1 2 2 1 1 2 0
]
nothing # hide

# To compute the persistent homology, simply feed the distance matrix to [`ripserer`](@ref).

result_icosa = ripserer(icosahedron; dim_max=2)

# ## Thresholding

# In our next example, we will show how to use thresholding to speed up computation. We
# start by defining a sampling function that generates ``n`` points from the square
# ``[-4,4]\times[-4,4]`` with a circular hole of radius 1 in the middle.

function cutout(n)
    points = NTuple{2,Float64}[]
    while length(points) < n
        x, y = (8rand() - 4, 8rand() - 4)
        if x^2 + y^2 > 1
            push!(points, (x, y))
        end
    end
    return points
end

# We sample 2000 points from this space.

cutout_2000 = cutout(2000)
scatter(cutout_2000; markersize=1, aspect_ratio=1, legend=false, title="Cutout")

# We calculate the persistent homology and time the calculation.

@time result_cut = ripserer(cutout_2000)
nothing # hide

#

plot(result_cut)

# Notice that while there are many 1-dimensional classes, one of them stands out. This class
# represents the hole in the middle of our square. Since the intervals are sorted by
# persistence, we know the last interval in the diagram will be the most persistent.

most_persistent = result_cut[2][end]

# Notice the death time of this interval is around 1.83 and that no intervals occur after
# that time. This means that we could stop computing when we reach this time and the result
# should not change. Let's try it out.

@time result_cut_thresh_2 = ripserer(cutout_2000; threshold=2)
nothing # hide

#

plot(result_cut_thresh_2; title="Persistence Diagram, threshold=2")

# Indeed, the result is exactly the same, but it took less than a third of the time to
# compute.

@assert result_cut_thresh_2 == result_cut # hide
result_cut_thresh_2 == result_cut

# If we pick a threshold that is too low, we still detect the interval, but its death time
# becomes infinite.

@time result_cut_thresh_1 = ripserer(cutout_2000; threshold=1)
nothing # hide

#

result_cut_thresh_1[2][end]

# ## Persistence Diagrams

# The result of a computation is returned as a vector of
# [`PersistenceDiagram`](@ref)s. Let's take a closer look at one of those.

diagram = result_cut[2]

# The diagram is a structure that acts as a vector of
# [`PersistenceInterval`](@ref)s. As such, you can use standard Julia
# functions on the diagram.

# For example, to extract the last three intervals by birth time, you can do
# something like this.

sort(diagram; by=birth, rev=true)[1:3]

# To find the [`persistence`](@ref)s of all the intervals, you can use broadcasting.

persistence.(diagram)

# Unlike regular vectors, a [`PersistenceDiagram`](@ref) has additional metadata attached
# to it. To see all metadata, use
# [`propertynames`](https://docs.julialang.org/en/v1/base/base/#Base.propertynames).

propertynames(diagram)

# You can access the properties with the dot syntax.

diagram.field

# The attributes `dim` and `threshold` are given special treatment and can be extracted
# with appropriately named functions.

dim(diagram), threshold(diagram)

# Now, let's take a closer look at one of the intervals.

interval = diagram[end]

# An interval is very similar to a tuple of two `Float64`s, but also has some metadata
# associated with it.

interval[1], interval[2]

# [`birth`](@ref), [`death`](@ref), [`persistence`](@ref), and [`midlife`](@ref) can be used
# to query commonly used values.

birth(interval), death(interval), persistence(interval), midlife(interval)

# Accessing metadata works in a similar manner as with diagrams.

propertynames(interval)

#

interval.birth_simplex

#

interval.death_simplex

# ## Simplices

# In the previous section, we saw each interval has an associated [`birth_simplex`](@ref)
# and [`death_simplex`](@ref). These values are of the type [`Simplex`](@ref). Let's take a
# closer look at simplices.

simplex = interval.death_simplex

# [`Simplex`](@ref) is an internal data structure that uses some tricks to increase
# efficiency. For example, if we were to
# [`dump`](https://docs.julialang.org/en/v1/base/io-network/#Base.dump) it, we notice the
# vertices are not actually stored in the simplex itself.

dump(simplex)

# To access the vertices, we use [`vertices`](@ref).

vertices(simplex)

# Other useful attributes a simplex has are [`index`](@ref), [`dim`](@ref), and
# [`birth`](@ref).

index(simplex), dim(simplex), birth(simplex)

# A few additional notes on simplex properties.

# * A `D`-dimensional simplex is of type `Simplex{D}` and has `D + 1` vertices.
# * [`vertices`](@ref) are always sorted in descending order.
# * [`index`](@ref) and [`dim`](@ref) can be used to uniquely identify a given simplex.
# * [`birth`](@ref) determines when a simplex is added to a filtration.

# ## Conclusion

# This concludes the basic usage of Ripserer. For more detailed information, please check
# out the [API](@ref) page, as well as other examples.


text:
using Dates: Hour, Minute, Second, Millisecond, days, hour, minute, second, millisecond

"""
    timezone(::ZonedDateTime) -> TimeZone

Returns the `TimeZone` used by the `ZonedDateTime`.
"""
timezone(zdt::ZonedDateTime) = zdt.timezone

Dates.days(zdt::ZonedDateTime) = days(DateTime(zdt))

for period in (:Hour, :Minute, :Second, :Millisecond)
    accessor = Symbol(lowercase(string(period)))
    @eval begin
        Dates.$accessor(zdt::ZonedDateTime) = $accessor(DateTime(zdt))
        Dates.$period(zdt::ZonedDateTime) = $period($accessor(zdt))
    end
end

Base.eps(::ZonedDateTime) = Millisecond(1)

ext: jl


text:
using ErrorFreeTransforms
using FactCheck

include("add.jl")
include("mul.jl")
include("fma.jl")

exitstatus()


text:
# This file was generated, do not modify it. # hide
pipe.knn.K = 2
pipe.hot.drop_last = true;


text:
mutable struct ref x end

k = 0
c = ref(0)
while (k <= 13)
	c.x = c.x + 1
	x = 7
	while (c.x % 3 != 0)
		y = x
		z = y * 9
		c.x = c.x + 1
	end
	k = k + 1
end
println(c.x)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
#
# oidata.jl --
#
# Implement OI-FITS data types.
#
#------------------------------------------------------------------------------
#
# This file is part of OIFITS.jl which is licensed under the MIT "Expat"
# License:
#
# Copyright (C) 2015: Éric Thiébaut.
#
#------------------------------------------------------------------------------

###############
# BASIC TYPES #
###############

# Conversion to real type (no-op. if already of the correct type).
to_real(x::Cdouble) = x
to_real(x::Array{Cdouble}) = x
to_real{T<:Real}(x::Array{T}) = convert(Array{Cdouble}, x)
to_real(x::Real) = convert(Cdouble, x)

to_integer(x::Int) = x
to_integer(x::Array{Int}) = x
to_integer{T<:Integer}(x::Array{T}) = convert(Array{Int}, x)
to_integer(x::Int) = convert(Int, x)

# OI-FITS files stores the following 4 different data types:
const _DTYPE_LOGICAL = 1 # for format letter 'L'
const _DTYPE_INTEGER = 2 # for format letters 'I' or 'J'
const _DTYPE_REAL    = 3 # for format letters 'D' or 'E'
const _DTYPE_STRING  = 4 # for format letter 'A'

# The following dictionary is used for quick conversion of FITS format
# letter to data type.
const _DATATYPES = Dict('l' =>  _DTYPE_LOGICAL,
                        'L' =>  _DTYPE_LOGICAL,
                        'i' =>  _DTYPE_INTEGER,
                        'I' =>  _DTYPE_INTEGER,
                        'j' =>  _DTYPE_INTEGER,
                        'J' =>  _DTYPE_INTEGER,
                        'e' =>  _DTYPE_REAL,
                        'E' =>  _DTYPE_REAL,
                        'd' =>  _DTYPE_REAL,
                        'D' =>  _DTYPE_REAL,
                        'a' =>  _DTYPE_STRING,
                        'A' =>  _DTYPE_STRING)

is_logical(::Any) = false
is_logical(::Bool) = true
is_logical(::Array{Bool}) = true

is_string(::Any) = false
is_string(::AbstractString) = true
is_string{T<:AbstractString}(::Array{T}) = true

is_integer(::Any) = false
is_integer(::Integer) = true
is_integer{T<:Integer}(::Array{T}) = true

is_real(::Any) = false
is_real(::Real) = true
is_real{T<:Real}(::Array{T}) = true

#################
# OI-FITS TYPES #
#################

typealias OIContents Dict{Symbol,Any}

abstract OIDataBlock

type OITarget <: OIDataBlock
    owner
    contents::OIContents
    OITarget(contents::OIContents) = new(nothing, contents)
end

type OIArray <: OIDataBlock
    owner
    contents::OIContents
    OIArray(contents::OIContents) = new(nothing, contents)
end

type OIWavelength <: OIDataBlock
    owner
    contents::OIContents
    OIWavelength(contents::OIContents) = new(nothing, contents)
end

type OISpectrum <: OIDataBlock
    owner
    contents::OIContents
    OISpectrum(contents::OIContents) = new(nothing, contents)
end

type OIVis <: OIDataBlock
    owner
    arr::Union{OIArray,Void}
    ins::Union{OIWavelength,Void}
    contents::OIContents
    OIVis(contents::OIContents) = new(nothing, nothing, nothing, contents)
end

type OIVis2 <: OIDataBlock
    owner
    arr::Union{OIArray,Void}
    ins::Union{OIWavelength,Void}
    contents::OIContents
    OIVis2(contents::OIContents) = new(nothing, nothing, nothing, contents)
end

type OIT3 <: OIDataBlock
    owner
    arr::Union{OIArray,Void}
    ins::Union{OIWavelength,Void}
    contents::OIContents
    OIT3(contents::OIContents) = new(nothing, nothing, nothing, contents)
end

# Correspondance between OI-FITS data-block names and Julia types.
_DATABLOCKS = Dict{ASCIIString,DataType}("OI_TARGET"     => OITarget,
                                         "OI_WAVELENGTH" => OIWavelength,
                                         "OI_ARRAY"      => OIArray,
                                         "OI_SPECTRUM"   => OISpectrum,
                                         "OI_VIS"        => OIVis,
                                         "OI_VIS2"       => OIVis2,
                                         "OI_T3"         => OIT3)
_EXTNAMES = Dict{DataType,ASCIIString}()

for (key, val) in _DATABLOCKS
    _EXTNAMES[val] = key
end

get_dbname(db::OIDataBlock) = _EXTNAMES[typeof(db)]

# OIData is any OI-FITS data-block which contains interferometric data.
typealias OIData Union{OIVis,OIVis2,OIT3}

# OIDataBlock can be indexed by the name (either as a string or as a
# symbol) of the field.
getindex(db::OIDataBlock, key::Symbol) = get(db.contents, key, nothing)
getindex(db::OIDataBlock, key::AbstractString) = getindex(db, symbol(key))
haskey(db::OIDataBlock, key::Symbol) = haskey(db.contents, key)
haskey(db::OIDataBlock, key::AbstractString) = haskey(db.contents, symbol(key))
keys(db::OIDataBlock) = keys(db.contents)

# OIDataBlock can be used as iterators.
start(db::OIDataBlock) = start(db.contents)
done(db::OIDataBlock, state) = done(db.contents, state)
next(db::OIDataBlock, state) = next(db.contents, state)

# OIMaster stores the contents of an OI-FITS file.  Data-blocks containing
# measurements (OI_VIS, OI_VIS2 and OI_T3) are stored into a vector and
# thus indexed by an integer.  Named data-blocks (OI_ARRAY and
# OI_WAVELENGTH) are indexed by their names (converted to upper case
# letters, with leading and trailing spaces stripped, multiple spaces
# replaced by a single ordinary space).
type OIMaster
    all::Vector{OIDataBlock}            # All data-blocks
    update_pending::Bool                # Update is needed?
    tgt::Union{OITarget,Void}
    arr::Dict{ASCIIString,OIArray}
    ins::Dict{ASCIIString,OIWavelength}
    function OIMaster()
        new(Array(OIDataBlock, 0),
            false,
            nothing,
            Dict{ASCIIString,OIArray}(),
            Dict{ASCIIString,OIWavelength}())
    end
end

function show(io::IO, db::OITarget)
    print(io, "OI_TARGET: ")
    if haskey(db, :target)
        tgt = db[:target]
        n = length(tgt)
        if n > 0
            print(io, "target=")
            for i in 1:n
                print(io, (i == 1 ? "[\"" : ", \""),
                      tgt[i], (i == n ? "\"]" : "\""))
            end
        end
    else
        n = 0
    end
    n == 0 && print(io, "<empty>")
end

function show(io::IO, db::OIArray)
    print(io, "OI_ARRAY: ")
    ntels = (haskey(db, :sta_index) ? length(db[:sta_index]) : 0)
    if haskey(db, :arrname)
        print(io, "arrname=\"", db[:arrname], "\", ", ntels, " telescope(s)")
    else
        print(io, ntels, " telescope(s)")
    end
end

function show(io::IO, db::OIWavelength)
    print(io, "OI_WAVELENGTH: ")
    if ! haskey(db, :insname)
        print(io, "<empty>")
    else
        print(io, "insname=\"", db[:insname], "\"")
        if haskey(db, :eff_wave)
            wave = db[:eff_wave]
            n = length(db[:eff_wave])
            if n < 1
                print(io, " with 0 spectral channels")
            elseif n == 1
                print(io, " with 1 spectral channel at ",
                      round(wave[1]*1e6, 3), " µm")
            else
                print(io, " with ", n, " spectral channels from ",
                      round(minimum(wave)*1e6, 3), " µm to ",
                      round(maximum(wave)*1e6, 3), " µm")
            end
        end
    end
end

function show(io::IO, db::OISpectrum)
    print(io, "OI_SPECTRUM")
end

function show(io::IO, db::OIDataBlock,
              name::AbstractString, nwaves::Integer, ntimes::Integer)
    print(io, name, ": ")
    if nwaves > 0 && ntimes > 0
        print(io, nwaves*ntimes, " measurements in ", nwaves,
              " spectral channel(s) and ", ntimes, " exposure(s)")
    else
        print(io, name, "<empty>")
    end
end

function show(io::IO, db::OIVis)
    if haskey(db, :visamp)
        dims = size(db[:visamp])
        nwaves = dims[1]
        ntimes = dims[2]
    else
        nwaves = 0
        ntimes = 0
    end
    show(io, db, "OI_VIS", nwaves, ntimes)
end

function show(io::IO, db::OIVis2)
    if haskey(db, :vis2data)
        dims = size(db[:vis2data])
        nwaves = dims[1]
        ntimes = dims[2]
    else
        nwaves = 0
        ntimes = 0
    end
    show(io, db, "OI_VIS2", nwaves, ntimes)
end

function show(io::IO, db::OIT3)
    if haskey(db, :t3amp)
        dims = size(db[:t3amp])
        nwaves = dims[1]
        ntimes = dims[2]
    else
        nwaves = 0
        ntimes = 0
    end
    show(io, db, "OI_T3", nwaves, ntimes)
end

function show(io::IO, master::OIMaster)
    print(io, "OI_MASTER: (", length(master.all), " data-block(s))")
    for db in master.all
        print(io, "\n    ")
        show(io, db)
    end
end

####################################
# PARSING OF DATABLOCK DEFINITIONS #
####################################

# OIFieldDef is used to store the definition of a keyword/column field.
type OIFieldDef
    name::ASCIIString    # keyword/column name as a string
    symb::Symbol         # keyword/column symbolic name
    keyword::Bool        # is keyword? (otherwise column)
    multiplier::Int      # multiplier; for keywords, 0 means optional and 1
                         # means required; for columns, a negative number
                         # means abs(n) times the number of spectral
                         # channels;
    dtype::Int           # data type
    units::ASCIIString   # units
    descr::ASCIIString   # description
end

# OIDataBlockDef is used to store the definition of data-block.
type OIDataBlockDef
    dbname::ASCIIString
    fields::Vector{Symbol}        # ordered field symbolic names
    spec::Dict{Symbol,OIFieldDef} # dictionary of field specifications
    function OIDataBlockDef(dbname::ASCIIString, vect::Vector{OIFieldDef})
        spec = Dict{Symbol,OIFieldDef}()
        fields = Array(Symbol, length(vect))
        for j in 1:length(vect)
            entry = vect[j]
            fields[j] = entry.symb
            spec[entry.symb] = entry
        end
        new(dbname, fields, spec)
    end
end

# OIFormatDef is used to store all the data-block definitions
# for a given revision number.
typealias OIFormatDef Dict{ASCIIString,OIDataBlockDef}

# _FORMATS array is indexed by the revision number.
_FORMATS = Array(OIFormatDef, 0)

# The default format version number is the highest registered one.
default_revision() = length(_FORMATS)

# _FIELDS is a dictionary indexed by the data-block name (e,g., "OI_VIS"), each
# entry stores a set of its fields.
_FIELDS = Dict{ASCIIString,Set{Symbol}}()

# get_def(db, revn) -- yields the OIDataBlockDef for datablock of type `db`
#                      (e.g., "OI_TARGET") in revison `revn` of OI-FITS
#                      standard.
#
function get_def(dbname::ASCIIString, revn::Integer)
    if revn < 0 || revn > length(_FORMATS)
        error("unsupported revision number: $revn")
    end
    if ! haskey(_FORMATS[revn], dbname)
        error("unknown data-block: $db")
    end
    _FORMATS[revn][dbname]
end

function get_def(db::OIDataBlock)
    get_def(get_dbname(db), get_revn(db))
end

function get_descr(db::OIDataBlock)
    name = get_dbname(db)
    revn = get_revn(db)
    return (name, revn, get_def(name, revn))
end

function add_def(dbname::ASCIIString, revn::Integer, tbl::Vector{ASCIIString})
    if ! startswith(dbname, "OI_") || dbname != uppercase(dbname) || contains(dbname, " ")
        error("invalid data-block name: \"$db\"")
    end
    if revn < 0 || revn > 2
        error("invalid revision number: $revn")
    end

    fields = get(_FIELDS, dbname, Set{Symbol}())
    def = Array(OIFieldDef, 0)
    keyword = true
    for j in 1:length(tbl)
        row = strip(tbl[j])
        m = match(r"^([^ ]+) +([^ ]+) +([^ ]+) +(.*)$", row)
        if m == nothing
            if match(r"^-+$", row) == nothing
                error("syntax error in OI_FITS definition: \"$row\"")
            end
            keyword = false
            continue
        end
        name = uppercase(m.captures[1])
        symb = name2symbol(name)
        dtype = get(_DATATYPES, m.captures[2][end], nothing)
        if dtype == nothing
            error("invalid format type in OI_FITS definition: \"$row\"")
        end
        multiplier = parse(Int, m.captures[2][1:end-1])
        units = m.captures[3]
        descr = m.captures[4]
        if keyword
            if ! (multiplier == 0 || multiplier == 1)
                error("multiplier must be 0 or 1 in OI_FITS definition: \"$row\"")
            end
        else
            if ! (multiplier == -1 || multiplier > 0)
                error("invalid multiplier in OI_FITS definition: \"$row\"")
            end
        end
        push!(def, OIFieldDef(name, symb, keyword, multiplier,
                              dtype, units, descr))
        push!(fields, symb)
    end

    # Insert the data-block definition in the global table.
    while revn > length(_FORMATS)
        push!(_FORMATS, OIFormatDef())
    end
    _FORMATS[revn][dbname] = OIDataBlockDef(dbname, def)
    _FIELDS[dbname] = fields
end

# This version takes care of converting the string into a symbol.
get_def(db::AbstractString, revn::Integer) =  get_def(symbol(db), revn)

# get_def(db) -- returns the definition for the default format version.
get_def(db) =  get_def(db, default_revision())

# Convert the name of an OI-FITS keyword/column into a valid symbol.
function name2symbol(name::AbstractString)
    key = lowercase(name)
    if key == "oi_revn"
        return :revn
    else
        return symbol(replace(key, r"[^a-z0-9_]", '_'))
    end
end

#########################
# BUILDING OF DATABLOCK #
#########################

# Build constructors for the various data-blocks types.
for (func, name) in ((:new_target,     "OI_TARGET"),
                     (:new_array,      "OI_ARRAY"),
                     (:new_wavelength, "OI_WAVELENGTH"),
                     (:new_spectrum,   "OI_SPECTRUM"),
                     (:new_vis,        "OI_VIS"),
                     (:new_vis2,       "OI_VIS2"),
                     (:new_t3,         "OI_T3"))
    @eval begin
        function $func(master::Union{OIMaster, Void}=nothing;
                       revn::Integer=default_revision(), args...)
            db = build_datablock($name, revn, args)
            if master != nothing
                attach!(master, db)
            end
            return db
        end
    end
end

function build_datablock(dbname::ASCIIString, revn::Integer, args)
    def = get_def(dbname, revn)
    contents = OIContents()
    nrows = -1     # number of measurements
    ncols = -1     # number of spectral channels
    nerrs = 0      # number of errors so far
    for (field, value) in args
        # Check whether this field exists.
        spec = get(def.spec, field, nothing)
        if spec == nothing
            error("data-block $dbname has no field \"$field\"")
        end

        # Check value type.
        if spec.dtype == _DTYPE_LOGICAL
            if ! is_logical(value)
                error("expecting boolean value for field \"$field\" in $dbname")
            end
        elseif spec.dtype == _DTYPE_INTEGER
            if ! is_integer(value)
                error("expecting integer value for field \"$field\" in $dbname")
            end
            value = to_integer(value)
        elseif spec.dtype == _DTYPE_REAL
            if ! is_real(value)
                error("expecting real value for field \"$field\" in $dbname")
            end
            value = to_real(value)
        elseif spec.dtype == _DTYPE_STRING
            if ! is_string(value)
                error("expecting string value for field \"$field\" in $dbname")
            end
        else
            error("*** CORRUPTED WORKSPACE ***")
        end

        # Check value dimensions.
        dims = (isa(value, Array) ? size(value) : ())
        rank = length(dims)
        if spec.keyword
            if rank != 0
                error("expecting a scalar value for field \"$field\" in $dbname")
            end
        else
            if spec.multiplier < 0 && rank == 1
                # Fix data arrays with a single spectral channel.
                dims = (1, length(value))
                value = reshape(value, dims)
                rank = length(dims)
            end
            n = (spec.multiplier == 1 || spec.dtype == _DTYPE_STRING ? 1 : 2)
            if rank != n
                error("expecting a $n-D array value for field \"$field\" in $dbname")
            end
            if nrows == -1
                nrows = dims[end]
            elseif nrows != dims[end]
                error("incompatible number of rows for value of field \"$field\" in $dbname")
            end
            if spec.multiplier < 0
                if ncols == -1
                    ncols = dims[2]
                elseif ncols != dims[2]
                    error("incompatible number of columns for value of field \"$field\" in $dbname")
                end
            end
        end

        # Store value.
        contents[field] = value
    end

    # Check that all mandatory fields have been given.
    nerrs = 0
    for field in def.fields
        spec = def.spec[field]
        if ! haskey(contents, field) && spec.multiplier != 0
            warn("missing value for field \"$field\" in $dbname")
            nerrs += 1
        end
    end
    if nerrs > 0
        error("some fields are missing in $dbname")
    end
    return _DATABLOCKS[dbname](contents)
end

###################################################
# MANAGING THE DATABLOCK CONTAINER (THE "MASTER") #
###################################################

# Letter case and trailing spaces are insignificant according to FITS
# conventions.
fixname(name::AbstractString) = uppercase(rstrip(name))

# Make OIMaster an iterator.
start(master::OIMaster) = start(update!(master).all)
done(master::OIMaster, state) = done(master.all, state)
next(master::OIMaster, state) = next(master.all, state)

function select(master::OIMaster, args::AbstractString...)
    datablocks = Array(OIDataBlock, 0)
    for db in master
        if get_dbname(db) ∈ args
            push!(datablocks, db)
        end
    end
    return datablocks
end

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_target(mst)
```
yields the `OITarget` data-block of `mst` if any, `nothing` otherwise.

Assuming `tgt` is an instance of `OITarget`, then:
```
    get_target(tgt)
```
yields the "TARGET" column of `tgt` which is an array of target names.
"""
get_target(master::OIMaster) = update!(master).tgt

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_array(mst, arrname)
```
yields the `OIArray` data-block of `mst` whose name is `arrname` if any,
`nothing` otherwise.

Assuming `db` is an instance of a sub-type of `OIDataBlock`, then:
```
    get_array(db)
```
yields corresponding instance of `OIArray` if any, `nothing` otherwise.
"""
function get_array(master::OIMaster, arrname::AbstractString)
    get(update!(master).arr, fixname(arrname), nothing)
end
get_array(db::OIDataBlock) = nothing
get_array(db::Union{OIVis,OIVis2,OIT3}) = db.arr
get_array(db::OIArray) = db

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_instrument(mst, insname)
```
yields the `OIWavelength` data-block of `mst` whose name is `insname` if any,
`nothing` otherwise.

Assuming `db` is an instance of a sub-type of `OIDataBlock`, then:
```
    get_instrument(db)
```
yields corresponding instance of `OIWavelength` if any, `nothing` otherwise.
"""
function get_instrument(master::OIMaster, insname::AbstractString)
    get(update!(master).ins, fixname(insname), nothing)
end
get_instrument(db::OIDataBlock) = nothing
get_instrument(db::Union{OIVis,OIVis2,OIT3}) = db.ins
get_instrument(db::OIWavelength) = db

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_targets(mst)
```
yields the names of the targets defined in `mst`.
"""
function get_targets(master::OIMaster)
    tgt = get_target(master)
    return tgt == nothing ? Array(ASCIIString, 0) : get_target(tgt)
end

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_arrays(mst)
```
yields the names of the interferometric arrays defined in `mst`.
"""
get_arrays(master::OIMaster) = collect(keys(update!(master).arr))

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_instruments(mst)
```
yields the names of the instruments defined in `mst`.
"""
get_instruments(master::OIMaster) = collect(keys(update!(master).ins))

function new_master(datablocks::OIDataBlock...)
    master = OIMaster()
    for db in datablocks
        attach!(master, db)
    end
    return update!(master)
end

function new_master(datablocks::Array{OIDataBlock})
    master = OIMaster()
    for db in datablocks
        attach!(master, db)
    end
    return update!(master)
end

function attach!(master::OIMaster, db::OIDataBlock)
    db.owner == nothing || error("data-block already attached")
    if isa(db, OITarget)
        if master.tgt != nothing
            error("only one OI_TARGET data-block can be attached")
        end
        master.tgt = db
    elseif isa(db, OIWavelength)
        insname = fixname(db[:insname])
        if haskey(master.ins, insname)
            error("master already have an OI_WAVELENGTH data-block with INSNAME=\"$insname\"")
        end
        master.ins[insname] = db
    elseif isa(db, OIArray)
        arrname = fixname(db[:arrname])
        if haskey(master.arr, arrname)
            error("master already have an OI_ARRAY data-block with ARRNAME=\"$arrname\"")
        end
        master.arr[arrname] = db
    end
    push!(master.all, db)
    master.update_pending = true
    db.owner = master
    nothing
end

function update!(master::OIMaster)
    if master.update_pending
        if master.tgt == nothing
            error("missing mandatory OI_TARGET data-block")
        end
        for db in master.all
            if isa(db, OIData)
                insname = fixname(db[:insname])
                arrname = fixname(db[:arrname])
                db.ins = get(master.ins, insname, nothing)
                db.arr = get(master.arr, arrname, nothing)
                if db.ins == nothing
                    error("OI_WAVELENGTH data-block with INSNAME=\"$insname\" not found in master")
                end
                if db.arr == nothing
                    warn("OI_ARRAY data-block with ARRNAME=\"$arrname\" not found in master")
                end
            end
        end
        master.update_pending = false
    end
    return master
end

"""
### Clone a data-block

This method clones an existing data-block.  The array data are shared between
the clones but not the links (owner, target, instrument, etc.).  Only the
fields defined by OI-FITS standard are cloned.  This method is intended when
a data-block from a given `OIMaster` instance is to be attached to another
`OIMaster` instance.
"""
function clone(db::OIDataBlock)
    (name, revn, defn) = get_descr(db)
    data = Dict{Symbol,Any}()
    for (key, val) in db.contents
        if haskey(defn.spec, key)
            data[key] = val
        end
    end
    return build_datablock(name, revn, data)
end

"""
### Select data for a given target

The method `select_target` selects a subset of data corresponding to a given
target.   The general syntax is:
```
    out = select_target(inp, tgt)
```
where `inp` is the input data (can be an instance of `OIMaster` or of any
`OIDataBlock` sub-types), `tgt` is the target number or name.

The result `out` may share part of its contents with the input data `inp`.

The result may be `nothing` if the input contains no data for the given target.
"""
select_target(db::OIDataBlock, target::Integer) = clone(db)

function select_target(db::OITarget, target::Integer)
    (name, revn, defn) = get_descr(db)
    k = findfirst(id -> id == target, get_target_id(db))
    k == 0 && return
    data = Dict{Symbol,Any}()
    for (key, val) in db.contents
        spec = get(defn.spec, key, nothing)
        spec != nothing || continue
        if spec.keyword
            data[key] = db[key]
        else
            data[key] = db[key][k:k]
        end
    end
    return build_datablock(name, revn, data)
end

function select_target(db::OIData, target::Integer)
    (name, revn, defn) = get_descr(db)
    target_id = get_target_id(db)
    sel = find(id -> id == target, target_id)
    length(sel) > 0 || return
    data = Dict{Symbol,Any}()
    cpy = (length(sel) == length(target_id)) # just copy?
    for (key, val) in db.contents
        spec = get(defn.spec, key, nothing)
        spec != nothing || continue
        if cpy || spec.keyword
            data[key] = db[key]
        else
            # Copy a sub-array corresponding to the selection.  (The target is
            # specified for each last index.)
            src = db[key]
            @assert(size(src)[end] == length(target_id))
            rank = ndims(src)
            dims = ntuple(i -> i == rank ? length(sel) : size(src, i), rank)
            dst = Array(eltype(src), dims)
            data[key] = dst
            if rank == 1
                for j in 1:length(sel)
                    dst[j] = src[sel[j]]
                end
            elseif rank == 2
                for j in 1:length(sel)
                    dst[:,j] = src[:,sel[j]]
                end
            else
                error("unexpected rank $rank")
            end
        end
    end
    return build_datablock(name, revn, data)
end

function select_target(master::OIMaster, target::Integer)
    result = OIMaster()
    for db in master
        sel = select_target(db, target)
        sel != nothing && attach!(result, sel)
    end
    return update!(result)
end

function select_target(master::OIMaster, target::ASCIIString)
    db = get_target(master)
    if db != nothing
        k = findfirst(name -> name == target, get_target(db))
        k != 0 && return select_target(master, get_target_id(db)[k])
    end
    error("target name not found")
end

"""
### Select data at given wavelengths

The method `select_wavelength` selects a subset of data on the basis of their
wavelength.   The general syntax is:
```
    out = select_wavelength(inp, sel)
```
where `inp` is the input data (can be an instance of `OIMaster` or of any
`OIDataBlock` sub-types), `sel` is a predicate function which takes a
wavelength (in meters) as argument and returns true if this wavelength is to be
selected and false otherwise.

An alternative is:
```
    out = select_wavelength(inp, wmin, wmax)
```
to select wavelengths `w` such that `wmin ≤ w ≤ wmax`.  The wavelength
bounds are in meters.

The result `out` may share part of its contents with the input data `inp`.

The result may be `nothing` if the input contains no data at selected
wavelengths.
"""
function select_wavelength(inp::Union{OIMaster,OIDataBlock},
                           wavemin::Real, wavemax::Real)
    const wmin = convert(Cdouble, wavemin)
    const wmax = convert(Cdouble, wavemax)
    select_wavelength(inp, w -> wmin ≤ w ≤ wmax)
end

select_wavelength(db::OIDataBlock, selector::Function) = clone(db)

function select_wavelength(db::Union{OIWavelength,OIVis,OIVis2,OIT3,OISpectrum},
                           selector::Function)
    (name, revn, defn) = get_descr(db)
    wave = get_eff_wave(db)
    sel = find(selector, wave)
    length(sel) > 0 || return
    data = Dict{Symbol,Any}()
    cpy = (length(sel) == length(wave)) # just copy?
    iswave = (typeof(db) == OIWavelength)
    for (key, val) in db.contents
        spec = get(defn.spec, key, nothing)
        spec != nothing || continue
        if cpy || spec.keyword || (spec.multiplier >= 0 && ! iswave)
            data[key] = db[key]
        else
            # Copy a sub-array corresponding to the selection.  (The wavelength
            # corresponds to the first index.)
            src = db[key]
            @assert(size(src, 1) == length(wave))
            rank = ndims(src)
            dims = ntuple(i -> i == 1 ? length(sel) : size(src, i), rank)
            dst = Array(eltype(src), dims)
            data[key] = dst
            if rank == 1
                for j in 1:length(sel)
                    dst[j] = src[sel[j]]
                end
            elseif rank == 2
                for j in 1:length(sel)
                    dst[j,:] = src[sel[j],:]
                end
            else
                error("unexpected rank $rank")
            end
        end
    end
    return build_datablock(name, revn, data)
end

function select_wavelength(master::OIMaster, selector::Function)
    result = OIMaster()
    for db in master
        sel = select_wavelength(db, selector)
        sel != nothing && attach!(result, sel)
    end
    return update!(result)
end


text:
abstract AbstractGrid{T,N} <: AbstractArray{T,N}

abstract AbstractGrid1d{T<:Number,N} <: AbstractGrid{T,N}
abstract AbstractGrid2d{T<:Number,N} <: AbstractGrid{Point2{T},N}
abstract AbstractGrid3d{T<:Number,N} <: AbstractGrid{Point3{T},N}

Base.length(g::AbstractGrid) = length(g.c)
Base.eltype{T}(g::AbstractGrid{T}) = T
Base.ndims(g::AbstractGrid) = ndims(g.c)
Base.size(g::AbstractGrid) = size(g.c)
Base.size(g::AbstractGrid, i::Integer) = size(g.c, i)
Base.getindex(g::AbstractGrid, i...) = g.c[i...]
function Base.setindex!(g::AbstractGrid, value, i...)
    g.c[i...] = value
end
Base.transpose(g::AbstractGrid) = transpose(g.c)
Base.ctranspose(g::AbstractGrid) = ctranspose(g.c)
Base.eachindex(g::AbstractGrid) = eachindex(g.c)
Base.linearindexing{T<:AbstractGrid}(::Type{T}) = Base.LinearFast()
# Base.convert(T::Type{Array},g::AbstractGrid) = convert(T,g.c)

integrate(g::AbstractGrid, f) = integrate(g, f(g))
integrate(g::AbstractGrid, f::AbstractVector) = vecdot(g.w, f)

function integrate(g::AbstractGrid, f::AbstractArray)
    n,m = length(f),length(g)
    shp = (m,Int(n/m))
    At_mul_B(reshape(f, shp), vec(g.w))
end

# integrate(g::AbstractGrid, f::AbstractArray) = sum(g.w .* f, 1:ndims(g.w)) |> squeezeall

# integrate{T}(g::AbstractGrid{T,1}, f::AbstractMatrix) = At_mul_B(f, g.w)

∫ = integrate

function Base.writemime{T}(io::IO, mime::MIME"text/plain", g::AbstractGrid{T,1})
    write(io, summary(g), "\n")
    write(io, "[coords weights] =\n")
    writemime(io, mime, [g.c g.w])
end

function Base.writemime{T}(io::IO, mime::MIME"text/plain", g::AbstractGrid{T,2})
    write(io, summary(g), "\n")
    write(io, "coords =\n")
    writemime(io, mime, g.c)
end

function setposition(g::AbstractGrid3d, R::Point3{Float64})
    for i in 1:length(g)
        g[i] += R - g.R
    end
end

type LogarithmicGrid{T,A<:AbstractVector} <: AbstractGrid1d{T,1}
    c::A  # coords
    w::A  # weights
    dx::T  # spacing in log space
    function LogarithmicGrid(c::AbstractVector{T},w::AbstractVector{T},dx::T)
        @assert size(c) == size(w)
        new(c,w,dx)
    end
end

LogarithmicGrid{T}(c::AbstractVector{T},w::AbstractVector{T},dx::T) =
    LogarithmicGrid{eltype(c),typeof(c)}(c,w,dx)

function loggrid(start::Real, stop::Real, n::Int)
    x = collect(linspace(start, stop, n))
    dx = x[2] - x[1]
    r = exp(x)
    w = r.*dx
    LogarithmicGrid(r, w, dx)
end

# returns a logarithmic grid over same range but with different number of points
loggrid(r::LogarithmicGrid, n::Int) = loggrid(log(r[1]),log(r[end]),n)

type UniformGrid{T,A<:AbstractVector} <: AbstractGrid1d{T,1}
    c::A  # coordinates
    w::A  # weights
    dx::T  # spacing
    function UniformGrid(c::AbstractVector{T},w::AbstractVector{T},dx::T)
        @assert size(c) == size(w)
        new(c,w,dx)
    end
end

UniformGrid{T}(c::AbstractVector{T},w::AbstractVector{T},dx::T) =
    UniformGrid{eltype(c),typeof(c)}(c,w,dx)

function lingrid(start::Real, stop::Real, n::Int)
    x = collect(linspace(start, stop, n))
    dx = x[2]-x[1]
    UniformGrid(x, fill(dx, length(x)), dx)
end

function simpsgrid(start::Real, stop::Real, n::Int)
    @assert isodd(n) "Simpson's rule requires an odd number of grid points"
    x = collect(linspace(start, stop, n))
    dx = x[2] - x[1]
    w = fill(dx/3, n)
    w[2:2:end] .*= 4
    w[3:2:end-1] .*= 2
    UniformGrid(x, w, dx)
end

type AtomicRadialGrid{T,A<:AbstractVector} <: AbstractGrid1d{T,1}
    c::A  # coords
    w::A  # weights
    rmax::T
    function AtomicRadialGrid(c::AbstractVector{T},w::AbstractVector{T},rmax::T)
        @assert size(c) == size(w)
        new(c,w,rmax)
    end
end

AtomicRadialGrid{T}(c::AbstractVector{T},w::AbstractVector{T},rmax::T) =
    AtomicRadialGrid{eltype(c),typeof(c)}(c,w,rmax)

function atomicradgrid(n::Int, rmax::Real)
    s = collect(1:n)
    r = rmax * log(1 - (s/(n+1.)).^2) ./ log(1 - (n/(n+1))^2)
    w = -2rmax*s ./ ( (1+n-s) .* (1+n+s) .* log(1 - (n/(n+1.))^2) )
    AtomicRadialGrid(r, w, rmax)
end

loggrid(r::AtomicRadialGrid, n::Int) = loggrid(log(r[1]),log(r[end]),n)

type AngularGrid{T,A<:AbstractVector,B<:AbstractVector} <: AbstractGrid3d{T,1}
    c::A  # coords
    w::B  # weights
    function AngularGrid(c::AbstractVector{Point3{T}},w::AbstractVector{T})
        @assert size(c) == size(w)
        new(c,w)
    end
end

AngularGrid{T}(c::AbstractVector{Point3{T}},w::AbstractVector{T}) =
    AngularGrid{eltype(w),typeof(c),typeof(w)}(c,w)

typealias AngularGridf AngularGrid{Float64,Vector{Tuple{Float64,Float64,Float64}},Vector{Float64}}

lebgrid(n::Int) = lebedev[n]::AngularGridf

# function expand(grid, basis, f)
#     S = basis.data' * (basis .* x.w)
#     (integrate(x, basis.*h)' / S)'
# end

# function expand(grid, basis, f)
#     S = Symmetric(LinAlg.BLAS.syrk('U', 'T', basis.data .* sqrt(grid.w)))
#     (integrate(grid, basis.*f)' / S)'
# end

# LinAlg.BLAS.herk('U', 'C', 1.0, mat(basis.data .* sqrt(g.w)))

function qrgramschmidt(grid, basis)
    Q,R = qr(basis .* sqrt(grid.w));
    Q ./= sqrt(grid.w)
    Q
end

max_stars_repo_name: johncsnyder/Element


text:
export readaln, readseeds, writealn

"""
    readseeds(file::AbstractString,
              nodes1::AbstractVector,
              nodes2::AbstractVector) -> Matrix{Int} : n x 2

Outputs n x 2 matrix of node indices
associates with nodes1 and nodes2
"""
function readseeds(fd::IO,
                   nodes1::AbstractVector,nodes2::AbstractVector)
    aln = readdlm(fd,String)
    d1 = indexmap(nodes1)
    d2 = indexmap(nodes2)
    paln = zeros(Int,size(aln))
    for i = 1:size(paln,1)
        paln[i,1] = d1[aln[i,1]]
        paln[i,2] = d2[aln[i,2]]
    end
    paln
end
readseeds(file::AbstractString, args...) =
    open(fd -> readseeds(fd, args...), file, "r")

"""
    readaln(file::AbstractString, nodes1::Vector,
             nodes2::Vector, flip=false)

Read alignment file for pairwise network alignment.  Each line will
contain a node pair, with the first node from nodes1, and the second
node from nodes2.  Returns permutation from nodes1 to nodes2
corresponding to the node pairs (so need length(nodes1) <=
length(nodes2)) If flip=true, then returns permutation from nodes2 to
nodes1, (so need length(nodes2) <= length(nodes1)) where first node in
each line is from nodes2, and the second node is from nodes1.
"""
function readaln(fd::IO, nodes1::AbstractVector,
                 nodes2::AbstractVector, flip=false)
    aln = readdlm(fd,String)
    if flip
        aln2perm(flipdim(aln,2),nodes2,nodes1)
    else
        aln2perm(aln,nodes1,nodes2)
    end
end
readaln(file::AbstractString, args...) =
    open(fd -> readaln(fd, args...), file, "r")

import NetalignMeasures: aln2perm
"""
Given aln :: n x 2 matrix of node names, 1st column from nodes1,
second column from nodes2, convert to a permutation from nodes1 to
nodes2 If fill=true and there are nodes in nodes2 that are not
considered, map to those nodes randomly, but all nodes from nodes1 need
to be considered.
"""
function aln2perm(aln::AbstractMatrix,nodes1::AbstractVector,nodes2::AbstractVector,fill=true)
    m = length(nodes1)
    n = length(nodes2)
    if size(aln,1)!=m error("Need aln size to match") end
    d1 = indexmap(nodes1)
    d2 = indexmap(nodes2)
    p = zeros(Int, ifelse(fill,n,m))
    for i = 1:m
        p[d1[aln[i,1]]] = d2[aln[i,2]]
    end
    if fill
        mapped_to = zeros(Int,n)
        for i = 1:n
            if p[i] > 0
                mapped_to[p[i]] = 1
            end
        end
        j = 1
        for i = 1:n
            if mapped_to[i] == 0
                mapped_to[j] = i
                j += 1
            end
        end
        rp = randperm(n-m)
        for i = 1:(n-m)
            p[m+i] = mapped_to[rp[i]]
        end
    end
    p
end

"""
    writealn(fd::IO, nodes1::AbstractVector, nodes2::AbstractVector)
    writealn(file::AbstractString, nodes1::AbstractVector, nodes2::AbstractVector)

Write alignment to file
"""
function writealn(fd::IO, nodes1::AbstractVector, nodes2::AbstractVector)
    for i = 1:length(f)
        println(fd, nodes1[i], " ", nodes2[f[i]])
    end
end
writealn(file::AbstractString, args...) =
    open(fd -> writealn(fd,args...), file,"w")


text:
using Test
using Fjage

# start fjåge

println("Starting fjåge...")
cp = replace(read(`find lib -name '*.jar'`, String), "\n" => ":")
master = open(`java -cp $cp org.arl.fjage.shell.GroovyBoot etc/initrc.groovy`)
sleep(2)

# tests

println("Starting tests...")
try

  @testset "Fjage" begin

    gw = Gateway("localhost", 5081)
    @testset "Gateway" begin
      @test typeof(gw) <: Gateway
      @test typeof(gw.agentID) <: AgentID
    end

    shell = agentforservice(gw, "org.arl.fjage.shell.Services.SHELL")
    @testset "agentforservice (+)" begin
      @test typeof(shell) <: AgentID
      @test shell.name == "shell"
      @test shell.istopic == false
    end

    @testset "agentforservice (-)" begin
      dummy = agentforservice(gw, "org.arl.dummy.Services.DUMMY")
      @test dummy == nothing
    end

    @testset "agentsforservice (+)" begin
      alist = agentsforservice(gw, "org.arl.fjage.shell.Services.SHELL")
      @test typeof(alist) <: Array
      @test length(alist) == 1
      @test typeof(alist[1]) <: AgentID
      @test alist[1].name == "shell"
      @test alist[1].istopic == false
    end

    @testset "agentsforservice (-)" begin
      alist = agentsforservice(gw, "org.arl.dummy.Services.DUMMY")
      @test typeof(alist) <: Array
      @test length(alist) == 0
    end

    ShellExecReq = MessageClass("org.arl.fjage.shell.ShellExecReq")
    @testset "MessageClass" begin
      @test ShellExecReq <: Message
    end

    @testset "send & receive (gw)" begin
      send(gw, ShellExecReq(recipient=shell, cmd="ps"))
      rsp = receive(gw, 1000)
      @test typeof(rsp) <: Message
      @test rsp.performative == "AGREE"
    end

    @testset "send & receive (aid)" begin
      send(shell, ShellExecReq(cmd="ps"))
      rsp = receive(gw, 1000)
      @test typeof(rsp) <: Message
      @test rsp.performative == "AGREE"
    end

    @testset "request (gw)" begin
      rsp = request(gw, ShellExecReq(recipient=shell, cmd="ps"))
      @test typeof(rsp) <: Message
      @test rsp.performative == "AGREE"
    end

    @testset "request (aid)" begin
      rsp = request(shell, ShellExecReq(cmd="ps"))
      @test typeof(rsp) <: Message
      @test rsp.performative == "AGREE"
    end

    @testset "<< (aid, +)" begin
      rsp = shell << ShellExecReq(cmd="ps")
      @test typeof(rsp) <: Message
      @test rsp.performative == "AGREE"
    end

    dummy = agent(gw, "dummy")
    @testset "agent" begin
      @test typeof(dummy) <: AgentID
    end

    @testset "<< (aid, -)" begin
      rsp = dummy << ShellExecReq(cmd="ps")
      @test rsp == nothing
    end

    @testset "flush" begin
      send(gw, ShellExecReq(recipient=shell, cmd="ps"))
      sleep(1)
      flush(gw)
      rsp = receive(gw, 1000)
      @test rsp == nothing
    end

    ntf = topic(gw, "broadcast")
    @testset "topic" begin
      @test typeof(ntf) <: AgentID
    end

    @testset "subscribe (-)" begin
      send(ntf, ShellExecReq(cmd="ps"))
      msg = receive(gw, 1000)
      @test msg == nothing
    end

    @testset "subscribe (+)" begin
      subscribe(gw, ntf)
      send(ntf, ShellExecReq(cmd="ps"))
      msg = receive(gw, 1000)
      @test typeof(msg) <: ShellExecReq
    end

    @testset "unsubscribe" begin
      unsubscribe(gw, ntf)
      send(ntf, ShellExecReq(cmd="ps"))
      msg = receive(gw, 1000)
      @test msg == nothing
    end

    close(gw)

  end

finally

# stop fjåge

  println("Stopping fjåge...")
  kill(master)

end


text:
using Statistics

function split(X, y::AbstractVector; dims=1, ratio_train=0.8, kwargs...)
    n = length(y)
    size(X, dims) == n || throw(DimensionMismatch("..."))

    n_train = round(Int, ratio_train*n)
    i_rand = randperm(n)
    i_train = i_rand[1:n_train]
    i_test = i_rand[n_train+1:end]

    return selectdim(X, dims, i_train), y[i_train], selectdim(X, dims, i_test), y[i_test]
end

function normalize(X_train, X_test; dims=1, kwargs...)
    col_mean = mean(X_train; dims)
    col_std = std(X_train; dims)

    return (X_train .- col_mean) ./ col_std, (X_test .- col_mean) ./ col_std
end

function prepare_data(X, y; do_normal=true, do_onehot=true, kwargs...)
    X_train, y_train, X_test, y_test = split(X, y; kwargs...)

    if do_normal
        X_train, X_test = normalize(X_train, X_test; kwargs...)
    end

    classes = unique(y)

    if do_onehot
        y_train = onehot(y_train, classes)
        y_test = onehot(y_test, classes)
    end

    return X_train, y_train, X_test, y_test, classes
end

function onehot(y, classes)
    y_onehot = falses(length(classes), length(y))
    for (i, class) in enumerate(classes)
        y_onehot[i, y .== class] .= 1
    end
    return y_onehot
end

onecold(y, classes) = [classes[argmax(y_col)] for y_col in eachcol(y)]

accuracy(x, y) = mean(onecold(m(x), classes) .== onecold(y, classes))

max_stars_repo_name: petrmvala/JuliaCourse

max_issues_repo_name: petrmvala/JuliaCourse

max_forks_repo_name: petrmvala/JuliaCourse


text:
using Revise
using PyPlot
using AdFem
using MAT

function calc_residual_and_jacobian(nu, uvnext, uv)
    unext, vnext = uvnext[1:mmesh.ndof], uvnext[mmesh.ndof+1:end]
    u, v = uv[1:mmesh.ndof], uv[mmesh.ndof+1:end]
    ugauss = eval_grad_on_gauss_pts1(unext, mmesh)
    vgauss = eval_grad_on_gauss_pts1(vnext, mmesh)
    ux, uy = ugauss[:,1], ugauss[:,2]
    vx, vy = vgauss[:,1], vgauss[:,2]
    ug = dof_to_gauss_points(unext, mmesh)
    vg = dof_to_gauss_points(vnext, mmesh)
    uo = dof_to_gauss_points(u, mmesh)
    vo = dof_to_gauss_points(v, mmesh)
    F = compute_fem_source_term((ug-uo)/Δt + ug * ux + vg * uy, (vg-vo)/Δt + ug * vx + vg * vy,  mmesh)
    G = [compute_fem_laplace_term1(unext, nu, mmesh); compute_fem_laplace_term1(vnext, nu, mmesh)]
    r = F - G

    F11 = compute_fem_mass_matrix1(1/Δt + ux, mmesh) + compute_fem_advection_matrix1(ug, vg, mmesh)
    F12 = compute_fem_mass_matrix1(uy, mmesh)
    F22 = compute_fem_mass_matrix1(1/Δt + vy, mmesh) + compute_fem_advection_matrix1(ug, vg, mmesh)
    F21 = compute_fem_mass_matrix1(vx, mmesh)
    DF = [F11 F12; F21 F22]

    DG = [compute_fem_laplace_matrix1(nu, mmesh) spzero(mmesh.ndof)
        spzero(mmesh.ndof) compute_fem_laplace_matrix1(nu, mmesh)]
    J = DF - DG

    r, J
end

function calc_next_step(θ, uv)
    function f(θ, uvnext)
        r, J = calc_residual_and_jacobian(θ, uvnext, uv)
        J, r = impose_Dirichlet_boundary_conditions(J, r, bdnode, zeros(length(bdnode)))
        r, J
    end
    uvnext = newton_raphson_with_grad(f, uv, θ)
end

"""
θ -- viscosity field (must be a vector of size ngauss)
NT -- total number of time steps
u0 -- initial condition
"""
function solve_burgers(u0, NT, θ)
    function condition(i, uarr)
        i <= NT
    end
    function body(i, uarr)
        uv = read(uarr, i)
        uvnext = calc_next_step(θ, uv)
        i+1, write(uarr, i+1, uvnext)
    end
    i = constant(1, dtype = Int32)
    uarr = TensorArray(NT+1)
    uarr = write(uarr, 1, u0)
    _, u = while_loop(condition, body, [i, uarr])
    set_shape(stack(u), NT+1, length(u0))
end


text:
using MocosSim: tanh_modulation

@testset "InfectionModulations" begin
  @testset "ReferenceFunction" begin
    @test tanh_modulation(-Inf, 5, 1, -5, 7) == -5
    @test tanh_modulation(Inf, 5, 1, -5, 7) == 7
    @test tanh_modulation(5, 5, 1, -5, 7) == (-5+7)/2

    @test tanh_modulation(-Inf, 5, 1, 3, 3) == 3
    @test tanh_modulation(0, 5, 1, 3, 3) == 3
    @test tanh_modulation(+Inf, 5, 1, 3, 3) == 3

    @testset "Randomized testing" begin
      rng = Random.MersenneTwister(13)
      for i in 1:100
        loc = rand(rng) * 200 - 100
        scale = rand(rng) * 10
        lhs = rand(rng) * 20 - 10
        rhs = rand(rng) * 20 - 10
        @test tanh_modulation(-Inf, loc, scale, lhs, rhs) ≈ lhs
        @test tanh_modulation(+Inf, loc, scale, lhs, rhs) ≈ rhs
        @test tanh_modulation(loc, loc, scale, lhs, rhs) ≈ (lhs + rhs)/2
      end
    end
  end

  struct TanhMockState <: MocosSim.AbstractSimState
    rng::MersenneTwister
  end

  MocosSim.time(::TanhMockState) = 10
  MocosSim.numdetected(::TanhMockState) = 0
  MocosSim.numdead(::TanhMockState) = 0

  struct TanhMockParams <: MocosSim.AbstractSimParams end

@testset "TanhModulation" begin
  modulation = MocosSim.TanhModulation(
      weight_detected=0,
      weight_deaths=0,
      weight_days=1,
      loc=10.0,
      scale=1.0,
      limit_value=0.0)

  mock_state = TanhMockState(MersenneTwister(13))
  strain = MocosSim.ChineseStrain

  hits = 0
  num_samples = 100
  for i in 1:num_samples
    household_infection = Event(Val(MocosSim.TransmissionEvent), 0.0, 0, 0, MocosSim.HouseholdContact, strain)
    constant_infection = Event(Val(MocosSim.TransmissionEvent), 0.0, 0, 0, MocosSim.ConstantKernelContact, strain)

    @test modulation(mock_state, TanhMockParams(), household_infection) == true

    hits += modulation(mock_state, TanhMockParams(), constant_infection)
  end

  @test hits == num_samples/2

end

@testset "IncreasingTanhModulation" begin
  modulation = MocosSim.IncreasingTanhModulation(
      weight_detected=0,
      weight_deaths=0,
      weight_days=1,
      loc=10.0,
      scale=1.0,
      initial_value=0.0)

  mock_state = TanhMockState(MersenneTwister(13))
  strain = MocosSim.ChineseStrain

  hits = 0
  num_samples = 100
  for i in 1:num_samples
    household_infection = Event(Val(MocosSim.TransmissionEvent), 0.0, 0, 0, MocosSim.HouseholdContact, strain)
    constant_infection = Event(Val(MocosSim.TransmissionEvent), 0.0, 0, 0, MocosSim.ConstantKernelContact, strain)

    @test modulation(mock_state, TanhMockParams(), household_infection) == true

    hits += modulation(mock_state, TanhMockParams(), constant_infection)
  end

  @test hits == num_samples/2

end

end


text: Dict{String, Any}("a" => Dict{String, Any}())


text:
using DiffEqBase
using Test
using Random

mutable struct TestDiffEqOperator{T} <: DiffEqBase.AbstractDiffEqLinearOperator{T}
    m::Int
    n::Int
end

TestDiffEqOperator(A::AbstractMatrix{T}) where {T} =
    TestDiffEqOperator{T}(size(A)...)

Base.size(A::TestDiffEqOperator) = (A.m, A.n)

A = TestDiffEqOperator([0 0; 0 1])
B = TestDiffEqOperator([0 0 0; 0 1 0; 0 0 2])

@test_throws ErrorException AffineDiffEqOperator{Int64}((A,B),())

@testset "DiffEq linear operators" begin
    Random.seed!(0)
    M = rand(2,2); A = DiffEqArrayOperator(M)
    b = rand(2)
    u = rand(2)
    p = rand(1)
    t = rand()
    As_list = [(A,), (A, A)]#, (A, α)]
    bs_list = [(), (b,), (2b,), (b, 2b)]
    @testset "combinations of A's and b's" for As in As_list, bs in bs_list
        L = AffineDiffEqOperator{Float64}(As, bs, zeros(2))
        mysum = sum(A*u for A in As)
        for b in bs; mysum .+= b; end
        @test L(u,p,t) == mysum
    end
end


text:
@testset "points.jl" begin
    p1 = Point(80.0, 190.0)
    p2 = Point(40.0, -10.0)

    mp = midpoint(p1, p2)
    @test mp.ϕ ≈ 69.58472222222 atol = 0.001
    @test mp.λ ≈ -15.62638888888 atol = 0.001
    @test midpoint(Arc(p1, p2)).ϕ ≈ 69.58472222222 atol = 0.001

    ip = intermediate_point(p1, p2, 0.5)
    @test ip.ϕ ≈ 69.58472222222 atol = 0.001
    @test ip.λ ≈ -15.62638888888 atol = 0.001

    p3 = Point(0.0, -10.0)
    p4 = Point(0.0, 0.0)
    p5 = Point(0.0, 10.0)
    lss = Arcs([p3, p4, p5])
    ip_lss = intermediate_point(lss, 0.6)
    @test ip_lss.ϕ ≈ 0.0 atol = 0.01
    @test ip_lss.λ ≈ 2.0 atol = 0.01
    @test intermediate_point(lss, 0.0).λ ≈ -10.0 atol = 0.01
    @test intermediate_point(lss, 1.0).λ ≈ 10.0 atol = 0.01
    @test midpoint(lss).λ ≈ 0.0 atol = 0.01

    ad = 89.9321605919
    @test SphericalGeometry.normalize(destination_point(p1, ad, 45.0).λ, -180.0, 180.0) ≈ -35.44694444444 atol = 0.001
    @test destination_point(Line(p1, 45.0), ad).ϕ ≈ 7.12027777777 atol = 0.001

    intp = intersection_point(p4, p2, 45.0, 20.0)
    @test intp.ϕ ≈ 27.43166666666 atol = 0.001
    @test intp.λ ≈ 148.73138888888 atol = 0.001
    intp2 = intersection_point(Line(p4, 45.0), Line(p2, 20.0))
    @test intp2.ϕ ≈ 27.43166666666 atol = 0.001
    @test intp2.λ ≈ 148.73138888888 atol = 0.001

    p1ni = Point(50.0, 10.0)
    p2ni = Point(50.0, 30.0)
    p3ni = Point(40.0, 20.0)
    p4ni = Point(60.0, 20.0)

    @test intersection_point(p1ni, p2ni, p3ni, p4ni).ϕ ≈ 50.43250 atol = 0.01
    @test intersection_point(p1ni, p2ni, p3ni, p4ni).λ ≈ 20.0 atol = 0.0001
    @test intersection_point(p1ni, p2ni, p3ni, p2ni).λ ≈ 30.0 atol = 0.0001
    @test intersection_point(Arc(p1ni, p2ni), Arc(p3ni, p2ni)).λ ≈ 30.0 atol = 0.0001

    pa11 = Point(20.0, 10.0)
    pa12 = Point(90.0, 60.0)
    pa21 = Point(10.0, 50.0)
    pa22 = Point(80.0, 5.0)

    @test intersection_point(pa11, pa12, pa21, pa22).λ ≈ 10.0 atol = 0.1
    @test intersection_point(pa11, pa12, pa21, pa22).ϕ ≈ 79.06854 atol = 0.01

    p4ni2 = Point(49.0, 20.0)
    @test isnan(intersection_point(p1ni, p2ni, p3ni, p4ni2).λ)

    # p1a = Point(10.0374, 0.0)
    # p2a = Point(-90.0, 0.0)
    # p3a = Point(0.0, 0.0)
    # @test isinf(intersection_point(p3a, p2a, p1a, p3a).ϕ)

    p1ip1 = Point(5.0, -10.0)
    p2ip1 = Point(5.0, -5.0)
    p3ip1 = Point(5.0, 5.0)
    p4ip1 = Point(5.0, 6.0)
    p5ip1 = Point(5.0, 11.0)
    p6ip1 = Point(-1.0, 5.0)
    p1ip2 = Point(0.0, 0.0)
    p2ip2 = Point(0.0, 10.0)
    p3ip2 = Point(10.0, 10.0)
    p4ip2 = Point(10.0, 0.0)
    p5ip2 = Point(1.0, 0.0)

    @test length(intersection_points(Arcs([p1ip1, p2ip1, p3ip1, p4ip1, p5ip1,
    p6ip1]), Arcs([p1ip2, p2ip2, p3ip2, p4ip2, p5ip2]))) == 4
    @test length(intersection_points(Arcs([p1ip1, p2ip1, p3ip1, p4ip1, p5ip1,
    p6ip1]), Polygon(p1ip2, [p1ip2, p2ip2, p3ip2, p4ip2, p5ip2]))) == 4
    @test length(intersection_points(Polygon(p1ip2, [p1ip1, p2ip1, p3ip1, p4ip1, p5ip1,
    p6ip1]), Polygon(p1ip2, [p1ip2, p2ip2, p3ip2, p4ip2, p5ip2]))) == 6

    #concave polygon
    pc1 = Point(64.0, 27.0)
    pc2 = Point(47.0, 0.0)
    pc3 = Point(42.0, 25.0)
    pc4 = intermediate_point(pc2,pc3,0.33)
    pc5 = Point(52.0, 16.0)
    pc6 = intermediate_point(pc2,pc3,0.67)
    pc7 = Point(55.0, -7.0)
    pc8 = intermediate_point(pc2,pc1,0.33)
    pc9 = Point(62.0, 4.0)
    pc10 = intermediate_point(pc2,pc1,0.67)
    pc11 = Point(56.0, 21.0)

    polygon_concave = Polygon(pc11,
    [pc1,pc3,pc6,pc5,pc4,pc2,pc7,pc8,pc9,pc10])

    @test length(intersection_points(Arc(Point(63.0, 27.0), Point(65.0, 27.0)), polygon_concave)) == 1
    @test length(intersection_points(Arc(Point(63.0, 28.0), Point(65.0, 28.0)), polygon_concave)) == 0
    @test length(intersection_points(Arc(Point(62.0, 4.0), Point(55.0, -7.0)), polygon_concave)) == 2
    @test length(intersection_points(Arc(Point(56.0, 21.0), Point(55.0, 21.0)), polygon_concave)) == 0
    @test length(intersection_points(Arc(Point(64.0, 26.0), Point(47.0, -1.0)), polygon_concave)) == 4
    @test length(intersection_points(Arc(Point(56.0, 21.0), Point(52.0, 16.0)), polygon_concave)) == 1
    @test length(intersection_points(Arc(Point(52.0, 16.0), Point(56.0, 21.0)), polygon_concave)) == 1
    @test length(intersection_points(Arc(Point(56.0, 21.0), Point(56.0, 30.0)), polygon_concave)) == 1

    pt1 = Point(20.0, -5.0)
    pt2 = Point(21.0, 5.0)
    pt3 = Point(22.0, 15.0)
    pa = Point(40.0, 5.0)
    pb = Point(30.0, 5.0)
    pc = Point(10.0, 5.0)
    pd = Point(0.0, 5.0)
    arcs123 = Arcs([pt1, pt2, pt3])

    z = intersection_points(arcs123, Arcs([pt2,pc,pd]))

    @test length(intersection_points(arcs123, Arcs([pt1,pc,pd]))) == 1
    @test length(intersection_points(arcs123, Arcs([pt2,pc,pd]))) == 1
    @test length(intersection_points(arcs123, Arcs([pt3,pc,pd]))) == 1
    @test length(intersection_points(arcs123, Arcs([pb,pt1,pc]))) == 1
    @test length(intersection_points(arcs123, Arcs([pb,pt2,pc]))) == 1
    @test length(intersection_points(arcs123, Arcs([pb,pt3,pc]))) == 1
    @test length(intersection_points(arcs123, Arcs([pa,pb,pt1]))) == 1
    @test length(intersection_points(arcs123, Arcs([pa,pb,pt2]))) == 1
    @test length(intersection_points(arcs123, Arcs([pa,pb,pt3]))) == 1

    @test length(intersection_points(arcs123, Arcs([pb,Point(10.0, 4.0),pd]))) == 1
    @test length(intersection_points(arcs123, Arcs([pa,pb,pc]))) == 1

    @test !ison(p1ip1, Arc(p2ip1, p3ip1))
    @test ison(Point(0.0, 5.0), Arc(p1ip2, p2ip2))
    @test ison(Point(0.0, 5.0), Arcs([p1ip2, p2ip2, p3ip2]))
    @test ison(Point(0.0, 5.0), Polygon(p5ip2, [p1ip2, p2ip2, p3ip2]))

    p1i = Point(0.0, 0.0)
    p2i = Point(0.0, 10.0)
    p3i = Point(10.0, 10.0)
    p4i = Point(5.0, 5.0)
    p5i = Point(10.0, 0.0)
    px_inside = Point(7.0, 8.0)

    p6i_out = Point(15.0, 0.0)
    p7i_out = Point(6.0, 5.0)
    p8i_out = Point(6.0, -2.0)
    p9i_in = Point(6.0, 2.0)
    p10i_in = Point(1.0, 6.0)
    p11i_in = Point(1.0, 1.0)

    poly_inside = Polygon(px_inside, [p1i, p2i, p3i, p4i, p5i])
    @test !isinside(p6i_out, poly_inside)
    @test !isinside(p7i_out, poly_inside)
    @test !isinside(p8i_out, poly_inside)
    @test isinside(p9i_in, poly_inside)
    @test isinside(p10i_in, poly_inside)

    @test isinside(Arc(p9i_in, p10i_in), poly_inside)
    @test !isinside(Arc(p8i_out, p10i_in), poly_inside)
    @test !isinside(Arc(p10i_in, p8i_out), poly_inside)
    @test !isinside(Arc(p6i_out, p8i_out), poly_inside)

    @test isinside(Arcs([p9i_in, p10i_in]), poly_inside)
    @test !isinside(Arcs([p8i_out, p10i_in]), poly_inside)
    @test !isinside(Arcs([p10i_in, p8i_out]), poly_inside)
    @test !isinside(Arcs([p6i_out, p8i_out]), poly_inside)

    poly2 = Polygon(Point(2.0, 2.0), [p9i_in, p11i_in, p10i_in])
    @test isinside(poly2, poly_inside)

    poly3 = Polygon(px_inside, [p1i, p2i, p5i, p3i])
    arcs3 = Arcs([p1i, p3i, p2i, p5i])
    @test !isselfintersecting(poly2)
    @test isselfintersecting(poly3)
    @test isselfintersecting(arcs3)
    @test !isselfintersecting(Arcs([p1i, p2i, p3i, p4i, p5i]))

    # https://www.onboardintelligence.com/GC for generating test (sphere)
    pz1 = Point(50.0, 70.0)
    pz2 = Point(50.0, -70.0)
    azi = azimuth(pz1, pz2)
    @test highest_latitude(pz1, azi) ≈ 73.9871 atol = 0.0001
    @test lowest_latitude(pz1, azi) ≈ -73.9871 atol = 0.0001
    @test highest_latitude_point(pz1, pz2).ϕ ≈ 73.9871 atol = 0.0001
    @test lowest_latitude_point(pz1, pz2).ϕ ≈ 50.0 atol = 0.001
    @test highest_latitude_point(pz1, Point(60.0, 70.0)).ϕ ≈ 60.0 atol = 0.001
    @test highest_latitude_point(Point(60.0, 70.0), pz1).ϕ ≈ 60.0 atol = 0.001
    @test lowest_latitude_point(pz1, Point(60.0, 70.0)).ϕ ≈ 50.0 atol = 0.001
    @test lowest_latitude_point(Point(60.0, 70.0), pz1).ϕ ≈ 50.0 atol = 0.001
    @test highest_latitude_point(pz1, Point(60.0, 80.0)).ϕ ≈ 60.0 atol = 0.0001

    azi2 = azimuth(Point(30.0, 89.0), Point(50.0, -60.0))
    @test highest_latitude_point(Point(30.0, 89.0), azi2).ϕ ≈ 73.2626 atol=0.0001
    @test highest_latitude_point(Point(30.0, 89.0), azi2).λ ≈ 8.9986 atol=0.0001
    @test lowest_latitude_point(Point(30.0, 89.0), azi2).λ ≈ -171.0014 atol=0.0001

    @test highest_latitude_point(pz1, azi).ϕ ≈ 73.9871 atol = 0.0001
    @test lowest_latitude_point(pz1, azi).ϕ ≈ -73.9871 atol = 0.0001
    @test highest_latitude_point(pz1, azi).λ ≈ 0.0 atol = 0.0001

    @test highest_latitude_point(Point(30.0, 90.0), 0.0).ϕ ≈ 90.0 atol = 0.1

end


text:
# Copyright 2021 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import FirstOrderLp
using Plots
using LaTeXStrings

include("primal_dual_hybrid_gradient.jl")
include("utils.jl")

function generate_problem(alpha::Float64, beta::Float64)
  matrix_two_variables = [
    -1.0 -2.0
    -3.0 -1.0
    1.0 1.0
  ]
  matrix_three_variables = [
    -1.0 -2.0 0.0
    -3.0 -1.0 0.0
    1.0 1.0 0.0
  ]
  variable_lower_bounds = alpha == 0.0 ? [-Inf, -Inf] : [-Inf, -Inf, -Inf]
  variable_upper_bounds = alpha == 0.0 ? [Inf, Inf] : [Inf, Inf, Inf]
  objective_vector = alpha == 0.0 ? [1.0, 1.0] : [1.0, 1.0, -1.0]
  constraint_matrix =
    alpha == 0.0 ? matrix_two_variables : matrix_three_variables
  right_hand_side = [-2.0, -2.0, beta]
  return FirstOrderLp.linear_programming_problem(
    variable_lower_bounds,
    variable_upper_bounds,
    objective_vector,
    0.0, # objective_constant
    constraint_matrix,
    right_hand_side,
    0, #num_equalities
  )
end

function generate_component_wise_plots(
  primal_iterates,
  dual_iterates,
  results_dir::String,
  prefix::String,
)
  gr()
  fntsm = font("serif-roman", pointsize = 18)
  fntlg = font("serif-roman", pointsize = 22)
  default(
    titlefont = fntlg,
    guidefont = fntlg,
    tickfont = fntsm,
    legendfont = fntsm,
  )
  colors = [1, 2, 4]
  # Plot primal
  plot(
    [primal_iterates[j][1] for j in 1:length(primal_iterates)],
    label = L"x_0",
    line = (3, :solid),
    color = colors[1],
  )
  plot!(
    [primal_iterates[j][2] for j in 1:length(primal_iterates)],
    label = L"x_1",
    line = (3, :solid),
    color = colors[2],
  )
  if length(primal_iterates[1]) > 2
    plot!(
      [primal_iterates[j][3] for j in 1:length(primal_iterates)],
      label = L"x_2",
      line = (3, :solid),
      color = colors[3],
    )
  end
  xaxis!("Iteration count")
  yaxis!("Value of the variables")
  savefig(joinpath(results_dir, prefix * "p.pdf"))

  # Plot dual
  plot(
    [dual_iterates[j][1] for j in 1:length(dual_iterates)],
    label = L"y_0",
    line = (3, :solid),
    color = colors[1],
  )
  plot!(
    [dual_iterates[j][2] for j in 1:length(dual_iterates)],
    label = L"y_1",
    line = (3, :solid),
    color = colors[2],
  )
  if length(dual_iterates[1]) > 2
    plot!(
      [dual_iterates[j][3] for j in 1:length(dual_iterates)],
      label = L"y_2",
      line = (3, :solid),
      color = colors[3],
    )
  end
  xaxis!("Iteration count")
  yaxis!("Value of the variables")
  savefig(joinpath(results_dir, prefix * "d.pdf"))
end

function main()
  params = PdhgParameters(
    nothing, # step_size (forces the solver to use a provably correct step size)
    1.0, # primal_weight
    1000, # printing frequency
    false, # verbose
    1000, # iteration limit
    nothing, # two_dimensional_subspace
    0.0, # eps_primal_infasible
    0.0, # eps_dual_infeasible
    true, # store_iterates
    nothing, # reference_point
  )
  both_feasible_problem = generate_problem(0.0, 1.0)
  both_infeasible_problem = generate_problem(1.0, 2.0)
  dual_unbounded_problem = generate_problem(0.0, 2.0)
  primal_unbounded_problem = generate_problem(1.0, 1.0)
  output = optimize(params, both_feasible_problem)
  generate_component_wise_plots(
    output.primal_iterates,
    output.dual_iterates,
    ARGS[1],
    "pfdf",
  )
  output = optimize(params, both_infeasible_problem)
  generate_component_wise_plots(
    output.primal_iterates,
    output.dual_iterates,
    ARGS[1],
    "pidi",
  )
  output = optimize(params, dual_unbounded_problem)
  generate_component_wise_plots(
    output.primal_iterates,
    output.dual_iterates,
    ARGS[1],
    "pidf",
  )
  output = optimize(params, primal_unbounded_problem)
  generate_component_wise_plots(
    output.primal_iterates,
    output.dual_iterates,
    ARGS[1],
    "pfdi",
  )
end

main()


text:
# Exercise 2-1
# 5
# x = 5
# x + 1
# The output is empty when running in script.
# With println, the output is 5 followed by error from trying to print x = 5.

# Exercise 2-2
# 1. 42 is not a valid variable name
# 2. Multiple assignment x = y = 1 is valid.
# 3. In the REPL, ; suppresses output
# 4. For integers, it converts it to a float. Used for calling vectorised functions.
# 5. x y is invalid. 5x is valid.

# Exercise 2-3
println("Volume of sphere with radius 5 = $(4/3 * pi * 5^3)")
println("Total cost of 60 books = \$$(round(0.6*24.95 * 60 + 3 + 59 * 0.75, digits=2))")
t1 = 8*60 + 15
t2 = 7*60 + 12
total = (2*t1 + 3*t2) / 60
start = 6*60 + 52
hrs, mins = divrem(start + total, 60)
println("Time to reach home = $(floor(Int, hrs)):$(floor(Int, mins)) am")


text:
using GeometricIntegrators
using Weave

GeometricIntegrators.set_config(:verbosity, 0)

Weave.set_chunk_defaults!(:echo => false, :results => "raw")

include("../src/lotka-volterra-2d-symmetric.jl")

weave("../weave/lotka-volterra-2d-symmetric-vprk-gauss.jmd",
         out_path = "src/lotka-volterra-2d-symmetric",
         doctype = "github",
         mod = LotkaVolterra2dSymmetricSRK)

weave("../weave/lotka-volterra-2d-symmetric-vprk-lobatto-ab.jmd",
         out_path = "src/lotka-volterra-2d-symmetric",
         doctype = "github",
         mod = LotkaVolterra2dSymmetricSRK)

weave("../weave/lotka-volterra-2d-symmetric-vprk-lobatto-ba.jmd",
         out_path = "src/lotka-volterra-2d-symmetric",
         doctype = "github",
         mod = LotkaVolterra2dSymmetricSRK)

weave("../weave/lotka-volterra-2d-symmetric-vprk-radau.jmd",
         out_path = "src/lotka-volterra-2d-symmetric",
         doctype = "github",
         mod = LotkaVolterra2dSymmetricSRK)

weave("../weave/lotka-volterra-2d-symmetric-firk.jmd",
         out_path = "src/lotka-volterra-2d-symmetric",
         doctype = "github",
         mod = LotkaVolterra2dSymmetricSRK)

weave("../weave/lotka-volterra-2d-symmetric-srk.jmd",
         out_path = "src/lotka-volterra-2d-symmetric",
         doctype = "github",
         mod = LotkaVolterra2dSymmetricSRK)


text:
# ---------------------------------------------------------------------------- #
#
#   porder.jl
#
#   Abstract type for polynomial order
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    Porder

Porder abstract type:
Overarching abstract type for polynomial order types.
"""
abstract type Porder

end

### P = 1
struct P1 <: Porder
  p::Int64
end
function P1(  )
  P1(1)
end
function PGdef( P::P1 )
  return PG3()
end

### P = 2
struct P2 <: Porder
  p::Int64
end
function P2(  )
  P2(2)
end
function PGdef( P::P2 )
  return PG6()
end

### P = 3
struct P3 <: Porder
  p::Int64
end
function P3(  )
  P3(3)
end
function PGdef( P::P3 )
  return PG9()
end

### P = 4
struct P4 <: Porder
  p::Int64
end
function P4(  )
  P4(4)
end
function PGdef( P::P4 )
  return PG12()
end

# ### P = 5
# struct P5 <: Porder
#   p::Int64
# end
#
# function P5(  )
#   P5(5)
# end
#
# function PGdef( P::P5 )
#   return PG15()
# end


text:
@info "Testing discriminant_parameters.jl"

### Model Parameters

@testset "set_dimensions!(Θ, y, X, dims)" begin
    nₘ = [400; 500; 600]
    n = sum(nₘ)
    m = length(nₘ)
    p = 10

    for T in (Float32, Float64)
        X, y, M = random_data(T, nₘ, p)
        Xt = copy(transpose(X))

        for (dims, X_test) in [(2,X), (1,Xt)]
            Θ_test = DA.DiscriminantParameters{T}()

            Θ_res = DA.set_dimensions!(Θ_test, y, X_test, dims)

            @test Θ_test === Θ_res
            @test Θ_test.dims == dims
            @test Θ_test.m == m
            @test Θ_test.p == p
        end
    end
end

@testset "set_gamma!(Θ, γ)" begin
    for T in (Float32, Float64)
        Θ_test = DA.DiscriminantParameters{T}()

        for γ in range(zero(T), stop=one(T), length=5)
            Θ_res = DA.set_gamma!(Θ_test, γ)

            @test Θ_test === Θ_res
            @test Θ_test.γ == γ
        end

        for γ in [zero(T)-eps(zero(T)), one(T)+eps(one(T))]
            @test_throws DomainError DA.set_gamma!(Θ_test, γ)
        end
    end
end

@testset "set_statistics!(Θ, y, X, centroids)" begin
    nₘ = [400; 500; 600]
    n = sum(nₘ)
    m = length(nₘ)
    p = 10

    for T in (Float32, Float64)
        X, y, M = random_data(T, nₘ, p)
        Xt = copy(transpose(X))
        Mt = copy(transpose(M))

        Θ_test = DA.DiscriminantParameters{T}()
        Θ_test.m = m
        Θ_test.p = p

        Θ_test.dims = 2
        for M_test in [nothing, perturb(M)]
            Θ_res = DA.set_statistics!(Θ_test, y, X, M_test)

            @test Θ_res === Θ_test
            @test isapprox(Θ_test.M, M_test === nothing ? M : M_test)
            @test Θ_test.nₘ == nₘ
        end

        Θ_test.dims = 1
        for Mt_test in [nothing, perturb(Mt)]
            Θ_res = DA.set_statistics!(Θ_test, y, Xt, Mt_test)

            @test Θ_res === Θ_test
            @test isapprox(Θ_test.M, Mt_test === nothing ? Mt : Mt_test)
            @test Θ_test.nₘ == nₘ
        end
    end
end

@testset "set_priors!(Θ, π)" begin
    nₘ = [400; 500; 600]
    n = sum(nₘ)
    m = length(nₘ)

    for T in (Float32, Float64)
        π = convert(Vector{T},nₘ)/n

        Θ_test = DA.DiscriminantParameters{T}()
        Θ_test.m = m
        Θ_test.nₘ = nₘ

        for π_test in [nothing, ones(T,m)/m]
            Θ_res = DA.set_priors!(Θ_test, π_test)

            @test Θ_test === Θ_res
            @test isapprox(Θ_test.π, π_test === nothing ? π : π_test)
        end
    end
end


text:
# Tests that every option is documented
@testset "Documentation" begin
    @testset "Options" begin
        cwd = pwd()
        dir = pathof(ConstraintSolver)[1:(end - 20)]
        cd(dir)
        cd("../docs/src")
        options_md = readlines("options.md")
        all_options = fieldnames(CS.SolverOptions)
        found_all = true
        for option in all_options
            found = false
            option_str = String(option)
            for line in options_md
                if startswith(line, "## `$option_str`")
                    found = true
                    break
                end
            end
            if !found && option != :no_prune
                found_all = false
                @error "Option $option is not documented"
            end
        end
        @test found_all

        cd(cwd)
    end
end


text:
function test_measures()
    m1 = LebesgueMeasure()
    @test !isdiscrete(m1)
    @test iscontinuous(m1)
    @test !isnormalized(m1)
    @test domaintype(m1) == Float64
    @test DomainIntegrals.unsafe_weight(m1, 0.4) == 1

    m2 = UnitLebesgueMeasure()
    @test !isdiscrete(m2)
    @test iscontinuous(m2)
    @test isnormalized(m2)
    @test domaintype(m2) == Float64
    @test support(m2) == UnitInterval()
    @test weight(m2, 0.4) == 1
    @test weight(m2, 1.4) == 0

    m3 = LegendreMeasure()
    @test !isdiscrete(m3)
    @test iscontinuous(m3)
    @test !isnormalized(m3)
    @test domaintype(m3) == Float64
    @test support(m3) == ChebyshevInterval()
    @test weight(m3, 0.4) == 1
    @test weight(m3, -0.4) == 1
    @test weight(m3, 1.4) == 0
    @test weight(m3, -1.4) == 0

    x = 0.5
    m4 = DiracMeasure(x)
    @test isdiscrete(m4)
    @test !iscontinuous(m4)
    @test isnormalized(m4)
    @test domaintype(m4) == Float64
    @test support(m4) == Point(x)
    @test point(m4) == x
    @test weight(m4, x) == Inf
    @test weight(m4, x+1) == 0
    @test weight(m4, big(x)) == Inf
    @test weight(m4, big(x+1)) == 0

    m5 = GaussianMeasure{SVector{2,Float64}}()
    @test !isdiscrete(m5)
    @test iscontinuous(m5)
    @test isnormalized(m5)
    @test domaintype(m5) == SVector{2,Float64}
    @test support(m5) == FullSpace{SVector{2,Float64}}()
    @test weight(m5, SVector(0.0,0.0)) ≈ 1/(2pi)
    @test weight(m5, SVector(0,0)) ≈ 1/(2pi)
    @test weight(m5, SVector(big(0),big(0))) ≈ 1/(2pi)
    @test weight(m5, SVector(big(0.0),big(0.0))) ≈ 1/(2pi)

    m6 = LaguerreMeasure(0.0)
    @test !isdiscrete(m6)
    @test iscontinuous(m6)
    @test isnormalized(m6)
    @test !isnormalized(LaguerreMeasure(0.1))
    @test domaintype(m6) == Float64
    @test support(m6) == HalfLine()
    @test weight(m6, 0.4) == exp(-0.4)
    @test weight(m6, -0.4) == 0
    @test weight(m6, big(0.4)) == exp(-big(0.4))
end


text:
using DiffEqBase, OrdinaryDiffEq, ApproxFun, DiffEqApproxFun
using Base.Test

S=Fourier()
u0=Fun(θ->cos(cos(θ-0.1))-cos(cos(0-0.1)),S)
c=Fun(cos,S)

prob = ODEProblem((u,p,t)->u''+(c+1)*u',u0,(0.,1.))

@time sol=solve(prob,Euler(),dt=1/1000)

tstops =  [0.0,0.00630546,0.0196241,0.0365809,0.0580536,0.0840865,0.1151,
           0.151529,0.195094,0.24658,0.292971,0.316851,0.338594,0.36075,
           0.384077,0.408548,0.433682,0.458922,0.483926,0.508619,0.533097,
           0.557498,0.581922,0.606414,0.630967,0.655551,0.68014,0.704717,
           0.729278,0.753828,0.778375,0.802923,0.827475,0.85203,0.876586,
           0.901142,0.925698,0.950253,0.974808,0.999362,1.0]

@time sol=solve(prob,Tsit5(),dt=1/40, adaptive=false, tstops=tstops)

#@time sol=solve(prob,Tsit5())
#@time sol=solve(prob,Rosenbrock23(),dt=1/40, adaptive=false, tstops=tstops)


text:
## core stream types ##

# the first argument to any IO MUST be a POINTER (to a JL_STREAM) or using show on it will cause memory corruption

# Generic IO functions

## byte-order mark, ntoh & hton ##

const ENDIAN_BOM = reinterpret(Uint32,uint8([1:4]))[1]

if ENDIAN_BOM == 0x01020304
    ntoh(x) = x
    hton(x) = x
    ltoh(x) = bswap(x)
    htol(x) = bswap(x)
elseif ENDIAN_BOM == 0x04030201
    ntoh(x) = bswap(x)
    hton(x) = bswap(x)
    ltoh(x) = x
    htol(x) = x
else
    error("seriously? what is this machine?")
end

isreadonly(s) = isreadable(s) && !iswritable(s)

## binary I/O ##

# all subtypes should implement this
write(s::IO, x::Uint8) = error(typeof(s)," does not support byte I/O")

write(io::IO, x) = throw(MethodError(write, (io, x)))
write(io::IO, xs...) = for x in xs write(io, x) end

if ENDIAN_BOM == 0x01020304
    function write(s::IO, x::Integer)
        sz = sizeof(x)
        for n = sz:-1:1
            write(s, uint8((x>>>((n-1)<<3))))
        end
        sz
    end
else
    function write(s::IO, x::Integer)
        sz = sizeof(x)
        for n = 1:sz
            write(s, uint8((x>>>((n-1)<<3))))
        end
        sz
    end
end

write(s::IO, x::Bool)    = write(s, uint8(x))
write(s::IO, x::Float16) = write(s, reinterpret(Int16,x))
write(s::IO, x::Float32) = write(s, reinterpret(Int32,x))
write(s::IO, x::Float64) = write(s, reinterpret(Int64,x))

function write(s::IO, a::AbstractArray)
    nb = 0
    for i = 1:length(a)
        nb += write(s, a[i])
    end
    nb
end

function write(s::IO, c::Char)
    if c < 0x80
        write(s, uint8(c))
        return 1
    elseif c < 0x800
        write(s, uint8(( c >> 6          ) | 0xC0))
        write(s, uint8(( c        & 0x3F ) | 0x80))
        return 2
    elseif c < 0x10000
        write(s, uint8(( c >> 12         ) | 0xE0))
        write(s, uint8(((c >> 6)  & 0x3F ) | 0x80))
        write(s, uint8(( c        & 0x3F ) | 0x80))
        return 3
    elseif c < 0x110000
        write(s, uint8(( c >> 18         ) | 0xF0))
        write(s, uint8(((c >> 12) & 0x3F ) | 0x80))
        write(s, uint8(((c >> 6)  & 0x3F ) | 0x80))
        write(s, uint8(( c        & 0x3F ) | 0x80))
        return 4
    else
        return write(s, '\ufffd')
    end
end

function write(s::IO, p::Ptr, n::Integer)
    for i=1:n
        write(s, unsafe_load(p, i))
    end
    n
end

function write(io::IO, s::Symbol)
    pname = convert(Ptr{Uint8}, s)
    write(io, pname, int(ccall(:strlen, Csize_t, (Ptr{Uint8},), pname)))
end

# all subtypes should implement this
read(s::IO, x::Type{Uint8}) = error(typeof(s)," does not support byte I/O")

function read{T <: Integer}(s::IO, ::Type{T})
    x = zero(T)
    for n = 1:sizeof(x)
        x |= (convert(T,read(s,Uint8))<<((n-1)<<3))
    end
    return x
end

read(s::IO, ::Type{Bool})    = (read(s,Uint8)!=0)
read(s::IO, ::Type{Float16}) = box(Float16,unbox(Int16,read(s,Int16)))
read(s::IO, ::Type{Float32}) = box(Float32,unbox(Int32,read(s,Int32)))
read(s::IO, ::Type{Float64}) = box(Float64,unbox(Int64,read(s,Int64)))

read{T}(s::IO, t::Type{T}, d1::Int, dims::Int...) =
    read(s, t, tuple(d1,dims...))
read{T}(s::IO, t::Type{T}, d1::Integer, dims::Integer...) =
    read(s, t, map(int,tuple(d1,dims...)))

read{T}(s::IO, ::Type{T}, dims::Dims) = read!(s, Array(T, dims))

function read!{T}(s::IO, a::Array{T})
    for i = 1:length(a)
        a[i] = read(s, T)
    end
    return a
end

function read(s::IO, ::Type{Char})
    ch = read(s, Uint8)
    if ch < 0x80
        return char(ch)
    end

    # mimic utf8.next function
    trailing = Base.utf8_trailing[ch+1]
    c::Uint32 = 0
    for j = 1:trailing
        c += ch
        c <<= 6
        ch = read(s, Uint8)
    end
    c += ch
    c -= Base.utf8_offset[trailing+1]
    char(c)
end

function readuntil(s::IO, delim::Char)
    if delim < 0x80
        data = readuntil(s, uint8(delim))
        enc = byte_string_classify(data)
        return (enc==1) ? ASCIIString(data) : UTF8String(data)
    end
    out = IOBuffer()
    while !eof(s)
        c = read(s, Char)
        write(out, c)
        if c == delim
            break
        end
    end
    takebuf_string(out)
end

function readuntil{T}(s::IO, delim::T)
    out = T[]
    while !eof(s)
        c = read(s, T)
        push!(out, c)
        if c == delim
            break
        end
    end
    out
end

# based on code by Glen Hertz
function readuntil(s::IO, t::String)
    l = length(t)
    if l == 0
        return ""
    end
    if l > 40
        warn("readuntil(IO,String) will perform poorly with a long string")
    end
    out = IOBuffer()
    m = Array(Char, l)  # last part of stream to match
    t = collect(t)
    i = 0
    while !eof(s)
        i += 1
        c = read(s, Char)
        write(out, c)
        if i <= l
            m[i] = c
        else
            # shift to last part of s
            for j = 2:l
                m[j-1] = m[j]
            end
            m[l] = c
        end
        if i >= l && m == t
            break
        end
    end
    return takebuf_string(out)
end

readline(s::IO) = readuntil(s, '\n')
readchomp(x) = chomp!(readall(x))

# read up to nb bytes into nb, returning # bytes read
function readbytes!(s::IO, b::AbstractArray{Uint8}, nb=length(b))
    olb = lb = length(b)
    nr = 0
    while nr < nb && !eof(s)
        a = read(s, Uint8)
        nr += 1
        if nr > lb
            lb = nr * 2
            resize!(b, lb)
        end
        b[nr] = a
    end
    if lb > olb
        resize!(b, nr) # shrink to just contain input data if was resized
    end
    return nr
end

# read up to nb bytes from s, returning a Vector{Uint8} of bytes read.
function readbytes(s::IO, nb=typemax(Int))
    b = Array(Uint8, nb == typemax(Int) ? 1024 : nb)
    nr = readbytes!(s, b, nb)
    resize!(b, nr)
end

function readall(s::IO)
    b = readbytes(s)
    return is_valid_ascii(b) ? ASCIIString(b) : UTF8String(b)
end
readall(filename::String) = open(readall, filename)

## high-level iterator interfaces ##

type EachLine
    stream::IO
    ondone::Function
    EachLine(stream) = EachLine(stream, ()->nothing)
    EachLine(stream, ondone) = new(stream, ondone)
end
eachline(stream::IO) = EachLine(stream)

start(itr::EachLine) = nothing
function done(itr::EachLine, nada)
    if !eof(itr.stream)
        return false
    end
    itr.ondone()
    true
end
next(itr::EachLine, nada) = (readline(itr.stream), nothing)
eltype(itr::EachLine) = ByteString

readlines(s=STDIN) = collect(eachline(s))

# IOStream Marking

# Note that these functions expect that io.mark exists for
# the concrete IO type.  This may not be true for IO types
# not in base.

function mark(io::IO)
    io.mark = position(io)
end

function unmark(io::IO)
    !ismarked(io) && return false
    io.mark = -1
    return true
end

function reset(io::IO)
    !ismarked(io) && error(io, " not marked")
    m = io.mark
    seek(io, m)
    io.mark = -1 # must be after seek, or seek may fail
    return m
end

ismarked(io::IO) = io.mark >= 0

max_stars_repo_name: Aerlinger/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
#jl #! format: off
#[PSSE 240 Bus Case system with Renewables](https://www.nrel.gov/grid/test-case-repository.html)

# **Originally Contributed by**: José Daniel Lara

# # Introduction

# This tutorial will introduce the industry models of Renewable Energy the comparisons between
# DiffEq Integration techniques for comparison. We show the uses of Sundials and OrdinaryDiffEq
# to obtain the transient response of a system to a perturbation.

using PowerSimulationsDynamics
using PowerSystemCaseBuilder
using PowerSystems
using Sundials
using Plots
using OrdinaryDiffEq
PSD = PowerSimulationsDynamics

# ## Load the system
# The system data is provided through PowerSystemCaseBuilder.
sys = build_system(PSSETestSystems, "psse_240_case_renewable_sys")

# ## Build the simulation and initialize the problem

# The next step is to create the simulation structure. This will create the indexing
# of our system that will be used to formulate the differential-algebraic system of
# equations. To do so, it is required to specify the perturbation that will occur in
# the system. In this case, we will use a ResidualModel formulation, for more details
# about the formulation checkout the [Models Section](https://nrel-siip.github.io/PowerSimulationsDynamics.jl/stable/models/)
# in `PowerSimulationsDynamics.jl` documentation

using Logging
sim_ida = PSD.Simulation(
    ResidualModel,
    sys, #system
    pwd(),
    (0.0, 20.0), #time span
    BranchTrip(1.0, Line, "CORONADO    -1101-PALOVRDE    -1401-i_10");
    console_level = Logging.Info,
)

# ## Run the simulation using Sundials
# We will now run the simulation using Sundials.jl solver IDA() by specifying the maximum
# dt we want for the simulation. In our experience with this solver, solution times are faster
# when supplying information about the maximum time step than the tolerances as we can see in the example
PSD.execute!(sim_ida, IDA(), dtmax = 0.01)

# ## Read the results and plot a system variable
# After the simulation is completed, we can extract the results and make plots as desired.
# In this case, we will plot the voltage magnitude at the bus at which the line was connected.
res_ida = read_results(sim_ida)
v1101_ida = get_voltage_magnitude_series(res_ida, 1101);
plot(v1101_ida)

# ## Run the simulation using Rodas4()
# In this case, we will use a MassMatrixModel formulation, for more details
# about the formulation checkout the [Models Section](https://nrel-siip.github.io/PowerSimulationsDynamics.jl/stable/models/)
# in `PowerSimulationsDynamics.jl` documentation

sim_rodas = PSD.Simulation(
    MassMatrixModel,
    sys, #system
    pwd(),
    (0.0, 20.0), #time span
    BranchTrip(1.0, Line, "CORONADO    -1101-PALOVRDE    -1401-i_10");
    console_level = Logging.Info,
)

# We will now run the simulation using OrdinaryDiffEq.jl solver Rodas4() by specifying the
# tolerance we want for the simulation. In our experience with this solver, solution times are faster
# when supplying information about the atol and rtol values as we can see in the example. The solver will also
# work with a specified dtmax but take a significantly longer time to solve.
# When using OrdinaryDiffEq.jl solvers always pass the option `initializealg = NoInit()` to avoid
# unnecessary re-initialization of the algebraic equations.
PSD.execute!(
    sim_rodas,
    Rodas4(),
    saveat = 0.01,
    atol = 1e-10,
    rtol = 1e-10,
    initializealg = NoInit(),
)

# ## Read the results
# After the simulation is completed, we can extract the results and make plots as desired.
# In this case, we will plot the voltage magnitude at the bus at which the line was connected.
res_rodas = read_results(sim_rodas)

# ## Compare the results
# After the simulation is completed, we can extract the results and make plots as desired.
# In this case, we will plot the voltage magnitude at the bus at which the line was connected.
# For both of the solution techniques.
v1101 = get_voltage_magnitude_series(res_rodas, 1101);
plot(v1101, label = "RODAS4")
plot!(v1101_ida, label = "IDA")


text:
export PointMassConstraint

"""
Description:

    Constraints the marginal of the connected variable to a point-mass.
    Implementation according to (Senoz et al. "Variational Message Passing
    and Local Constraint Manipulation in Factor Graphs", Entropy, 2021).

Interfaces:

    1. out

Construction:

    PointMassConstraint(out; id=:my_node)
"""
mutable struct PointMassConstraint <: DeltaFactor
    id::Symbol
    interfaces::Array{Interface,1}
    i::Dict{Symbol, Interface}

    function PointMassConstraint(out; id=ForneyLab.generateId(PointMassConstraint))
        @ensureVariables(out)
        self = new(id, Vector{Interface}(undef, 1), Dict{Symbol,Interface}())
        ForneyLab.addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)

        return self
    end
end

slug(::Type{PointMassConstraint}) = "δ"

# A breaker message is required if interface is partnered with a point-mass constraint
requiresBreaker(interface::Interface, partner_interface::Interface, partner_node::PointMassConstraint) = true

breakerParameters(interface::Interface, partner_interface::Interface, partner_node::PointMassConstraint) = (Message{GaussianMeanVariance, Univariate}, ()) # Univariate only

isDeltaConstraint(::PointMassConstraint) = true


text:
# # Custom halftoning with SDFs
# Ordered dithering & halftoning algorithms typically use predefined threshold matrices
# with unique values ranging from `1//n` to `(n-1)/n`, where `n` is the number of elements
# in the matrix.
#
# In this example, we will be generating these threshold matrices
# from signed distance functions (SDFs). A good resource for those is Inigo Quilez'
# blogpost ["2D distance functions"](https://www.iquilezles.org/www/articles/distfunctions2d/distfunctions2d.htm).
#
# Here we reimplement his function `sdStar5`:
using Images
using DitherPunk
using ColorSchemes
using TestImages

function star_sdf(x, y; r=0.5, rf=2.0)
    k1 = [0.809016994375, -0.587785252292]
    k2 = [-k1[1], k1[2]]
    p = [abs(x), y]
    p -= 2.0 * maximum([k1 ⋅ p, 0.0]) * k1
    p -= 2.0 * maximum([k2 ⋅ p, 0.0]) * k2
    p = [abs(p[1]), p[2] - r]
    ba = rf * [-k1[2], k1[1]] - [0, 1]
    h = clamp((p ⋅ ba) / (ba ⋅ ba), 0.0, r)
    return norm(p - ba * h) * sign(p[2] * ba[1] - p[1] * ba[2])
end;

# This SDF can be visualized by its contour:
#
# ![](https://user-images.githubusercontent.com/20258504/144755939-19b13e29-e50b-4744-bdeb-234f33790ca8.png)
#
# To convert an SDF to an `OrderedDither` algorithm, use this function:
function sdf2halftone(sdf, n)
    rg = range(-1, 1; length=n)
    A = [sdf(x, y) for y in rg, x in rg]
    p = sortperm(reshape(-A, :))
    B = Vector{Int}(undef, n^2)
    B[p] .= 1:(n^2)
    return OrderedDither(reshape(B, size(A)...)//(n^2 + 1))
end;

# Let's try it on a test image:
img = testimage("fabio_gray_512")
alg = sdf2halftone(star_sdf, 7)
dither(img, alg)

# We can also change the size of the matrix.
# Large halftoning matrices look great on high resolution images.
alg = sdf2halftone(star_sdf, 15)
dither(img, alg)


text:
## Load transfer scheme
#==========================================================================================#

# Sum adjacent values
adjacent_adder(x1, x2) = @views [ [ x1[1] ]; x1[2:end] .+ x2[1:end-1]; [ x2[end] ] ]

# Compute moments for each section with local beam nodes as origins
section_moment(vlm_ac, fem_pts, half_vlm_force) = @. (vlm_ac - fem_pts) × half_vlm_force
section_moments(vlm_acs, fem_pts, half_vlm_forces) = sum(x -> section_moment(x[1], fem_pts, x[2]), zip(eachrow(vlm_acs), eachrow(half_vlm_forces)))

function compute_loads(vlm_acs, vlm_forces, fem_mesh)
    # Forces
    sec_forces   = vec(sum(vlm_forces, dims = 1)) / 2
    beam_forces  = adjacent_adder(sec_forces / 2, sec_forces / 2)

    # Moments
    M_ins        = @views section_moments(vlm_acs, fem_mesh[1:end-1], vlm_forces / 2)
    M_outs       = @views section_moments(vlm_acs, fem_mesh[2:end],   vlm_forces / 2)
    beam_moments = adjacent_adder(M_ins, M_outs)

    # Concatenate forces and moments into loads array
    [ combinedimsview(beam_forces); combinedimsview(beam_moments) ]
end

# Generate load vector for FEM system
fem_load_vector(vlm_acs, vlm_forces, fem_mesh) = [ zeros(6); vec(compute_loads(vlm_acs, vlm_forces, fem_mesh)) ]

## Displacement transfer scheme
#==========================================================================================#

# Build cross product as an antisymmetric bilinear form.
# (This is just a fancy way of saying antisymmetric matrix in 3 dimensions. The cross product, more generally, is actually the exterior product in higher dimensions.)
rotation_matrix(Ωx, Ωy, Ωz) = @SMatrix [  0  -Ωz  Ωy ;
                                          Ωz  0  -Ωx ;
                                         -Ωy  Ωx  0  ]

rotation_matrix(θs) = rotation_matrix.(θs[1,:], θs[2,:], θs[3,:])

# Transfer states by summing the displacements including rotations.
transfer_displacement(xyz, dx, rot, r) = xyz + dx + rot * (xyz - r)
transfer_displacements(dxs, Ts, chord_mesh, fem_mesh) = permutedims(combinedimsview(map(xyz -> transfer_displacement.(xyz, dxs, Ts, fem_mesh), eachrow(chord_mesh))))

mesh_translation(δs) = @views SVector.(δs[1,:], δs[2,:], δs[3,:])
mesh_rotation(δs)    = @views rotation_matrix(δs[4:6,:])

# Make new horseshoes
function new_horseshoes(dxs, Ts, chord_mesh, camber_mesh, fem_mesh)
    new_chord_mesh = transfer_displacements(dxs, Ts, chord_mesh, fem_mesh)
    new_camber_mesh = transfer_displacements(dxs, Ts, camber_mesh, fem_mesh)
    Horseshoe.(make_panels(new_chord_mesh), normal_vector.(make_panels(new_camber_mesh)))
end


text:
using Mimi
using Test
using DataFrames
using VegaLite
using Electron

import Mimi:
    dataframe_or_scalar, _spec_for_item, menu_item_list, getdataframe, dim_names

@defcomp MyComp begin
    a = Parameter(index=[time, regions])
    b = Parameter(index=[time])
    c = Parameter(index=[regions])
    d = Parameter()
    e = Parameter(index=[four])
    f::Array{Float64, 2} = Parameter()

    x = Variable(index=[time, regions])

    function run_timestep(p, v, d, t)
        for r in d.regions
            v.x[t, r] = 0
        end
    end
end

m = Model()
set_dimension!(m, :time, 2000:2100)
set_dimension!(m, :regions, 3)
set_dimension!(m, :four, 4)

add_comp!(m, MyComp)
set_param!(m, :MyComp, :a, ones(101,3))
set_param!(m, :MyComp, :b, 1:101)
set_param!(m, :MyComp, :c, [4,5,6])
set_param!(m, :MyComp, :d, .5)
set_param!(m, :MyComp, :e, [1,2,3,4])
set_param!(m, :MyComp, :f, [1.0 2.0; 3.0 4.0])

run(m)

# 1.  dataframe helper functions
@test typeof(dataframe_or_scalar(m, :MyComp, :a)) == DataFrame
@test typeof(dataframe_or_scalar(m, :MyComp, :d)) == Float64

#2.  Specs and menu
items = [:a, :b, :c, :d, :e, :f, :x]
for item in items
    static_spec = _spec_for_item(m, :MyComp, item; interactive = false)
    interactive_spec = _spec_for_item(m, :MyComp, item)
    if length(dim_names(m, :MyComp, item)) == 0
        name =  string(:MyComp, " : ", item, " = ", m[:MyComp, item])
    else
        name = string(:MyComp, " : ", item)
    end
    @test static_spec["name"] == interactive_spec["name"] == name
end

s = menu_item_list(m)
@test typeof(s) == Array{Any, 1}
@test length(s) == 7

#3.  explore(m::Model, title = "Electron")
w = explore(m, title = "Testing Window")
@test typeof(w) == Electron.Window
close(w)

#4.  Mim.plot(m::Model, comp_name::Symbol, datum_name::Symbol;
#       dim_name::Union{Nothing, Symbol} = nothing)
items = [:a, :b, :c, :d, :e, :f, :x]
for item in items
    p = Mimi.plot(m, :MyComp, item)
    @test typeof(p) == VegaLite.VLSpec{:plot}
end

#6.  errors and warnings
@defcomp MyComp2 begin

    a = Parameter(index = [time, regions, four])
    x = Variable(index=[time, regions])

    function run_timestep(p, v, d, t)
        for r in d.regions
            v.x[t, r] = rand(10)[1]
        end
    end
end

m2 = Model()
set_dimension!(m2, :time, 2000:2100)
set_dimension!(m2, :regions, 3)
set_dimension!(m2, :four, 4)

add_comp!(m2, MyComp2)
set_param!(m2, :MyComp2, :a, ones(101, 3, 4))

run(m2)

#spec creation for MyComp.a should warn because over 2 indexed dimensions
@test_logs (:warn, "MyComp2.a has > 2 indexed dimensions, not yet implemented in explorer") explore(m2)
@test_logs (:warn, "MyComp2.a has > 2 indexed dimensions, not yet implemented in explorer") _spec_for_item(m2, :MyComp2, :a)