text:
"""
    Dims_create!(nnodes::Integer, [ndims::Integer], dims)

Create a division of `nnodes` processes in a Cartesian grid.

# External links
$(_doc_external("MPI_Dims_create"))
"""
function Dims_create!(nnodes::Integer, ndims::Integer, dims::MPIBuffertype{T}) where {T <: Integer}
    # int MPI_Dims_create(int nnodes, int ndims, int dims[])
    @mpichk ccall((:MPI_Dims_create, libmpi), Cint,
                  (Cint, Cint, Ptr{Cint}), nnodes, ndims, dims)
end

function Dims_create!(nnodes::Integer, dims::AbstractArray{T,N}) where {T<:Integer, N}
    cdims = Cint.(dims[:])
    ndims = length(cdims)
    Dims_create!(nnodes, ndims, cdims)
    dims[:] .= cdims
end

"""
    comm_cart = Cart_create(comm_old::Comm, [ndims::Integer], dims, periods, reorder)

Create new MPI communicator with Cartesian structure from an existent communicator.

# External links
$(_doc_external("MPI_Cart_create"))
"""
function Cart_create(comm_old::Comm, ndims::Integer, dims::MPIBuffertype{Cint}, periods::MPIBuffertype{Cint}, reorder)
    comm_cart = Comm()
    # int MPI_Cart_create(MPI_Comm comm_old, int ndims, const int dims[],
    #                     const int periods[], int reorder, MPI_Comm *comm_cart)
    @mpichk ccall((:MPI_Cart_create, libmpi), Cint,
                  (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{MPI_Comm}),
                  comm_old, ndims, dims, periods, reorder, comm_cart)
    if comm_cart != COMM_NULL
        finalizer(free, comm_cart)
    end
    comm_cart
end

function Cart_create(comm_old::Comm, dims::AbstractArray{T,N}, periods::Array{T,N}, reorder) where T <: Integer where N
    cdims    = Cint.(dims[:])
    cperiods = Cint.(periods[:])
    ndims    = length(cdims)
    Cart_create(comm_old, ndims, cdims, cperiods, reorder)
end

"""
    rank = Cart_rank(comm::Comm, coords)

Determine process rank in communicator `comm` with Cartesian structure.
The `coords` array specifies the Cartesian coordinates of the process.

# External links
$(_doc_external("MPI_Cart_rank"))
"""
function Cart_rank(comm::Comm, coords::MPIBuffertype{Cint})
    rank = Ref{Cint}()
    # int MPI_Cart_rank(MPI_Comm comm, const int coords[], int *rank)
    @mpichk ccall((:MPI_Cart_rank, libmpi), Cint,
                  (MPI_Comm, Ptr{Cint}, Ptr{Cint}),
                  comm, coords, rank)
    Int(rank[])
end

function Cart_rank(comm::Comm, coords::AbstractArray{T}) where {T <: Integer}
    ccoords = Cint.(coords[:])
    Cart_rank(comm, ccoords)
end

"""
    dims, periods, coords = Cart_get(comm::Comm)

Obtain information on the Cartesian topology of dimension `N` underlying the
communicator `comm`. This is specified by two `Cint` arrays of `N` elements
for the number of processes and periodicity properties along each Cartesian dimension.
A third `Cint` array is returned, containing the Cartesian coordinates of the calling process.

# External links
$(_doc_external("MPI_Cart_get"))
"""
function Cart_get(comm::Comm)
    maxdims = Cartdim_get(comm)
    # preallocate with nontrivial values
    dims    = Cint[-1 for i = 1:maxdims]
    periods = Cint[-1 for i = 1:maxdims]
    coords  = Cint[-1 for i = 1:maxdims]
    # int MPI_Cart_get(MPI_Comm comm, int maxdims, int dims[], int periods[], int coords[])
    @mpichk ccall((:MPI_Cart_get, libmpi), Cint,
                  (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
                  comm, maxdims, dims, periods, coords)
    return dims, periods, coords
end

"""
    ndims = Cartdim_get(comm::Comm)

Return number of dimensions of the Cartesian topology associated with the communicator `comm`.

# External links
$(_doc_external("MPI_Cartdim_get"))
"""
function Cartdim_get(comm::Comm)
    dims = Ref{Cint}()
    # int MPI_Cartdim_get(MPI_Comm comm, int *ndims)
    @mpichk ccall((:MPI_Cartdim_get, libmpi), Cint,
                  (MPI_Comm, Ptr{Cint}),
                  comm, dims)
    return Int(dims[])
end

"""
    Cart_coords!(comm::Comm, rank::Integer, coords)

Determine coordinates of a given process in Cartesian communicator.

# External links
$(_doc_external("MPI_Cart_coords"))
"""
function Cart_coords!(comm::Comm, rank::Integer, coords::MPIBuffertype{Cint})
    maxdims = Cartdim_get(comm)
    # int MPI_Cart_coords(MPI_Comm comm, int rank, int maxdims, int coords[])
    @mpichk ccall((:MPI_Cart_coords, libmpi), Cint,
                  (MPI_Comm, Cint, Cint, Ptr{Cint}),
                  comm, rank, maxdims, coords)
end

"""
    Cart_coords(comm::Comm)

Determine coordinates of local process in Cartesian communicator.

See also [`Cart_coords!`](@ref).
"""
function Cart_coords(comm::Comm)
    maxdims = Cartdim_get(comm)
    ccoords = Vector{Cint}(undef, maxdims)
    rank    = Comm_rank(comm)
    Cart_coords!(comm, rank, ccoords)
    Int.(ccoords)
end

"""
    rank_source, rank_dest = Cart_shift(comm::Comm, direction::Integer, disp::Integer)

Return the source and destination ranks associated to a shift along a given
direction.

# External links
$(_doc_external("MPI_Cart_shift"))
"""
function Cart_shift(comm::Comm, direction::Integer, disp::Integer)
    rank_source = Ref{Cint}()
    rank_dest   = Ref{Cint}()
    # int MPI_Cart_shift(MPI_Comm comm, int direction, int disp,
    #                    int *rank_source, int *rank_dest)
    @mpichk ccall((:MPI_Cart_shift, libmpi), Cint,
                  (MPI_Comm, Cint, Cint, Ptr{Cint}, Ptr{Cint}),
                  comm, direction, disp, rank_source, rank_dest)
    Int(rank_source[]), Int(rank_dest[])
end

"""
    comm_sub = Cart_sub(comm::Comm, remain_dims)

Create lower-dimensional Cartesian communicator from existent Cartesian
topology.

`remain_dims` should be a boolean vector specifying the dimensions that should
be kept in the generated subgrid.

# External links
$(_doc_external("MPI_Cart_sub"))
"""
function Cart_sub(comm::Comm, remain_dims::MPIBuffertype{Cint})
    comm_sub = Comm()
    # int MPI_Cart_sub(MPI_Comm comm, const int remain_dims[], MPI_Comm *comm_new)
    @mpichk ccall((:MPI_Cart_sub, libmpi), Cint,
                  (MPI_Comm, Ptr{Cint}, Ptr{MPI_Comm}),
                  comm, remain_dims, comm_sub)
    if comm_sub != COMM_NULL
        finalizer(free, comm_sub)
    end
    comm_sub
end

function Cart_sub(comm::Comm, remain_dims)
    cremain_dims = [Cint(dim) for dim in remain_dims]
    Cart_sub(comm, cremain_dims)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# See https://github.com/jonniedie/ConcreteStructs.jl/issues/4

# %%
using ConcreteStructs

Base.@kwdef @concrete struct Problem g; y0; v0; tspan end
Problem(;
    g = 9.80665,
    y0 = 0.0,
    v0 = 30.0,
    tspan = (0.0, 8.0)
) = Problem(g, y0, v0, tspan)

prob1 = Problem()

# %%
prob2 = Problem(v0 = 100.0)

# %%
using MonteCarloMeasurements
prob3 = Problem(y0 = 0.0 ± 0.0, v0 = 30.0 ± 1.0)

# %%

ext: jl

max_stars_repo_path: 0014/Base.@kwdef @concrete.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0014/Base.@kwdef @concrete.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0014/Base.@kwdef @concrete.jl

max_forks_repo_name: genkuroki/public


text:
# helper function
check_dim_mul(A, B) = size(A, 2) == size(B, 1)

struct CompositeMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
    maps::As # stored in order of application to vector
    function CompositeMap{T, As}(maps::As) where {T, As}
        N = length(maps)
        for n in 2:N
            check_dim_mul(maps[n], maps[n-1]) || throw(DimensionMismatch("CompositeMap"))
        end
        for n in 1:N
            promote_type(T, eltype(maps[n])) == T || throw(InexactError())
        end
        new{T, As}(maps)
    end
end
CompositeMap{T}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} = CompositeMap{T, As}(maps)

# basic methods
Base.size(A::CompositeMap) = (size(A.maps[end], 1), size(A.maps[1], 2))
Base.isreal(A::CompositeMap) = all(isreal, A.maps) # sufficient but not necessary

# the following rules are sufficient but not necessary
const RealOrComplex = Union{Real,Complex}
for (f, _f, g) in ((:issymmetric, :_issymmetric, :transpose),
                    (:ishermitian, :_ishermitian, :adjoint),
                    (:isposdef, :_isposdef, :adjoint))
    @eval begin
        LinearAlgebra.$f(A::CompositeMap) = $_f(A.maps)
        $_f(maps::Tuple{}) = true
        $_f(maps::Tuple{<:LinearMap}) = $f(maps[1])
        function $_f(maps::Tuple{Vararg{<:LinearMap}}) # length(maps) >= 2
            if maps[1] isa UniformScalingMap{<:RealOrComplex}
                return $f(maps[1]) && $_f(Base.tail(maps))
            elseif maps[end] isa UniformScalingMap{<:RealOrComplex}
                return $f(maps[end]) && $_f(Base.front(maps))
            else
                return maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps)))
            end
        end
    end
end

# scalar multiplication and division
function Base.:(*)(α::Number, A::LinearMap)
    T = promote_type(typeof(α), eltype(A))
    return CompositeMap{T}(tuple(A, UniformScalingMap(α, size(A, 1))))
end
function Base.:(*)(α::Number, A::CompositeMap)
    T = promote_type(typeof(α), eltype(A))
    Alast = last(A.maps)
    if Alast isa UniformScalingMap
        return CompositeMap{T}(tuple(Base.front(A.maps)..., α * Alast))
    else
        return CompositeMap{T}(tuple(A.maps..., UniformScalingMap(α, size(A, 1))))
    end
end
function Base.:(*)(A::LinearMap, α::Number)
    T = promote_type(typeof(α), eltype(A))
    return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A))
end
function Base.:(*)(A::CompositeMap, α::Number)
    T = promote_type(typeof(α), eltype(A))
    Afirst = first(A.maps)
    if Afirst isa UniformScalingMap
        return CompositeMap{T}(tuple(Afirst * α, Base.tail(A.maps)...))
    else
        return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A.maps...))
    end
end
Base.:(\)(α::Number, A::LinearMap) = inv(α) * A
Base.:(/)(A::LinearMap, α::Number) = A * inv(α)
Base.:(-)(A::LinearMap) = -1 * A

# composition of linear maps
"""
    A::LinearMap * B::LinearMap

Construct a `CompositeMap <: LinearMap`, a (lazy) representation of the product
of the two operators. Products of `LinearMap`/`CompositeMap` objects and
`LinearMap`/`CompositeMap` objects are reduced to a single `CompositeMap`.
In products of `LinearMap`s and `AbstractMatrix`/`UniformScaling` objects, the
latter get promoted to `LinearMap`s automatically.

# Examples
```jldoctest; setup=(using LinearAlgebra, LinearMaps)
julia> CS = LinearMap{Int}(cumsum, 3)::LinearMaps.FunctionMap;

julia> LinearMap(ones(Int, 3, 3)) * CS * I * rand(3, 3);
```
"""
function Base.:(*)(A₁::LinearMap, A₂::LinearMap)
    check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*"))
    T = promote_type(eltype(A₁), eltype(A₂))
    return CompositeMap{T}(tuple(A₂, A₁))
end
function Base.:(*)(A₁::LinearMap, A₂::CompositeMap)
    check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*"))
    T = promote_type(eltype(A₁), eltype(A₂))
    return CompositeMap{T}(tuple(A₂.maps..., A₁))
end
function Base.:(*)(A₁::CompositeMap, A₂::LinearMap)
    check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*"))
    T = promote_type(eltype(A₁), eltype(A₂))
    return CompositeMap{T}(tuple(A₂, A₁.maps...))
end
function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap)
    check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*"))
    T = promote_type(eltype(A₁), eltype(A₂))
    return CompositeMap{T}(tuple(A₂.maps..., A₁.maps...))
end
Base.:(*)(A₁::LinearMap, A₂::UniformScaling) = A₁ * A₂.λ
Base.:(*)(A₁::UniformScaling, A₂::LinearMap) = A₁.λ * A₂

# special transposition behavior
LinearAlgebra.transpose(A::CompositeMap{T}) where {T} = CompositeMap{T}(map(transpose, reverse(A.maps)))
LinearAlgebra.adjoint(A::CompositeMap{T}) where {T}   = CompositeMap{T}(map(adjoint, reverse(A.maps)))

# comparison of LinearCombination objects
Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps)

