text:
export TUDataset

using DataDeps
# using ..MLDatasets: datafile, datadir
using DelimitedFiles: readdlm

function __init__tudataset()
    DEPNAME = "TUDataset"
    LINK = "https://www.chrsmrrs.com/graphkerneldatasets"
    DOCS = "https://chrsmrrs.github.io/datasets/docs/home/"
    DATA = "PROTEINS.zip"

    register(DataDep(
        DEPNAME,
        """
        Dataset: The $DEPNAME dataset.
        Website: $LINK)
        """,
        "$LINK/$DATA",
        # "81de017067dc045ebdb8ffd5c0e69a209973ffdb1fe2d5b434e94d3614f3f5c7",  # if checksum omitted, will be generated by DataDeps
        post_fetch_method = unpack
    ))
end

struct TUDataset
    num_nodes::Int
    num_edges::Int
    num_graphs::Int
    source::Vector{Int}
    target::Vector{Int}
    graph_indicator
    node_labels::Union{Nothing, Vector{Int}}
    edge_labels::Union{Nothing, Vector{Int}}
    graph_labels
    node_attributes
    edge_attributes
    graph_attributes
end

"""
    TUDataset(name; dir=nothing)

A variety of graph benchmark datasets, *.e.g.* "QM9", "IMDB-BINARY",
"REDDIT-BINARY" or "PROTEINS", collected from the [TU Dortmund University](https://chrsmrrs.github.io/datasets/).
Retrieve from TUDataset collection the dataset `name`, where `name`
is any of the datasets available [here](https://chrsmrrs.github.io/datasets/docs/datasets/).

A `TUDataset` object can be indexed to retrieve a specific graph or a subset of graphs.

# Internal fields

```
num_nodes           # total number of nodes (considering all graphs)
num_edges           # total number of edges (considering all graphs)
num_graphs          # total number of graphs
source              # vector of edges' source vectors
target              # vector of edges' target vectors
graph_indicator     # graph to which a node belongs too
node_labels
edge_labels
graph_labels
node_attributes
edge_attributes
graph_attributes
```

See [here](https://chrsmrrs.github.io/datasets/docs/format/) for an in-depth
description of the format.

# Usage Example

```julia
using MLDatasets: TUDataset
using LightGraphs: SimpleGraph, add_edge!

data = TUDataset("PROTEINS")

# Access first graph
d1 = data[1]

# Create a LightGraphs' graph
g = SimpleGraph(d1.num_nodes)
for (s, t) in zip(d1.source, d1.target)
    add_edge!(g, s, t)
end

# Node features
X = d1.node_attributes # (nfeatures x nnodes) matrix
```
"""
function TUDataset(name; dir=nothing)
    d = datadir_tudataset(name, dir)
    # See here for the file format: https://chrsmrrs.github.io/datasets/docs/format/

    st = readdlm(joinpath(d, "$(name)_A.txt"), ',', Int)
    # Check that the first node is labeled 1.
    # TODO this will fail if the first node is isolated
    @assert minimum(st) == 1
    source, target = st[:,1], st[:,2]

    graph_indicator = readdlm(joinpath(d, "$(name)_graph_indicator.txt"), Int) |> vec
    @assert all(sort(unique(graph_indicator)) .== 1:length(unique(graph_indicator)))

    num_nodes = length(graph_indicator)
    num_edges = length(source)
    num_graphs = length(unique(graph_indicator))

    # LOAD OPTIONAL FILES IF EXIST

    node_labels = isfile(joinpath(d, "$(name)_node_labels.txt")) ?
                    readdlm(joinpath(d, "$(name)_node_labels.txt"), Int) |> vec :
                    nothing
    edge_labels = isfile(joinpath(d, "$(name)_edge_labels.txt")) ?
                    readdlm(joinpath(d, "$(name)_edge_labels.txt"), Int) |> vec :
                    nothing
    graph_labels = isfile(joinpath(d, "$(name)_graph_labels.txt")) ?
                    readdlm(joinpath(d, "$(name)_graph_labels.txt"), Int) |> vec :
                    nothing

    node_attributes = isfile(joinpath(d, "$(name)_node_attributes.txt")) ?
                        readdlm(joinpath(d, "$(name)_node_attributes.txt"), ',', Float32)' |> collect :
                        nothing
    edge_attributes = isfile(joinpath(d, "$(name)_edge_attributes.txt")) ?
                        readdlm(joinpath(d, "$(name)_edge_attributes.txt"), ',', Float32)' |> collect :
                        nothing
    graph_attributes = isfile(joinpath(d, "$(name)_graph_attributes.txt")) ?
                        readdlm(joinpath(d, "$(name)_graph_attributes.txt"), ',', Float32)' |> collect :
                        nothing

    TUDataset( num_nodes, num_edges, num_graphs,
                source, target,
                graph_indicator,
                node_labels,
                edge_labels,
                graph_labels,
                node_attributes,
                edge_attributes,
                graph_attributes)
end

function datadir_tudataset(name, dir = nothing)
    dir = isnothing(dir) ? datadep"TUDataset" : dir
    LINK = "https://www.chrsmrrs.com/graphkerneldatasets/$name.zip"
    d  = joinpath(dir, name)
    if !isdir(d)
        DataDeps.fetch_default(LINK, dir)
        currdir = pwd()
        cd(dir) # Needed since `unpack` extracts in working dir
        DataDeps.unpack(joinpath(dir, "$name.zip"))
        cd(currdir)
    end
    @assert isdir(d)
    return d
end

Base.getindex(data::TUDataset, i::Int) = getindex(data, [i])

function Base.getindex(data::TUDataset, i::AbstractVector{Int})
    node_mask = data.graph_indicator .∈ Ref(i)

    nodes = (1:data.num_nodes)[node_mask]
    node_labels = isnothing(data.node_labels) ? nothing : data.node_labels[node_mask]
    nodemap = Dict(v => vnew for (vnew, v) in enumerate(nodes))

    graphmap = Dict(i => inew for (inew, i) in enumerate(i))
    graph_indicator = [graphmap[i] for i in data.graph_indicator[node_mask]]

    edge_mask = data.source .∈ Ref(nodes)
    source = [nodemap[i] for i in data.source[edge_mask]]
    target = [nodemap[i] for i in data.target[edge_mask]]
    edge_labels = isnothing(data.edge_labels) ? nothing : data.edge_labels[edge_mask]

    graph_labels = isnothing(data.graph_labels) ? nothing : data.graph_labels[i]

    node_attributes = isnothing(data.node_attributes) ? nothing : data.node_attributes[:,node_mask]
    edge_attributes = isnothing(data.edge_attributes) ? nothing : data.edge_attributes[:,edge_mask]
    graph_attributes = isnothing(data.graph_attributes) ? nothing : data.graph_attributes[:,i]

    num_nodes = length(graph_indicator)
    num_edges = length(source)
    num_graphs = length(i)

    TUDataset(num_nodes, num_edges, num_graphs,
            source, target,
            graph_indicator,
            node_labels,
            edge_labels,
            graph_labels,
            node_attributes,
            edge_attributes,
            graph_attributes)
end


text:
include(joinpath(@__DIR__, "config.jl"))
include(joinpath(@__DIR__, "..", "subspace_util.jl"))

using OneClassActiveLearning, SVDD, JuMP, Gurobi, MLKernels

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)

data_dirs = ["HeartDisease", "Stamps", "Pima", "Cardiotocography", "SpamBase", "Annthyroid", "PageBlocks"]
data_info = "03-baselines"

### learning scenario ###
initial_pool_strategy = [("Pu", Dict())]
split_strategy = [("Sf", Dict())]
classify_precision = SVDD.OPT_PRECISION

#### models ####

models = [Dict(:type => :SSAD, :param => Dict{Symbol, Any}(:κ => 0.1)),
          Dict(:type => :SVDDneg, :param => Dict{Symbol, Any}())]

init_strategies = [SimpleCombinedStrategy(SVDD.WangGammaStrategy(solver), BoundedTaxErrorEstimate(0.05, 0.02, 0.98))]

subspace_generators = [1]

#### oracle ####
oracle_param = Dict{Symbol, Any}(
    :type => PoolOracle,
    :param => Dict{Symbol, Any}()
)

#### query strategies ####

query_strategies = [
    Dict(:type => :RandomOutlierPQs, :param => Dict{Symbol, Any}()),
    Dict(:type => :HighConfidencePQs, :param => Dict{Symbol, Any}()),
    Dict(:type => :DecisionBoundaryPQs, :param => Dict{Symbol, Any}())]

num_al_iterations = 50
num_resamples_initial_pool = 1
num_resamples_subspaces = 1
data_output_root = joinpath(data_root, "output", "evaluation-part1")

ext: jl

max_stars_repo_name: holtri/subsvdd-evaluation

max_issues_repo_name: holtri/subsvdd-evaluation

max_forks_repo_name: holtri/subsvdd-evaluation


text:
#############################################################################
# ECOS.jl
# Wrapper around the ECOS solver https://github.com/ifa-ethz/ecos
# See http://github.com/JuliaOpt/ECOS.jl
#############################################################################
# test/runtests.jl
# Test the ECOS.jl solver wrapper
#############################################################################

tests = ["direct.jl",
         "options.jl",
         "mpb_linear.jl"]

println("Running tests:")

for curtest in tests
    println(" Test: $(curtest)")
    include(curtest)
end

# Run the conic interface test from MathProgBase.jl
import ECOS
include(joinpath(Pkg.dir("MathProgBase"),"test","conicinterface.jl"))
coniclineartest(ECOS.ECOSSolver(), duals=true)
conicSOCtest(ECOS.ECOSSolver(), duals=true)
conicEXPtest(ECOS.ECOSSolver(), duals=true)

include(joinpath(Pkg.dir("MathProgBase"),"test","quadprog.jl"))
socptest(ECOS.ECOSSolver())


text:
const MISCELLANEOUS_GPU_RELATED_REPOS = MirrorUpdater.Types.SrcDestPair[
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaGPU/docker",
        destination_repo_name = "JuliaGPU-docker",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/JuliaGPU/gitlab-ci",
        destination_repo_name = "JuliaGPU-gitlab-ci",
        ),
    ]


text:
using Revise
using JDF
using CSV, DataFrames
using Serialization:serialize,deserialize

@time a = gf()

@time metadatas = savejdf("c:/data/a_pre_compress.jdf", a);

# before = Base.summarysize(a)
#
@time type_compress!(a; verbose = true)
#
# after = Base.summarysize(a)
#
# before/after
# 2+2

@time savejdf("c:/data/a.jdf", a);

using Revise
using DataFrames, JDF
@time a2 = loadjdf("c:/data/a.jdf");

@time a2 = sloadjdf("c:/data/a.jdf");

type_compress!(a2)

@time savejdf("c:/data/a_p.jdf", a2)
@time ssavejdf("c:/data/a_p.jdf", a2)

#a2 = ploadjdf("c:/data/a_p.jdf")

all(names(a) .== names(a2))
all(skipmissing([all(a2[!,name] .== Array(a[!,name])) for name in names(a2)]))


text:
module QuBase

    using Compat

    ####################
    # String Constants #
    ####################
        const otimes = "\u2297"
        const vdots ="\u205E"

    ##################
    # Abstract Types #
    ##################
        abstract AbstractStructure
        abstract Orthogonal <: AbstractStructure
        abstract Orthonormal <: Orthogonal
        abstract AbstractQuantum{S<:AbstractStructure}

        # Various constructor methods in this repo allow an argument
        # of type Type{BypassFlag} to be passed in in order to
        # circumvent value precalculation/checking. This is useful for
        # conversion methods and the like, where you know the input
        # has already been vetted elsewhere. Don't use this unless
        # you're sure of what you're doing, and don't export this.
        abstract BypassFlag

    #############
    # Functions #
    #############
        # This should be the only `structure`
        # method that needs to be defined for
        # type *instances*
        structure{S}(::AbstractQuantum{S}) = S
        structure(::DataType) = AbstractStructure

        # ...and all relevant singleton types
        # should have it defined as well:
        structure{S}(::Type{AbstractQuantum{S}}) = S

        # an n-arity form of the tensor
        # product, reduction is done via
        # the binary definition of tensor()
        # defined in the files included above.
        tensor(s...) = reduce(tensor, s)

    ######################
    # Include Statements #
    ######################
        include("bases/bases.jl")
        include("arrays/arrays.jl")

    export AbstractStructure,
        AbstractQuantum,
        Orthogonal,
        Orthonormal,
        structure,
        tensor
end


text:
load_mesh = load
save_mesh = save
function create_group_df(data_folder)
    fname1= joinpath(data_folder,"cor_gc_gda.csv")
    fname2= joinpath(data_folder,"gda.tmdf")

    cor_gc_gda=CSV.read(fname1)
    gda = CSV.read(fname2)[!,[:cellid,:trackid,:x,:y,:z,:t]]
    n_sample =size(gda,1)
    rename!(cor_gc_gda,[:cellid,:t,:gc_cellid])
    # df = join(gda,cor_gc_gda,on=[:t,:cellid],kind=:inner)
    df = innerjoin(gda,cor_gc_gda,on=[:t,:cellid])

    dropmissing!(df)
    df[!,:t] = convert.(Float64,df[!,:t])
    return df
end


text:
using Printf