# multiplication with vectors
function A_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector)
    # no size checking, will be done by individual maps
    N = length(A.maps)
    if N==1
        A_mul_B!(y, A.maps[1], x)
    else
        T = eltype(y)
        dest = Array{T}(undef, size(A.maps[1], 1))
        A_mul_B!(dest, A.maps[1], x)
        source = dest
        if N>2
            dest = Array{T}(undef, size(A.maps[2], 1))
        end
        for n in 2:N-1
            try
                resize!(dest, size(A.maps[n], 1))
            catch err
                if err == ErrorException("cannot resize array with shared data")
                    dest = Array{T}(undef, size(A.maps[n], 1))
                else
                    rethrow(err)
                end
            end
            A_mul_B!(dest, A.maps[n], source)
            dest, source = source, dest # alternate dest and source
        end
        A_mul_B!(y, A.maps[N], source)
    end
    return y
end

At_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)

Ac_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)


text:
@testset "Plotting" begin
  Random.seed!(123)
  z₁ = randn(10000)
  z₂ = z₁ + randn(10000)
  d = georef((z₁=z₁,z₂=z₂), CartesianGrid(100,100))

  agri = georef(CSV.File(joinpath(datadir,"agriculture.csv")), (:x,:y))

  if visualtests
    @test_reference "data/agriculture.png" plot(agri, (:band4,:crop), ms=0.2)
  end

  @testset "HScatter" begin
    if visualtests
      sdata = georef(CSV.File(joinpath(datadir,"samples2D.tsv")), (:x,:y))
      p0 = hscatter(sdata, :value, lag=0)
      p1 = hscatter(sdata, :value, lag=1)
      p2 = hscatter(sdata, :value, lag=2)
      p3 = hscatter(sdata, :value, lag=3)
      plt = plot(p0, p1, p2, p3, layout=(2,2), size=(600,600))
      @test_reference "data/hscatter.png" plt
    end
  end
end


text:
const NORMAL = :normal
const LOGNORMAL = :lognormal
const SIGMA = :sigma
const MEAN = :mean

# CSV methods

function add_measurements!(
  scenario::Scenario,
  vector::AbstractVector;
  subset::AbstractVector{P} = Pair{Symbol, Symbol}[]
) where P <: Pair{Symbol, Symbol}
  selected_rows = _subset(vector, subset)

  for row in selected_rows
    _add_measurement!(scenario, row)
  end
end

function add_measurements!(
  platform::Platform,
  vector::AbstractVector;
  subset::AbstractVector{P} =  Pair{Symbol, Symbol}[]
) where P <: Pair{Symbol, Symbol}
  selected_rows = _subset(vector, subset)

  # we will store here lost scenario names
  lost_refs = Symbol[]

  for row in selected_rows
    scenario_ref = row[:scenario]

    if !haskey(platform.scenarios, scenario_ref)
      push!(lost_refs, scenario_ref)
    else
      scenario = platform.scenarios[scenario_ref]
      _add_measurement!(scenario, row)
    end
  end

  if length(lost_refs) > 0
    @warn "Lost scenario names: $(unique(lost_refs)). Some measurement points haven't been added."
  end

  return nothing
end

# DataFrame methods
"""
    add_measurements!(
      scenario::Scenario,
      df::DataFrame;
      kwargs...
    )

Adds measurements to `Scenario`

Arguments:

- `scenario` : simulation scenario of type [`Scenario`](@ref)
- `df` : `DataFrame` with measurements, typically obtained with [`read_measurements`](@ref) function
- `subset` : subset of measurements which will be added to the `Scenario`. Default `Pair{Symbol, Symbol}[]` adds all measurements from the `df`
"""
function add_measurements!(
  scenario::Scenario,
  df::DataFrame;
  kwargs...
)
  add_measurements!(scenario, eachrow(df); kwargs...)
end

"""
    add_measurements!(
      platform::Platform,
      df::DataFrame;
      kwargs...
    )

Adds measurements to `Scenario`

Arguments:

- `platform` : platform of [`Platform`](@ref) type
- `df` : `DataFrame` with measurements, typically obtained with [`read_measurements`](@ref) function
- `subset` : subset of measurements which will be added to the `Scenario`. Default `Pair{Symbol, Symbol}[]` adds all measurements from the `df`
"""
function add_measurements!(
  platform::Platform,
  df::DataFrame;
  kwargs...
)
  add_measurements!(platform, eachrow(df); kwargs...)
end

# private function to add one measurement row into scenario

function _add_measurement!(scenario::Scenario, row::Any) # maybe not any
  _t = row[:t]
  _val = row[:measurement]
  _scope = haskey(row, :scope) ? row[:scope] : missing

  _type = haskey(row, Symbol("prob.type")) ? row[Symbol("prob.type")] : missing
  type = ismissing(_type) ? NORMAL : _type

  if type in [NORMAL, LOGNORMAL]
    _mean = typed(row[Symbol("prob.$MEAN")])
    _sigma = typed(row[Symbol("prob.$SIGMA")])

    point = (type == LOGNORMAL) ? LogNormalMeasurementPoint(_t, _val, _scope, _mean, _sigma) : NormalMeasurementPoint(_t, _val, _scope, _mean, _sigma)
  else
    error("Distribution value $type is wrong or not supported. Supported distributions are: $Normal, $Lognormal")
  end

  push!(scenario.measurements, point)
end

# helper to read from csv and xlsx

function read_measurements_csv(filepath::String; kwargs...)
  df = DataFrame(CSV.File(
    filepath;
    typemap = Dict(Int64=>Float64, Int32=>Float64),
    types = Dict(:t=>Float64, :measurement=>Float64, :scope=>Symbol, :scenario=>Symbol, Symbol("prob.type")=>Symbol),
    kwargs...)
  )
  assert_measurements(df)

  return df
end

function read_measurements_xlsx(filepath::String, sheet=1; kwargs...)
  df = DataFrame(XLSX.readtable(filepath, sheet,infer_eltypes=true)...)
  assert_measurements(df)

  df[!,:t] .= Float64.(df[!,:t])
  df[!,:measurement] .= Float64.(df[!,:measurement])
  hasproperty(df, :scope) && (df[!,:scope] .= Symbol.(df[!,:scope]))
  df[!,:scenario] .= Symbol.(df[!,:scenario])
  hasproperty(df, Symbol("prob.type")) && (df[!,Symbol("prob.type")] .= Symbol.(df[!,Symbol("prob.type")]))

  return df
end

"""
    read_measurements(filepath::String, sheet=1; kwargs...)

Reads table file with measurements to `DataFrame`

Arguments:

- `filepath` : path to table file. Supports ".csv" and ".xlsx" files
- `sheet` : number of sheet in case of ".xlsx" file. Default is `1`
- kwargs : other arguments supported by `CSV.File`
"""
function read_measurements(filepath::String, sheet=1; kwargs...)
  ext = splitext(filepath)[end]

  if ext == ".csv"
    df = read_measurements_csv(filepath; kwargs...)
  elseif ext == ".xlsx"
    df = read_measurements_xlsx(filepath, sheet)
  else
    error("Extension $ext is not supported.")
  end
  return df
end

function assert_measurements(df)
  names_df = names(df)
  for f in ["t", "measurement", "scenario"]
    @assert f ∈ names_df "Required column name $f is not found in measurements table."
  end
  return nothing
end

#################################### GET MEASUREMENTS ###########################

# ! Current assumption is that error dist is normal or lognormal

measurements_as_table(sr::SimResult) = measurements_as_table(sr.scenario)

function measurements_as_table(scn::Scenario)
  meas = measurements(scn)
  @assert !isempty(meas) "Scenario doesn't contain measurements."

  lm = length(meas)
  _t = Vector{Float64}(undef,lm)
  _mea = []
  _scope = Vector{Symbol}(undef,lm)
  _mu = []
  _sigma = []
  _type = Vector{Symbol}(undef,lm)

  for (i,m) in enumerate(meas)
    _t[i] = m.t
    push!(_mea, m.val)
    _scope[i] = m.scope
    push!(_mu, m.μ)
    push!(_sigma, m.σ)
    _type[i] = typeof(m) <: NormalMeasurementPoint ? NORMAL : LOGNORMAL
  end
  df = DataFrame(t =_t, measurement = _mea, scope = _scope)
  df[!,"prob.mean"] = _mu
  df[!,"prob.sigma"] = _sigma
  df[!,"prob.type"] = _type
  return df
end


text:
function BFS(G,s)
    u = 1:size(G.AdjMatrix,1)
    for v in u
        if v != s
            G.V[v].color = "W"
            G.V[v].d = typemax(Int)
            G.V[v].pred = -1
        end
    end

    G.V[s].color = "G"
    G.V[s].d = 0
    G.V[s].prd = -1

    Q = Int[]
    push!(Q,s)
    while !isempty(Q)
        k = Q[1]
        deleteat!(Q,1)
        for m in findall(G.AdjMatrix[k,:] .== 1)
            if G.V[m].color == "W"
                G.V[m].color == "G"
                G.V[m].d = G.V[k].d + 1
                G.V[m].pred = k
                push!(Q,m)
            end
        end
        G.V[k].color = "B"
    end
end

function getPath(G,idStart,idEnd)
    path = Int[]
    temp = idEnd
    while temp != idStart
        path = [temp; path]
        temp = G.V[temp].pred
    end
    return [idStart;path]
end


text:
module BlockTests

using Blox
using FactCheck

facts("Blox") do
  @fact 1 => 1

  x = collect(1:4)
  y = collect(5:6)
  vw = view(x, y)
  @fact length(vw) => length(x) + length(y)
  @fact size(vw) => (length(vw),)

  for i in 1:6
    @fact vw[i] => i
  end
  @fact_throws vw[0]
  @fact_throws vw[7]

  vw[6] = 555
  @fact vw[6] => 555

  # default vals at the end
  vw = view(x, y; len=10)
  for i in 1:5
    @fact vw[i] => i
  end
  for i in 7:10
    @fact vw[i] => 0
  end

  vw = view(x, x; def=1, len=10)
  @fact vw[3] => vw[7]
  @fact vw[9] => 1
  vw[1] = 33
  @fact vw[1] => 33
  @fact vw[1] => vw[5]

  # staggered blocks
  defval = 1
  vw = view(x, y; def=defval, len=12, indices=[2,8])
  @fact vw[1] => defval
  @fact vw[2] => 33
  @fact vw[6] => defval
  @fact vw[8] => 5
  @fact vw[10] => defval
  @fact_throws vw[7] = 10

end

FactCheck.exitstatus()
end


text:
module WebDriver

using PyCall

import Base: size, close, quit, parent, get

include("defs.jl")
include("generated.jl")
include("lib.jl")

end # module


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# ╔═╡ 6dd76507-be07-4c3a-9747-99eca74d3519
begin
	filter!(path -> path != "@v#.#", LOAD_PATH)
	using Pkg
	Pkg.activate(Base.current_project())
	Pkg.instantiate()
	using PlutoUI, StatisticsWithJulia
	PlutoUI.TableOfContents(title = "初识 Julia")
end

# ╔═╡ 3baed4bd-1c0c-4d44-b7ab-f7c0627fda1b
begin
	Pkg.add("AbstractTrees") # command + /: 切换为注释形式
	using AbstractTrees
	AbstractTrees.children(x::Type) = subtypes(x)
	with_terminal(() -> print_tree(Number))
end

# ╔═╡ 2edcc72c-bfa4-4623-a34f-220c0df9c4d3
begin
	Pkg.add("HTTP")
	Pkg.add("JSON")
	using HTTP, JSON
end

# ╔═╡ d507db6a-f869-4dc3-b7ef-fd3518fc1660
Pkg.add.(["Plots", "LaTeXStrings", "Measures"])

# ╔═╡ 2fdb3281-7229-4f18-b9d9-bc5e3086d903
using Statistics

# ╔═╡ f924796f-52f3-4b4b-aa47-61a179874647
using Plots, LaTeXStrings, Measures

# ╔═╡ ec3f7b5c-cee9-432f-b33c-a08e16c0171a
StatisticsWithJulia.header(title = @__FILE__, author="")

# ╔═╡ 9d98fa7a-3654-40a0-8a96-75c3fca007ae
md"""
## 本章讲义

- [讲义下载](http://likan.oss-cn-beijing.aliyuncs.com/StatisticsWithJulia/01-%E5%88%9D%E8%AF%86Julia.pdf)

"""

# ╔═╡ 42da8f74-1b44-11ec-32f2-276bdcef1c69
md"""
## 1.1. 你好世界
"""

# ╔═╡ e71aacbb-efff-4b93-b5c4-cbdfa88981ad
md"""
- `command/control + m`: 把当下单元格切换到markdown格式
"""

# ╔═╡ 0660d564-9a26-417d-9a2c-42ef181187c2
with_terminal() do
    println("你好， 世界!")
end

# ╔═╡ 77c5ddd3-ddfc-4a57-ab48-e9db765057d3
with_terminal(() -> println("你好， 世界!"))

# ╔═╡ 8eefb8aa-ae85-4231-abf5-9948cdc92ec1
md"""
`do f() end` 等价于 `（） -> f()`
"""

# ╔═╡ d22fc28f-c591-4db2-8275-54d93fda85e9
f(x) = BigInt(x)^BigInt(x)

# ╔═╡ 4ece5fc8-18d5-408a-851b-31ce58015e72
map(x -> 3x, 1:10)

# ╔═╡ c266f6e2-4e8e-4870-ac38-d4d0e8b9d042
map(1:10) do x
	3x
end

# ╔═╡ 8f3277e9-7905-4beb-9ef1-c3e18b2f449a
map(x -> rand(), 1:10)

# ╔═╡ f6e99ac9-2183-4ea2-991e-f3422b1d1f20
map(__ -> rand(), 1:10)

# ╔═╡ cd0ce1b4-9bc9-4218-95e9-e576abeb1af1
# map(() -> rand(), 1:10) # 不合语法

# ╔═╡ 37fec978-f817-4aa8-85fc-4912eaeb402e
helloArray = ["你好", "👋", "世界"]

# ╔═╡ adc183db-ff23-43de-8cf6-a961d13c2f01
with_terminal() do
	for i ∈ eachindex(helloArray)
		println(helloArray[i])
	end
end

# ╔═╡ 4f1ce95f-2f2f-4817-9d5f-348184a1b800
[i^2 for i in 1:10 for __ in 1:10] # 一维数组

# ╔═╡ dec64158-9161-4b37-b19f-432a5635719e
[i^2 for i in 1:10, __ in 1:10]   # 矩阵

# ╔═╡ 5a401f10-09d2-4939-87fe-621f469a68ee
Square = [i^2 for i in 1:10]

# ╔═╡ c2c7c916-1ba8-4e14-a25b-e15faf676f8a
with_terminal() do
	for i ∈ eachindex(Square)
			println(Square[i])
	end
end

# ╔═╡ 93b0afb2-e9f3-434d-b21a-2ccad6a12d90
[rand() for __ in 1:100]

# ╔═╡ 8bcc15c7-d66a-40cc-9f33-50c6de62fadc
md"""
### 运行速度
"""

# ╔═╡ 428c124a-1052-45c8-a49b-b1047267d644
md"""
- LaTeX Syntax

```math
y = \frac{\displaystyle\sum_{i=1}^\infty}{\epsilon\times\sigma}

```
"""

# ╔═╡ 863624c7-f378-467a-bed2-94182bae73b7
with_terminal() do
	@time begin
		global data = Float64[]
		for __ in 1:10^6
			group = Float64[]
			for __ in 1:5*10^2
				push!(group, rand())
			end
			push!(data, mean(group))
		end
	end
end

# ╔═╡ 0e8b6a32-ba23-4e37-b78c-8ffb3779a43e
"98%的平均值位于 $(quantile(data, 0.01)) 和 $(quantile(data, 0.99)) 之间。"

# ╔═╡ 9e35d6a1-df94-4a27-9e04-234548495780
with_terminal() do
	@time global data2 = [mean(rand(5 * 10 ^ 2)) for __ in 1 : 10 ^ 6]
end

# ╔═╡ 1b45a9bf-56c7-4b7b-bdd5-b4fcaa572ec0
"98%的平均值位于 $(quantile(data2, 0.01)) 和 $(quantile(data2, 0.99)) 之间。"

# ╔═╡ 60b263eb-f786-40c4-b60b-2daabdacba7b
md"""
### 类型和多重派发
"""

# ╔═╡ b0db7a9c-2498-46c5-bf3c-280d3083ba88
# @edit 1 + 2

# ╔═╡ 412ee0c1-26ad-4b9f-9d6e-c552d8431e76
# @edit 1.0 + 2

# ╔═╡ aa5e313d-19d9-4300-a65e-f547f7e571d0
mutable struct zhan
	dt::Int64
end

# ╔═╡ 0ba25fa7-72cb-4d6e-bfaf-eddcf12b9f9c
x = zhan(4)

# ╔═╡ d5251a23-b71e-4faf-9f5a-f6ab1d2ffb7d
y = zhan(8)

# ╔═╡ b717d872-0d16-4556-a6c8-1ccc3c7b76c3
Base.:+(x::zhan, y::zhan) = x.dt * 100 + y.dt * 200

# ╔═╡ a0322171-cff9-48bb-9f73-792c8ea4798f
function g(x)
	x + 1000
end

# ╔═╡ 3711c338-748f-448e-a06d-d3c29094a14b
map(g, 1:10)

# ╔═╡ f77b7762-29ee-4472-a4ee-a5c1a6dfb369
sqrt.(Square) .+ 100

# ╔═╡ e9d5ee6c-de1a-41ee-850e-a6f5471537ea
@. sqrt(Square) + 100

# ╔═╡ 90fe2c2b-d62d-40b5-9a34-0eb1e9836250
[10 + x for x in 1:100]

# ╔═╡ 72b3ed37-3f3f-4a7d-9036-dadb06c3e1e4
x + y

# ╔═╡ 8c6e966d-9185-4a44-8cd9-e69528cee54b
# methods(+)

# ╔═╡ 1e9be5f5-ab15-4a29-8fd6-d20fee084ad5
abstract type LIKANZHAN <: Signed end #抽象类型

# ╔═╡ 59b884e5-a55a-42f6-9ed4-404d55eae5e9
md"""
## 1.2. 几个例子
"""

# ╔═╡ 3cbd7219-7fb0-4c6e-a2ad-724e97b67d94
function bubbleSort!(a)
	n = length(a)
	for i ∈ 1:n-1
		for j ∈ 1:n-i
			if a[j] > a[j+1]
				a[j], a[j+1] = a[j+1], a[j] # Julian Way
			end
		end
	end
	return a
end

# ╔═╡ 93a86467-cf02-4d3c-b127-79247a588250
data3 = [65, 66, 55, 44, 20, 34, 78, 80]

# ╔═╡ 0423223e-90c5-4d86-8370-66ddee6ba496
with_terminal() do
	@time bubbleSort!(rand(10^5))
end

# ╔═╡ 721c67fd-6d40-4171-89c4-21a1c0064880
with_terminal() do
	@time sort!(rand(10^5))
end

# ╔═╡ 5a4dd40b-d917-4634-a4d4-255c63d68d3f
md"""
- 网络接口和 JSON
"""

# ╔═╡ 20097ff6-c516-454c-978f-4c92e6c72d52
data4 = HTTP.request("GET","https://ocw.mit.edu/ans7870/6/6.006/s08/lecturenotes/files/t8.shakespeare.txt");

# ╔═╡ fbbc0558-7916-4308-b8be-4428f47ca4e5
shakespeare = String(data4)

# ╔═╡ 290d3b7f-ee12-476c-a4f2-189671c84398
# shakespeareWords = split(shakespeare)

# ╔═╡ cc728038-b20f-4b93-97a7-883d5734a0b2
# jsonWords = HTTP.request("GET", "https://raw.githubusercontent.com/"*"h-Klok/StatsWithJuliaBook/master/data/jsonCode.json")

# ╔═╡ 2eca9295-370d-41e9-95a4-b56ee2b0e1dc
"julia" * " is good"

# ╔═╡ b4c5a72d-5ce1-4e21-a04f-8344d1023648
string("Julia", " is Good", ". Python", " is slow.")

# ╔═╡ 4f11cc1a-a78e-4d6c-8bd2-857f85f48aa8
ls = 1:7

# ╔═╡ 50039da8-2c71-4a1d-8cc1-2d4bd8126988
typeof(ls)

# ╔═╡ efb2b80c-af6c-4f6e-bcba-e57c7c39cab3
collect(ls)

# ╔═╡ 1ca893dd-2040-47df-8fa0-c3cb4dd0d7dc
md"""
- 数据存储的常见格式：数组
"""

# ╔═╡ dbbcd955-e0c9-4a6c-81bf-fe5c5ad2f018
md"""
- 数据存储的另外格式： 词典
"""

# ╔═╡ 5565b0a4-792b-40f6-82a2-f8f642432129
dct = Dict("x" => collect(ls), "y" => collect(ls .+ 1) )

# ╔═╡ 6d4175a4-af6e-4a08-bcd2-0b1966ae9035
getindex(dct, "x")

# ╔═╡ c4a58e43-2ea8-48fa-a99d-09f83b1b8ba8
md"""
- 数据存储的另外格式： 元组
"""

# ╔═╡ e0e5ef4a-fc35-4d9b-9a34-a9439d87051d
vec = [1, 2, 3]; typeof(vec)

# ╔═╡ c9d041cb-6a29-449a-8909-b529e8c562a6
yz = (1, 2, 4); typeof(yz)

# ╔═╡ c77241fb-bd9e-4270-bd04-e8b5bcb2406b
md"""
## 1.3. 做图初步
"""

# ╔═╡ 87296377-fe6d-4225-b114-a8c3b8bea1bc
with_terminal() do
	Pkg.status()
end

# ╔═╡ aa61efe4-b3a2-4379-a748-f37cd7f9c6fb
f(x, y) = x^2 + y^2

# ╔═╡ 9d9f4d96-7755-4fd1-851a-092929ad1ac7
f(6)

# ╔═╡ 8b0db1dd-5d79-42e5-9763-ca4ef7259060
map(f, 1:10)

# ╔═╡ 182ab07b-a3e6-47f0-8de5-d59113027c43
f0(x) = f(x, 0)

# ╔═╡ 6262451f-4dcf-4f84-8be6-2f8fbf2631bd
f2(x) = f(x, 2)

# ╔═╡ 3d664567-d879-41f7-9e1a-cdfc08e4e7c9
xVals, yVals = -5:0.1:5, -5:0.1:5

# ╔═╡ 5e03a17e-861e-4749-b18e-00e293d24f28
plot([f0, f2], xVals) # command + /

# ╔═╡ 4d5c7119-f2bf-44ed-bc20-5597a795fbcd
plot(xVals, [f0.(xVals), f2.(xVals)],
	color = [:blue :green],
	xlims = (-6, 6),
	ylims = (-5, 30),
	xlabel = "Range of X",
	ylabel = L"f(x, \cdot)",
	label = [L"f(x, 0)" L"f(x, 2)"],
	legends = :bottomleft
)

# ╔═╡ af2e86c1-4c58-4577-9ae9-313f8ae78475
p1 = annotate!(0, -0.2, text("(0, 0) the mainimum \n of f(x, 0)", :left, :top, 10))

# ╔═╡ 9461292e-3658-41e2-8349-45790fa7c9e0
z = [f(x, y) for y in yVals, x in xVals];

# ╔═╡ 7e03b8d1-1ddb-435f-b9ba-b72b0c44da18
p2 = surface(xVals, yVals, z,
	color = cgrad([:blue, :red]),
	legend = :none,
	ylabel = "y",
	zlabel = L"f(x, y)"
)

# ╔═╡ e54420c1-1ff3-4f38-b9b8-0e3a18127c9a
M = z[1:10, 1:10];

# ╔═╡ 458dfaeb-0370-4e61-90a2-f5b54a0fab42
p3 = heatmap(M,
	color = cgrad([:blue, :red]),
	ylabel = "y",
	xticks = ([1:10; ], xVals),
	yticks = ([1:10; ], yVals)
)

# ╔═╡ b7729c93-781a-46ec-b057-65c4e9e15a04
plot(p1, p2, p3, layout = (1, 3), size = (1200, 400), xlabel = "x", margin = 5mm)

# ╔═╡ 543de112-675c-4368-957d-f37647ad85bd
 L"\frac{\delta}{\gamma\times \cdots}"

# ╔═╡ 37f6486e-e03c-4ddb-bca6-d13157d69822
plot([f0, f2], 1:0.5:10,
	seriestype = [:scatter, :line],
	markersize = 6,
	linewidth = 3, linecolor = :green,
	legend = false)

# ╔═╡ 68d89124-d6e9-4609-b12a-d3b508733680
# surface() # 等价于 plot(.., seriestype = :surface,...,)

# ╔═╡ 6003cd52-89ea-4d93-ac00-43ca395357b8
md"""
### 1.3.1 希尔顿序列
"""

# ╔═╡ 92c42555-8550-4c1b-9b1a-502813c40237
md"""
- 条件句
"""

# ╔═╡ 5c3f1925-ebfb-4bef-8425-635262927bf1
function ff1(x)
	if x % 3 == 0
		return "$x 除以 3 的余数是 0"
	elseif x % 3 == 1
		return "$x 除以 3 的余数是 1"
	else
		return "$x 除以 3 的余数是 2"
	end
end

# ╔═╡ ec02b4a7-529b-4d37-ba4b-bf146a33ee0f
[ff1(x) for x in 1:5]

# ╔═╡ 954f166b-e8db-4cf7-a0a2-75ab3b3acbe0
function ff2(x)
	if x % 3 == 0
		return "x 除3的余数是 0"
	end
end

# ╔═╡ e97dcedb-4b53-4ce8-a0e9-ed20b77d8b2d
[ff2(x) for x in 1:5]

# ╔═╡ f662c9f4-8004-43d6-ad1e-872ad3d24e48
md"""
- 类型稳定性： 除号可能会改变数据的类型
"""

# ╔═╡ c72335d5-39e3-4ab9-a398-e3960dcafe0c
xd = Int(2)

# ╔═╡ 8ebae10c-3148-43fe-9b14-f07934c3b2ea
typeof(xd)

# ╔═╡ e49a1426-bcf5-453a-8753-2bb8387e67eb
yd = xd/2

# ╔═╡ 46bd26e6-c1bf-4fa3-890a-2d424865c2dc
typeof(yd)

# ╔═╡ 55e79482-e30e-4eea-b083-cac2a9667348
md"""
- 希尔顿序列
"""