function print_header(ws::COSMO.Workspace)

	m, n = ws.p.model_size
	settings = ws.settings
	settings.scaling > 0 ? scaling_status = "on" : scaling_status = "off"
	nnz_in_P = count(!iszero,ws.p.P) - count(!iszero,diag(ws.p.P)) + n
	nnz_in_M = 2*count(!iszero,ws.p.A) + nnz_in_P + m
	println("-"^66 * "\n" * " "^10 * "COSMO v0.5.0 - A Quadratic Objective Conic Solver\n" * " "^25 * "Michael Garstka\n"  * " "^16 * "University of Oxford, 2017 - 2019\n" * "-"^66 * "\n")

	println("Problem:  x ∈ R^{$(n)},\n          constraints: A ∈ R^{$(m)x$(n)} ($(count(!iszero, ws.p.A)) nnz),\n          matrix size to factor: $(n + m)x$(n + m) ($(nnz_in_M) nnz)")
	for (iii, set) in enumerate(sort(ws.p.C.sets, by = x -> -x.dim))
		set_name = split(string(typeof(set).name), ".")[end]
		if isa(set, PsdCone) || isa(set, PsdConeTriangle) || isa(set, PsdConeTriangleLOBPCG)
			set_dim = string(set.sqrt_dim, "x", set.sqrt_dim)
		else
			set_dim = string(set.dim)
		end
		iii == 1 ? println("Sets:"*" "^5*"$(set_name) of dim: $(set_dim)") : println(" "^10*"$(set_name) of dim: $(set_dim)")
		if iii > 5
			println(" "^10*"... and $(length(ws.p.C.sets)-5) more")
			break
		end
	end
	if ws.settings.decompose && ws.ci.num_decomposable > 0
		println("Decomp:   Num of original PSD cones: $(ws.ci.num_psd_cones)\n" * " "^10*"Num decomposable PSD cones: $(ws.ci.num_decomposable)\n" * " "^10*"Num PSD cones after decomposition: $(ws.ci.num_decom_psd_cones)")
	end
	println("Settings: ϵ_abs = $(@sprintf("%.1e", settings.eps_abs)), ϵ_rel = $(@sprintf("%.1e", settings.eps_rel)),\n" * " "^10 * "ϵ_prim_inf = $(@sprintf("%.1e", settings.eps_prim_inf)), ϵ_dual_inf = $(@sprintf("%.1e", settings.eps_dual_inf)),\n" * " "^10 * "ρ = $(settings.rho), σ = $(settings.sigma), α = $(settings.alpha),\n" * " "^10 * "max_iter = $(settings.max_iter),\n" * " "^10 * "scaling iter = $(settings.scaling) ($(scaling_status)),\n" * " "^10 * "check termination every $(settings.check_termination) iter,\n" * " "^10 * "check infeasibility every $(settings.check_infeasibility) iter,\n" *	 " "^10 * "KKT system solver: $(print_lin_sys(settings.kkt_solver.ObjectType))")

	println("Setup Time: $(round.(ws.times.setup_time*1000; digits=2))ms\n")
	if settings.psd_projector == PsdConeTriangleLOBPCG || (isa(settings.psd_projector, OptionsFactory) && settings.psd_projector.ObjectType == PsdConeTriangleLOBPCG)
		println("Iter:\tObjective:\tPrimal Res:\tDual Res:\tRho:\t\tTime(s):\tSubspaces:\tLOBPCG avg it:\tExact avg it:")
	else
		println("Iter:\tObjective:\tPrimal Res:\tDual Res:\tRho:\t\tTime(s):")
	end
	flush(stdout)
	nothing
end

function print_iteration(ws::COSMO.Workspace, iter::Int64, cost::Float64, r_prim::Float64, r_dual::Float64, running_time::Float64)
	settings = ws.settings
	if mod(iter, 1) == 0 || iter == 1 || iter == 2 || iter == settings.max_iter
		if mod(iter, settings.check_termination) == 0
			lobpcg_dims_str, lobpcg_iterations, exact_projections = extract_lobpcg_statistics(ws)
			if isnan(lobpcg_iterations)
				@printf("%d\t%.4e\t%.4e\t%.4e\t%.4e\t%.2e\n", iter, cost, r_prim, r_dual, ws.ρ, running_time)
			else
				@printf("%d\t%.4e\t%.4e\t%.4e\t%.4e\t%.2e\t%s\t%.2f\t\t%.2f\n", iter, cost, r_prim, r_dual, ws.ρ, running_time,
					lobpcg_dims_str, lobpcg_iterations, exact_projections)
			end
		else
			@printf("%d\t%.4e\t ---\t\t\t---\t%.2e\n", iter, cost, running_time)
		end
	end
	flush(stdout)
	nothing
end

function extract_lobpcg_statistics(ws; max_printed_dimensions=3)
	lobpcg_iterations, exact_projections = 0, 0
	lobpcg_dims = Int[]
	sets_counter = 0
	for set in ws.p.C.sets
		if isa(set, COSMO.PsdConeTriangleLOBPCG)
			if length(lobpcg_dims) < max_printed_dimensions
				lobpcg_dim = set.subspace_dim_history[end]
				if set.is_subspace_positive
					append!(lobpcg_dims, lobpcg_dim)
				else
					append!(lobpcg_dims, -lobpcg_dim)
				end
			end

			lobpcg_iterations += set.lobpcg_iterations
			exact_projections += set.exact_projections
			set.lobpcg_iterations = 0
			set.exact_projections = 0
			sets_counter += 1
		end
	end
	lobpcg_dims_str = string(lobpcg_dims)
	if length(lobpcg_dims) < sets_counter
		lobpcg_dims_str = string(lobpcg_dims_str[1:end-1], ", ...]")
	end
	lobpcg_dims_str = rpad(lobpcg_dims_str, 12, " ")
	return lobpcg_dims_str, lobpcg_iterations/sets_counter, exact_projections/sets_counter
end

function print_result(status::Symbol, iter::Int64, cost::Float64, rt::Float64)
	println("\n" * "-"^66 * "\n>>> Results\nStatus: $(status)\nIterations: $(iter)\nOptimal objective: $(round.(cost; digits = 4))\nRuntime: $(round.(rt; digits = 3))s ($(round.(rt * 1000; digits = 2))ms)\n")
	nothing
end

function print_lin_sys(s::Type{<:AbstractKKTSolver})
	lin_sys_dict = Dict("QdldlKKTSolver" => "QDLDL", "CholmodKKTSolver" => "CHOLMOD", "PardisoDirectKKTSolver" => "Pardiso (direct)", "PardisoIndirectKKTSolver" => "Pardiso (indirect)",  "MKLPardisoKKTSolver" => "MKL Pardiso")
	return get(lin_sys_dict, string(s), string(s))
end


text:
##### Momentum tendencies

using CLIMAParameters.Planet: Omega

#####
##### First order fluxes
#####

function flux(::Advect{Momentum}, atmos, args)
    @unpack state = args
    return state.ρu .* (state.ρu / state.ρ)'
end