# ╔═╡ 9537df9b-7fc8-4b32-bc7c-d8ee8730d73a
function hailLength1(x::Int)
	n = 0
	while x != 1
		if x % 2 == 0
			x = Int(x / 2)
		else
			x = 3x + 1
		end
		n += 1 # 等价于 n = n + 1, *=, /=
	end
	return n
end

# ╔═╡ 6434a9f1-aa8a-491f-ac7e-c58b3a4b2638
function hailLength(x::Int)
	n = 0
	while x != 1
		x = x % 2 == 0 ? Int(x / 2) : 3x + 1
		n += 1
	end
	return n
end

# ╔═╡ 4abd57c6-1a0c-4e7c-a195-d07ddc299b2c
md"""
注意： `condition ? res1 : res2` 等价于 `if condition res1 else res2`
"""

# ╔═╡ aa2c895f-7219-4b98-bf50-c96d676ead09
x2 = 13

# ╔═╡ 5fb204e1-d150-4a3c-ac19-6fa5b67b072b
if x2 % 2 == 1 "奇数" else "偶数" end

# ╔═╡ c0df29b8-368d-40d1-bad0-73493531903a
lengths = [hailLength(x0) for x0 in 2:10^7];

# ╔═╡ 51046c66-6e14-420e-8a5e-5a58580841a6
# histogram(lengths)

# ╔═╡ 3659f664-3fc3-48bd-94e3-f57ae84e7301
plot(lengths, seriestype = :histogram,
	bins = 1000,
	# fill = (:green, true) # 元组
	fill = true, color = :green,
	normed = true, legend = false,
	linealpha = 0.2
)

# ╔═╡ 9392cdfb-ceef-40ac-a67b-8c7ca33de123
md"""
### 1.3.2 动图
"""

# ╔═╡ 20605ee6-f245-4e77-b014-a8e42a520e8d
n = 10

# ╔═╡ c69f8278-9627-458b-8204-91ae08efbf92
points = [exp(2 * pi * im * k / n) for k in 1:n];

# ╔═╡ a5f4d797-db11-4d6a-8906-b4af6804a936
plot(points, seriestype = :scatter, legend = false);

# ╔═╡ 46f62dad-820b-42d3-8cdd-ab02cacabc47
xPts, yPts = real.(points), imag.(points);

# ╔═╡ 758082fe-038e-4090-b658-111db3ee1f13
edges  = [(v, u) for v in 1:n for u in (v+1):n];

# ╔═╡ 931e4862-c109-4b17-9048-7826f3dd0c6f
animate = Animation();

# ╔═╡ 3017c6db-acee-422d-bb97-692495972631
plot(xPts, yPts, seriestype = :scatter, legend = false) # 等价于 scatter(xPts, yPts)

# ╔═╡ 404716dd-7a05-4f3b-b336-cab7fc679ed1
for i in eachindex(edges)
	u, v = edges[i][1], edges[i][2]
	xpoints = [xPts[u], xPts[v]]
	ypoints = [yPts[u], yPts[v]]
	plot!(xpoints, ypoints, line = :red, legend = false)
	frame(animate)
end

# ╔═╡ 43b3b163-05a7-4c06-a948-bf5c001b7af0
gif(animate)

# ╔═╡ 4c0169d7-b280-4326-aa79-9c979280077a
md"""
## 1.4. 随机数
"""

# ╔═╡ 694d2ceb-233f-49ea-9416-9e8a8f272d38
md"""
## 1.6. 嵌入其他语言
"""

# ╔═╡ 72a095fa-c5cf-4eb1-909d-d8b49af0bdc9
StatisticsWithJulia.list_notebooks(@__FILE__)

# ╔═╡ 8a20fa4c-338c-48c9-a12a-02e6d67153df
StatisticsWithJulia.footer() ## FFFFF

# ╔═╡ Cell order:
# ╟─ec3f7b5c-cee9-432f-b33c-a08e16c0171a
# ╟─6dd76507-be07-4c3a-9747-99eca74d3519
# ╟─9d98fa7a-3654-40a0-8a96-75c3fca007ae
# ╟─42da8f74-1b44-11ec-32f2-276bdcef1c69
# ╟─e71aacbb-efff-4b93-b5c4-cbdfa88981ad
# ╠═0660d564-9a26-417d-9a2c-42ef181187c2
# ╠═77c5ddd3-ddfc-4a57-ab48-e9db765057d3
# ╟─8eefb8aa-ae85-4231-abf5-9948cdc92ec1
# ╠═d22fc28f-c591-4db2-8275-54d93fda85e9
# ╠═9d9f4d96-7755-4fd1-851a-092929ad1ac7
# ╠═4ece5fc8-18d5-408a-851b-31ce58015e72
# ╠═8b0db1dd-5d79-42e5-9763-ca4ef7259060
# ╠═a0322171-cff9-48bb-9f73-792c8ea4798f
# ╠═3711c338-748f-448e-a06d-d3c29094a14b
# ╠═c266f6e2-4e8e-4870-ac38-d4d0e8b9d042
# ╠═8f3277e9-7905-4beb-9ef1-c3e18b2f449a
# ╠═f6e99ac9-2183-4ea2-991e-f3422b1d1f20
# ╠═cd0ce1b4-9bc9-4218-95e9-e576abeb1af1
# ╠═37fec978-f817-4aa8-85fc-4912eaeb402e
# ╠═adc183db-ff23-43de-8cf6-a961d13c2f01
# ╠═4f1ce95f-2f2f-4817-9d5f-348184a1b800
# ╠═dec64158-9161-4b37-b19f-432a5635719e
# ╠═5a401f10-09d2-4939-87fe-621f469a68ee
# ╠═c2c7c916-1ba8-4e14-a25b-e15faf676f8a
# ╠═2fdb3281-7229-4f18-b9d9-bc5e3086d903
# ╠═f77b7762-29ee-4472-a4ee-a5c1a6dfb369
# ╠═e9d5ee6c-de1a-41ee-850e-a6f5471537ea
# ╠═93b0afb2-e9f3-434d-b21a-2ccad6a12d90
# ╠═90fe2c2b-d62d-40b5-9a34-0eb1e9836250
# ╟─8bcc15c7-d66a-40cc-9f33-50c6de62fadc
# ╠═428c124a-1052-45c8-a49b-b1047267d644
# ╠═863624c7-f378-467a-bed2-94182bae73b7
# ╠═0e8b6a32-ba23-4e37-b78c-8ffb3779a43e
# ╠═9e35d6a1-df94-4a27-9e04-234548495780
# ╠═1b45a9bf-56c7-4b7b-bdd5-b4fcaa572ec0
# ╟─60b263eb-f786-40c4-b60b-2daabdacba7b
# ╠═b0db7a9c-2498-46c5-bf3c-280d3083ba88
# ╠═412ee0c1-26ad-4b9f-9d6e-c552d8431e76
# ╠═aa5e313d-19d9-4300-a65e-f547f7e571d0
# ╠═0ba25fa7-72cb-4d6e-bfaf-eddcf12b9f9c
# ╠═d5251a23-b71e-4faf-9f5a-f6ab1d2ffb7d
# ╠═b717d872-0d16-4556-a6c8-1ccc3c7b76c3
# ╠═72b3ed37-3f3f-4a7d-9036-dadb06c3e1e4
# ╠═8c6e966d-9185-4a44-8cd9-e69528cee54b
# ╠═1e9be5f5-ab15-4a29-8fd6-d20fee084ad5
# ╠═3baed4bd-1c0c-4d44-b7ab-f7c0627fda1b
# ╟─59b884e5-a55a-42f6-9ed4-404d55eae5e9
# ╠═3cbd7219-7fb0-4c6e-a2ad-724e97b67d94
# ╠═93a86467-cf02-4d3c-b127-79247a588250
# ╠═0423223e-90c5-4d86-8370-66ddee6ba496
# ╠═721c67fd-6d40-4171-89c4-21a1c0064880
# ╠═5a4dd40b-d917-4634-a4d4-255c63d68d3f
# ╠═2edcc72c-bfa4-4623-a34f-220c0df9c4d3
# ╠═20097ff6-c516-454c-978f-4c92e6c72d52
# ╠═fbbc0558-7916-4308-b8be-4428f47ca4e5
# ╠═290d3b7f-ee12-476c-a4f2-189671c84398
# ╠═cc728038-b20f-4b93-97a7-883d5734a0b2
# ╠═2eca9295-370d-41e9-95a4-b56ee2b0e1dc
# ╠═b4c5a72d-5ce1-4e21-a04f-8344d1023648
# ╠═4f11cc1a-a78e-4d6c-8bd2-857f85f48aa8
# ╠═50039da8-2c71-4a1d-8cc1-2d4bd8126988
# ╠═efb2b80c-af6c-4f6e-bcba-e57c7c39cab3
# ╟─1ca893dd-2040-47df-8fa0-c3cb4dd0d7dc
# ╟─dbbcd955-e0c9-4a6c-81bf-fe5c5ad2f018
# ╠═5565b0a4-792b-40f6-82a2-f8f642432129
# ╠═6d4175a4-af6e-4a08-bcd2-0b1966ae9035
# ╟─c4a58e43-2ea8-48fa-a99d-09f83b1b8ba8
# ╠═e0e5ef4a-fc35-4d9b-9a34-a9439d87051d
# ╠═c9d041cb-6a29-449a-8909-b529e8c562a6
# ╟─c77241fb-bd9e-4270-bd04-e8b5bcb2406b
# ╠═d507db6a-f869-4dc3-b7ef-fd3518fc1660
# ╠═f924796f-52f3-4b4b-aa47-61a179874647
# ╠═87296377-fe6d-4225-b114-a8c3b8bea1bc
# ╠═aa61efe4-b3a2-4379-a748-f37cd7f9c6fb
# ╠═182ab07b-a3e6-47f0-8de5-d59113027c43
# ╠═6262451f-4dcf-4f84-8be6-2f8fbf2631bd
# ╠═3d664567-d879-41f7-9e1a-cdfc08e4e7c9
# ╠═5e03a17e-861e-4749-b18e-00e293d24f28
# ╠═4d5c7119-f2bf-44ed-bc20-5597a795fbcd
# ╠═af2e86c1-4c58-4577-9ae9-313f8ae78475
# ╠═9461292e-3658-41e2-8349-45790fa7c9e0
# ╠═7e03b8d1-1ddb-435f-b9ba-b72b0c44da18
# ╠═e54420c1-1ff3-4f38-b9b8-0e3a18127c9a
# ╠═458dfaeb-0370-4e61-90a2-f5b54a0fab42
# ╠═b7729c93-781a-46ec-b057-65c4e9e15a04
# ╠═543de112-675c-4368-957d-f37647ad85bd
# ╠═37f6486e-e03c-4ddb-bca6-d13157d69822
# ╠═68d89124-d6e9-4609-b12a-d3b508733680
# ╟─6003cd52-89ea-4d93-ac00-43ca395357b8
# ╟─92c42555-8550-4c1b-9b1a-502813c40237
# ╠═5c3f1925-ebfb-4bef-8425-635262927bf1
# ╠═ec02b4a7-529b-4d37-ba4b-bf146a33ee0f
# ╠═954f166b-e8db-4cf7-a0a2-75ab3b3acbe0
# ╠═e97dcedb-4b53-4ce8-a0e9-ed20b77d8b2d
# ╟─f662c9f4-8004-43d6-ad1e-872ad3d24e48
# ╠═c72335d5-39e3-4ab9-a398-e3960dcafe0c
# ╠═8ebae10c-3148-43fe-9b14-f07934c3b2ea
# ╠═e49a1426-bcf5-453a-8753-2bb8387e67eb
# ╠═46bd26e6-c1bf-4fa3-890a-2d424865c2dc
# ╟─55e79482-e30e-4eea-b083-cac2a9667348
# ╠═9537df9b-7fc8-4b32-bc7c-d8ee8730d73a
# ╠═6434a9f1-aa8a-491f-ac7e-c58b3a4b2638
# ╠═4abd57c6-1a0c-4e7c-a195-d07ddc299b2c
# ╠═aa2c895f-7219-4b98-bf50-c96d676ead09
# ╠═5fb204e1-d150-4a3c-ac19-6fa5b67b072b
# ╠═c0df29b8-368d-40d1-bad0-73493531903a
# ╠═51046c66-6e14-420e-8a5e-5a58580841a6
# ╠═3659f664-3fc3-48bd-94e3-f57ae84e7301
# ╟─9392cdfb-ceef-40ac-a67b-8c7ca33de123
# ╠═20605ee6-f245-4e77-b014-a8e42a520e8d
# ╠═c69f8278-9627-458b-8204-91ae08efbf92
# ╠═a5f4d797-db11-4d6a-8906-b4af6804a936
# ╠═46f62dad-820b-42d3-8cdd-ab02cacabc47
# ╠═758082fe-038e-4090-b658-111db3ee1f13
# ╠═931e4862-c109-4b17-9048-7826f3dd0c6f
# ╠═3017c6db-acee-422d-bb97-692495972631
# ╠═404716dd-7a05-4f3b-b336-cab7fc679ed1
# ╠═43b3b163-05a7-4c06-a948-bf5c001b7af0
# ╟─4c0169d7-b280-4326-aa79-9c979280077a
# ╟─694d2ceb-233f-49ea-9416-9e8a8f272d38
# ╟─72a095fa-c5cf-4eb1-909d-d8b49af0bdc9
# ╟─8a20fa4c-338c-48c9-a12a-02e6d67153df


text:
include("solver_preamble.jl")
import ECOS
factory = with_optimizer(ECOS.Optimizer, verbose=false)
config = MOI.Test.Config(atol=1e-5, rtol=1e-5)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SOC" begin
    Tests.soc_test(factory, config, [
    # K = 30: Test Failed at /home/blegat/.julia/dev/SumOfSquares/test/Tests/options_pricing.jl:31
    #   Expression: JuMP.primal_status(model) == MOI.FEASIBLE_POINT
    #    Evaluated: MathOptInterface.NEARLY_FEASIBLE_POINT == MathOptInterface.FEASIBLE_POINT
    # K = 30: Test Failed at /home/blegat/.julia/dev/SumOfSquares/test/Tests/options_pricing.jl:32
    #   Expression: ≈(JuMP.objective_value(model), expected, atol=atol, rtol=rtol)
    #    Evaluated: 21.511440549896246 ≈ 21.51 (atol=1.0e-5, rtol=1.0e-5)
        "sdsos_options_pricing"])
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Empirikos

@test NormalSample(1.0, 2.0) < NormalSample(2.0, 1.0)
@test NormalSample(1.0, 2.0) < NormalSample(1.0, 3.0)
@test NormalSample(1.0, 2.0) <= NormalSample(1.0, 2.0)

@test BinomialSample(1, 3) < BinomialSample(2, 2)
@test BinomialSample(2, 2) < BinomialSample(2, 3)
@test BinomialSample(3, 4) > BinomialSample(2, 3)
@test BinomialSample(3, 3) <= BinomialSample(3, 3)

@test PoissonSample(1, 3) < PoissonSample(2, 2)
@test PoissonSample(2, 2) < PoissonSample(2, 3)
@test PoissonSample(3, 4) > PoissonSample(2, 3)
@test PoissonSample(3, 3.1) <= PoissonSample(3, 3.1)

Zs_comp = compound( [NormalSample(1,2), NormalSample(2,1),NormalSample(1,1), NormalSample(0.5,2)])
@test Zs_comp[1] < Zs_comp[2]
@test Zs_comp[1] == Zs_comp[3]

Zs = PoissonSample.(sample(1:100, 1000))

Zs_summary = summarize(Zs)

Zs_unique = collect(keys(Zs_summary))
@test issorted(Zs_unique)

@test collect(values(Zs_summary))[2] == sum(Zs .== Zs_unique[2])
@test collect(keys(Zs_summary.store)) == collect(keys(Zs_summary))

@test issorted([response(k) for (k, v) in Zs_summary.store])


text: # This file should contain file-parser for reading in data to a specific struct


text:
export iauHd2ae
"""
Equatorial to horizon coordinates:  transform hour angle and
declination to azimuth and altitude.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   ha       double       hour angle (local)
   dec      double       declination
   phi      double       site latitude

Returned:
   *az      double       azimuth
   *el      double       altitude (informally, elevation)

Notes:

   1.  All the arguments are angles in radians.

   2.  Azimuth is returned in the range 0-2pi;  north is zero, and east
   is +pi/2.  Altitude is returned in the range +/- pi/2.

   3.  The latitude phi is pi/2 minus the angle between the Earth's
   rotation axis and the adopted zenith.  In many applications it
   will be sufficient to use the published geodetic latitude of the
   site.  In very precise (sub-arcsecond) applications, phi can be
   corrected for polar motion.

   4.  The returned azimuth az is with respect to the rotational north
   pole, as opposed to the ITRS pole, and for sub-arcsecond
   accuracy will need to be adjusted for polar motion if it is to
   be with respect to north on a map of the Earth's surface.

   5.  Should the user wish to work with respect to the astronomical
   zenith rather than the geodetic zenith, phi will need to be
   adjusted for deflection of the vertical (often tens of
   arcseconds), and the zero point of the hour angle ha will also
   be affected.

   6.  The transformation is the same as Vh = Rz(pi)*Ry(pi/2-phi)*Ve,
   where Vh and Ve are lefthanded unit vectors in the (az,el) and
   (ha,dec) systems respectively and Ry and Rz are rotations about
   first the y-axis and then the z-axis.  (n.b. Rz(pi) simply
   reverses the signs of the x and y components.)  For efficiency,
   the algorithm is written out rather than calling other utility
   functions.  For applications that require even greater
   efficiency, additional savings are possible if constant terms
   such as functions of latitude are computed once and for all.

   7.  Again for efficiency, no range checking of arguments is carried
   out.

Last revision:   2017 September 12

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauHd2ae (double ha, double dec, double phi,
#                double *az, double *el)

function iauHd2ae(ha::Real, dec::Real, phi::Real)
   ref_az = Ref{Float64}(0.0)
   ref_el = Ref{Float64}(0.0)

   ccall((:iauHd2ae, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble,
         Ref{Cdouble}, Ref{Cdouble}),
          convert(Float64, ha), convert(Float64, dec),
          convert(Float64, phi),
          ref_az, ref_el)

   return ref_az[], ref_el[]
end


text:
mutable struct ncmpop_data
    n::Int # number of all variables
    m::Int # number of all constraints
    numeq::Int # number of equality constraints
    d::Int # relaxation order
    supp # support data
    coe # coefficient data
    obj # "eigen" or "trace"
    ksupp # extending support at the k-th step
    basis # monomial bses
    cql # number of cliques
    cliques # cliques of variables
    cliquesize # numbers of cliques
    J # constraints associated to each clique
    ncc # constraints associated to no clique
    blocks # block structure
    cl # numbers of blocks
    blocksize # sizes of blocks
    sb # sizes of different blocks
    numb # numbers of different blocks
end

function cs_nctssos_first(f, x; d=0, CS="MF", minimize=false, TS="block", merge=false, md=3,
    QUIET=false, obj="eigen", solve=true)
    println("***************************NCTSSOS***************************")
    println("NCTSSOS is launching...")
    n,supp,coe = poly_info(f, x)
    if d == 0
        d = ceil(Int, maxdegree(f)/2)
    end
    opt,data = cs_nctssos_first(supp, coe, n, d=d, CS=CS, minimize=minimize, TS=TS, merge=merge, md=md, QUIET=QUIET, obj=obj, solve=solve)
    return opt,data
end

function cs_nctssos_first(supp::Vector{Vector{UInt16}}, coe, n::Int; d=0, CS="MF",
    minimize=false, TS="block", merge=false, md=3, QUIET=false, obj="eigen", solve=true)
    if obj == "trace"
        supp,coe = cyclic_canon(supp, coe)
    else
        supp,coe = sym_canon(supp, coe)
    end
    time = @elapsed begin
    cliques,cql,cliquesize = clique_decomp(n, supp)
    end
    if CS != false && QUIET == false
        mc = maximum(cliquesize)
        println("Obtained the variable cliques in $time seconds. The maximal size of cliques is $mc.")
    end
    if TS != false && QUIET == false
        println("Starting to compute the block structure...")
    end
    time = @elapsed begin
    blocks,cl,blocksize,sb,numb,basis,_ = get_blocks_mix(d, supp, cliques, cql, cliquesize, TS=TS, merge=merge, md=md, obj=obj)
    end
    if TS != false && QUIET == false
        mb = maximum(maximum.(sb))
        println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
    end
    opt,ksupp = blockupop_mix(n, supp, coe, basis, cliques, cql, cliquesize, blocks, cl, blocksize, obj=obj, solve=solve, QUIET=QUIET)
    data = ncmpop_data(n, 0, 0, d, supp, coe, obj, ksupp, basis, cql, cliques, cliquesize, [], [], blocks, cl, blocksize, sb, numb)
    return opt,data
end

"""
    opt,data = cs_nctssos_first(pop, x, d; numeq=0, CS="MF", TS="block", merge=false, md=3,
    QUIET=false, obj="eigen", solve=true)

Compute the first step of the CS-NCTSSOS hierarchy for constrained noncommutative polynomial
optimization with relaxation order `d`. Return the optimum and other auxiliary data.

# Arguments
- `pop`: the vector of the objective function, inequality constraints, and equality constraints.
- `x`: the set of noncommuting variables.
- `d`: the relaxation order of the moment-SOHS hierarchy.
- `numeq`: the number of equality constraints.
"""
function cs_nctssos_first(pop, x, d; numeq=0, CS="MF", minimize=false, assign="first", TS="block", merge=false, md=3, QUIET=false, obj="eigen", solve=true)
    n,supp,coe = polys_info(pop, x)
    opt,data = cs_nctssos_first(supp, coe, n, d, numeq=numeq, CS=CS, minimize=minimize, assign=assign, TS=TS, QUIET=QUIET, obj=obj, solve=solve)
    return opt,data
end

"""
    opt,data = cs_nctssos_first(supp::Vector{Vector{Vector{UInt16}}}, coe, n::Int, d::Int; numeq=0,
    CS="MF", TS="block", merge=false, md=3, QUIET=false, obj="eigen", solve=true)

Compute the first step of the CS-NCTSSOS hierarchy for constrained noncommutative polynomial optimization
with relaxation order `d`. Here the polynomial optimization problem is defined by `supp` and `coe`,
corresponding to the supports and coeffients of `pop` respectively. Return the optimum and other auxiliary data.

# Arguments
- `supp`: the supports of the polynomial optimization problem.
- `coe`: the coeffients of the polynomial optimization problem.
- `d`: the relaxation order of the moment-SOHS hierarchy.
- `numeq`: the number of equality constraints.
"""
function cs_nctssos_first(supp::Vector{Vector{Vector{UInt16}}}, coe, n::Int, d::Int; numeq=0, CS="MF",
    minimize=false, assign="first", TS="block", merge=false, md=3, QUIET=false, obj="eigen", solve=true)
    println("***************************NCTSSOS***************************")
    println("NCTSSOS is launching...")
    m = length(supp)-1
    dg = [maximum(length.(supp[i])) for i=2:m+1]
    if obj == "trace"
        supp[1],coe[1] = cyclic_canon(supp[1], coe[1])
    else
        supp[1],coe[1] = sym_canon(supp[1], coe[1])
    end
    time = @elapsed begin
    cliques,cql,cliquesize = clique_decomp(n, m, d, dg, supp, alg=CS, minimize=minimize)
    end
    if CS != false && QUIET == false
        mc = maximum(cliquesize)
        println("Obtained the variable cliques in $time seconds. The maximal size of cliques is $mc.")
    end
    if TS != false && QUIET == false
        println("Starting to compute the block structure...")
    end
    J,ncc = assign_constraint(m, supp, cliques, cql, cliquesize, assign=assign)
    time = @elapsed begin
    blocks,cl,blocksize,sb,numb,basis,status = get_cblocks_mix(d, dg, J, m, supp, cliques, cql, cliquesize, TS=TS, obj=obj)
    end
    if TS != false && QUIET == false
        mb = maximum(maximum.(sb))
        println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
    end
    opt,ksupp = blockcpop_mix(n, m, supp, coe, basis, cliques, cql, cliquesize, J, ncc, blocks, cl, blocksize, numeq=numeq, QUIET=QUIET, obj=obj, solve=solve)
    data = ncmpop_data(n, m, numeq, d, supp, coe, obj, ksupp, basis, cql, cliques, cliquesize, J, ncc, blocks, cl, blocksize, sb, numb)
    return opt,data
end

"""
    opt,data = cs_nctssos_higher!(data; TS="block", QUIET=false, merge=false, md=3, solve=true)