function flux(::PressureGradient{Momentum}, atmos, args)
    @unpack state, aux = args
    @unpack ts = args.precomputed
    pad = (state.ρu .* (state.ρu / state.ρ)') * 0
    if atmos.ref_state isa HydrostaticState
        return pad + (air_pressure(ts) - aux.ref_state.p) * I
    else
        return pad + air_pressure(ts) * I
    end
end

#####
##### Second order fluxes
#####

struct ViscousStress{PV <: Momentum} <: TendencyDef{Flux{SecondOrder}, PV} end
function flux(::ViscousStress{Momentum}, atmos, args)
    @unpack state, aux, t, diffusive = args
    pad = (state.ρu .* (state.ρu / state.ρ)') * 0
    ν, D_t, τ = turbulence_tensors(atmos, state, diffusive, aux, t)
    return pad + τ * state.ρ
end

function flux(::MoistureDiffusion{Momentum}, atmos, args)
    @unpack state, aux, t, diffusive = args
    ν, D_t, τ = turbulence_tensors(atmos, state, diffusive, aux, t)
    d_q_tot = (-D_t) .* diffusive.moisture.∇q_tot
    return d_q_tot .* state.ρu'
end

#####
##### Sources
#####

export Gravity
struct Gravity{PV <: Momentum} <: TendencyDef{Source, PV} end
Gravity() = Gravity{Momentum}()
function source(s::Gravity{Momentum}, m, args)
    @unpack state, aux = args
    if m.ref_state isa HydrostaticState
        return -(state.ρ - aux.ref_state.ρ) * aux.orientation.∇Φ
    else
        return -state.ρ * aux.orientation.∇Φ
    end
end

export Coriolis
struct Coriolis{PV <: Momentum} <: TendencyDef{Source, PV} end
Coriolis() = Coriolis{Momentum}()
function source(s::Coriolis{Momentum}, m, args)
    @unpack state = args
    FT = eltype(state)
    _Omega::FT = Omega(m.param_set)
    # note: this assumes a SphericalOrientation
    return -SVector(0, 0, 2 * _Omega) × state.ρu
end

export GeostrophicForcing
struct GeostrophicForcing{PV <: Momentum, FT} <: TendencyDef{Source, PV}
    f_coriolis::FT
    u_geostrophic::FT
    v_geostrophic::FT
end
GeostrophicForcing(::Type{FT}, args...) where {FT} =
    GeostrophicForcing{Momentum, FT}(args...)
function source(s::GeostrophicForcing{Momentum}, m, args)
    @unpack state, aux = args
    u_geo = SVector(s.u_geostrophic, s.v_geostrophic, 0)
    ẑ = vertical_unit_vector(m, aux)
    fkvector = s.f_coriolis * ẑ
    return -fkvector × (state.ρu .- state.ρ * u_geo)
end

export RayleighSponge
"""
    RayleighSponge{PV <: Momentum, FT} <: TendencyDef{Source, PV}

Rayleigh Damping (Linear Relaxation) for top wall momentum components
Assumes laterally periodic boundary conditions for LES flows. Momentum components
are relaxed to reference values (zero velocities) at the top boundary.
"""
struct RayleighSponge{PV <: Momentum, FT} <: TendencyDef{Source, PV}
    "Maximum domain altitude (m)"
    z_max::FT
    "Altitude at with sponge starts (m)"
    z_sponge::FT
    "Sponge Strength 0 ⩽ α_max ⩽ 1"
    α_max::FT
    "Relaxation velocity components"
    u_relaxation::SVector{3, FT}
    "Sponge exponent"
    γ::FT
end
RayleighSponge(::Type{FT}, args...) where {FT} =
    RayleighSponge{Momentum, FT}(args...)
function source(s::RayleighSponge{Momentum}, m, args)
    @unpack state, aux = args
    z = altitude(m, aux)
    if z >= s.z_sponge
        r = (z - s.z_sponge) / (s.z_max - s.z_sponge)
        β_sponge = s.α_max * sinpi(r / 2)^s.γ
        return -β_sponge * (state.ρu .- state.ρ * s.u_relaxation)
    else
        FT = eltype(state)
        return SVector{3, FT}(0, 0, 0)
    end
end

max_stars_repo_name: jm-c/CLIMA

max_issues_repo_name: jm-c/CLIMA

max_forks_repo_name: jm-c/CLIMA


text:
include("types.jl")
include("accessors.jl")
include("rendering.jl")


text:
export Kronecker

immutable Kronecker <: GraphGenerator
	scale::Int64
	edgefactor::Int64
    a::Float64
    b::Float64
    c::Float64

    function Kronecker(
        scale::Int64,
        edgefactor::Int64;
        a::Float64=0.57,
        b::Float64=0.19,
        c::Float64 = 0.19
        )
        new(scale, edgefactor, a, b, c)
    end
end

function generate(generator::Kronecker)
    #Ported from https://gitorious.org/graph500/graph500?p=graph500:graph500.git;a=blob;f=octave/kronecker_generator.m;h=064b6d4f42194cc95ec4ab375dcbee7c3243b233;hb=01a32a9c16994c3cf30bbf4caabf9be128379416.
    n = 2 ^ generator.scale
    m = generator.edgefactor * n
    ij = ones(Int64, 2, m)
    ab = generator.a + generator.b
    cnorm = generator.c/(1-ab)
    anorm = generator.a/ab

    for ib in 1:generator.scale
        ii_bit = rand(1, m) .> ab
        jj_bit = rand(1, m) .> ( cnorm * ii_bit + anorm * !(ii_bit) )
        ij = ij .+ 2^(ib-1) .* [ii_bit; jj_bit]
    end

    p = randperm(n)
    ij = p[ij]

    p = randperm(m)
    ij = ij[:, p]

    return ij
end


text:
@info "JLD2.jl is available - including functionality to read / write JLD2 files"

import .JLD2

function read_table(::TableIOInterface.JLD2Format, filename:: AbstractString, tablename:: AbstractString; kwargs...)
    return JLD2.jldopen(filename, "r") do file
        file[tablename]
    end
end

function read_table(::TableIOInterface.JLD2Format, filename:: AbstractString; kwargs...)
    return JLD2.jldopen(filename, "r") do file
        table_list = keys(file) |> sort
        file[first(table_list)]
    end
end

function write_table!(::TableIOInterface.JLD2Format, filename:: AbstractString, tablename:: AbstractString, table; kwargs...)
    _checktable(table)
    _checktablename(tablename)
    JLD2.jldopen(filename, "a+") do file
        file[tablename] = table
    end
    nothing
end

function list_tables(::TableIOInterface.JLD2Format, filename:: AbstractString)
    return JLD2.jldopen(filename, "r") do file
        keys(file) |> sort
    end
end


text:
"""
    _get_titles_from_CSV(csv::String) -> Vector{String}

Extract the values of the `Title` column from the CSV file.

# Arguments
- `csv::String`: $(CSV_ARGUMENT)

# Returns
- `Vector{String}`: the titles

# Throws
- Exceptions from: [`get_all`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights
HEADER = PDFHighlights.Internal.Constants.HEADER

_file, io = mktemp()
println(io, HEADER, '\\n', ",\\"Girl, Interrupted\\",,,1")
flush(io)
file = _file * ".csv"
mv(_file, file)

PDFHighlights.Internal.CSV._get_titles_from_CSV(file) == ["Girl, Interrupted"]

# output

true
```
"""
function _get_titles_from_CSV(csv::String)::Vector{String}
    _, titles, _, _, _ = get_all(csv)
    return titles
end


text:
"""
    oneline(s)

Remove newlines, replacing them with spaces when it seems appropriate.
"""
function oneline(v::AbstractString)
    R = split(v, '\n')
    buf = IOBuffer()
    for (i, r) in enumerate(R)
        if i != lastindex(R) && (r[end] == '-' || occursin(' ', r))
            print(buf, r, " ")
        else
            print(buf, r)
        end
    end
    String(take!(buf))
end
oneline(v::Any) = v

function findbreak(v, t)
    i = t == 59 ?
        findlast('\n', v) :
        1
    j = i+t-1
    j >= lastindex(v) && return nothing
    x = findlast(c -> c == ' ' || c == '-', v[1:j])
    if isnothing(x)
        return j:j+1
    elseif x == ' '
        return x-1:x+1
    else
        return x:x+1
    end
end

"""
    multiline(v, s::Symbol)

Return a `String` with newlines and spaces added so that it conforms to the GenBank line width.
"""
function multiline(v, s)
    v = string(v)
    s = string(s)
    if length(v) + length(s) > 55 && !occursin('\n', v)
        b = findbreak(v, 55 - length(s))
        v = v[1:b.start] * "\n" * v[b.stop:end]
        b = findbreak(v, 59)
        while !isnothing(b) && b.start < length(v)
            v = v[1:b.start] * "\n" * v[b.stop:end]
            b = findbreak(v, 59)
        end
    end
    return v
end

"""
    readgbk(input)

Parse GenBank-formatted file, returning a `Vector{GenBank.Record}`. File names ending in ".gz" are assumed to be gzipped and are decompressed.
"""
function readgbk(input)
    collect(open(GenBank.Reader, input))
end

"""
    readgff(input)

Parse GFF3-formatted file, returning a `Vector{Record}`. File names ending in ".gz" are assumed to be gzipped and are decompressed.
"""
function readgff(input)
    collect(open(GFF.Reader, input))
end


text:
using Symbolics
using SatisfiabilityInterface

vars = @variables x, y, z, w

constraints = [
    x ∈ 1:2
    y ∈ 2:5
    z ∈ 3:7
    z == x * y
]

prob = DiscreteCSP(constraints)

solve(prob)

all_solutions(prob)

###

vars = @variables x, y, z, w

constraints = [
    x ∈ 1:2
    y ∈ 2:5
    z ∈ 3:7
    w ∈ 4:6
    z == x + y
    w == z + x
]

prob = DiscreteCSP(constraints)
solve(prob)
all_solutions(prob)

###

constraints = [
    x ∈ 1:2
    y ∈ 2:5
    z ∈ 3:7
    z == 2x
]

all_solutions(DiscreteCSP(constraints))

###

vars = @variables x, y, z

# constraints = [
#     x ∈ 1:4
#     y ∈ 2:5
#     z == x + y
#     z ≤ 4
# ]

### Subset sum Problem

# https://en.wikipedia.org/wiki/Subset_sum_problem

# Given a set of integers like {−7, −3, −2, 5, 8},
# do a subset sum to 0?

s = [-7, -3, -2, 5, 8]
n = length(s)
x = [Num(Variable(:x, i)) for i in 1:n]

constraints = [
    [x[i] ∈ (0:1) for i in 1:n]
    sum(x[i] * s[i] for i in 1:n) == 0
    sum(x) ≠ 0   # at least one is not 0
    ]

prob = DiscreteCSP(constraints)

all_solutions(prob)

#####

@variables x, y, z

constraints = [
    x ∈ 1:2
    y ∈ 2:5
    z ∈ 3:7
    z == 2x
]

prob = DiscreteCSP(constraints)
all_solutions(prob)

## Example from https://github.com/newptcai/BeeEncoder.jl

function make_vector(name::Symbol, range)
    return Num.(Variable.(name, range))
end

@variables x, y, z, w
b = make_vector(:b, 1:4)
# @variables b[1:4]

¬(x::Num) = 1 - x   # assumes a Boolean variable

constraints = [
    x ∈ 0:5
    y ∈ -4:9
    z ∈ -5:10
    w ∈ 0:10

    [bb ∈ 0:1 for bb ∈ b]

    x + y == z

    b[1] == ¬b[2]
    b[2] == 1

    ¬b[1] + b[2] + ¬b[3] + b[4] == w
]

prob = DiscreteCSP(constraints)
solve(prob)

# @time length(all_solutions(prob))   # 1184 solutions in 38 seconds

###

constraints = [
    x ∈ -5:5
    y ∈ -4:9
    x - y^2 == 3
    # sum(b) == w

]

prob = DiscreteCSP(constraints)
all_solutions(prob)

###

@variables x, y

constraints = [
    x ∈ -5:5
    y ∈ -4:9
    10x + y == 9
    # sum(b) == w
]

prob = DiscreteCSP(constraints)
all_solutions(prob)


text:
abstract type Pat end
struct StarPat <: Pat end
struct CstrPat <: Pat
	cstr::Symbol
	args::Vector{Pat}
end
struct VarPat <: Pat
	var::Symbol
end
struct OrPat <: Pat
	left::Pat
	right::Pat
end

iswildcard(::StarPat) = true
iswildcard(::Pat) = false

heads(::StarPat) = Set{Symbol}()
heads(c::CstrPat) = Set([c.cstr])
heads(o::OrPat) = union(heads(o.left), heads(o.right))

Base.:(==)(a::StarPat, b::StarPat) = true
Base.:(==)(a::CstrPat, b::CstrPat) = a.cstr == b.cstr && all(a.args .== b.args)
Base.:(==)(a::VarPat, b::VarPat) = a.var == b.var
Base.:(==)(a::OrPat, b::OrPat) = a.left == b.left && a.right == b.right

extract_pattern(sym::Any) =
	if sym == :(_)
		StarPat()
	else
		VarPat(sym)
	end
function extract_pattern(ex::Expr)
	if ex.head == :call
		if ex.args[1] == :(|)
			return OrPat(extract_pattern(ex.args[2]), extract_pattern(ex.args[3]))
		else
			arg_patterns = extract_pattern.(ex.args[2:end])
			nargs = length(arg_patterns)
			return CstrPat(ex.args[1], arg_patterns)
		end
	else
		throw("Invalid match expression $ex")
	end
end

function extract_case(expr::Expr)
	@assert(expr.head == :call)
	@assert(expr.args[1] == :(=>))
	@assert(length(expr.args) == 3)
	if expr.args[2].head == :tuple
		exprs = extract_pattern.(expr.args[2].args)
	else
		exprs = Pat[extract_pattern(expr.args[2])]
	end
	term = expr.args[3]
	return Pair{Vector{Pat}, Any}(exprs, term)
end
function extract_patterns(expr)
	if expr.head == :block
		return Pair{Vector{Pat}, Any}[extract_case(t) for t in expr.args if t isa Expr]
	elseif expr.head == :call && expr.args[1] == :(=>)
		return Pair{Vector{Pat}, Any}[extract_case(expr)]
	end
end


text:
using BEE
using Suppressor
using Test

function capture_render(c)
    @capture_out render(c)
end

function capture_render()
    @capture_out render()
end

@testset "BEE.jl" begin
    @testset "simple" begin
        BEE.reset()

        example="""new_int(w, 0, 10)
        new_int(x, 0, 5)
        new_int(z, -5, 10)
        new_int(y, -4, 9)
        new_bool(x1)
        new_bool(x4)
        new_bool(x2)
        new_bool(x3)
        int_plus(x, y, z)
        bool_eq(x1, -x2)
        bool_eq(x2, true)
        bool_array_sum_eq([-x1, x2, -x3, x4], w)
        solve satisfy
        """
        ret = @capture_out include("../example/simple-example.jl")
        @test ret == example

        @test "new_int(w, 0, 10)\n" == capture_render(w)

        c = x+y == z
        @test "int_plus(x, y, z)\n" == capture_render(c)

        c = xl[1] == -xl[2]
        @test "bool_eq(x1, -x2)\n" == capture_render(c)

        BEE.reset()
        @constrain xl[1] == -xl[2]
        @test "bool_eq(x1, -x2)\nsolve satisfy\n" == capture_render()

        BEE.reset()
        constrain(xl[1] == -xl[2])
        @test "bool_eq(x1, -x2)\nsolve satisfy\n" == capture_render()

    end

    @testset "Declaring Variable" begin
        BEE.reset()

        @beebool x1
        @test "new_bool(x1)\n" == capture_render(x1)

        x2 = beebool("x2")
        @test "new_bool(x2)\n" == capture_render(x2)

        x3 = beebool(:x3)
        @test "new_bool(x3)\n" == capture_render(x3)

        yl = @beebool y1 y2 y3
        for i in 1:3
            @test "new_bool(y$i)\n" == capture_render(yl[i])
        end

        zl = @beebool z[1:10]
        for i in 1:10
            @test "new_bool(z$i)\n" == capture_render(zl[i])
        end

        (a, b, cl, d) = @beebool a b c[1:10] d
        for i in 1:10
            @test "new_bool(c$i)\n" == capture_render(cl[i])
        end

        BEE.reset()

        @beeint xx 3 55

        @test "new_int(xx, 3, 55)\n" == capture_render(xx)

        il = @beeint i[1:10] 3 7
        for i in 1:10
            @test "new_int(i$i, 3, 7)\n" == capture_render(il[i])
        end
    end

    @testset "Boolean statements" begin
        BEE.reset()

        @beebool x1
        @beebool x2
        x3 = BeeBool("x3")

        c = x1 == x2
        @test "bool_eq(x1, x2)\n" == capture_render(c)

        c = -x1 == x2
        @test "bool_eq(-x1, x2)\n" == capture_render(c)

        c = true == and([x1, -x2, x3])
        @test "bool_array_and([x1, -x2, x3])\n" == capture_render(c)

        c = and([x1, x2]) == -x3
        @test "bool_array_and_reif([x1, x2], -x3)\n" == capture_render(c)

        c = -x3 == BEE.xor(-x1, x2)
        @test "bool_xor_reif(-x1, x2, -x3)\n" == capture_render(c)
    end

    @testset "Integer statements" begin
        BEE.reset()

        @beeint x1 3 7
        @beeint x2 4 6
        x3 = beeint("x3", 10, 15)

        c = (x1 < x2) == true
        @test "int_lt(x1, x2)\n" == capture_render(c)

        c = true == (x1 < x2)
        @test "int_lt(x1, x2)\n" == capture_render(c)

        c = false == (x1 < x2)
        @test "int_lt_reif(x1, x2, false)\n" == capture_render(c)

        c = (x1 < x2) == false
        @test "int_lt_reif(x1, x2, false)\n" == capture_render(c)

        c = alldiff([x1, x2, x3]) == true
        @test "int_array_allDiff([x1, x2, x3])\n" == capture_render(c)

        c = true == alldiff([x1, x2, x3])
        @test "int_array_allDiff([x1, x2, x3])\n" == capture_render(c)

        c = alldiff([x1, x2]) == x3
        @test "int_array_allDiff_reif([x1, x2], x3)\n" == capture_render(c)

        c = x3 == alldiff([x1, x2])
        @test "int_array_allDiff_reif([x1, x2], x3)\n" == capture_render(c)

        c = x1 * x2 == x3
        @test "int_times(x1, x2, x3)\n" == capture_render(c)

        c = x3 == x1 * x2
        @test "int_times(x1, x2, x3)\n" == capture_render(c)

        c = x1 + x2 == 33
        @test "int_plus(x1, x2, 33)\n" == capture_render(c)

        c = 33 == x1 + x2
        @test "int_plus(x1, x2, 33)\n" == capture_render(c)

        c = max([x1, x2]) == x3
        @test "int_array_max([x1, x2], x3)\n" == capture_render(c)

        c = x3 == max([x1, x2])
        @test "int_array_max([x1, x2], x3)\n" == capture_render(c)
    end

    @testset "Cardinality" begin
        BEE.reset()

        xl = [beebool("x$i") for i in 1:3]

        il = [beeint("i$i", 3, 5) for i in 1:3]

        c = sum(il) == il[1]
        @test "int_array_sum_eq([i1, i2, i3], i1)\n" == capture_render(c)

        c = sum(xl) >= il[1]
        @test "bool_array_sum_geq([x1, x2, x3], i1)\n" == capture_render(c)

        c = sum(xl) > il[1]
        @test "bool_array_sum_gt([x1, x2, x3], i1)\n" == capture_render(c)
    end

    @testset "Boolean arrays relation" begin
        BEE.reset()

        xl = [beebool("x$i") for i in 1:3]

        il = [beebool("i$i") for i in 1:3]

        @beebool b

        c = (xl == il) == true
        @test "bool_arrays_eq([x1, x2, x3], [i1, i2, i3])\n" == capture_render(c)

        c = (xl != il) == true
        @test "bool_arrays_neq([x1, x2, x3], [i1, i2, i3])\n" == capture_render(c)

        c = (xl <= il) == true
        @test "bool_arrays_lex([x1, x2, x3], [i1, i2, i3])\n" == capture_render(c)

        c = (xl < il) == true
        @test "bool_arrays_lexLt([x1, x2, x3], [i1, i2, i3])\n" == capture_render(c)

        c = (xl <= il) == b
        @test "bool_arrays_lex_reif([x1, x2, x3], [i1, i2, i3], b)\n" == capture_render(c)

        c = (xl < il) == b
        @test "bool_arrays_lexLt_reif([x1, x2, x3], [i1, i2, i3], b)\n" == capture_render(c)
    end
end


text:
export yet_anthother_a, yet_anthother_b,yet_anthother_c

function yet_anthother_a()
	1
end

yet_anthother_b= 5

yet_anthother_c(x)=x

type yet_anthother_d <: AbstractString
	s
end

max_stars_repo_name: oxinabox/julia-vim-omnicomplete

max_issues_repo_name: oxinabox/julia-vim-omnicomplete

max_forks_repo_name: oxinabox/julia-vim-omnicomplete


text:
# В данном файле собраны служебные функции, необходимые для:
# 1. обработки входных аргументов вызова скрипта;
# 2. чтения больших файлов в заданном формате с возможным пропуском закомментированных или ненужных строк  заголовка;
# 3. вывод числа типа Int с добавлением необходимого числа нулей в начале;
# 4. функция логирования событий.

# Функции обработки входящего аргумента скрипта, задающего набор степеней полиномов для вывода результатов аппроксимации.
# Данная функция обрабатывает параметр, заключенный в квадратные скобки
# и на основании указанных данных формирует массив элементов
# Пример:
# на вход подается строка "[1,2,3,4,5:5:20,30,40]"
# на выходе формируется массив Array{Int}[1,2,3,4,5,10,15,20,30,40]

function make_numberlist(arg_string)

    arg_split = split(strip(arg_string, ['[',']']), ",")

    arg_deg = Any[]
    for arg in arg_split
        s_arg = split(arg, ":")
        if length(s_arg) >1
            if length(s_arg) == 2
                st = parse(Int, s_arg[1])
                en = parse(Int, s_arg[2])
                push!(arg_deg, st:en)
            end
            if length(s_arg) == 3
                st = parse(Int, s_arg[1])
                sp = parse(Int, s_arg[2])
                en = parse(Int, s_arg[3])
                push!(arg_deg, st:sp:en)
            end
        else
            push!(arg_deg, parse(Int, arg))
        end
    end
    degs = Int[]
    for els in arg_deg
        if typeof(els) <: Integer
            push!(degs, els)
        else
            for el_in_collect in els
                push!(degs, el_in_collect)
            end
        end
    end
    return degs
end

# Данная функция предназначена для формирования набора окон для обработки
# На вход функции подается строка аргументов с параметрами окна и максимальный размер данных, на которые формируются окна
# Строка аргументов может быть задана тремя способами:
# M:N - одним окном;
# M:W:N - окнами с перекрытием в одну точку;
# M:W:S:N - окнами с заданным перекрытием.
# При этом всюду выше:
# M - начальная точка первого окна,
# N - последняя точка последнего окна (в качестве последней точки
# можно указать конец файла, для этого используется ключевое слово end),
# W - размер окна,
# S - перекрытие соседних окон.
# !!!ВАЖНОЕ ЗАМЕЧАНИЕ - 1!!!
# При формировании окон последнее окно формируется так, чтобы
# последняя точка окна совпадала с последней заданной точкой N,
# при этом перекрытие с предпоследним окном может быть больше, чем указано в исходных параметрах
# !!!ВАЖНОЕ ЗАМЕЧАНИЕ - 2!!!
# Если размер окна задан больше, чем количество доступных данных, то скрипт выдает ошибку в лог-файл и завершает работу
# **********************************************************# Примеры
# "1:10", 20 => [(1,10)]
# "1:10:19", 18 => [(1, 10), (9, 18)]
# "1:10:19", 20 => [(1, 10), (10, 19)]
# "1:10:20", 20 => [(1, 10), (10, 19), (11, 20)]
# "1:10:3:20", 20 => [(1, 10), (8, 17), (11, 20)]

function make_windows(arg_string, data_length)
    windows_array = Any[]
    arg = split(arg_string, ':')
    if length(arg) == 2 # Вариант, когда используется 1 окно
        start_index = parse(Int, arg[1])
        if arg[2] == "end"
            end_index = data_length
        else
            end_index =  parse(Int, arg[2])
            if end_index>data_length
                return windows_array
            end
        end
        windows_array = [(start_index, end_index)]
    elseif length(arg) == 3 # окна c перекрытием в 1 точку
        start_index = parse(Int, arg[1])
        step_size = parse(Int, arg[2])
        if arg[3] == "end"
            end_index = data_length
        else
            end_index =  parse(Int, arg[3])
            if end_index>data_length
                return windows_array
            end
        end
        for start_window = start_index:step_size-1:end_index - step_size
            push!(windows_array, (start_window, start_window + step_size-1))
        end
        if length(windows_array) ==0
            return windows_array
        end
        if windows_array[end][2] < end_index
            push!(windows_array, (end_index - step_size+1, end_index))
        end
    elseif length(arg) == 4 # окна с перекрытием
        start_index = parse(Int, arg[1])
        step_size = parse(Int, arg[2])
        overlap = parse(Int, arg[3])
        if arg[4] == "end"
            end_index = data_length
        else
            end_index =  parse(Int, arg[4])
            if end_index>data_length
                println("return")
                return windows_array
            end
        end
        println("read param: start$start_index - step$step_size - overlap$overlap - end$end_index")
        for start_window = start_index:(step_size - overlap):(end_index - step_size)
            push!(windows_array, (start_window, start_window + step_size-1))
        end
        if length(windows_array) ==0
            return windows_array
        end
        if windows_array[end][2] < end_index
            push!(windows_array, (end_index - step_size + 1, end_index))
        end
    else
        windows_array = [(1, data_length)]
    end
    return windows_array
end

# функция для чтения данных
function readbigfile(filename::AbstractString,           # имя входного файла
    Columns::Array{Int,1},     # номера стобцов, которые необходимо считать. Если указать 0 как один из номеров для считывания, то будет сформирован столбец с номерами строк (пропущенные при чтении строки не считаются)
    T::Type                    # тип данных, в который будут сконвертированы считанные столбцы
    ;
    skip_first_row = 0        # число строк в начале файла, которые необходимо пропустить
    )::Array{T,2}             # функция возвращает двумерный массив

    # Открываем файл для чтения
    in_file = readlines(filename)
    skip_row = 0
    for line in in_file
        if skip_first_row > skip_row
            skip_row=skip_row+1
            continue
        end
        # пропускаем строки, начинающиеся (В ПЕРВОЙ ПОЗИЦИИ!) c #
        if line[1] == '#'
            skip_row = skip_row+1
            continue
        end
        # if length(split(line)) < length(Columns)
        #     skip_row = skip_row+1
        #     continue
        # end
    end
    println("skipped rows: $skip_row")
    # Создаем нулевой двумерный массив типа T
    out_array = zeros(T, (length(in_file) - skip_row, length(Columns)))
    # Построчно обрабатываем файл
    num_out_array_row = 1
    for num_str = 1+skip_first_row:length(in_file)
        if in_file[num_str][1] == '#'
            continue
        end
#  пропускаем строки с недостаточным числом столбцов
        # if length(split(in_file[num_str])) < length(Columns)
        #     continue
        # end
        value_str = split(in_file[num_str])
        for num_col = 1:length(Columns)
            if Columns[num_col] == 0
                out_array[num_out_array_row,num_col] = num_out_array_row
            else
                out_array[num_out_array_row,num_col] = parse(T, value_str[Columns[num_col]])
            end
        end
        num_out_array_row = num_out_array_row + 1
    end
    # Возвращаем считанный массив
    return out_array
end

# Данная функция принимает на вход два числа:
# 1. целое положительное число, которое необходимо напечатать
# 2. минимальное число символов в выводимом числе.
# Если символов в числе меньше указанного, перед числом добавляются нули так, чтобы число с нулями стало необходимой минимальной длины.
# Если число длиннее, выводится просто указанное число
function print_with_zero(num::Integer, min_length_num::Integer)
    str_num = "$num"
    zero_in_begin = "0"^maximum([0,(min_length_num-length(str_num))])
    return "$(zero_in_begin)$(str_num)"
end

# Функция дописывает сообщение, указанное вторым аргументом функции в файл с именем указанным как первый аргумент функции,
# при этом указывается время записи
function add_log_messange(log_name, message)
    out_file = open(log_name, "a")
    println(out_file, "[$(now())] $(message)")
    flush(out_file)
    close(out_file)
end

# Функция для вывода двумерных и одномерных массивов с указанием разделителя между элементами
# На вход подается входной массив и могут быть указаны разделители между элементами строки (col_delims)
# и между строками (row_delims),
# если не указывать - по умолчанию между элементами табуляция и между строками перенос строки
# также может быть указан поток вывода, по умолчанию STDOUT
function println_array(InputArray::Array;col_delims = "\t", row_delims="\n", out_stream = STDOUT)
    dims = ndims(InputArray)
    if dims == 1
        println(out_stream,join(map(x->"$(x)",InputArray),col_delims))
    elseif dims == 2
        num_rows = size(InputArray)[1]
        for ind_row = 1:num_rows
            print(out_stream,join(map(x->"$(x)",InputArray[ind_row, :]),col_delims))
            print(out_stream,row_delims)
        end
        println(out_stream)
    else
        println("warning: Dims > 2")
    end
end

max_stars_repo_name: sptsarev/high-deg-polynomial-fitting

max_issues_repo_name: sptsarev/high-deg-polynomial-fitting

max_forks_repo_name: sptsarev/high-deg-polynomial-fitting


text:
module DiametricalClustering
using StatsBase
function initialAssignment(M::Matrix, k::Integer)
    (m,n) = size(M)
    local indexes =  sample(1:k,m)
    return indexes
end
function clusteringError(G::Matrix{Float32},vectors::Vector{Vector{Float32}}, assignments, k::Integer)
    local cerror = 0
    for j=1:k
        local cluster = find(assignments.==j)
        for c in cluster
            cerror = cerror + dot(G[c,:], vectors[j])^2
        end
    end
    return  cerror
end

function initialCentroids(M::Matrix{<:AbstractFloat}, k::Integer)
    (d,n) =size(M)
    Mu = zeros(d,k)
    vectorMean = vec(mean(M,2))
    Mu0 = vectorMean/norm(vectorMean)
    perturbation = 0.1
    for i=1:k
        randomVector = rand(d)-0.5
        randomNorm   = perturbation*rand()
        smallRandVec =  randomNorm*randomVector/norm(randomVector)
        v            =  Mu0 + smallRandVec;
        Mu[:,i]      =  v/norm(v)
    end
    return Mu
end
function clusterize(X::Matrix{<:AbstractFloat}, k::Integer; maxIterations::Integer = 5000, threshold = 0.00000001)
    local mu = initialCentroids(X,k)
    local assignments = nothing
    local diff    = 1;
    local epsilon = 0.001;
    local value   = 100;
    local iteration = 1;
    local oldvalue = 5
    while (value-oldvalue > epsilon && iteration<200)
        iteration = iteration+1;

        oldvalue  = value;
        # assign points to nearest cluster
        simMat          =  (mu'*X).^2;
        (simax,indexes) =  findmax(simMat,1);

        assignments   = ind2sub(size(simMat),vec(indexes))[1]

        #compute objective function value
        value         = sum(simax);
        # compute cluster centroids
        for h=1:k
            if length(find(assignments.==h)) > 0
                clustData = X[:,find(assignments.==h)];
                (svd,_)   = svds(clustData',nsv=1);
                V = svd[:V]
                mu[:,h]   = V';
            end
        end
    end
    return (assignments, mu)
end
end


text:
"""
    main_ptdf
A main file to generate a PTDF file

Parameters:
    input_path : path to input data directory describing a grid
                (not the pscopf_ files but branches.txt and buses.txt)
"""

using Dates

root_path = dirname(@__DIR__)
push!(LOAD_PATH, root_path);
cd(root_path)
include(joinpath(root_path, "src", "PTDF.jl"));

MATPOWER_NETWORKS = [
    "case118"
    ,"case1354pegase"
    ,"case13659pegase"
    ,"case14"
    ,"case145"
    ,"case1888rte"
    ,"case1951rte"
    ,"case2383wp"
    ,"case24_ieee_rts"
    ,"case2736sp"
    ,"case2737sop"
    ,"case2746wop"
    ,"case2746wp"
    ,"case2848rte"
    ,"case2868rte"
    ,"case2869pegase"
    ,"case30"
    ,"case300"
    ,"case3012wp"
    ,"case30pwl"
    ,"case30Q"
    ,"case3120sp"
    ,"case3375wp"
    ,"case39"
    ,"case4gs"
    ,"case5"
    ,"case57"
    ,"case6468rte"
    ,"case6470rte"
    ,"case6495rte"
    ,"case6515rte"
    ,"case6ww"
    ,"case89pegase"
    ,"case9"
    ,"case9241pegase"
    ,"case9Q"
    ,"case9target"
];

function compute_ptdf(input_path::String)
    output_path = joinpath(input_path, "pscopf_ptdf.txt")
    ref_bus_num = 1
    distributed = true

    #########################
    # EXECUTION
    #########################
    network = PTDF.read_network(input_path)
    ptdf = PTDF.compute_ptdf(network, ref_bus_num)
    if distributed
        ptdf = PTDF.distribute_slack(ptdf);
        # coeffs = Dict([ "poste_1_0" => .2,
        #                 "poste_2_0" => .8])
        # ptdf = PTDF.distribute_slack(ptdf, coeffs, network);
    end
    PTDF.write_PTDF(output_path, network, ptdf, distributed, ref_bus_num)

end

#########################
# INPUT & PARAMS
#########################

# input_path = ( length(ARGS) > 0 ? ARGS[1] : joinpath(@__DIR__, "..", "data_matpower", "case1354pegase") )
for input_path in MATPOWER_NETWORKS
    compute_ptdf(joinpath(@__DIR__, "..", "data_matpower", input_path));
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using LinearAlgebra
using StatsBase

"""
    shift_to_non_negative(matrix::Array)

Returns a matrix in which values are non-negative. This is done by finding the
minimal value in the input matrix and adding its absolute value to the matrix
elements.
"""
function shift_to_non_negative(matrix::Array)

    min_val = findmin(matrix)[1]
    if min_val < 0
        return matrix .-= min_val
    else
        return matrix
    end
end

"""
    normalize_to_01(matrix::Array; use_factor=false, norm_factor=256)

Returns a matrix which values are in range [0, 1]. If 'use_factor' is set to
'true' then values are normalized to 'norm_factor' (by default set to 256).

If the values in the input matrix are below 0, then they are shifted so that only positive numbers are in
the matrix (the values are normalized to new maximal value or norm_factor).
"""
function normalize_to_01(matrix::Array; use_factor = false, norm_factor = 256)
    normalized_matrix = copy(matrix)

    min_val = findmin(normalized_matrix)[1]
    if min_val < 0
        normalized_matrix .+= abs(min_val)
    else
        normalized_matrix .-= abs(min_val)
    end

    max_val = findmax(normalized_matrix)[1]

    if use_factor
        if max_val > norm_factor
            @warn "Maximal values exceed \'norm_factor\'."
        end
        normalized_matrix = normalized_matrix ./ norm_factor
    else
        normalized_matrix = normalized_matrix ./ max_val
    end

    return normalized_matrix
end

# function symmetrize_image(image)
"""
    function diagonal_symmetrize(image::Matrix; below_over_upper::Bool=false)

Takes an 'image' in the form of a matrix and return a copy which is symmetric
with respect to diagonal- values above diagonal are copied over values below the
diagonal. This can be inverted by setting 'below_over_upper=true'.

If the input matrix is not square, then square matrix is created by taking
matrix of k times 'k' elements, 'k=min(r,c)' where 'r' is number of rows and 'c'
is number of columns.
"""
function diagonal_symmetrize(image::Matrix; below_over_upper::Bool = false)
    w, h = size(image)
    mat_size = findmin([w, h])[1]

    img = copy(image[1:mat_size, 1:mat_size])

    # Get all cartesian indices from input matrix
    matrix_indices = CartesianIndices((1:mat_size, 1:mat_size))
    # Filter out indices below diagonal
    if below_over_upper
        matrix_indices = findall(x -> x[1] > x[2], matrix_indices)
    else
        matrix_indices = findall(x -> x[2] > x[1], matrix_indices)
    end

    # how many elements are above diagonal
    repetition_number = Int(ceil((mat_size * (mat_size - 1)) / 2))

    # Substitute elements
    for k = 1:repetition_number
        # n_pos = matrix_indices[k]
        mat_ind = matrix_indices[k]
        # ordered_matrix[mat_ind] = k
        img[mat_ind[2], mat_ind[1]] = img[mat_ind]
    end

    try
        checksquare(img)
    catch err
        if isa(err, DimensionMismatch)
            @error "Resulting matrix is not a square matrix"
            throw(err)
        end
    end
    # issymmetric(Float64.(img))
    return img
end

# =====
# matrix ordering

"""
    function get_ordered_matrix(in_matrix::Matrix;
                                assign_same_values::Bool = false,
                                force_symmetry::Bool = false,
                                small_dist_grouping::Bool = false,
                                min_dist::Number = 1e-16,
                                total_dist_groups::Int = 0,
                                ordering_start::Int=1)

Takes a @input_matrix and returns ordered form of this matrix.
The ordered form is a matrix which elements represent ordering from smallest to
highest values in @input_matrix.

If @input_matrix is symmetric, then ordering happens only with upper diagonal.
Lower diagonal is symetrically copied from values above diagonal.

By default, if there is a geoup of entriess with the same value, they all are
assigned with the same ordering number. This can be changed with
@assign_same_values parameter.

Symetry ordering can be froced with @force_symmetry parameter.

By setting 'small_dist_grouping' to true, all the values that difference is
lower than 'min_dist', will be assigned with the same order number.

# Examples
```julia-repl
julia> a = [0 11 12;
            11 0 13;
            12 13 0];
julia> get_ordered_matrix(a)
3×3 Array{Int64,2}:
 0  1  2
 1  0  3
 2  3  0
```

```julia-repl
julia> b = [38 37 36 30;
            37 34 30 32;
            36 30 31 30;
            30 32 30 29]
julia> get_ordered_matrix(b; assign_same_values=false)
4×4 Array{Int64,2}:
0  6  5  2
6  0  1  4
5  1  0  3
2  4  3  0

julia> get_ordered_matrix(b; assign_same_values=true)
4×4 Array{Int64,2}:
0  4  3  1
4  0  1  2
3  1  0  1
1  2  1  0
```
"""
function get_ordered_matrix(in_matrix::Matrix;
                            assign_same_values::Bool = false,
                            force_symmetry::Bool = false,
                            small_dist_grouping::Bool = false,
                            min_dist::Number = 1e-16,
                            total_dist_groups::Int = 0,
                            ordering_start::Int=1)

    # TODO Symmetry must be forced for matrix in which there are NaN elements- needs
    #   to be further investigated
    # TODO not working for negative only values
    # TODO check for square matrix

    # ==
    mat_size = size(in_matrix)
    ord_mat = zeros(Int, mat_size)

    # how many elements are above diagonal
    if issymmetric(in_matrix) || force_symmetry
        matrix_indices =
            generate_indices(mat_size, symmetry_order = true, include_diagonal = false)
        do_symmetry = true
    else
        matrix_indices = generate_indices(mat_size, symmetry_order = false)
        do_symmetry = false
    end
    total_elements = length(matrix_indices)

    # Collect vector of indices
    all_ind_collected = arr_to_vec(matrix_indices)

    # Sort indices vector according to inpu array
    # TODO Cant this be done with sortperm? in_matrix > UpperTriangular |> sortperm
    index_sorting = sort_indices_by_values(in_matrix, all_ind_collected)

    ordering_number = ordering_start
    for k = 1:total_elements
        # global ordering_number
        next_sorted_pos = index_sorting[k]
        mat_ind = matrix_indices[next_sorted_pos]

        if assign_same_values && k != 1
            prev_sorted_pos = index_sorting[k-1]
            prev_mat_ind = matrix_indices[prev_sorted_pos]

            cond1 = in_matrix[prev_mat_ind] == in_matrix[mat_ind]
            cond2 = small_dist_grouping
            cond3 = abs(in_matrix[prev_mat_ind] - in_matrix[mat_ind]) < min_dist

            if cond1 || (cond2 && cond3)
                ordering_number -= 1
            end
        end
        set_values!(ord_mat, mat_ind, ordering_number; do_symmetry = do_symmetry)
        ordering_number += 1
        # else
        #     set_values!(ord_mat, mat_ind, ordering_number; do_symmetry=do_symmetry)
        #     ordering_number+=1
        # end
    end

    return ord_mat
end

# TODO this one has to be specified for 3 dim matrix
function get_ordered_matrix(input_array::Array{Any,3}; do_slices = true, dims = 0)
    arr_size = size(input_array)
    out_arr = zeros(Int, arr_size)

    if do_slices
        # param check
        if dims > length(arr_size)
            throw(DomainError("Given dimension is greater than total size of array."))
        elseif dims > 0
            throw(DomainError("Given dimension must be positive value."))
        elseif dims <= 3
            throw(DomainError("Given dimension must be lower than 3."))
        end

        for dim = 1:arr_size[dims]
            if dims == 1
                out_arr[dim, :, :] = get_ordered_matrix(input_array[dim, :, :])
            elseif dims == 2
                out_arr[:, dim, :] = get_ordered_matrix(input_array[:, dim, :])
            elseif dims == 3
                out_arr[:, :, dim] = get_ordered_matrix(input_array[:, :, dim])
            end
        end
    else
        out_arr = get_ordered_matrix(input_array)
    end
end

function get_ordered_matrix(input_array::Array)
    out_array = copy(input_array)
    arr_size = size(input_array)
    total_elements = length(input_array)

    # Collect vector of indices
    all_ind_collected = collect(reshape(generate_indices(arr_size), (length(input_array))))

    # Sort indices vector according to inpu array
    index_sorting = sort!(
        [1:total_elements;],
        by = i -> (input_array[all_ind_collected][i], all_ind_collected[i]),
    )

    for k = 1:total_elements
        target = index_sorting[k]
        out_array[target] = k
    end

    return out_array
end

# Care must be taken so that values from 'input_matrix' are within distance
# groups, otherwise error is thrown.
"""
    function group_distances(input_matrix::Array, total_dist_groups::Int)

Takes a matrix and rearranges values into 'total_dist_groups' number of groups.
Every group is assigned with number value from range '<0,1>'.
"""
function group_distances(input_matrix::Array, total_dist_groups::Int)

    normed_matrix = normalize_to_01(input_matrix)
    target_matrix = copy(normed_matrix)

    h, w = size(input_matrix)

    if h * w < total_dist_groups
        throw(DomainError("Total number of groups exceed total number of entries in input matrix"))
    end

    total_borders = total_dist_groups + 1

    range_val = collect(range(0, 1, length = total_borders))

    for k = 2:total_borders
        indices = findall(x -> x >= range_val[k-1] && x <= range_val[k], normed_matrix)
        target_matrix[indices] .= range_val[k]
    end
    unique(target_matrix)

    # Sets last range to values smaller than unity, just in case this might cause trobules
    # normed_matrix[normed_matrix .> range_val[end-1]] .= 0.99
    return target_matrix
end

"""
    generate_indices(matrix_size::Tuple; symmetry_order::Bool=false, include_diagonal::Bool=true)

Return all the possible indices of the matrix of size 'matrix_size'.
'matrix_size' may be a tuple or a series of integer arguments corresponding to
the lengths in each dimension.

If 'symetry_order' is set to'true', then only indices of values below diagonal
are returned.
"""
function generate_indices(matrix_size::Tuple;
                            symmetry_order::Bool = false,
                            include_diagonal::Bool = true)

    # Get all cartesian indices from input matrix
    matrix_indices = CartesianIndices(matrix_size)
    # Filter out indices below diagonal
    if symmetry_order
        matrix_indices = findall(x -> x[1] <= x[2], matrix_indices)
    else
        matrix_indices = findall(x -> true, matrix_indices)
    end

    if !include_diagonal
        filter!(x -> x[1] != x[2], matrix_indices)
    end

    return matrix_indices
end

"""
    generate_indices(matrix_size::Int; symmetry_order::Bool=false, include_diagonal::Bool=true)

Generate indices for a matrix of given dimensions. 'generate_indices' is a
series of integer arguments corresponding to the lengths in each dimension.
"""
function generate_indices(matrix_size::Int;
                            symmetry_order::Bool = false,
                            include_diagonal::Bool = true)
    return generate_indices(
        (matrix_size, matrix_size);
        symmetry_order = symmetry_order,
        include_diagonal = include_diagonal,
    )
end

"""
    arr_to_vec(some_array::Array)

Takes an array and reshapes it into a vector.
"""
function arr_to_vec(some_array::Array)
    return collect(reshape(some_array, length(some_array)))
end

function cartesianInd_to_vec(some_array::CartesianIndices)
    return collect(reshape(some_array, length(some_array)))
end

"""
    sort_indices_by_values(values_matrix::T, index_vector) where {T<:VecOrMat}

Sorts the 'index_vector' according to corresponding values in the 'values_matrix'
and returns a Vector of intigers which is an list of ordering of
'sorted index_vector'.
"""
function sort_indices_by_values(values_matrix::T, index_vector) where {T<:VecOrMat}
    if !isa(index_vector, Vector)
        throw(TypeError(
            :sort_indices_by_values,
            "\'index_vector\' must be a vector, otherwise an ordering list can no be created!",
            Vector,
            typeof(index_vector),
        ))
    end
    total_elements = length(index_vector)
    return sort!(
        [1:total_elements;],
        by = i -> (values_matrix[index_vector][i], index_vector[i]),
    )
end

"""
    set_values!(input_matrix::Matrix, position::CartesianIndex, target_value::Number; do_symmetry=false)

Assigns 'target_value' to indices at 'input_matrix[position[1], position[2]]'.
If 'do_symmetry' is set to 'true', then the 'target_value' is also assigned at
position 'input_matrix[position[2], position[1]]'.
"""
function set_values!(input_matrix::Matrix,
                        position::CartesianIndex,
                        target_value::Number;
                        do_symmetry::Bool = false)
    input_matrix[position[1], position[2]] = target_value
    if do_symmetry
        input_matrix[position[2], position[1]] = target_value
    end
    return input_matrix
end

# matrix ordering
# =====

function get_high_dim_ordered_matrix(input_matrix)
    matrix_size = size(input_matrix)
    ordered_matrix_3D = zeros(Int, matrix_size)

    for slice = 1:matrix_size[1]
        ordered_matrix_3D[slice, :, :] = get_ordered_matrix(input_matrix[slice, :, :])
    end
    return ordered_matrix_3D
end

"""
    reduce_arrs_to_min_len(arrs)

Takes vector of vectors of different length and returns array of arrays which
are of the same length. Length in the output is the shortest vector length from
the input- values above this size are discarded.
"""
function reduce_arrs_to_min_len(arrs::Array)
    @debug "Argument specific"
    new_arr = copy(arrs)

    simulation = size(new_arr, 1)
    min_size = Inf
    for m = 1:simulation
        @debug "Simulation number" m
        current_size = size(new_arr[m], 1)
        @debug "Current size: " current_size
        if convert(Float64, current_size) < min_size
            min_size = current_size
            @debug "min size changed to: " min_size
        end
    end
    # min_size = Int.(min_size)
    @debug "Concatenating"
    for m = 1:simulation
        new_arr[m] = new_arr[m][1:min_size, :]
    end
    min_size = Inf
    return new_arr
end

"""
    increase_arrs_to_max_len(arrs)

Takes vector of vectors of different length and returns array of arrays which
are of the same length. Length in the output is the longest vector length from
the input- values above this size are discarded.
"""
function increase_arrs_to_max_len(arrs)
    new_arr = copy(arrs)

    simulation = size(new_arr, 1)
    max_size = 0
    for m = 1:simulation
        @debug "Simulation number" m
        current_size = size(new_arr[m], 1)
        @debug "Current size: " current_size
        if convert(Float64, current_size) > max_size
            max_size = current_size
            @debug "min size changed to: " max_size
        end
    end
    # max_size = Int.(max_size)
    @debug "Concatenating"
    for m = 1:simulation
        correct_len_arr = zeros(Int, max_size, 3)
        correct_len_arr[1:size(arrs[m], 1), :] = new_arr[m][:, :]
        new_arr[m] = correct_len_arr
    end
    # min_size = Inf
    return new_arr
end


text:
const ci_cache = GPUCompiler.CodeCache()

struct oneAPICompilerParams <: AbstractCompilerParams end

oneAPICompilerJob = CompilerJob{SPIRVCompilerTarget,oneAPICompilerParams}

GPUCompiler.runtime_module(::oneAPICompilerJob) = oneAPI

# TODO: eagerly lower these using the translator API
GPUCompiler.isintrinsic(job::oneAPICompilerJob, fn::String) =
    invoke(GPUCompiler.isintrinsic,
           Tuple{CompilerJob{SPIRVCompilerTarget}, typeof(fn)},
           job, fn) ||
    in(fn, opencl_builtins)

function GPUCompiler.finish_module!(job::oneAPICompilerJob, mod::LLVM.Module)
    invoke(GPUCompiler.finish_module!,
           Tuple{CompilerJob{SPIRVCompilerTarget}, typeof(mod)},
           job, mod)
    ctx = LLVM.context(mod)

    # OpenCL 2.0
    push!(metadata(mod)["opencl.ocl.version"],
          MDNode([ConstantInt(Int32(2); ctx),
                  ConstantInt(Int32(0); ctx)]; ctx))

    # SPIR-V 1.5
    push!(metadata(mod)["opencl.spirv.version"],
          MDNode([ConstantInt(Int32(1); ctx),
                  ConstantInt(Int32(5); ctx)]; ctx))
end

GPUCompiler.ci_cache(::oneAPICompilerJob) = ci_cache

GPUCompiler.method_table(::oneAPICompilerJob) = method_table


text:
using Peacock
using PyPlot
#using Peacock.FDFD

function change_cell(kappa, Nkx, Nky, d1, d2, bands)
#kappa = 0im
mu1 = [14 kappa 0;-kappa 14 0;0 0 1]
# Permittivity
function epf(x,y)
    # equation of a circle with radius 0.2a
    s = 1/2
    r = 0.1/sqrt(2)
    if x^2+y^2 <= r^2 || (x+s)^2+(y+s)^2 < r^2 || (x+s)^2+(y-s)^2 < r^2 || (x-s)^2+(y+s)^2 < r^2 || (x-s)^2+(y-s)^2 < r^2
        # dielectric inside the circle
        return 15
    else
        # air outside the circle
        return 1
    end
end
function muf(x,y)
    # equation of a circle with radius 0.2a
    s = 1/2
    r = 0.1/sqrt(2)
    if x^2+y^2 <= r^2 || (x+s)^2+(y+s)^2 < r^2 || (x+s)^2+(y-s)^2 < r^2 || (x-s)^2+(y+s)^2 < r^2 || (x-s)^2+(y-s)^2 < r^2
        # dielectric inside the circle
        return mu1
    else
        # air outside the circle
        return 1
    end
end

a1 = [1, 0]  # first lattice vector
a2 = [0, 1]  # second lattice vector
#a2 = 0.5*[sqrt(3), 1]

#d1 = 0.005/4  # resolution along first lattice vector
#d2 = 0.005/4  # resolution along second lattice vector

G = BrillouinZoneCoordinate(  0,   0, "Γ")
X = BrillouinZoneCoordinate(1/2,   0, "X")
M = BrillouinZoneCoordinate(1/2, 1/2, "M")
ks = [G,X,M,G]

geometry_FDFD = Geometry(epf, muf, a1, a2, d1, d2, TM)
#solver_FDFD = Peacock.FDFD.Solver(geometry_FDFD, [2*d1, 2*d2])
#figure(size(3,4)); plot_band_diagram(solver_FDFD, ks, TM, bands=1:6, dk=0.4, frequency_scale=1/2pi)

Chern_number(epf, muf, a1, a2, d1, d2, TM, Nkx=Nkx, Nky=Nky, bands=bands)
end
#=
geometry_PWE = Geometry(epf, muf, a1, a2, d1, d2)
solver_PWE = Solver(geometry_PWE, 7)
figure(); plot_band_diagram(solver_PWE, ks, TM, bands=1:4, dk=0.4, frequency_scale=1/2pi)
=#


text:
using NativeSVG

dr =
    Drawing(width = "5cm", height = "4cm", viewBox = "0 0 500 400") do
        title() do
            str("Example cubic01- cubic Bézier commands in path data")
        end
        desc() do
            str("""Picture showing a simple example of path data
                   using both a "C" and an "S" command,
                   along with annotations showing the control points
                   and end points""")
        end
        style(type = "text/css") do
            cdata(""".Border { fill:none; stroke:blue; stroke-width:1 }
                     .Connect { fill:none; stroke:#888888; stroke-width:2 }
                     .SamplePath { fill:none; stroke:red; stroke-width:5 }
                     .EndPoint { fill:none; stroke:#888888; stroke-width:2 }
                     .CtlPoint { fill:#888888; stroke:none }
                     .AutoCtlPoint { fill:none; stroke:blue; stroke-width:4 }
                     .Label { font-size:22; font-family:Verdana }""")
        end
        rect(class = "Border", x = "1", y = "1", width = "498", height = "398")
        polyline(class = "Connect", points = "100,200 100,100")
        polyline(class = "Connect", points = "250,100 250,200")
        polyline(class = "Connect", points = "250,200 250,300")
        polyline(class = "Connect", points = "400,300 400,200")
        path(class = "SamplePath", d = "M100,200 C100,100 250,100 250,200
                                        S400,300 400,200")
        circle(class = "EndPoint", cx = "100", cy = "200", r = "10")
        circle(class = "EndPoint", cx = "250", cy = "200", r = "10")
        circle(class = "EndPoint", cx = "400", cy = "200", r = "10")
        circle(class = "CtlPoint", cx = "100", cy = "100", r = "10")
        circle(class = "CtlPoint", cx = "250", cy = "100", r = "10")
        circle(class = "CtlPoint", cx = "400", cy = "300", r = "10")
        circle(class = "AutoCtlPoint", cx = "250", cy = "300", r = "9")
        text(class = "Label", x = "25", y = "70") do
            str("M100,200 C100,100 250,100 250,200")
        end
        text(
            class = "Label",
            x = "325",
            y = "350",
            style = "text-anchor:middle"
        ) do
            str("S400,300 400,200")
        end
    end
display(dr)

dr =
    Drawing(width = "12cm", height = "6cm", viewBox = "0 0 1200 600") do
        title() do
            str("Example quad01 - quadratic Bézier commands in path data")
        end
        desc() do
            str("""Picture showing a "Q" a "T" command,
                   along with annotations showing the control points
                   and end points""")
        end
        rect(
            x = "1",
            y = "1",
            width = "1198",
            height = "598",
            fill = "none",
            stroke = "blue",
            stroke_width = "1"
        )
        path(
            d = "M200,300 Q400,50 600,300 T1000,300",
            fill = "none",
            stroke = "red",
            stroke_width = "5"
        )
    # End points
        g(fill = "black") do
            circle(cx = "200", cy = "300", r = "10")
            circle(cx = "600", cy = "300", r = "10")
            circle(cx = "1000", cy = "300", r = "10")
        end
    # Control points and lines from end points to control points
        g(fill = "#888888") do
            circle(cx = "400", cy = "50", r = "10")
            circle(cx = "800", cy = "550", r = "10")
        end
        path(
            d = "M200,300 L400,50 L600,300 L800,550 L1000,300",
            fill = "none",
            stroke = "#888888",
            stroke_width = "2"
        )
    end
display(dr)

dr =
    Drawing(width = "12cm", height = "5.25cm", viewBox = "0 0 1200 400") do
        title() do
            str("Example arcs01 - arc commands in path data")
        end
        desc() do
            str("""Picture of a pie chart with two pie wedges and
                   a picture of a line with arc blips""")
        end
        rect(
            x = "1",
            y = "1",
            width = "1198",
            height = "398",
            fill = "none",
            stroke = "blue",
            stroke_width = "1"
        )

        path(
            d = "M300,200 h-150 a150,150 0 1,0 150,-150 z",
            fill = "red",
            stroke = "blue",
            stroke_width = "5"
        )
        path(
            d = "M275,175 v-150 a150,150 0 0,0 -150,150 z",
            fill = "yellow",
            stroke = "blue",
            stroke_width = "5"
        )

        path(
            d = """M600,350 l 50,-25
                   a25,25 -30 0,1 50,-25 l 50,-25
                   a25,50 -30 0,1 50,-25 l 50,-25
                   a25,75 -30 0,1 50,-25 l 50,-25
                   a25,100 -30 0,1 50,-25 l 50,-25""",
            fill = "none",
            stroke = "red",
            stroke_width = "5"
        )
    end

max_stars_repo_name: BenLauwens/NativeSVG

max_issues_repo_name: BenLauwens/NativeSVG

max_forks_repo_name: BenLauwens/NativeSVG


text:
function build_blockchain_data_url(fname::String; kwargs...)
    if fname == "available_coin_list"
        url = "https://min-api.cryptocompare.com/data/blockchain/list?"
    elseif fname == "latest"
        url = "https://min-api.cryptocompare.com/data/blockchain/latest?"
    elseif fname == "historical_daily"
        url = "https://min-api.cryptocompare.com/data/blockchain/histo/day?"
    elseif fname == "mining_calculator_data"
        url = "https://min-api.cryptocompare.com/data/blockchain/mining/calculator?"
    end

    # collect url parts in list
    url_parts = []

    for (key, value) in kwargs
        # conversion
        if key == :fsym
            push!(url_parts,"fsym=$(value)")
        elseif key == :fsyms
            push!(url_parts,"fsyms=$(join(value,","))")
        elseif key == :tsyms
            push!(url_parts,"tsyms=$(join(value,","))")

        # limit
        elseif key==:limit
            push!(url_parts, "limit=$(value)")
        # toTs
        elseif key==:toTs && value ≠ ""
            push!(url_parts, "toTs=$(value)")

        # extraParams
        elseif key == :extraParams && value ≠ ""
            push!(url_parts,"extraParams=$(value)")
        # sign
        elseif key == :sign && value==true
            push!(url_parts,"sign=$(value)")

        # api_key
        elseif key == :api_key && value!=""
            push!(url_parts, "api_key=$(value)")
        end

    end

    # put url together
    url = url * join(url_parts,"&")

    return url

end

function build_exchange_vol_url(fname::String; kwargs...)
    if fname in ["day", "hour", "minute"]
        url = "https://min-api.cryptocompare.com/data/exchange/histo$(fname)?"
    end

    # collect url parts in list
    url_parts = []

    for (key, value) in kwargs
        # conversion

        if key == :tsym
            push!(url_parts,"tsym=$(value)")
        # exchange
        elseif key == :e && value!= "all"
            push!(url_parts, "e=$(value)")
        elseif key == :e && value == "all"
            push!(url_parts, "e=CCCAGG")

        # aggregate
        elseif key==:aggregate && value>1
            push!(url_parts, "aggregate=$(value)")
        # aggregate Predictable Time Periods
        elseif key==:aggregatePredictableTimePeriods && value==false
            push!(url_parts, "aggregatePredictableTimePeriods=$(value)")
        # limit
        elseif key==:limit
            push!(url_parts, "limit=$(value)")
        # toTs
        elseif key==:toTs && value!=""
            push!(url_parts, "toTs=$(value)")
        # extraParams
        elseif key == :extraParams && value ≠ ""
            push!(url_parts,"extraParams=$(value)")
        # sign
        elseif key == :sign && value==true
            push!(url_parts,"sign=$(value)")
        end

    end

    # put url together
    url = url * join(url_parts,"&")

    return url
end

function build_symbol_vol_url(fname::String; kwargs...)
    if fname in ["day", "hour"]
        url = "https://min-api.cryptocompare.com/data/symbol/histo$(fname)?"
    end

    # collect url parts in list
    url_parts = []

    for (key, value) in kwargs
        # conversion
        if key == :tsym
            push!(url_parts,"tsym=$(value)")
        elseif key == :fsym
            push!(url_parts,"fsym=$(value)")

        # aggregate
        elseif key==:aggregate && value>1
            push!(url_parts, "aggregate=$(value)")
        # limit
        elseif key==:limit
            push!(url_parts, "limit=$(value)")
        # toTs
        elseif key==:toTs && value!=""
            push!(url_parts, "toTs=$(value)")
        # extraParams
        elseif key == :extraParams && value ≠ ""
            push!(url_parts,"extraParams=$(value)")
        # sign
        elseif key == :sign && value==true
            push!(url_parts,"sign=$(value)")
        end

    end

    # put url together
    url = url * join(url_parts,"&")

    return url
end

function build_symbol_vol_single_exchange_url(fname::String; kwargs...)
    if fname in ["day", "hour"]
        url = "https://min-api.cryptocompare.com/data/exchange/symbol/histo$(fname)?"
    end

    # collect url parts in list
    url_parts = []

    for (key, value) in kwargs
        # conversion
        if key == :tsym
            push!(url_parts,"tsym=$(value)")
        elseif key == :fsym
            push!(url_parts,"fsym=$(value)")
        # exchange
        elseif key == :e && value!= "all"
            push!(url_parts, "e=$(value)")
        elseif key == :e && value == "all"
            push!(url_parts, "e=CCCAGG")
        # aggregate
        elseif key==:aggregate && value>1
            push!(url_parts, "aggregate=$(value)")
        # limit
        elseif key==:limit
            push!(url_parts, "limit=$(value)")
        # toTs
        elseif key==:toTs && value!=""
            push!(url_parts, "toTs=$(value)")
        # extraParams
        elseif key == :extraParams && value ≠ ""
            push!(url_parts,"extraParams=$(value)")
        # sign
        elseif key == :sign && value==true
            push!(url_parts,"sign=$(value)")
        end

    end

    # put url together
    url = url * join(url_parts,"&")

    return url
end

function build_historical_data_url(fname::String; kwargs...)
    if fname in ["day", "hour", "minute"]
        url = "https://min-api.cryptocompare.com/data/v2/histo$(fname)?"
    end

    # collect url parts in list
    url_parts = []

    for (key, value) in kwargs
        # conversion
        if key == :try_conversion && value==false
            push!(url_parts, "tryConversion=false")
        elseif key == :fsym
            push!(url_parts,"fsym=$(value)")
        elseif key == :fsyms
            push!(url_parts,"fsyms=$(join(value,","))")
        elseif key == :tsym
            push!(url_parts,"tsym=$(value)")
        elseif key == :tsyms
            push!(url_parts,"tsyms=$(join(value,","))")

        # exchange
        elseif key == :e && value!= "all"
            push!(url_parts, "e=$(value)")
        elseif key == :e && value == "all"
            push!(url_parts, "e=CCCAGG")

        # aggregate
        elseif key==:aggregate && value>1
            push!(url_parts, "aggregate=$(value)")
        # aggregate Predictable Time Periods
        elseif key==:aggregatePredictableTimePeriods && value==false
            push!(url_parts, "aggregatePredictableTimePeriods=$(value)")
        # limit
        elseif key==:limit
            push!(url_parts, "limit=$(value)")
        # allData
        elseif key==:allData && value==true
            push!(url_parts, "allData=$(value)")
        # toTs
        elseif key==:toTs && value ≠ ""
            push!(url_parts, "toTs=$(value)")
        # explainPath
        elseif key==:explainPath && value==true
            push!(url_parts, "explainPath=$(value)")
        # extraParams
        elseif key == :extraParams && value ≠ ""
            push!(url_parts,"extraParams=$(value)")
        # sign
        elseif key == :sign && value==true
            push!(url_parts,"sign=$(value)")
        end

    end

    # put url together
    url = url * join(url_parts,"&")

    return url

end

function build_price_url(fname::String; kwargs...)
    if fname=="single_symbol_price"
        url = "https://min-api.cryptocompare.com/data/price?"
    elseif fname == "multiple_symbol_price"
        url = "https://min-api.cryptocompare.com/data/pricemulti?"
    elseif fname == "multiple_symbol_full"
        url = "https://min-api.cryptocompare.com/data/pricemultifull?"
    elseif fname == "generateAvg"
        url = "https://min-api.cryptocompare.com/data/generateAvg?"
    end

    # collect url parts in list
    url_parts = []

    for (key, value) in kwargs
        if key == :fsym
            push!(url_parts,"fsym=$(value)")
        elseif key == :fsyms
            push!(url_parts,"fsyms=$(join(value,","))")
        elseif key == :tsym
            push!(url_parts,"tsym=$(value)")
        elseif key == :tsyms
            push!(url_parts,"tsyms=$(join(value,","))")

        # exchange
        elseif key == :e && value!= "all"
            push!(url_parts, "e=$(value)")
        elseif key == :e && value == "all"
            push!(url_parts, "e=cccagg_or_exchange")

        # conversion
        elseif key == :try_conversion && value==false
            push!(url_parts, "tryConversion=false")

        elseif key == :relaxedValidation && value==false
            push!(url_parts, "relaxedValidation=false")

        elseif key == :extraParams && value ≠ ""
            push!(url_parts,"extraParams=$(value)")
        elseif key == :sign && value==true
            push!(url_parts,"sign=$(value)")
        end

    end

    # put url together
    url = url * join(url_parts,"&")

    return url

end

function build_news_url(fname::String; kwargs...)
    # latest_news_articles_url = "https://min-api.cryptocompare.com/data/v2/news/?lang=EN"
    # latest_news_feeds_url = "https://min-api.cryptocompare.com/data/news/feeds"
    # news_articles_categories_url = "https://min-api.cryptocompare.com/data/news/categories"
    # list_news_feeds_and_categories_ulr = "https://min-api.cryptocompare.com/data/news/feedsandcategories"

    # keys2 = [feeds::String, categories::String, excludeCategories::String, ITs::Int64, lang::String, sortOrder::String]
    # keys = [extraParams::String, sign::Bool]
    if fname == "latest_news_articles"
        url = "https://min-api.cryptocompare.com/data/v2/news/?"
    elseif fname in ["feeds", "categories", "feedsandcategories"]
        url = "https://min-api.cryptocompare.com/data/news/$(fname)?"
    end

    url_parts = []

    for (key,value) in kwargs
        if key == :feeds
            if size(value)[1]≥1
                push!(url_parts,"feeds=$(join(value,","))")
            end
        elseif key == :categories
            if size(value)[1]≥1
                push!(url_parts,"categories=$(join(value,","))")
            end
        elseif key == :excludeCategories
            if size(value)[1]≥1
                push!(url_parts,"excludeCategories=$(join(value,","))")
            end
        elseif key == :lTs
            push!(url_parts,"lTs=$(value)")
        elseif key == :lang
            push!(url_parts,"lang=$(value)")
        elseif key == :sortOrder
            push!(url_parts,"sortOrder=$(value)")
        elseif key == :extraParams
            if value ≠ ""
                push!(url_parts,"extraParams=$(value)")
            end
        elseif key == :sign
            push!(url_parts,"sign=$(value)")
        end
    end

    # put url together
    url = url * join(url_parts,"&")
    #return url_parts
    return url
end

function build_url2(fname::String; kwargs...)

    if fname in ["price", "coinlist", "coinsnapshot", "miningcontracts","miningequipment"]
        url = "https://www.cryptocompare.com/api/data/$(fname)?"
    elseif fname in ["exchanges", "volumes", "pairs"]
        url = "https://min-api.cryptocompare.com/data/top/$(fname)?"
    elseif fname in ["minute", "hour", "day"]
        url = "https://min-api.cryptocompare.com/data/histo$(fname)?"
    else
        url = "https://min-api.cryptocompare.com/data/$(fname)?"
    end

    # collect url parts in list
    url_parts = []

    for (key, value) in kwargs
        if key == :fsym
            append!(url_parts,"fsym=$(value)")
        elseif key == :fsyms
            append!(url_parts,"fsyms=$(join(value,","))")
        elseif key == :tsym
            append!(url_parts,"tsym=$(value)")
        elseif key == :tsyms
            append!(url_parts,"tsyms=$(join(value,","))")

        # exchange
        elseif key == :e && value!= "all"
            append!(url_parts, "e=$(value)")
        elseif key == :e && value == "all"
            append!(url_parts, "e=CCCAGG")

        # conversion
        elseif key == :try_conversion && value==false
            append!(url_parts, "tryConversion=false")

        # markets
        elseif key == :markets && value != "all"
            append!(url_parts, "e=$(join(value, ","))")
        elseif key == :markets && value == "all"
            append!(url_parts, "e=CCCAGG")

        elseif key == :avg_type && value!= "HourVWAP"
            append!(url_parts, "avgType=$(value)")
        elseif key == :utc_hour_diff && value!=0
            append!(url_parts, "UTCHourDiff=$(value)")
        # elseif key == :to_ts && value==true
        #     timestamp = to_ts_args_to_timestamp(value)
        #     if timestamp:
        #         url_parts.append("toTs={}".format(timestamp))
        elseif key==:aggregate && value!= 1
            append!(url_parts, "aggregate=$(value)")

        elseif key == :limit
            append!(url_parts, "limit=$(value)")
        end

    end

    # put url together
    url = url * join(url_parts,"&")

    return url
end


text:
@traitfn function mean_f(
    model::TGP, kernel_matrices
) where {TGP <: AbstractGPModel; !IsMultiOutput{TGP}}
    return mean_f.(model.f, kernel_matrices)
end

@traitfn function var_f(
    model::TGP, kernel_matrices
) where {TGP <: AbstractGPModel; !IsMultiOutput{TGP}}
    return var_f.(model.f, kernel_matrices)
end

## return the expectation gradient given μ ##
@traitfn function ∇E_μ(m::TGP, y, state) where {TGP <: AbstractGPModel; !IsMultiOutput{TGP}}
    return ∇E_μ(likelihood(m), opt(inference(m)), y, state)
end

## return the expectation gradient given diag(Σ) ##
@traitfn function ∇E_Σ(m::TGP, y, state) where {TGP <: AbstractGPModel; !IsMultiOutput{TGP}}
    return ∇E_Σ(likelihood(m), opt(inference(m)), y, state)
end

##
@traitfn function mean_f(
    model::TGP, kernel_matrices
) where {T,TGP<:AbstractGPModel{T};IsMultiOutput{TGP}}
    μ_q = mean_f.(model.f, kernel_matrices)
    return ntuple(n_output(model)) do i
        ntuple(model.nf_per_task[i]) do j
            sum(model.A[i][j] .* μ_q)
        end
    end
end

##
@traitfn function var_f(
    model::TGP, kernel_matrices
) where {T,TGP<:AbstractGPModel{T};IsMultiOutput{TGP}}
    Σ_q = var_f.(model.f, kernel_matrices)
    return ntuple(n_output(model)) do i
        ntuple(model.nf_per_task[i]) do j
            sum(model.A[i][j] .^ 2 .* Σ_q)
        end
    end
end

## return the linear sum of the expectation gradient given μ ##
@traitfn function ∇E_μ(
    m::TGP, y, state
) where {T,TGP<:AbstractGPModel{T};IsMultiOutput{TGP}}
    ∇ = [zeros(T, batchsize(inference(m))) for i in 1:n_latent(m)]
    ∇Eμs = ∇E_μ.(likelihood(m), Ref(opt(inference(m))), y, state.local_vars)
    ∇EΣs = ∇E_Σ.(likelihood(m), Ref(opt(inference(m))), y, state.local_vars)
    μ_f = mean_f.(m.f, state.kernel_matrices)
    for t in 1:n_output(m)
        for j in 1:m.nf_per_task[t]
            for q in 1:n_latent(m)
                ∇[q] .+=
                    m.A[t][j][q] * (
                        ∇Eμs[t][j] -
                        2 * ∇EΣs[t][j] .*
                        sum(m.A[t][j][qq] * μ_f[qq] for qq in 1:n_latent(m) if qq != q)
                    )
            end
        end
    end
    return ∇
end

## return the linear sum of the expectation gradient given diag(Σ) ##
@traitfn function ∇E_Σ(
    m::TGP, y, state
) where {T,TGP<:AbstractGPModel{T};IsMultiOutput{TGP}}
    ∇ = [zeros(T, batchsize(inference(m))) for _ in 1:n_latent(m)]
    ∇Es = ∇E_Σ.(likelihood(m), Ref(opt(inference(m))), y, state.local_vars)
    for t in 1:n_output(m)
        for j in 1:m.nf_per_task[t]
            for q in 1:n_latent(m)
                ∇[q] .+= m.A[t][j][q]^2 * ∇Es[t][j]
            end
        end
    end
    return ∇
end

##
function update_A!(m::TGP, state, ys) where {TGP<:AbstractGPModel}
    if !isnothing(m.A_opt)
        local_vars = state.local_vars
        μ_f = mean_f.(m.f, state.kernel_matrices) # κμ || μ
        Σ_f = var_f.(m.f, state.kernel_matrices) # Diag(K̃ + κΣκ) || Diag(Σ)
        ∇Eμ = ∇E_μ.(likelihood(m), Ref(opt(inference(m))), ys, local_vars)
        ∇EΣ = ∇E_Σ.(likelihood(m), Ref(opt(inference(m))), ys, local_vars)
        for t in 1:n_output(m)
            for j in 1:m.nf_per_task[t]
                ∇A = zero(m.A[t][j])
                for q in 1:n_latent(m)
                    x1 =
                        dot(∇Eμ[t][j], μ_f[q]) -
                        2 * dot(
                            ∇EΣ[t][j],
                            μ_f[q] .*
                            sum(m.A[t][j][qq] * μ_f[qq] for qq in 1:n_latent(m) if qq != q),
                        )
                    x2 = dot(∇EΣ[t][j], abs2.(μ_f[q]) + Σ_f[q])
                    # new_A[t,j,q] = x1/(2*x2)
                    ∇A[q] = x1 - 2 * m.A[t][j][q] * x2
                end
                state_A = state.A_state[t][j]
                state_A, ΔA = Optimisers.apply(m.A_opt, state_A, m.A[t][j], ∇A)
                m.A[t][j] .+= ΔA
                m.A[t][j] /= sqrt(sum(abs2, m.A[t][j])) # Projection on the unit circle
                state.A_state[t][j] = state_A
            end
        end
    end
    return state
end


text:
# calls to official packages ...
using DataFrames
using DifferentialEquations
using CSV
using Convex
using Optim
using GLPK
using Plots
using UUIDs
using Pluto
using PlutoUI


text:
import Base: +, -, ∘, real, imag, abs

# Set it to negative of itself
function (-)(p_in::ScalarData)
  p = deepcopy(p_in)
  @. p.data = -p.data
  return p
end

function (-)(p_in::VectorData)
  p = deepcopy(p_in)
  @. p.u = -p.u
  @. p.v = -p.v
  return p
end

# Add and subtract the same type
function (-)(p1::T,p2::T) where {T<:ScalarData}
  return T(p1.data .- p2.data)
end

function (+)(p1::T,p2::T) where {T<:ScalarData}
  return T(p1.data .+ p2.data)
end

function (-)(p1::T,p2::T) where {T<:VectorData}
  return T(p1.u - p2.u, p1.v - p2.v)
end

function (+)(p1::T,p2::T) where {T<:VectorData}
  return T(p1.u + p2.u, p1.v + p2.v)
end

# Multiply and divide by a constant
function (*)(p::T,c::Number) where {T<:ScalarData}
  return T(c*p.data)
end

function (/)(p::T,c::Number) where {T<:ScalarData}
  return T(p.data ./ c)
end

function (*)(p::T,c::Number) where {T<:VectorData}
  return T(c*p.u,c*p.v)
end

(*)(c::Number,p::T) where {T<:PointData} = *(p,c)

function (/)(p::T,c::Number) where {T<:VectorData}
  return T(p.u / c, p.v / c)
end

function (*)(p1::T,p2::T) where {T<:ScalarData}
  return T(p1.data .* p2.data)
end

function (*)(p1::T,p2::T) where {T<:VectorData}
  return T(p1.u * p2.u, p1.v * p2.v)
end

"""
    (+)(X::VectorData,a::Tuple{T,T}) where {T<:Number} -> VectorData
    (-)(X::VectorData,a::Tuple{T,T}) where {T<:Number} -> VectorData

Adds or subtracts the tuple `a` component by component to each element of `X`. All data in `a` are
converted to Float64. Can also switch the arguments.

# Example

```jldoctest
julia> f = VectorData(5);

julia> f + (2,3)
5 points of vector-valued Float64 data
5×2 Array{Float64,2}:
 2.0  3.0
 2.0  3.0
 2.0  3.0
 2.0  3.0
 2.0  3.0
```
"""
function (+)(A::VectorData,a::Tuple{T,T}) where {T <: Number}
   B = deepcopy(A)
    u, v = a
    B.u .+= u
    B.v .+= v
    return B
end

function (-)(A::VectorData,a::Tuple{T,T}) where {T <: Number}
   B = deepcopy(A)
    u, v = a
    B.u .-= u
    B.v .-= v
    return B
end

### Hadamard products

"""
    product!(out::PointData,p::PointData,q::PointData)

Compute the Hadamard (i.e. element by element) product of point data
data `p` and `q` and return the result in `out`. Note that `p` and `q`
can be of mixed type (scalar, vector, tensor), as long as one of them is
a scalar. Also, `out` must have element type that is consistent with the promoted
type of `p` and `q`.

# Example
```jldoctest
julia> fcs = ScalarData(5,dtype=ComplexF64);

julia> fill!(fcs,2im)
5 points of scalar-valued Complex{Float64} data
5-element Array{Complex{Float64},1}:
 0.0 + 2.0im
 0.0 + 2.0im
 0.0 + 2.0im
 0.0 + 2.0im
 0.0 + 2.0im

julia> frt = TensorData(fcs,dtype=Float64);

julia> fill!(frt,1.0)
5 points of tensor-valued Float64 data dudx, dudy, dvdx, dvdy
5×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> out = similar(frt,element_type=ComplexF64);

julia> product!(out,frt,fcs)
5 points of tensor-valued Complex{Float64} data dudx, dudy, dvdx, dvdy
5×4 Array{Complex{Float64},2}:
 0.0+2.0im  0.0+2.0im  0.0+2.0im  0.0+2.0im
 0.0+2.0im  0.0+2.0im  0.0+2.0im  0.0+2.0im
 0.0+2.0im  0.0+2.0im  0.0+2.0im  0.0+2.0im
 0.0+2.0im  0.0+2.0im  0.0+2.0im  0.0+2.0im
 0.0+2.0im  0.0+2.0im  0.0+2.0im  0.0+2.0im
```
"""
function product!(out::ScalarData{N},
                  p::ScalarData{N},
                  q::ScalarData{N}) where {N}

    @inbounds for x in 1:N
        out.data[x] = p.data[x] * q.data[x]
    end
    out
end

function product!(out::VectorData{N},
                  p::VectorData{N},
                  q::VectorData{N}) where {N}

    @inbounds for x in 1:N
        out.u[x] = p.u[x] * q.u[x]
        out.v[x] = p.v[x] * q.v[x]
    end
    out
end

function product!(out::TensorData{N},
                  p::TensorData{N},
                  q::TensorData{N}) where {N}

    @inbounds for x in 1:N
        out.dudx[x] = p.dudx[x] * q.dudx[x]
        out.dudy[x] = p.dudy[x] * q.dudy[x]
        out.dvdx[x] = p.dvdx[x] * q.dvdx[x]
        out.dvdy[x] = p.dvdy[x] * q.dvdy[x]
    end
    out
end

function product!(out::VectorData{N},
                  p::ScalarData{N},
                  q::VectorData{N}) where {N}

    @inbounds for x in 1:N
        out.u[x] = p[x] * q.u[x]
        out.v[x] = p[x] * q.v[x]
    end
    out
end

function product!(out::TensorData{N},
                  p::ScalarData{N},
                  q::TensorData{N}) where {N}

    @inbounds for x in 1:N
        out.dudx[x] = p[x] * q.dudx[x]
        out.dudy[x] = p[x] * q.dudy[x]
        out.dvdx[x] = p[x] * q.dvdx[x]
        out.dvdy[x] = p[x] * q.dvdy[x]
    end
    out
end

"""
    product(p::PointData,q::PointData) -> PointData
    (∘)(p::PointData,q::PointData) -> PointData

Compute the Hadamard (i.e. element by element) product of point data
data `p` and `q`. Works similarly to `product!`.
"""
function product end

for f in (:ScalarData, :VectorData, :TensorData)
    @eval function product(p::$f{N}, q::$f{N}) where {N}
        product!($f(N,dtype=promote_type(eltype(p),eltype(q))), p, q)
    end

    @eval (∘)(p::$f{N}, q::$f{N}) where {N} = product(p,q)
end

for f in (:VectorData, :TensorData)
    @eval product!(out::$f{N},p::$f{N},q::ScalarData{N}) where {N} = product!(out,q,p)

    @eval function product(p::$f{N},q::ScalarData{N}) where {N}
        product!($f(N,dtype=promote_type(eltype(p),eltype(q))), p, q)
    end

    @eval product(q::ScalarData{N},p::$f{N}) where {N} = product(p,q)

    @eval (∘)(p::$f{N}, q::ScalarData{N}) where {N} = product(p,q)

    @eval (∘)(q::ScalarData{N}, p::$f{N}) where {N} = product(p,q)

end

### Operations between tuples and vectors

a::(Tuple{T,T} where {T}) + A::VectorData = A + a
a::(Tuple{T,T} where {T}) - A::VectorData = (B = A - a; @. B.u = -B.u; @. B.v = -B.v; return B)

"""
    cross(a::Number/ScalarData,A::VectorData) -> VectorData
    ×(a::Number/ScalarData,A::VectorData) -> VectorData

Compute the cross product between the scalar `a` (treated as an out-of-plane component of a vector)
and the planar vector data `A`.
"""
function cross(a::Union{Number,ScalarData},A::VectorData)
    B = deepcopy(A)
    @. B.u = -a*A.v
    @. B.v = a*A.u
    return B
end

### Vector operations

"""
    cross(A::VectorData,B::VectorData) -> ScalarData
    ×(A::VectorData,A::VectorData) -> ScalarData

Compute the cross product between the vector point data `A` and `B`
and return the result as scalar data (treated as an out-of-plane
component of a vector).
"""
function cross(A::VectorData{N},B::VectorData{N}) where N
    C = ScalarData(N,dtype=promote_type(eltype(A),eltype(B)))
    @. C = A.u*B.v - A.v*B.u
    return C
end

"""
    dot(A::Tuple{T,T},B::VectorData) where {T<:Number} -> ScalarData
    ⋅(A::Tuple{T,T},B::VectorData) where {T<:Number} -> ScalarData

Computes the dot product between the tuple `v` and the elements of a tensor `B` on
a set of points and returns scalar data on the same set of points.
"""
function dot(A::Tuple{T,T},B::VectorData) where {T<:Number}
    C = ScalarData(B)
    x, y = A
    @. C.data = x*B.u + y*B.v
    return C
end

"""
    dot(A::Tuple{T,T},B::TensorData) where {T<:Number} -> VectorData
    ⋅(A::Tuple{T,T},B::TensorData) where {T<:Number} -> VectorData

Computes the dot product between the tuple `A` and the elements of a tensor `B` on
a set of points and returns vector data on the same set of points.
"""
function dot(A::Tuple{T,T},B::TensorData) where {T<:Number}
    C = VectorData(B)
    x, y = A
    @. C.u = x*B.dudx + y*B.dudy
    @. C.v = x*B.dvdx + y*B.dvdy
    return C
end

### Operations on complex point data

for f in (:real, :imag, :abs)
  @eval function $f(A::PointData{N,T}) where {N,T <: ComplexF64}
      Acopy = similar(A,element_type=Float64)
      Acopy .= broadcast($f,A)
      return Acopy
  end
end


text:
using Test

@testset "creation" "creation of the FourVectors, Interfaces" begin
    include("creation.jl") end

@testset "transfrormations" "boost and rotations on example" begin
    include("transformations.jl") end


text:
using Distributions

# Uncertain parameters defined in "Projections and Uncertainties about Climate
# Change in an Era of Minimal Climate Policies" (Nordhaus, 2018) EViews file
# provided by the AEJ : https://www.aeaweb.org/articles?id=10.1257/pol.20170046

dice_default_mcs = @defsim begin
    climatedynamics.t2xco2 = LogNormal(1.106, 0.2646)
	totalfactorproductivity.ga0 = Normal(0.076, 0.056)
	emissions.gsigma1 = Normal(-0.0152, 0.0032)
	co2cycle.mueq = LogNormal(5.851, 0.2649)
	damages.a2 = Normal(0.00227, 0.001135)

    # save(damages.DAMAGES)
	# save(grosseconomy.YGROSS)
	# save(emissions.E)
	# save(co2cycle.MU)
	# save(damages.DAMFRAC)
end

function getmcs()
    return deepcopy(dice_default_mcs)
end


text:
mutable struct AutoAbstolAffect{T}
    current_max::T
end

function (p::AutoAbstolAffect)(integrator)
    p.current_max = max(p.current_max, maximum(integrator.u))
    integrator.opts.abstol = p.current_max * integrator.opts.reltol
    SciMLBase.u_modified!(integrator, false)
end

function AutoAbstol(save = true, init_current_max = 1e-9)
    affect! = AutoAbstolAffect(init_current_max)
    condtion = (u, t, integrator) -> true
    save_positions = (save, false)
    return SciMLBase.DiscreteCallback(condtion, affect!, save_positions = save_positions)
end


text:
# Section 6.1 : Truck without trailer of
# [RMT13] Rungger, Matthias and Mazo Jr, Manuel and Tabuada, Paulo
# Specification-guided controller synthesis for linear systems and safe linear-time temporal logic
# Proceedings of the 16th international conference on Hybrid systems: computation and control, 2013

N = 2

using FillArrays
using HybridSystems
A = LightAutomaton(N)
add_transition!(A, 1, 2, 1)
add_transition!(A, 2, 2, 1)
##add_transition!(A, 2, 3, 1)
##add_transition!(A, 3, 4, 1)
##add_transition!(A, 4, 4, 1)
#add_transition!(A, 1, 5, 1)
#add_transition!(A, 5, 6, 1)
#add_transition!(A, 6, 7, 1)

using Polyhedra
using CDDLib
const U = 4
const va = 15.6
const vb = 24.5
P0 = polyhedron(SimpleHRepresentation([-1. 0; 1 0; 0 -1; 0 1], [0., 35, U, U]), CDDLibrary())
Pa = polyhedron(SimpleHRepresentation([1. 0], [va]), CDDLibrary())
Pb = polyhedron(SimpleHRepresentation([1. 0], [vb]), CDDLibrary())

is = DiscreteIdentitySystem(2)
#s = DiscreteLinearControlSystem([1. 0; 0 1], reshape([1.; 0], 2, 1), U)
const h = 0.4
s = DiscreteLinearControlSystem([1 h; 0 0], reshape([0; h], 2, 1)) # TODO [0; 1] should be equivalent and more correct

sw = AutonomousSwitching()

using SemialgebraicSets
fs = FullSpace()

M = LightGraphs.ne(A.G)

S = Fill(is, N)
Gu = Fill(fs, M)
Re = Fill(s, M)
Sw = Fill(sw, N)

#I = [P0, P0, P0, P0 ∩ Pa]
I = [P0, P0 ∩ Pa]

hs = HybridSystem(A, S, I, Gu, Re, Sw)


text:
using LinearAlgebra
using Plots

function flower(x, a=1, b=1, c=4)
    return a * norm(x) + b * sin(c * atan(x[2], x[1]));
end

contour(-3:0.1:3, -3:0.1:3, (x, y)->flower([x, y]))