Compute higher steps of the CS-NCTSSOS hierarchy.
Return the optimum and other auxiliary data.
"""
function cs_nctssos_higher!(data::ncmpop_data; TS="block", QUIET=false, merge=false, md=3, solve=true)
    n = data.n
    m = data.m
    numeq = data.numeq
    d = data.d
    supp = data.supp
    coe = data.coe
    obj = data.obj
    ksupp = data.ksupp
    basis = data.basis
    cql = data.cql
    cliques = data.cliques
    cliquesize = data.cliquesize
    J = data.J
    ncc = data.ncc
    blocks = data.blocks
    cl = data.cl
    blocksize = data.blocksize
    sb = data.sb
    numb = data.numb
    if m == 0
        time = @elapsed begin
        blocks,cl,blocksize,sb,numb,basis,status = get_blocks_mix(d, supp, cliques, cql, cliquesize, basis=basis, sb=sb, numb=numb, TS=TS, merge=merge, md=md, obj=obj)
        end
        if status == 1
            if QUIET == false
                mb = maximum(maximum.(sb))
                println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
            end
            opt,ksupp = blockupop_mix(n, supp, coe, basis, cliques, cql, cliquesize, blocks, cl, blocksize, obj=obj, solve=solve, QUIET=QUIET)
        end
    else
        time = @elapsed begin
        blocks,cl,blocksize,sb,numb,basis,status = get_cblocks_mix(d, [], J, m, supp, cliques, cql, cliquesize, ksupp=ksupp, basis=basis, blocks=blocks, cl=cl, blocksize=blocksize, sb=sb, numb=numb, TS=TS, merge=merge, md=md, obj=obj)
        end
        if status==1
            if QUIET == false
                mb = maximum(maximum.(sb))
                println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
            end
            opt,ksupp = blockcpop_mix(n, m, supp, coe, basis, cliques, cql, cliquesize, J, ncc, blocks, cl, blocksize, numeq=numeq, QUIET=QUIET, obj=obj, solve=solve)
        end
    end
    if status == 0
        opt = nothing
        println("No higher CS-NCTSSOS hierarchy!")
    end
    data.ksupp = ksupp
    data.blocks = blocks
    data.cl = cl
    data.blocksize = blocksize
    data.sb = sb
    data.numb = numb
    return opt,data
end

function blockupop_mix(n, supp, coe, basis, cliques, cql, cliquesize, blocks, cl, blocksize; QUIET=false, obj="eigen", solve=true)
    ksupp = Vector{UInt16}[]
    for i = 1:cql, j = 1:cl[i], k = 1:blocksize[i][j], r = k:blocksize[i][j]
        @inbounds bi = [basis[i][blocks[i][j][k]][end:-1:1]; basis[i][blocks[i][j][r]]]
        push!(ksupp, bi)
    end
    ksupp = reduce!.(ksupp, obj=obj)
    sort!(ksupp)
    unique!(ksupp)
    lksupp = length(ksupp)
    if QUIET == false
        println("There are $lksupp affine constraints.")
    end
    objv = nothing
    if solve == true
        if QUIET == false
            println("Assembling the SDP...")
        end
        model = Model(optimizer_with_attributes(Mosek.Optimizer))
        set_optimizer_attribute(model, MOI.Silent(), QUIET)
        time = @elapsed begin
        cons = [AffExpr(0) for i=1:lksupp]
        pos = Vector{Vector{Union{VariableRef,Symmetric{VariableRef}}}}(undef, cql)
        for i = 1:cql
            pos[i] = Vector{Union{VariableRef,Symmetric{VariableRef}}}(undef, cl[i])
            for k = 1:cl[i]
                if blocksize[i][k] == 1
                   pos[i][k] = @variable(model, lower_bound=0)
                   @inbounds bi = [basis[i][blocks[i][k][1]][end:-1:1]; basis[i][blocks[i][k][1]]]
                   bi = reduce!(bi, obj=obj)
                   Locb = ncbfind(ksupp, lksupp, bi)
                   @inbounds add_to_expression!(cons[Locb], pos[i][k])
                else
                   pos[i][k] = @variable(model, [1:blocksize[i][k], 1:blocksize[i][k]], PSD)
                   for j = 1:blocksize[i][k], r = j:blocksize[i][k]
                       @inbounds ind1 = blocks[i][k][j]
                       @inbounds ind2 = blocks[i][k][r]
                       @inbounds bi = [basis[i][ind1][end:-1:1]; basis[i][ind2]]
                       bi = reduce!(bi, obj=obj)
                       Locb = ncbfind(ksupp, lksupp, bi)
                       if j == r
                           @inbounds add_to_expression!(cons[Locb], pos[i][k][j,r])
                       else
                           @inbounds add_to_expression!(cons[Locb], 2, pos[i][k][j,r])
                       end
                   end
                end
            end
        end
        bc = zeros(lksupp)
        for i = 1:length(supp)
            Locb = ncbfind(ksupp, lksupp, supp[i])
            if Locb == 0
               @error "The monomial basis is not enough!"
               return nothing,nothing
            else
               bc[Locb] = coe[i]
            end
        end
        @variable(model, lower)
        cons[1] += lower
        @constraint(model, cons.==bc)
        @objective(model, Max, lower)
        end
        if QUIET == false
            println("SDP assembling time: $time seconds.")
            println("Solving the SDP...")
        end
        time=@elapsed begin
        optimize!(model)
        end
        if QUIET == false
            println("SDP solving time: $time seconds.")
        end
        status = termination_status(model)
        objv = objective_value(model)
        if status != MOI.OPTIMAL
           println("termination status: $status")
           status = primal_status(model)
           println("solution status: $status")
        end
        println("optimum = $objv")
    end
    return objv,ksupp
end

function blockcpop_mix(n, m, supp, coe, basis, cliques, cql, cliquesize, J, ncc, blocks, cl, blocksize; numeq=0, QUIET=false, obj="eigen", solve=true)
    ksupp = Vector{UInt16}[]
    for i = 1:cql
        for j = 1:cl[i][1], k = 1:blocksize[i][1][j], r = k:blocksize[i][1][j]
            @inbounds bi = [basis[i][1][blocks[i][1][j][k]][end:-1:1]; basis[i][1][blocks[i][1][j][r]]]
            push!(ksupp, bi)
        end
        for (j, w) in enumerate(J[i])
            for l = 1:cl[i][j+1], t = 1:blocksize[i][j+1][l], r = t:blocksize[i][j+1][l], s = 1:length(supp[w+1])
                ind1 = blocks[i][j+1][l][t]
                ind2 = blocks[i][j+1][l][r]
                @inbounds bi = [basis[i][j+1][ind1][end:-1:1]; supp[w+1][s]; basis[i][j+1][ind2]]
                push!(ksupp, bi)
            end
        end
    end
    for i ∈ ncc
        append!(ksupp, supp[i+1])
    end
    ksupp = reduce!.(ksupp, obj=obj)
    sort!(ksupp)
    unique!(ksupp)
    lksupp = length(ksupp)
    if QUIET == false
        println("There are $lksupp affine constraints.")
    end
    objv = nothing
    if solve == true
        if QUIET == false
            println("Assembling the SDP...")
        end
        model = Model(optimizer_with_attributes(Mosek.Optimizer))
        set_optimizer_attribute(model, MOI.Silent(), QUIET)
        time = @elapsed begin
        cons = [AffExpr(0) for i=1:lksupp]
        for i = 1:cql, l = 1:cl[i][1]
            if blocksize[i][1][l] == 1
               @inbounds pos = @variable(model, lower_bound=0)
               @inbounds bi = [basis[i][1][blocks[i][1][l][1]][end:-1:1]; basis[i][1][blocks[i][1][l][1]]]
               bi = reduce!(bi, obj=obj)
               Locb = ncbfind(ksupp, lksupp,bi)
               @inbounds add_to_expression!(cons[Locb], pos)
            else
               @inbounds bs = blocksize[i][1][l]
               @inbounds pos = @variable(model, [1:bs, 1:bs], PSD)
               for t = 1:bs, r = t:bs
                   @inbounds ind1 = blocks[i][1][l][t]
                   @inbounds ind2 = blocks[i][1][l][r]
                   @inbounds bi = [basis[i][1][ind1][end:-1:1]; basis[i][1][ind2]]
                   bi = reduce!(bi, obj=obj)
                   Locb = ncbfind(ksupp, lksupp, bi)
                   if t == r
                      @inbounds add_to_expression!(cons[Locb], pos[t,r])
                   else
                      @inbounds add_to_expression!(cons[Locb], 2, pos[t,r])
                   end
               end
            end
        end
        for k = 1:length(ncc)
            i = ncc[k]
            if i <= m-numeq
                pos = @variable(model, lower_bound=0)
            else
                pos = @variable(model)
            end
            for j = 1:length(supp[i+1])
                bi = reduce!(supp[i+1][j], obj=obj)
                Locb = ncbfind(ksupp, lksupp, bi)
                @inbounds add_to_expression!(cons[Locb], coe[i+1][j], pos)
            end
        end
        for i = 1:cql, (j, w) in enumerate(J[i])
            for l = 1:cl[i][j+1]
                bs = blocksize[i][j+1][l]
                if bs == 1
                    if j <= m-numeq
                        pos = @variable(model, lower_bound=0)
                    else
                        pos = @variable(model)
                    end
                    ind1 = blocks[i][j+1][l][1]
                    for s = 1:length(supp[w+1])
                        @inbounds bi = [basis[i][j+1][ind1][end:-1:1]; supp[w+1][s]; basis[i][j+1][ind1]]
                        bi = reduce!(bi, obj=obj)
                        Locb = ncbfind(ksupp, lksupp, bi)
                        @inbounds add_to_expression!(cons[Locb], coe[w+1][s], pos)
                    end
                else
                    if j <= m-numeq
                        pos = @variable(model, [1:bs, 1:bs], PSD)
                    else
                        pos = @variable(model, [1:bs, 1:bs], Symmetric)
                    end
                    for t = 1:bs, r = t:bs
                        ind1 = blocks[i][j+1][l][t]
                        ind2 = blocks[i][j+1][l][r]
                        for s = 1:length(supp[w+1])
                            @inbounds bi = [basis[i][j+1][ind1][end:-1:1]; supp[w+1][s]; basis[i][j+1][ind2]]
                            bi = reduce!(bi, obj=obj)
                            Locb = ncbfind(ksupp, lksupp, bi)
                            if t == r
                                @inbounds add_to_expression!(cons[Locb], coe[w+1][s], pos[t,r])
                            else
                                @inbounds add_to_expression!(cons[Locb], 2*coe[w+1][s], pos[t,r])
                            end
                        end
                    end
                end
            end
        end
        bc = zeros(lksupp)
        for i = 1:length(supp[1])
            Locb = ncbfind(ksupp, lksupp, supp[1][i])
            if Locb == 0
               @error "The monomial basis is not enough!"
               return nothing,nothing
            else
               bc[Locb] = coe[1][i]
            end
        end
        @variable(model, lower)
        cons[1] += lower
        @constraint(model, cons.==bc)
        @objective(model, Max, lower)
        end
        if QUIET == false
            println("SDP assembling time: $time seconds.")
            println("Solving the SDP...")
        end
        time=@elapsed begin
        optimize!(model)
        end
        if QUIET == false
            println("SDP solving time: $time seconds.")
        end
        status = termination_status(model)
        objv = objective_value(model)
        if status != MOI.OPTIMAL
           println("termination status: $status")
           status = primal_status(model)
           println("solution status: $status")
        end
        println("optimum = $objv")
    end
    return objv,ksupp
end

function get_blocks_mix(d, supp, cliques, cql, cliquesize; basis=[], sb=[], numb=[], TS="block", merge=false, md=3, obj="eigen")
    blocks = Vector{Vector{Vector{UInt16}}}(undef, cql)
    cl = Vector{UInt16}(undef, cql)
    blocksize = Vector{Vector{UInt16}}(undef, cql)
    status = ones(UInt8, cql)
    if isempty(basis)
        sb = Vector{Vector{UInt16}}(undef, cql)
        numb = Vector{Vector{UInt16}}(undef, cql)
        basis = Vector{Vector{Vector{UInt16}}}(undef, cql)
        flag = 1
    else
        flag = 0
    end
    for i = 1:cql
        nvar = cliquesize[i]
        ind = [issubset(supp[j], cliques[i]) for j=1:length(supp)]
        ksupp = copy(supp[ind])
        if flag == 1
            basis[i] = get_ncbasis(nvar, d, ind=cliques[i])
            if obj == "trace"
                append!(ksupp, [_cyclic_canon([basis[i][k][end:-1:1]; basis[i][k]]) for k=1:length(basis[i])])
            else
                append!(ksupp, [[basis[i][k][end:-1:1]; basis[i][k]] for k=1:length(basis[i])])
            end
            sort!(ksupp)
            unique!(ksupp)
            blocks[i],cl[i],blocksize[i],sb[i],numb[i],status[i] = get_ncblocks(ksupp, basis[i], TS=TS, obj=obj, QUIET=true, merge=merge, md=md)
        else
            blocks[i],cl[i],blocksize[i],sb[i],numb[i],status[i] = get_ncblocks(ksupp, basis[i], sb=sb[i], numb=numb[i], TS=TS, obj=obj, QUIET=true, merge=merge, md=md)
        end
    end
    return blocks,cl,blocksize,sb,numb,basis,maximum(status)
end

function get_cblocks_mix(d, dg, J, m, supp, cliques, cql, cliquesize; ksupp=[], basis=[], blocks=[], cl=[], blocksize=[], sb=[], numb=[], TS="block", merge=false, md=3, obj="eigen")
    if isempty(basis)
        blocks = Vector{Vector{Vector{Vector{UInt16}}}}(undef, cql)
        cl = Vector{Vector{UInt16}}(undef, cql)
        blocksize = Vector{Vector{Vector{UInt16}}}(undef, cql)
        sb = Vector{Vector{UInt16}}(undef, cql)
        numb = Vector{Vector{UInt16}}(undef, cql)
        basis = Vector{Vector{Vector{Vector{UInt16}}}}(undef, cql)
        ksupp = copy(supp[1])
        for i = 2:m+1
            append!(ksupp,  _sym_canon.(supp[i]))
        end
        sort!(ksupp)
        unique!(ksupp)
        flag = 1
    else
        flag = 0
    end
    status = ones(UInt8, cql)
    for i = 1:cql
        lc = length(J[i])
        nvar = cliquesize[i]
        ind = [issubset(ksupp[j], cliques[i]) for j=1:length(ksupp)]
        fsupp = copy(ksupp[ind])
        if flag == 1
            basis[i] = Vector{Vector{Vector{UInt16}}}(undef, lc+1)
            basis[i][1] = get_ncbasis(cliquesize[i], d, ind=cliques[i])
            for s = 1:lc
                basis[i][s+1] = get_ncbasis(nvar, d-ceil(Int, dg[J[i][s]]/2), ind=cliques[i])
            end
            blocks[i] = Vector{Vector{Vector{UInt16}}}(undef, lc+1)
            cl[i] = Vector{UInt16}(undef, lc+1)
            blocksize[i] = Vector{Vector{UInt16}}(undef, lc+1)
            sb[i] = Vector{UInt16}(undef, lc+1)
            numb[i] = Vector{UInt16}(undef, lc+1)
            if obj == "trace"
                append!(fsupp, [_cyclic_canon([basis[i][1][k][end:-1:1]; basis[i][1][k]]) for k=1:length(basis[i][1])])
            else
                append!(fsupp, [[basis[i][1][k][end:-1:1]; basis[i][1][k]] for k=1:length(basis[i][1])])
            end
            sort!(fsupp)
            unique!(fsupp)
            blocks[i],cl[i],blocksize[i],sb[i],numb[i],status[i] = get_nccblocks(lc, fsupp, supp[J[i].+1], basis[i], TS=TS, QUIET=true, merge=merge, md=md, obj=obj)
        else
            blocks[i],cl[i],blocksize[i],sb[i],numb[i],status[i] = get_nccblocks(lc, fsupp, supp[J[i].+1], basis[i], blocks=blocks[i], cl=cl[i], blocksize=blocksize[i], sb=sb[i], numb=numb[i], TS=TS, QUIET=true, merge=merge, md=md, obj=obj)
        end
    end
    return blocks,cl,blocksize,sb,numb,basis,maximum(status)
end

function assign_constraint(m, supp, cliques, cql, cliquesize; assign="first")
    J = [UInt16[] for i=1:cql]
    ncc = UInt16[]
    for i = 2:m+1
        rind = copy(supp[i][1])
        for j = 2:length(supp[i])
            append!(rind, supp[i][j])
        end
        rind = unique(rind)
        if assign == "first"
            ind = findfirst(k->issubset(rind, cliques[k]), 1:cql)
            if ind != nothing
                push!(J[ind], i-1)
            else
                push!(ncc, i-1)
            end
        else
            temp = UInt16[]
            for j = 1:cql
                if issubset(rind, cliques[j])
                    push!(temp,j)
                end
            end
            if !isempty(temp)
                if assign == "min"
                    push!(J[temp[argmin(cliquesize[temp])]], i-1)
                else
                    push!(J[temp[argmax(cliquesize[temp])]], i-1)
                end
            else
                push!(ncc, i-1)
            end
        end
    end
    return J,ncc
end

function clique_decomp(n::Int, supp::Vector{Vector{UInt16}}; alg="MF", minimize=false)
    if alg == false
        cliques = [UInt16[i for i=1:n]]
        cql = 1
        cliquesize = [n]
    else
        G = SimpleGraph(n)
        for j = 1:length(supp)
            add_clique!(G, unique(supp[j]))
        end
        if alg == "NC"
            cliques,cql,cliquesize = max_cliques(G)
        else
            cliques,cql,cliquesize = chordal_cliques!(G, method=alg, minimize=minimize)
        end
    end
    uc = unique(cliquesize)
    sizes = [sum(cliquesize.== i) for i in uc]
    println("------------------------------------------------------")
    println("The clique sizes of varibles:\n$uc\n$sizes")
    println("------------------------------------------------------")
    return cliques,cql,cliquesize
end

function clique_decomp(n::Int, m::Int, d::Int, dg::Vector{Int}, supp::Vector{Vector{Vector{UInt16}}}; alg="MF", minimize=false)
    if alg == false
        cliques = [UInt16[i for i=1:n]]
        cql = 1
        cliquesize=[n]
    else
        G = SimpleGraph(n)
        for i = 1:m+1
            if i == 1 || d == ceil(Int, dg[i-1]/2)
                for j = 1:length(supp[i])
                    add_clique!(G, unique(supp[i][j]))
                end
            else
                temp = copy(supp[i][1])
                for j = 2:length(supp[i])
                    append!(temp, supp[i][j])
                end
                add_clique!(G, unique(temp))
            end
        end
        if alg == "NC"
            cliques,cql,cliquesize = max_cliques(G)
        else
            cliques,cql,cliquesize = chordal_cliques!(G, method=alg, minimize=minimize)
        end
    end
    uc = unique(cliquesize)
    sizes = [sum(cliquesize.== i) for i in uc]
    println("------------------------------------------------------")
    println("The clique sizes of varibles:\n$uc\n$sizes")
    println("------------------------------------------------------")
    return cliques,cql,cliquesize
end


text:
module Test_selinv

include("../src/ldlt.jl")
include("../src/selinv.jl")

using Test
using Random
using SparseArrays
using LinearAlgebra

@testset "selinv_from_ldlt" begin
    function test(A,conj,adj)
        Ap,Ai,Av = A.colptr,A.rowval,A.nzval
        Fp,Fi = symbolic_ldlt(Ap,Ai)
        Fv = numeric_ldlt(Ap,Ai,Av,Fp,Fi; conj=conj)
        P = SparseMatrixCSC(size(A)...,Fp,Fi,ones(Bool,length(Fi)))
        Bv = selinv_from_ldlt(Fp,Fi,Fv; conj=conj)
        B = SparseMatrixCSC(size(A)...,Fp,Fi,Bv)
        Bref = inv(Matrix(A)) .* P
        @test B ≈ Bref
    end

    for T in (Float32,Float64,ComplexF32,ComplexF64)
        for (conj,adj) in ((conj,adjoint),(identity,transpose))
            @testset "n = 1" begin
                Random.seed!(1)
                Ap = [1,2]
                Ai = [1]
                Av = rand(T,1)
                A = SparseMatrixCSC(1,1,Ap,Ai,Av)
                A = 2I + A + adj(A)
                test(A,conj,adj)
            end

            @testset "n = 2" begin
                @testset "A = [1 0; 0 1]" begin
                    Random.seed!(1)
                    Ap = [1,2,3]
                    Ai = [1,2]
                    Av = rand(T,2)
                    A = SparseMatrixCSC(2,2,Ap,Ai,Av)
                    A = 2I + A + adj(A)
                    test(A,conj,adj)
                end
                @testset "A = [1 1; 1 0]" begin
                    Ap = [1,3,3]
                    Ai = [1,2]
                    Av = rand(T,2)
                    A = SparseMatrixCSC(2,2,Ap,Ai,Av)
                    A = 2I + A + adj(A)
                    test(A,conj,adj)
                end
                @testset "A = [1 1; 1 1]" begin
                    Ap = [1,3,4]
                    Ai = [1,2,2]
                    Av = rand(T,3)
                    A = SparseMatrixCSC(2,2,Ap,Ai,Av)
                    A = 2I + A + adj(A)
                    test(A,conj,adj)
                end
            end

            @testset "n in 3:100" begin
                @testset for seed = 1:100
                    Random.seed!(seed)
                    n = rand(3:100)
                    A = sprand(n,n,rand())
                    A = 2I + A + adj(A)
                    test(A,conj,adj)
                end
            end
        end
    end
end

end # module

max_stars_repo_name: ettersi/SparseFactorizations

max_issues_repo_name: ettersi/SparseFactorizations

max_forks_repo_name: ettersi/SparseFactorizations


text:
using Test: @testset, @test, @inferred
using RestrictedBoltzmannMachines: BinaryRBM, Binary, inputs_v_to_h, free_energy
using RestrictedBoltzmannMachines: visible, hidden, weights
using ConvolutionalRBMs: hankel_image, hankel_weight, hankel, output_size, out2in
using ConvolutionalRBMs: BinaryConvRBM, replicate

@testset "replicate" begin
    @test size(@inferred replicate(Binary(zeros(2,3)),5,7)) == (2,3,5,7)
end

@testset "hankel" begin
    channel_sz = (2,3)
    in_sz = (5,5)
    kernel_sz = (3,2)
    out_sz = output_size(kernel_sz, in_sz)
    batch_sz = (2,)
    hidden_sz = (2,3)
    w = randn(channel_sz..., kernel_sz..., hidden_sz...)
    v = randn(channel_sz..., in_sz..., batch_sz...)

    convrbm = BinaryConvRBM(zeros(channel_sz), zeros(hidden_sz), w)

    V = @inferred hankel_image(v, channel_sz, kernel_sz)
    @test size(V) == (channel_sz..., kernel_sz..., out_sz..., batch_sz...)
    for j in CartesianIndices(kernel_sz), k in CartesianIndices(out_sz)
        i = out2in(j, k)
        for b in CartesianIndices(batch_sz), c in CartesianIndices(channel_sz)
            @test V[c,j,k,b] == v[c,i,b]
        end
    end
    denserbm = BinaryRBM(zeros(channel_sz..., kernel_sz...), zeros(hidden_sz...), w)
    @test inputs_v_to_h(denserbm, V) ≈ inputs_v_to_h(convrbm, v)

    W = @inferred hankel_weight(w, channel_sz, in_sz)
    @test size(W) == (channel_sz..., in_sz..., hidden_sz..., out_sz...)
    for j in CartesianIndices(kernel_sz), k in CartesianIndices(out_sz)
        i = out2in(j, k)
        for μ in CartesianIndices(hidden_sz), c in CartesianIndices(channel_sz)
            @test W[c,i,μ,k] == w[c,j,μ]
        end
    end
    densrbm = BinaryRBM(zeros(channel_sz..., in_sz...), zeros(hidden_sz..., out_sz...), W)
    @test inputs_v_to_h(densrbm, v) ≈ inputs_v_to_h(convrbm, v)

    convrbm = BinaryConvRBM(randn(channel_sz...), randn(hidden_sz...), w)
    hankrbm = hankel(convrbm, in_sz)
    @test all(visible(hankrbm).θ .== visible(convrbm).θ)
    @test all(hidden(hankrbm).θ .== hidden(convrbm).θ)
    @test inputs_v_to_h(hankrbm, v) ≈ inputs_v_to_h(convrbm, v)
    @test free_energy(hankrbm, v) ≈ free_energy(convrbm, v)

    convrbm = BinaryConvRBM(randn(channel_sz...), randn(hidden_sz...), w; pad=3, dilation=(2,1), stride=2)
    hankrbm = hankel(convrbm, in_sz)
    @test all(visible(hankrbm).θ .== visible(convrbm).θ)
    @test all(hidden(hankrbm).θ .== hidden(convrbm).θ)
    @test inputs_v_to_h(hankrbm, v) ≈ inputs_v_to_h(convrbm, v)
    @test free_energy(hankrbm, v) ≈ free_energy(convrbm, v)
end

@testset "hankel 1D" begin
    C = 2
    N = 7
    J = 3
    K = only(output_size((J,), (N,)))
    M = 2
    w = randn(C, J, M)
    v = randn(C, N)

    convrbm = BinaryConvRBM(zeros(C), zeros(M), w)

    V = @inferred hankel_image(v, (C,), (J,))
    @test size(V) == (C, J, K)
    for j in 1:J, k in 1:K, c in 1:C
        i = out2in(CartesianIndex(j), CartesianIndex(k))
        @test V[c,j,k] == v[c,i]
    end
    denserbm = BinaryRBM(zeros(C, J), zeros(M), w)
    @test inputs_v_to_h(denserbm, V) ≈ inputs_v_to_h(convrbm, v)

    W = @inferred hankel_weight(w, (C,), (N,))
    @test size(W) == (C, N, M, K)
    for j in CartesianIndices(1:J), k in CartesianIndices(1:K)
        i = out2in(j, k)
        for μ in 1:M, c in 1:C
            @test W[c,i,μ,k] == w[c,j,μ]
        end
    end
    densrbm = BinaryRBM(zeros(C, N), zeros(M, K), W)
    @test inputs_v_to_h(densrbm, v) ≈ inputs_v_to_h(convrbm, v)

    convrbm = BinaryConvRBM(randn(C), randn(M), w)
    hankrbm = hankel(convrbm, (N,))
    @test all(visible(hankrbm).θ .== visible(convrbm).θ)
    @test all(hidden(hankrbm).θ .== hidden(convrbm).θ)
    @test inputs_v_to_h(hankrbm, v) ≈ inputs_v_to_h(convrbm, v)
    @test free_energy(hankrbm, v) ≈ free_energy(convrbm, v)

    convrbm = BinaryConvRBM(randn(C), randn(M), w; stride=3, pad=:same, dilation=2)
    hankrbm = hankel(convrbm, (N,))
    @test all(visible(hankrbm).θ .== visible(convrbm).θ)
    @test all(hidden(hankrbm).θ .== hidden(convrbm).θ)
    @test inputs_v_to_h(hankrbm, v) ≈ inputs_v_to_h(convrbm, v)
    @test free_energy(hankrbm, v) ≈ free_energy(convrbm, v)
end


text:
import DelayEmbeddings: AbstractDataset, minima, maxima
import StaticArrays: SVector, MVector

"""
    encode_as_bin(point, reference_point, edgelengths) → Vector{Int}

Encode a point into its integer bin labels relative to some `reference_point`
(always counting from lowest to highest magnitudes), given a set of box
`edgelengths` (one for each axis). The first bin on the positive side of
the reference point is indexed with 0, and the first bin on the negative
side of the reference point is indexed with -1.

See also: [`joint_visits`](@ref), [`marginal_visits`](@ref).

## Example

```julia
using Entropies

refpoint = [0, 0, 0]
steps = [0.2, 0.2, 0.3]
encode_as_bin(rand(3), refpoint, steps)
```
"""
function encode_as_bin end

function encode_as_bin(point::Vector{T}, reference_point, edgelengths) where {T <: Real}
    floor.(Int, (point .- reference_point) ./ edgelengths)
end

function encode_as_bin(point::SVector{dim, T}, reference_point, edgelengths) where {dim, T <: Real}
    floor.(Int, (point .- reference_point) ./ edgelengths)
end

function encode_as_bin(point::MVector{dim, T}, reference_point, edgelengths) where {dim, T <: Real}
    floor.(Int, (point .- reference_point) ./ edgelengths)
end

function encode_as_bin(points::Vector{T}, reference_point, edgelengths) where {T <: Union{Vector, SVector, MVector}}
    [encode_as_bin(points[i], reference_point, edgelengths) for i = 1:length(points)]
end

function encode_as_bin(points::AbstractDataset{dim, T}, reference_point, edgelengths) where {dim, T}
    [encode_as_bin(points[i], reference_point, edgelengths) for i = 1:length(points)]
end

"""
    joint_visits(points, binning_scheme::RectangularBinning) → Vector{Vector{Int}}

Determine which bins are visited by `points` given the rectangular binning
scheme `ϵ`. Bins are referenced relative to the axis minima, and are
encoded as integers, such that each box in the binning is assigned a
unique integer array (one element for each dimension).

For example, if a bin is visited three times, then the corresponding
integer array will appear three times in the array returned.

See also: [`marginal_visits`](@ref), [`encode_as_bin`](@ref).

# Example

```julia
using DelayEmbeddings, Entropies

pts = Dataset([rand(5) for i = 1:100]);
joint_visits(pts, RectangularBinning(0.2))
```
"""
function joint_visits(points, binning_scheme::RectangularBinning)
    axis_minima, box_edge_lengths = get_minima_and_edgelengths(points, binning_scheme)
    # encode points relative to axis minima, with boxes of fixed length
    encode_as_bin(points, axis_minima, box_edge_lengths)
end

"""
    marginal_visits(points, binning_scheme::RectangularBinning, dims) → Vector{Vector{Int}}

Determine which bins are visited by `points` given the rectangular binning
scheme `ϵ`, but only along the desired dimensions `dims`. Bins are referenced
relative to the axis minima, and are encoded as integers, such that each box
in the binning is assigned a unique integer array (one element for each
dimension in `dims`).

For example, if a bin is visited three times, then the corresponding
integer array will appear three times in the array returned.

See also: [`joint_visits`](@ref), [`encode_as_bin`](@ref).

# Example

```julia
using DelayEmbeddings, Entropies
pts = Dataset([rand(5) for i = 1:100]);

# Marginal visits along dimension 3 and 5
marginal_visits(pts, RectangularBinning(0.3), [3, 5])

# Marginal visits along dimension 2 through 5
marginal_visits(pts, RectangularBinning(0.3), 2:5)
```
"""
function marginal_visits(points, binning_scheme::RectangularBinning, dims)
    axis_minima, box_edge_lengths = get_minima_and_edgelengths(points, binning_scheme)
    dim = length(axis_minima)
    if length(dims) == 1
        dims = [dim]
    end
    if sort(collect(dims)) == sort(collect(1:dim))
        joint_visits(points, binning_scheme)
    else
        # encode point relative to axis minima, with boxes of fixed length
        [encode_as_bin(pt, axis_minima, box_edge_lengths)[dims] for pt in points]
    end
end

"""
    marginal_visits(joint_visits, dims) → Vector{Vector{Int}}

If joint visits have been precomputed using [`joint_visits`](@ref), marginal
visits can be returned directly without providing the binning again
using the `marginal_visits(joint_visits, dims)` signature.

See also: [`joint_visits`](@ref), [`encode_as_bin`](@ref).

# Example

```
using DelayEmbeddings, Entropies
pts = Dataset([rand(5) for i = 1:100]);

# First compute joint visits, then marginal visits along dimensions 1 and 4
jv = joint_visits(pts, RectangularBinning(0.2))
marginal_visits(jv, [1, 4])

# Marginals along dimension 2
marginal_visits(jv, 2)
```
"""
function marginal_visits(joint_visits, dims)
    [bin[dims] for bin in joint_visits]
end


text:
"""
macro
"""

export @model, @connect, @blk
export @simulate, @scope

mutable struct DefFunc
    name::Symbol
    params
    body
    DefFunc() = new(:f, [], [])
end

function _makeparams(params)::Expr
    args = map(params) do x
        if Meta.isexpr(x, :(=))
            Expr(:kw, x.args[1], Expr(:call, :signal, x.args[2]))
        elseif isa(x, Symbol)
            Expr(:kw, x, Expr(:call, :signal, Expr(:quote, x)))
        end
    end
    Expr(:parameters, args...)
end

function _deffunc(f::DefFunc)
    Expr(:function,
        Expr(:call, f.name, _makeparams(f.params)),
        Expr(:block, f.body...))
end

"""
int1 = Sim.IntegratorBlock(init = 1)
Sim.blk!(system, :int1, int1)
"""

macro model(name, desc)
    @assert Meta.isexpr(desc, :block)
    f = DefFunc()
    f.name = name
    s = gensym()
    push!(f.body, :(local $s = SystemBlock()))
    for line = desc.args
        if Meta.isexpr(line, :macrocall) && line.args[1] == Symbol("@parameter")
            append!(f.params, filter(x->isa(x, Union{Symbol,Expr}), line.args[2:end]))
        else
            push!(f.body, _modelexpand(s, line))
        end
    end
    push!(f.body, s)
    esc(_deffunc(f))
end

function _modelexpand(system, expr::Expr)
    if Meta.isexpr(expr, :macrocall) && expr.args[1] == Symbol("@blk")
        b = expr.args[3]
        v = expr.args[4:end]
        _blk(system, b, v)
    elseif Meta.isexpr(expr, :macrocall) && expr.args[1] == Symbol("@connect")
        _connect(expr.args[3])
    else
        args = Any[_expand(system, x) for x = expr.args]
        Expr(expr.head, args...)
    end
end

function _modelexpand(system, expr::Any)
    expr
end

"""
blk
"""

macro blk(system, blk, vars...)
    _blk(system, blk, vars)
end

function _blk(system, blk, vars)
    @assert Meta.isexpr(blk, :(=)) && isa(blk.args[1], Symbol) && Meta.isexpr(blk.args[2], :call)
    ex = []
    symbol = blk.args[1]
    qsymbol = Expr(:quote, symbol)
    blktype = blk.args[2].args[1]
    blkargs = blk.args[2].args[2:end]
    push!(ex, :($symbol = blk!($system, $qsymbol, $blktype($(blkargs...)))))
    for x = vars
        @assert Meta.isexpr(x, :call) && x.args[1] == :(:)
        if Meta.isexpr(x.args[2], :ref) && isa(x.args[2].args[1], Symbol)
            port = Expr(:ref, Expr(:., symbol, Expr(:quote, x.args[2].args[1])), x.args[2].args[2])
            var = x.args[3]
            qvar = Expr(:quote, var)
            push!(ex, :($var = assignvar!($system, $qvar, $port)))
        elseif isa(x.args[2], Symbol)
            port = Expr(:., symbol, Expr(:quote, x.args[2]))
            var = x.args[3]
            qvar = Expr(:quote, var)
            push!(ex, :($var = assignvar!($system, $qvar, $port)))
        end
    end
    Expr(:block, ex...)
end

"""
connect
"""

macro connect(con)
    _connect(con)
end

function _connect(con)
    if Meta.isexpr(con, :call) && con.args[1] == :(<=)
        src = _connectexpandoutport(con.args[3])
        dest = con.args[2]
        :(connect!($src, $dest))
    elseif Meta.isexpr(con, :call) && con.args[1] == :(=>)
        src = _connectexpandoutport(con.args[2])
        dest = con.args[3]
        :(connect!($src, $dest))
    else
        error("error")
    end
end

function _connectexpandoutport(x::Expr)
    if Meta.isexpr(x, :call)
        args = [_connectexpandoutport(u) for u = x.args[2:end]]
        Expr(:call, x.args[1], args...)
    else
        Expr(:call, :signal, x)
    end
#         args = [Expr(:call, :signal, u) for u = x.args[1:end]]
#         Expr(x.head, args...)
#     end
end

function _connectexpandoutport(x::Symbol)
    Expr(:call, :signal, x)
end

function _connectexpandoutport(x::Any)
    x
end

"""
simulate
"""

struct SimulationResult
    t
    graph
    odesolution
    f::Dict{Symbol,Any}
end

macro simulate(model, kargs...)
    vars = nothing
    tspan = nothing
    solver = nothing
    solverops = []
    for x = kargs
        if Meta.isexpr(x, :(=)) && x.args[1] == :scope
            vars = _parsescope(x)
        elseif Meta.isexpr(x, :(=)) && x.args[1] == :tspan
            tspan = _parsetspan(x)
        elseif Meta.isexpr(x, :(=)) && x.args[1] == :solver
            solver = x.args[2]
        elseif Meta.isexpr(x, :(=)) && x.args[1] == :solveroptions
            solverops = _parsesolver(x.args[2])
        end
    end
    if solver != nothing
        insert!(solverops, 1, solver)
    end
    @assert vars != nothing && tspan != nothing
    ex = quote
        n = length($vars)
        @assert n >= 1
        @assert all(x -> haskey($model.env, x), $vars)
        prob = odeproblem($model, tspan = $tspan)
        solution = DifferentialEquations.solve(prob, $(solverops...))
        vars = [scope($model, x, solution) for x = $vars]
        g = Plots.plot(vars, solution.t[1], solution.t[end], layout=(n,1), leg=false)
        result = SimulationResult(solution.t, g, solution, Dict{Symbol,Any}())
        for (i,s) = enumerate($vars)
            result.f[s] = vars[i]
        end
        result
    end
    esc(ex)
end

function _parsescope(scope)
    @assert Meta.isexpr(scope, :(=)) && scope.args[1] == :scope
    if Meta.isexpr(scope.args[2], :vect) || Meta.isexpr(scope.args[2], :tuple)
        args = [Expr(:quote, x) for x = scope.args[2].args]
        Expr(:vect, args...)
    else
        Expr(:vect, Exper(:quote, scope.args[2].args))
    end
end

function _parsetspan(tspan)
    @assert Meta.isexpr(tspan, :(=)) && tspan.args[1] == :tspan
    if (Meta.isexpr(tspan.args[2], :vect) || Meta.isexpr(tspan.args[2], :tuple)) && length(tspan.args[2].args) == 2
        Expr(:tuple, tspan.args[2].args[1], tspan.args[2].args[2])
    else
        Expr(:tuple, 0, tspan.args[2])
    end
end

function _parsesolver(ops)
    @assert Meta.isexpr(ops, :vect) || Meta.isexpr(ops, :tuple)
    map(ops.args) do x
        @assert Meta.isexpr(x, :(=))
        Expr(:kw, x.args[1], x.args[2])
    end
end

"""
scope
"""

macro scope(system, var, solution)
    esc(:(scope($system, $(Expr(:quote, var)), $solution)))
end


text:
struct Scaling{Dim,T} <: Transform
  vec::Point{Dim,T}
end
Scaling(vals::T...) where {T} = Scaling(Point{length(vals),T}(vals))
Scaling(vec::AbstractVector) where {Dim,T} = Scaling(Point{length(vec),eltype(vec)}(vec))
(s::Scaling)(p::AbstractVector) = s.vec .* p

Base.:∘(t1::Scaling, t2::Scaling) = Scaling(t1.vec .* t2.vec)

(≈)(x::Scaling, y::Scaling) = x.vec ≈ y.vec

inv(s::Scaling) = Scaling(inv.(s.vec))

const Scaled{O,Dim,T} = Transformed{O,Scaling{Dim,T}}
Scaled(obj::O, transf::Scaling{Dim,T}) where {O,Dim,T} = Scaled{O,Dim,T}(obj, transf)

Base.show(io::IO, s::Scaled{O,Dim,T}) where {O,Dim,T} = print(io, "Scaled{$O, $Dim, $T}($(s.obj), $(s.transf))")

struct UniformScaling{T} <: Transform
  factor::T
end
(s::UniformScaling)(p::AbstractVector) = s.factor .* p

Base.:∘(t1::UniformScaling, t2::UniformScaling) = UniformScaling(t1.factor * t2.factor)

inv(s::UniformScaling) = UniformScaling(inv(s.factor))


text:
using BinDeps
@BinDeps.setup

using Compat

version = "5.13"
url = "https://cache.julialang.org/ftp://ftp.atnf.csiro.au/pub/software/wcslib/wcslib-$version.tar.bz2"

wcs = library_dependency("libwcs")
provides(Sources, URI(url), wcs, unpacked_dir="wcslib-$version")
depsdir = BinDeps.depsdir(wcs)
builddir = joinpath(depsdir, "src/wcslib-$version")
prefix = joinpath(depsdir, "usr")
configopts = ["--disable-fortran", "--without-cfitsio", "--without-pgplot",
              "--disable-utils"]
if is_apple()
    libfilename = "libwcs.$version.dylib"
elseif is_unix()
    libfilename = "libwcs.so.$version"
end

provides(BuildProcess,
         (@build_steps begin
            GetSources(wcs)
            @build_steps begin
              ChangeDirectory(builddir)
              FileRule(joinpath(prefix,"lib",libfilename),
                       @build_steps begin
                         `./configure --prefix=$prefix $configopts`
                         `make install`
                         @static(if is_apple()
                            `ln -s $prefix/lib/libwcs.5.dylib $prefix/lib/libwcs.dylib`
                            end)
                       end)
            end
          end),
         wcs)

@BinDeps.install Dict(:libwcs => :libwcs)