text:
@doc raw"""
    device_timeseries_ub(canonical::Canonical,
                     ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}}
                     cons_name::Symbol,
                     var_name::Symbol)

Constructs upper bound for given variable and time series data and a multiplier.

# Constraint

```variable[name, t] <= ts_data[2][ix]*ts_data[3][ix][t] ```

# LaTeX

`` x_t \leq r^{val} r_t, \forall t ``

# Arguments
* canonical::Canonical : the canonical model built in PowerSimulations
* ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}}: timeseries data name (1), multiplier (2) and values (3)
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the variable
"""
function device_timeseries_ub(canonical::Canonical,
                              ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}},
                              cons_name::Symbol,
                              var_name::Symbol)

    time_steps = model_time_steps(canonical)
    names = (v[1] for v in ts_data)
    constraint = _add_cons_container!(canonical, cons_name, names, time_steps)
    variable = get_variable(canonical, var_name)

    for t in time_steps, data in ts_data
        name = data[1]
        forecast = data[4][t]
        multiplier = data[3]
        constraint[name, t] = JuMP.@constraint(canonical.JuMPmodel,
                                    variable[name, t] <= multiplier*forecast)

    end

    return

end

@doc raw"""
    device_timeseries_lb(canonical::Canonical,
                     ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}}
                     cons_name::Symbol,
                     var_name::Symbol)

Constructs lower bound for given variable subject to time series data and a multiplier.

# Constraint

``` ts_data[2][ix]*ts_data[3][ix][t] <= variable[name, t] ```

# LaTeX

`` r^{val} r_t \leq x_t, \forall t ``

where (ix, name) in enumerate(ts_data[1]).

# Arguments
* canonical::Canonical : the canonical model built in PowerSimulations
* ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}}: timeseries data name (1), multiplier (2) and values (3)
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the variable
"""
function device_timeseries_lb(canonical::Canonical,
                              ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}},
                              cons_name::Symbol,
                              var_name::Symbol)

    time_steps = model_time_steps(canonical)
    names = (v[1] for v in ts_data)
    constraint =_add_cons_container!(canonical, cons_name, names, time_steps)
    variable = get_variable(canonical, var_name)

    for t in time_steps, data in ts_data
        name = data[1]
        forecast = data[4][t]
        multiplier = data[3]
        constraint[name, t] = JuMP.@constraint(canonical.JuMPmodel,
                        multiplier*forecast <= variable[name, t])

    end

    return

end

#NOTE: there is a floating, unnamed lower bound constraint in this function. This may need to be changed.
@doc raw"""
    device_timeseries_param_ub(canonical::Canonical,
                                    ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}}
                                    cons_name::Symbol,
                                    param_reference::UpdateRef,
                                    var_name::Symbol)

Constructs upper bound for given variable using a parameter. The constraint is
    built with a time series data vector and a multiplier

# Constraint

``` variable[name, t] <= ts_data[2][ix]*param[name, t] ```

# LaTeX

`` x^{var}_t \leq r^{val} x^{param}_t, \forall t ``

# Arguments
* canonical::Canonical : the canonical model built in PowerSimulations
* ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}}: timeseries data name (1), multiplier (2) and values (3)
* cons_name::Symbol : name of the constraint
* param_reference::UpdateRef : UpdateRef to access the parameter
* var_name::Symbol : the name of the variable
"""
function device_timeseries_param_ub(canonical::Canonical,
                                    ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}},
                                    cons_name::Symbol,
                                    param_reference::UpdateRef,
                                    var_name::Symbol)

    time_steps = model_time_steps(canonical)
    ub_name = _middle_rename(cons_name, "_", "ub")
    variable = get_variable(canonical, var_name)
    names = (v[1] for v in ts_data)
    constraint = _add_cons_container!(canonical, ub_name, names, time_steps)
    param = _add_param_container!(canonical, param_reference, names, time_steps)

    for t in time_steps, data in ts_data
        name = data[1]
        forecast = data[4][t]
        multiplier = data[3]
        param[name, t] = PJ.add_parameter(canonical.JuMPmodel, forecast)
        constraint[name, t] = JuMP.@constraint(canonical.JuMPmodel,
                                        variable[name, t] <= multiplier*param[name, t])
    end

    return

end

@doc raw"""
    device_timeseries_param_lb(canonical::Canonical,
                                    ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}}
                                    cons_name::Symbol,
                                    param_reference::UpdateRef,
                                    var_name::Symbol)

Constructs lower bound for given variable using a parameter. The constraint is
    built with a time series data vector and a multiplier

# Constraint

``` ts_data[2][ix] * param[name, t] <= variable[name, t] ```

# LaTeX

`` r^{val} x^{param}_t \leq x^{var}_t, \forall t ``

# Arguments
* canonical::Canonical : the canonical model built in PowerSimulations
* ts_data::Tuple{Vector{String}, Vector{Vector{Float64}}} : timeseries data name (1) and values (2)
* cons_name::Symbol : name of the constraint
* param_reference::UpdateRef : UpdateRef to access the parameter
* var_name::Symbol : the name of the variable
"""
function device_timeseries_param_lb(canonical::Canonical,
                                    ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}},
                                    cons_name::Symbol,
                                    param_reference::UpdateRef,
                                    var_name::Symbol)

    time_steps = model_time_steps(canonical)
    variable = get_variable(canonical, var_name)
    lb_name = _middle_rename(cons_name, "_", "lb")
    names = (v[1] for v in ts_data)
    constraint =_add_cons_container!(canonical, lb_name, names, time_steps)
    param =_add_param_container!(canonical, param_reference, names, time_steps)

    for t in time_steps, data in ts_data
        name = data[1]
        forecast = data[4][t]
        multiplier = data[3]
        param[name, t] = PJ.add_parameter(canonical.JuMPmodel, forecast)
        constraint[name, t] = JuMP.@constraint(canonical.JuMPmodel, multiplier*param[name, t] <= variable[name, t])
    end

    return

end

@doc raw"""
    device_timeseries_ub_bin(canonical::Canonical,
                                    ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}}
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    binvar_name::Symbol)

Constructs upper bound for variable and time series or confines to 0 depending on binary variable.
    The upper bound is defined by a time series and a multiplier.

# Constraints

``` varcts[name, t] <= varbin[name, t]* ts_data[2][ix] * ts_data[3][ix][t] ```

where (ix, name) in enumerate(ts_data[1]).

# LaTeX

`` x^{cts}_t \leq r^{val} r_t x^{bin}_t, \forall t ``

# Arguments
* canonical::Canonical : the canonical model built in PowerSimulations
* ts_data::Tuple{Vector{String}, Vector{Vector{Float64}}} : timeseries data name (1) and values (2)
* cons_name::Symbol : name of the constraint
* var_name::Symbol :  name of the variable
* binvar_name::Symbol : name of binary variable
"""
function device_timeseries_ub_bin(canonical::Canonical,
                                    ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}},
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    binvar_name::Symbol)

    time_steps = model_time_steps(canonical)
    ub_name = _middle_rename(cons_name, "_", "ub")

    varcts = get_variable(canonical, var_name)
    varbin = get_variable(canonical, binvar_name)

    names = (v[1] for v in ts_data)
    con_ub =_add_cons_container!(canonical, ub_name, names, time_steps)

    for t in time_steps, data in ts_data
        name = data[1]
        forecast = data[4][t]
        multiplier = data[3]
        con_ub[name, t] = JuMP.@constraint(canonical.JuMPmodel, varcts[name, t] <= varbin[name, t]*multiplier*forecast)
    end

    return

end

@doc raw"""
    device_timeseries_ub_bigM(canonical::Canonical,
                                    ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}}
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    param_reference::UpdateRef,
                                    binvar_name::Symbol,
                                    M_value::Float64 = 1e6)

Constructs upper bound for variable and time series and a multiplier or confines to 0 depending on binary variable.
    Uses BigM constraint type to allow for parameter since ParameterJuMP doesn't support var*parameter

# Constraints

``` varcts[name, t] - ts_data[2][ix] * param[name, t] <= (1 - varbin[name, t])*M_value ```

``` varcts[name, t] <= varbin[name, t]*M_value ```

# LaTeX

`` x^{cts}_t - r^{val} x^{param}_t \leq M(1 - x^{bin}_t ), forall t ``

`` x^{cts}_t \leq M x^{bin}_t, \forall t ``

# Arguments
* canonical::Canonical : the canonical model built in PowerSimulations
* ts_data::Tuple{Vector{String}, Vector{Vector{Float64}}} : timeseries data name (1) and values (2)
* cons_name::Symbol : name of the constraint
* var_name::Symbol :  name of the variable
param_reference::UpdateRef : UpdateRef of access the parameters
* binvar_name::Symbol : name of binary variable
* M_value::Float64 : bigM
"""
function device_timeseries_ub_bigM(canonical::Canonical,
                                    ts_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}},
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    param_reference::UpdateRef,
                                    binvar_name::Symbol,
                                    M_value::Float64 = 1e6)

    time_steps = model_time_steps(canonical)
    ub_name = _middle_rename(cons_name, "_", "ub")
    key_status = _middle_rename(cons_name, "_", "status")

    varcts = get_variable(canonical, var_name)
    varbin = get_variable(canonical, binvar_name)
    names = (v[1] for v in ts_data)
    con_ub = _add_cons_container!(canonical, ub_name, names, time_steps)
    con_status =_add_cons_container!(canonical, key_status, names, time_steps)
    param =_add_param_container!(canonical, param_reference, names, time_steps)

    for t in time_steps, data in ts_data
        name = data[1]
        forecast = data[4][t]
        multiplier = data[3]
        param[name, t] = PJ.add_parameter(canonical.JuMPmodel, forecast)
        con_ub[name, t] = JuMP.@constraint(canonical.JuMPmodel, varcts[name, t] - param[name, t]*multiplier <= (1 - varbin[name, t])*M_value)
        con_status[name, t] =  JuMP.@constraint(canonical.JuMPmodel, varcts[name, t] <= varbin[name, t]*M_value)
    end

    return

end


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

name = "EarCut"
version = v"2.1.5"
# Collection of sources required to build Clipper
sources = [
    "https://github.com/SimonDanisch/EarCutBuilder.git" =>
    "cfa4233e26ac785a89954a72b4e2e84312b389c2",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd EarCutBuilder/
${CXX} -c -fPIC -std=c++11 cwrapper.cpp -o earcut.o
libdir="lib"
if [[ ${target} == *-mingw32 ]]; then     libdir="bin"; else     libdir="lib"; fi
mkdir ${prefix}/${libdir}
${CXX} $LDFLAGS -shared -o ${prefix}/${libdir}/earcut.${dlext} earcut.o
exit
"""

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

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "earcut", :earcut)
]

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

]

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

ext: jl


text:
__precompile__()

module ElasticNetVAR

	# Libraries
	using Distributed, Dates, Logging;
	using LinearAlgebra, Distributions, Statistics;
	using TSAnalysis;

	# Aliases for TSAnalysis
	compute_J1 = TSAnalysis.compute_J1;
	backwards_pass = TSAnalysis.backwards_pass;

	# Custom dependencies
	const local_path = dirname(@__FILE__);
	include("$local_path/types.jl");
	include("$local_path/methods.jl");
	include("$local_path/coordinate_descent.jl");
	include("$local_path/ecm.jl");
	include("$local_path/jackknife.jl");
	include("$local_path/validation.jl");

	# Export
	export JVector, JArray, EstimSettings, ValidationSettings, HyperGrid;
	export companion_form, ext_companion_form, no_combinations, rand_without_replacement;
	export coordinate_descent, ecm;
	export coordinate_descent, build_Γ;
	export block_jackknife, artificial_jackknife, optimal_d;
	export select_hyperparameters, fc_err, jackknife_err;
end


text:
"""
    SymbolicNode

SymbolicNode is the basic building block of the symbolic graph in MXNet.jl.

    (self :: SymbolicNode)(args :: SymbolicNode...)
    (self :: SymbolicNode)(; kwargs...)

Make a new node by composing `self` with `args`. Or the arguments
can be specified using keyword arguments.
"""
type SymbolicNode
  handle :: MX_SymbolHandle
end
function Base.unsafe_convert(::Type{MX_handle}, obj::SymbolicNode)
  Base.unsafe_convert(MX_handle, obj.handle)
end
Base.convert(t::Type{MX_handle}, obj::SymbolicNode) = Base.unsafe_convert(t, obj)
Base.cconvert(t::Type{MX_handle}, obj::SymbolicNode) = Base.unsafe_convert(t, obj)

"""
    deepcopy(self :: SymbolicNode)

Make a deep copy of a SymbolicNode.
"""
function Base.deepcopy(self :: SymbolicNode)
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCopy, (MX_handle, Ref{MX_handle}), self, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    copy(self :: SymbolicNode)

Make a copy of a SymbolicNode. The same as making a deep copy.
"""
function Base.copy(self :: SymbolicNode)
  Base.deepcopy(self)
end

@compat function (self::SymbolicNode)(args :: SymbolicNode...)
  s = deepcopy(self)
  _compose!(s, args...)
end
@compat function (self::SymbolicNode)(;kwargs...)
  s = deepcopy(self)
  _compose!(s; kwargs...)
end

macro _list_symbol_info(self, func_name)
  quote
    ref_sz    = Ref{MX_uint}(0)
    ref_names = Ref{char_pp}(0)
    @mxcall($func_name, (MX_handle, Ref{MX_uint}, Ref{char_pp}),
            $self, ref_sz, ref_names)
    narg = ref_sz[]
    names = unsafe_wrap(Array, ref_names[], narg)
    names = [Symbol(unsafe_wrap(String, x)) for x in names]
    return names
  end
end

"""
    list_arguments(self :: SymbolicNode)

List all the arguments of this node. The argument for a node contains both
the inputs and parameters. For example, a `FullyConnected` node will
have both data and weights in its arguments. A composed node (e.g. a MLP) will
list all the arguments for intermediate nodes.

Returns a list of symbols indicating the names of the arguments.
"""
function list_arguments(self :: SymbolicNode)
  @_list_symbol_info(self, :MXSymbolListArguments)
end

"""
    list_outputs(self :: SymbolicNode)

List all the outputs of this node.

Returns a list of symbols indicating the names of the outputs.
"""
function list_outputs(self :: SymbolicNode)
  @_list_symbol_info(self, :MXSymbolListOutputs)
end

"""
    list_auxiliary_states(self :: SymbolicNode)

List all auxiliary states in the symbool.

Auxiliary states are special states of symbols that do not corresponds to an argument,
and do not have gradient. But still be useful for the specific operations.
A common example of auxiliary state is the moving_mean and moving_variance in BatchNorm.
Most operators do not have Auxiliary states.

Returns a list of symbols indicating the names of the auxiliary states.
"""
function list_auxiliary_states(self :: SymbolicNode)
  @_list_symbol_info(self, :MXSymbolListAuxiliaryStates)
end

"""
    get_internals(self :: SymbolicNode)

Get a new grouped `SymbolicNode` whose output contains all the internal outputs of
this `SymbolicNode`.
"""
function get_internals(self :: SymbolicNode)
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolGetInternals, (MX_handle, Ref{MX_handle}), self, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    get_attr(self :: SymbolicNode, key :: Symbol)

Get attribute attached to this `SymbolicNode` belonging to key.

Returns the value belonging to key as a `Nullable`.
"""
function get_attr(self :: SymbolicNode, key :: Symbol)
  key_s = string(key)
  ref_out = Ref{Cstring}()
  ref_success = Ref{Cint}(-1)
  @mxcall(:MXSymbolGetAttr, (MX_handle, Cstring, Ref{Cstring}, Ref{Cint}),
          self, key_s, ref_out, ref_success)
  if ref_success[] == 1
    return Nullable{String}(unsafe_string(ref_out[]))
  else
    return Nullable{String}()
  end
end

"""
    list_attr(self :: SymbolicNode)

Get all attributes from a symbol.

Returns a dictionary of attributes.
"""
function list_attr(self :: SymbolicNode)
  ref_sz    = Ref{MX_uint}(0)
  ref_strings = Ref{char_pp}(0)
  @mxcall(:MXSymbolListAttrShallow, (MX_handle, Ref{MX_uint}, Ref{char_pp}),
            self, ref_sz, ref_strings)
  narg = 2*ref_sz[]
  strings = unsafe_wrap(Array, ref_strings[], narg)
  out = Dict{Symbol, String}()
  for i in 1:2:narg
    key = Symbol(unsafe_wrap(String, strings[i]))
    value = unsafe_string(strings[i+1]) # Creates a copy of string
    out[key] = value
  end
  return out
end

"""
    list_all_attr(self :: SymbolicNode)

Get all attributes from the symbol graph.

Returns a dictionary of attributes.
"""
function list_all_attr(self :: SymbolicNode)
  ref_sz    = Ref{MX_uint}(0)
  ref_strings = Ref{char_pp}(0)
  @mxcall(:MXSymbolListAttr, (MX_handle, Ref{MX_uint}, Ref{char_pp}),
            self, ref_sz, ref_strings)
  narg = 2*ref_sz[]
  strings = unsafe_wrap(Array, ref_strings[], narg)
  out = Dict{Symbol, String}()
  for i in 1:2:narg
    key = Symbol(unsafe_wrap(String, strings[i]))
    value = unsafe_string(strings[i+1])
    out[key] = value
  end
  return out
end

"""
    set_attr(self:: SymbolicNode, key :: Symbol, value :: AbstractString)

Set the attribute key to value for this `SymbolicNode`.

!!! note
    It is encouraged not to call this function directly, unless you know exactly what you are doing. The
    recommended way of setting attributes is when creating the `SymbolicNode`. Changing
    the attributes of a `SymbolicNode` that is already been used somewhere else might
    cause unexpected behavior and inconsistency.
"""
function set_attr(self :: SymbolicNode, key :: Symbol, value :: AbstractString)
  key_s = string(key)
  value_s = String(value)

  @mxcall(:MXSymbolSetAttr, (MX_handle, Cstring, Cstring), self, key_s, value_s)
end

"""
    grad(self :: SymbolicNode, wrt :: Vector{SymbolicNode})

Get the autodiff gradient of the current `SymbolicNode`. This function can
only be used if the current symbol is a loss function.

# Arguments:
* `self::SymbolicNode`: current node.
* `wrt::Vector{Symbol}`: the names of the arguments to the gradient.

Returns a gradient symbol of the corresponding gradient.
"""
function grad(self :: SymbolicNode, wrt :: Vector{Symbol})
  hdr_ref = Ref{MX_handle}(0)
  keys = String[string(key) for key in wrt]

  @mxcall(:MXSymbolGrad, (MX_handle, MX_uint, char_pp, Ptr{MX_handle}), self, length(keys), keys, hdr_ref)
  return SymbolicNode(MX_SymbolHandle(hdr_ref[]))
end

"""
    Variable(name :: Union{Symbol, AbstractString})

Create a symbolic variable with the given name. This is typically used as a placeholder.
For example, the data node, acting as the starting point of a network architecture.

# Arguments
* Dict{Symbol, AbstractString} attrs: The attributes associated with this `Variable`.
"""
function Variable(name :: Union{Symbol, AbstractString}; attrs = Dict())
  attrs = convert(Dict{Symbol, AbstractString}, attrs)
  hdr_ref = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCreateVariable, (char_p, Ref{MX_handle}), name, hdr_ref)
  node = SymbolicNode(MX_SymbolHandle(hdr_ref[]))
  for (k, v) in attrs
    set_attr(node, k, v)
  end
  node
end

"""
    Group(nodes :: SymbolicNode...)

Create a `SymbolicNode` by grouping nodes together.
"""
function Group(nodes :: SymbolicNode...)
  handles = MX_handle[nodes...]
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCreateGroup, (MX_uint, Ptr{MX_handle}, Ref{MX_handle}),
          length(handles), handles, ref_hdr)
  SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

function _build_shapes(shape_size::MX_uint, shape_ndim::Ptr{MX_uint}, shape_data::Ptr{Ptr{MX_uint}})
  shape_ndim = unsafe_wrap(Array, shape_ndim, shape_size)
  shape_data = unsafe_wrap(Array, shape_data, shape_size)
  shapes = map(1:shape_size) do i
    my_shape = unsafe_wrap(Array, shape_data[i], shape_ndim[i])
    tuple(flipdim(Int[my_shape...],1)...)
  end
  convert(Vector{Tuple}, shapes)
end

function _infer_shape(self, keys, indptr, sdata)
  ref_arg_shape_size = Ref{MX_uint}(0)
  ref_arg_shape_ndim = Ref{Ptr{MX_uint}}(0)
  ref_arg_shape_data = Ref{Ptr{Ptr{MX_uint}}}(0)
  ref_out_shape_size = Ref{MX_uint}(0)
  ref_out_shape_ndim = Ref{Ptr{MX_uint}}(0)
  ref_out_shape_data = Ref{Ptr{Ptr{MX_uint}}}(0)
  ref_aux_shape_size = Ref{MX_uint}(0)
  ref_aux_shape_ndim = Ref{Ptr{MX_uint}}(0)
  ref_aux_shape_data = Ref{Ptr{Ptr{MX_uint}}}(0)
  ref_complete       = Ref{Cint}(0)
  @mxcall(:MXSymbolInferShape,
          (MX_handle, MX_uint, char_pp, Ptr{MX_uint}, Ptr{MX_uint},
           Ref{MX_uint}, Ref{Ptr{MX_uint}}, Ref{Ptr{Ptr{MX_uint}}},
           Ref{MX_uint}, Ref{Ptr{MX_uint}}, Ref{Ptr{Ptr{MX_uint}}},
           Ref{MX_uint}, Ref{Ptr{MX_uint}}, Ref{Ptr{Ptr{MX_uint}}},
           Ref{Cint}),
          self, length(indptr)-1, keys, indptr, sdata,
          ref_arg_shape_size, ref_arg_shape_ndim, ref_arg_shape_data,
          ref_out_shape_size, ref_out_shape_ndim, ref_out_shape_data,
          ref_aux_shape_size, ref_aux_shape_ndim, ref_aux_shape_data,
          ref_complete)
  if ref_complete[] == 0
    return (nothing, nothing, nothing)
  else
    return (
      _build_shapes(ref_arg_shape_size[], ref_arg_shape_ndim[], ref_arg_shape_data[]),
      _build_shapes(ref_out_shape_size[], ref_out_shape_ndim[], ref_out_shape_data[]),
      _build_shapes(ref_aux_shape_size[], ref_aux_shape_ndim[], ref_aux_shape_data[])
    )
  end
end

"""
    infer_shape(self :: SymbolicNode, args...)
    infer_shape(self :: SymbolicNode; kwargs...)

Do shape inference according to the input shapes. The input shapes could be provided
as a list of shapes, which should specify the shapes of inputs in the same order as
the arguments returned by [`list_arguments`](@ref). Alternatively, the shape information
could be specified via keyword arguments.

Returns a 3-tuple containing shapes of all the arguments, shapes of all the outputs and
shapes of all the auxiliary variables. If shape inference failed due to incomplete
or incompatible inputs, the return value will be `(nothing, nothing, nothing)`.
"""
function infer_shape(self :: SymbolicNode; kwargs...)
  sdata  = MX_uint[]
  indptr = MX_uint[0]
  for (k,v) in kwargs
    append!(sdata, flipdim([v...],1))
    push!(indptr, length(sdata))
  end
  keys = AbstractString[string(x[1]) for x in kwargs]
  _infer_shape(self, keys, indptr, sdata)
end
function infer_shape(self :: SymbolicNode, args :: Union{Tuple, Void}...)
  sdata  = MX_uint[]
  indptr = MX_uint[0]
  for arg in args
    if isa(arg, Void); continue; end
    append!(sdata, flipdim([arg...],1))
    push!(indptr, length(sdata))
  end
  keys = Ptr{char_p}(0)
  _infer_shape(self, keys, indptr, sdata)
end

function _infer_type(self, keys, arg_type_data)
  ref_in_type_size  = Ref{MX_uint}()
  ref_in_type_data  = Ref{Ptr{Cint}}()
  ref_out_type_size = Ref{MX_uint}()
  ref_out_type_data = Ref{Ptr{Cint}}()
  ref_aux_type_size = Ref{MX_uint}()
  ref_aux_type_data = Ref{Ptr{Cint}}()
  ref_complete      = Ref{Cint}()

  @mxcall(:MXSymbolInferType,
          (MX_handle, MX_uint, char_pp, Ptr{Cint},
           Ref{MX_uint}, Ref{Ptr{Cint}},
           Ref{MX_uint}, Ref{Ptr{Cint}},
           Ref{MX_uint}, Ref{Ptr{Cint}},
           Ref{Cint}),
          self, length(arg_type_data)-1, keys, arg_type_data,
          ref_in_type_size, ref_in_type_data,
          ref_out_type_size, ref_out_type_data,
          ref_aux_type_size, ref_aux_type_data,
          ref_complete)

  if ref_complete[] == 0
    return (nothing, nothing, nothing)
  else
    in_type = unsafe_wrap(Array, ref_in_type_data[], ref_in_type_size[])
    out_type = unsafe_wrap(Array, ref_out_type_data[], ref_out_type_size[])
    aux_type = unsafe_wrap(Array, ref_aux_type_data[], ref_aux_type_size[])
    return ([fromTypeFlag(TypeFlag(t)) for t in in_type],
            [fromTypeFlag(TypeFlag(t)) for t in out_type],
            [fromTypeFlag(TypeFlag(t)) for t in aux_type])
  end
end

"""
    infer_type(self :: SymbolicNode; kwargs...)
    infer_type(self :: SymbolicNode, args...)

Do type inference according to the input types. The input types could be provided
as a list of types, which should specify the types of inputs in the same order as
the arguments returned by [`list_arguments`](@ref). Alternatively, the type information
could be specified via keyword arguments.

Returns a 3-tuple containing types of all the arguments, types of all the outputs and
types of all the auxiliary variables. If type inference failed due to incomplete
or incompatible inputs, the return value will be `(nothing, nothing, nothing)`.
"""
function infer_type(self :: SymbolicNode; kwargs...)
  types = Cint[toTypeFlag(x[2]) for x in kwargs]
  keys = AbstractString[string(x[1]) for x in kwargs]
  _infer_type(self, keys, types)
end

function infer_type(self :: SymbolicNode, args :: Union{Tuple, Void}...)
  types = Cint[]
  keys = Ptr{char_p}(0)

  for arg in args
    if isa(arg, Void); continue; end
    push!(types, toTypeFlag(arg))
  end
  _infer_type(self, keys, types)
end

"""
    getindex(self :: SymbolicNode, idx :: Union{Int, Base.Symbol, AbstractString})

Get a node representing the specified output of this node. The index could be
a symbol or string indicating the name of the output, or a 1-based integer
indicating the index, as in the list of [`list_outputs`](@ref).
"""
function Base.getindex(self :: SymbolicNode, idx :: Union{Base.Symbol, AbstractString})
  idx   = Symbol(idx)
  i_idx = find(idx .== list_outputs(self))
  @assert(length(i_idx) > 0, "Cannot find output with name '$idx'")
  @assert(length(i_idx) < 2, "Found duplicated output with name '$idx'")
  Base.getindex(self, i_idx[1])
end
function Base.getindex(self :: SymbolicNode, idx :: Int)
  ref_hdr = Ref{MX_handle}(0)
  # note Julia is 1-based, while MXNet is 0-based
  @mxcall(:MXSymbolGetOutput, (MX_handle, MX_uint, Ref{MX_handle}), self, idx-1, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

import Base: +, .+
function +(self :: SymbolicNode, args :: Union{SymbolicNode,Real}...)
  ret = self
  for arg in args
    if isa(arg, SymbolicNode)
      ret = _Plus(ret, arg)
    else
      ret = _PlusScalar(ret, scalar=MX_float(arg))
    end
  end
  ret
end
function .+(self :: SymbolicNode, args :: Union{SymbolicNode,Real}...)
  +(self, args...)
end
function +(s1 :: Real, self :: SymbolicNode, args :: Union{SymbolicNode,Real}...)
  +(self, s1, args...)
end
function .+(s1 :: Real, self :: SymbolicNode, args :: Union{SymbolicNode,Real}...)
  +(self, s1, args...)
end

import Base: -, .-
function -(self :: SymbolicNode, arg :: SymbolicNode)
  _Minus(self, arg)
end
function .-(self :: SymbolicNode, arg :: SymbolicNode)
  -(self, arg)
end
function -(self :: SymbolicNode, arg :: Real)
  _MinusScalar(self, scalar=MX_float(arg))
end
function .-(self :: SymbolicNode, arg :: Real)
  -(self, arg)
end

function -(arg :: Real, self :: SymbolicNode)
  _RMinusScalar(self, scalar=arg)
end
function .-(arg :: Real, self :: SymbolicNode)
  -(arg, self)
end

function -(self :: SymbolicNode)
  -(0, self)
end

import Base: .*, *
function .*(self :: SymbolicNode, args :: Union{SymbolicNode,Real}...)
  ret = self
  for arg in args
    if isa(arg, SymbolicNode)
      ret = _Mul(ret, arg)
    else
      ret = _MulScalar(ret, scalar=MX_float(arg))
    end
  end
  ret
end
function .*(arg :: Real, self :: SymbolicNode, args :: Union{SymbolicNode,Real}...)
  .*(self, arg, args...)
end
function *(arg :: Real, self :: SymbolicNode)
  _MulScalar(self, scalar=arg)
end
function *(self :: SymbolicNode, arg :: Real)
  *(arg, self)
end

import Base: ./, /
function ./(self :: SymbolicNode, arg :: SymbolicNode)
  _Div(self, arg)
end
function ./(self :: SymbolicNode, arg :: Real)
  _DivScalar(self, scalar=MX_float(arg))
end
function /(self :: SymbolicNode, arg :: Real)
  ./(self, arg)
end
function ./(arg :: Real, self :: SymbolicNode)
  _RDivScalar(self, scalar=arg)
end

import Base: .^, ^
function .^(self :: SymbolicNode, pow :: SymbolicNode)
  _Power(self, pow)
end
function .^(self :: SymbolicNode, pow :: AbstractFloat)
  _PowerScalar(self, scalar=pow)
end
function ^(self :: SymbolicNode, pow :: AbstractFloat)
  .^(self, pow)
end

function _compose!(node :: SymbolicNode; kwargs...)
  name     = char_p(0)
  arg_keys = AbstractString[]
  arg_vals = MX_handle[]

  for (k,v) in kwargs
    if k == :name
      name = string(v)
    else
      @assert(isa(v, SymbolicNode), "Compose expect `SymbolicNode` as arguments")
      push!(arg_keys, string(k))
      push!(arg_vals, v)
    end
  end

  @mxcall(:MXSymbolCompose,
          (MX_handle, char_p, MX_uint, Ptr{char_p}, Ptr{MX_handle}),
          node, name, length(arg_keys), arg_keys, arg_vals)
  return node
end
function _compose!(node :: SymbolicNode, args::SymbolicNode...)
  _compose!(node, char_p(0), args...)
end
function _compose!(node :: SymbolicNode, name :: Union{Base.Symbol, char_p}, args::SymbolicNode...)
  if isa(name, Base.Symbol); name = string(name); end
  arg_keys = Ptr{char_p}(0)
  arg_vals = MX_handle[args...]

  @mxcall(:MXSymbolCompose,
          (MX_handle, char_p, MX_uint, Ptr{char_p}, Ptr{MX_handle}),
          node, name, length(arg_vals), arg_keys, arg_vals)
  return node
end

"""
    to_json(self :: SymbolicNode)

Convert a `SymbolicNode` into a JSON string.
"""
function to_json(self :: SymbolicNode)
  ref_json = Ref{char_p}(0)
  @mxcall(:MXSymbolSaveToJSON, (MX_handle, Ref{char_p}), self, ref_json)
  return unsafe_string(ref_json[])
end

"""
    from_json(repr :: AbstractString, ::Type{SymbolicNode})

Load a `SymbolicNode` from a JSON string representation.
"""
function from_json(repr :: AbstractString, ::Type{SymbolicNode})
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCreateFromJSON, (char_p, Ref{MX_handle}), repr, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    load(filename :: AbstractString, ::Type{SymbolicNode})

Load a `SymbolicNode` from a JSON file.
"""
function load(filename :: AbstractString, ::Type{SymbolicNode})
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCreateFromFile, (char_p, Ref{MX_handle}), filename, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    save(filename :: AbstractString, node :: SymbolicNode)

Save a `SymbolicNode` to a JSON file.
"""
function save(filename :: AbstractString, node :: SymbolicNode)
  @mxcall(:MXSymbolSaveToFile, (MX_handle, char_p), node, filename)
end

################################################################################
# Atomic SymbolicNode functions dynamically imported from libmxnet
################################################################################
function _define_atomic_symbol_creator(hdr :: MX_handle)
  ref_name      = Ref{char_p}(0)
  ref_desc      = Ref{char_p}(0)
  ref_kv_nargs  = Ref{char_p}(0)
  ref_nargs     = Ref{MX_uint}(0)
  ref_arg_names = Ref{char_pp}(0)
  ref_arg_types = Ref{char_pp}(0)
  ref_arg_descs = Ref{char_pp}(0)
  ref_ret_type  = Ref{char_p}(0)

  @mxcall(:MXSymbolGetAtomicSymbolInfo,
          (MX_handle, Ref{char_p}, Ref{char_p}, Ref{MX_uint}, Ref{char_pp}, Ref{char_pp},
           Ref{char_pp}, Ref{char_p}, Ref{char_p}),
          hdr, ref_name, ref_desc, ref_nargs, ref_arg_names, ref_arg_types, ref_arg_descs,
          ref_kv_nargs, ref_ret_type)

  func_name_s = unsafe_wrap(String, ref_name[])
  func_name   = Symbol(func_name_s)
  kv_nargs_s  = unsafe_wrap(String, ref_kv_nargs[])
  kv_nargs    = Symbol(kv_nargs_s)

  signature = _format_signature(Int(ref_nargs[]), ref_arg_names)
  f_desc = "    " * func_name_s * "(" * signature * ")\n\n"
  f_desc *= unsafe_wrap(String, ref_desc[]) * "\n\n"
  if !isempty(kv_nargs_s)
    f_desc *= "This function support variable length positional `SymbolicNode` inputs.\n\n"
  end
  f_desc *= "# Arguments\n"
  f_desc *= _format_docstring(Int(ref_nargs[]), ref_arg_names, ref_arg_types, ref_arg_descs)
  f_desc *= "* `name::Symbol`: The name of the `SymbolicNode`. (e.g. `:my_symbol`), optional.\n"
  f_desc *= "* `attrs::Dict{Symbol, AbstractString}`: The attributes associated with this `SymbolicNode`.\n\n"
  f_desc *= "Returns `$(_format_typestring(unsafe_wrap(String, ref_ret_type[])))`."

  func_def = quote
  @doc $f_desc ->
  function $func_name(args::SymbolicNode...; kwargs...)
    idx = findfirst(x -> x[1] == :name, kwargs)
    if idx > 0
      name = kwargs[idx][2]
    else
      name = ""
    end

    # XXX: hacky way of solving the problem that the arguments of `dot` should be swapped
    # See https://github.com/dmlc/MXNet.jl/issues/55
    if $func_name_s == "dot"
      args = reverse(args)
    end

    # XXX: hacky way of solving the semantic difference of the axes parameter in Julia
    # and in libmxnet.
    # See https://github.com/dmlc/MXNet.jl/pull/123
    if $func_name_s == "transpose"
      kwargs = Any[key != :axes ? (key, arg) : (key, reverse(map(i->length(arg)-i, arg))) for (key, arg) in kwargs]
    end

    param_keys = String[]
    param_vals = String[]
    symbol_kws = Dict{Symbol, SymbolicNode}()
    attrs = Dict{Symbol, String}()

    $(if kv_nargs != Symbol("")
      quote
        if !in($kv_nargs_s, param_keys)
          push!(param_keys, $kv_nargs_s)
          push!(param_vals, string(length(args)))
        end
      end
    end)

    for (k,v) in kwargs
      if k == :name; continue; end
      if isa(v, SymbolicNode)
        symbol_kws[k] = v
      elseif k == :attrs
        if isa(v, Dict)
          attrs = convert(Dict{Symbol, AbstractString}, v)
        else
          throw(ArgumentError("attrs needs to be a Dictionary"))
        end
      else
        push!(param_keys, string(k))
        push!(param_vals, dump_mx_param(v))
      end
    end

    if length(args) != 0 && length(symbol_kws) != 0
      @assert(false, $func_name_s * " only accepts Symbols either as positional or keyword arguments, not both.")
    end
    $(if kv_nargs != Symbol("")
      quote
        if length(symbol_kws) > 0
          @assert(false, $func_name_s * " takes variable number of SymbolicNode arguments, " *
                         "please pass input Symbols via positional arguments, instead of keyword arguments.")
        end
      end
    end)

    local hdr = _get_symbol_creator($(QuoteNode(func_name)))

    # create the SymbolicNode
    ref_sym_hdr = Ref{MX_handle}()
    @mxcall(:MXSymbolCreateAtomicSymbol,
            (MX_handle, MX_uint, Ptr{char_p}, Ptr{char_p}, Ref{MX_handle}),
            hdr, length(param_keys), param_keys, param_vals, ref_sym_hdr)
    sym_hdr = ref_sym_hdr[]

    node = SymbolicNode(MX_SymbolHandle(sym_hdr))
    hint = lowercase($func_name_s)
    name = get!(DEFAULT_NAME_MANAGER, name, hint)

    # set attrs
    for (k, v) in attrs
      set_attr(node, k, v)
    end

    if length(args) != 0
      _compose!(node, name, args...)
    else
      _compose!(node; name=name, symbol_kws...)
    end

    return node
  end # function
  end # quote
  return func_def
end

function _get_atomic_symbol_creators()
  n_ref = Ref{MX_uint}(0)
  h_ref = Ref{Ptr{MX_handle}}(0)
  @mxcall(:MXSymbolListAtomicSymbolCreators, (Ref{MX_uint}, Ref{Ptr{MX_handle}}), n_ref, h_ref)

  return unsafe_wrap(Array, h_ref[], n_ref[])
end

function _get_atomic_symbol_name(handle :: MX_handle)
  name_r = Ref{char_p}(0)
  @mxcall(:MXSymbolGetAtomicSymbolName, (MX_handle, Ref{char_p}), handle, name_r)
  return unsafe_wrap(String, name_r[])
end

const _symbol_creator_cache = Dict{Symbol, MX_handle}()
function _populate_symbol_creator_cache!()
  empty!(_symbol_creator_cache)
  h_creators = _get_atomic_symbol_creators()
  for handle in h_creators
    name = Symbol(_get_atomic_symbol_name(handle))
    _symbol_creator_cache[name] = handle
  end
end

_get_symbol_creator(name :: Symbol) = _symbol_creator_cache[name]

macro _import_atomic_symbol_creators()
  h_creators = _get_atomic_symbol_creators()

  exprs = Expr[]
  for creator_hdr in h_creators
    expr = _define_atomic_symbol_creator(creator_hdr)
    push!(exprs, expr)
  end

  esc(quote
    $(exprs...)
  end)
end

@_import_atomic_symbol_creators()

################################################################################
# Utility macros to chain up symbols
################################################################################
macro chain(layers)
  exprs = []
  last_layer = nothing
  function _chain_layer(layer, last_layer)
    if isa(last_layer, Void)
      esc(layer)
    else
      @assert(isa(layer, Expr) && layer.head == :call, "Do not know how to chain up $layer")
      return Expr(:call, esc(layer.args[1]), last_layer, map(esc, layer.args[2:end])...)
    end
  end
  while true
    if layers.head == :(=>)
      new_layer = gensym()
      push!(exprs, :($new_layer = $(_chain_layer(layers.args[1], last_layer))))
      last_layer = new_layer
      layers = layers.args[2]
    else
      push!(exprs, _chain_layer(layers, last_layer))
      break
    end
  end
  return Expr(:block, exprs...)
end


text:
using GLM

tests = ["glmFit"]

println("Running tests:")

for t in tests
    tfile = joinpath("test", "$(t).jl")
    println(" * $(tfile) ...")
    include(tfile)
end


text:
using LinearAlgebra
using PROPACK
using Test
using Random

Random.seed!(1234)

for T ∈ (Float32, Float64, ComplexF32, ComplexF64)
  println("--- Test functions with ", T, " types ---")
  ϵ = eps(real(one(T)))

  # construct a matrix with chosen singular values
  m, n = 10, 5  # choose n ≤ m
  σ = collect(n:-1:1)  # in decreasing order because Propack always returns them so
  Σ = [diagm(0 => σ); zeros(T, m - n, m - n)]
  U, _ = qr(rand(T, m, m))
  V, _ = qr(rand(T, n, n))
  A = U * Σ * V'

  # compute a few leading singular values and vectors with lansvd
  k = 3
  tolin = sqrt(ϵ)
  U, s, V, bnd, nprod, nctprod = tsvd(A, k = k, tolin = tolin)
  @assert size(U) == (m, k)
  @assert length(s) == k
  @assert size(V) == (n, k)
  @assert length(bnd) == k
  println("computed singular values: ", s)
  println("exact singular values: ", σ[1:k])
  println("error bounds reported: ", bnd)
  println("number of matvecs and conjugate transposed matvecs: ", nprod, ", ", nctprod)
  for i = 1:k
    @assert bnd[i] ≤ max(16ϵ * s[1], tolin * s[i])
    @assert abs(s[i] - σ[i]) ≤ tolin * σ[i]
  end
  @assert norm(U * diagm(0 => s) * V' - A) ≤ 1.2 * σ[k + 1]

  # ensure tsvdvals returns the same singular value estimates
  s2, _, nprod, nctprod = tsvdvals(A, k = k, tolin = tolin)
  println("computed singular values: ", s2)
  println("number of matvecs and conjugate transposed matvecs: ", nprod, ", ", nctprod)
  @assert norm(s - s2) ≤ sqrt(ϵ) * norm(s)

  # compute the two smallest singular values with tsvd_irl
  # we use A' to skip zero singular values
  k = 2
  σs = σ[(n - k + 1):n]
  U, s, V, bnd, nprod, nctprod = tsvd_irl(A', k = k, tolin = tolin)
  @assert size(U) == (n, k)
  @assert length(s) == k
  @assert size(V) == (m, k)
  @assert length(bnd) == k
  println("computed smallest singular values: ", s)
  println("exact smallest singular values: ", σs)
  println("error bounds reported: ", bnd)
  println("number of matvecs and conjugate transposed matvecs: ", nprod, ", ", nctprod)
  for i = 1:k
    @assert bnd[i] ≤ max(16ϵ * s[1], tolin * s[i])
    @assert abs(s[i] - σs[i]) ≤ tolin * σs[i]
  end

  # ensure tsvdvals_irl returns the same singular value estimates
  s2, _, nprod, nctprod = tsvdvals_irl(A', k = k, tolin = tolin)
  println("computed smallest singular values: ", s2)
  println("number of matvecs and conjugate transposed matvecs: ", nprod, ", ", nctprod)
  @assert norm(s - s2) ≤ sqrt(ϵ) * norm(s)
  println()
end


text:
using MatrixMarket
using DataFrames
using CSV

include("./utils.jl")
include("../src/DirectMethod.jl")

const data = DataFrame(; N=Int64[], Time=Float64[], Error=Float64[], Space=Int64[])

foreach(readdir(DATA_DIR)) do filename
    print("Loading $filename... \n")
    A, t, m = @timed MatrixMarket.mmread(joinpath(DATA_DIR, filename))
    print("Loaded $filename [$t seconds, $m bytes] \n")

    T = solvematrix(A)
    push!(data, T)

    print("$T \n\n")
end

CSV.write(joinpath(OUTPUT_DIR, "$(getos()).csv"), data)

max_stars_repo_name: saiteki-kai/direct-methods

max_issues_repo_name: saiteki-kai/direct-methods

max_forks_repo_name: saiteki-kai/direct-methods


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

function delta_y_correction!(mapping, G, H)
    gsub = edgesubgraph(G, Set(keys(mapping)))
    hsub = edgesubgraph(H, Set(values(mapping)))
    revmap = Dict(v => k for (k, v) in mapping)
    for e in delta_y_edges(mapping, gsub, hsub)
        delete!(mapping, e)
    end
    for e in delta_y_edges(revmap, hsub, gsub)
        delete!(mapping, revmap[e])
    end
end

function delta_y_edges(mapping, gsub, hsub)
    dys = Int[]
    for gn in triangles(gsub)
        g_edges = edgeset(nodesubgraph(gsub, Set(gn)))
        h_edges = Set([mapping[e] for e in g_edges])
        if nodecount(edgesubgraph(hsub, h_edges)) != 3
            push!(dys, pop!(g_edges))
        end
    end
    return dys
end

function delta_y_mismatch(G, H, mapping)
    # Delta-Y check for edge induced subgraph isomorphism
    gsub = edgesubgraph(G, Set(keys(mapping)))
    hsub = edgesubgraph(H, Set(values(mapping)))
    revmap = Dict(v => k for (k, v) in mapping)
    return delta_y(gsub, hsub, mapping) || delta_y(hsub, gsub, revmap)
end

function delta_y(gsub, hsub, mapping)
    for gn in triangles(gsub)
        g_edges = edgeset(nodesubgraph(gsub, Set(gn)))
        h_edges = Set([mapping[e] for e in g_edges])
        if nodecount(edgesubgraph(hsub, h_edges)) != 3
            return true
        end
    end
    return false
end

function lgnodematcher(G::LineGraph, H::LineGraph,
                       nodematcher::Function, edgematcher::Function)
    return function (g, h)
        edgematcher(g, h) || return false
        ge = nodeattr(G, g)
        he = nodeattr(H, h)
        m1 = nodematcher(ge.n1, he.n1) && nodematcher(ge.n2, he.n2)
        m2 = nodematcher(ge.n1, he.n2) && nodematcher(ge.n2, he.n1)
        return m1 || m2
    end
end

function lgedgematcher(G::LineGraph, H::LineGraph, nodematcher::Function)
    return (g, h) -> nodematcher(edgeattr(G, g).node, edgeattr(H, h).node)
end


text:
@testset ExtendedTestSet "classical_chan_vese" begin
    @testset ExtendedTestSet "classical_chan_vese" begin
        ellipse = [
            -1024 -1024 -1024 200 200 200 200 -1024 -1024 -1024
            -1024 200 200 200 200 200 200 200 200 -1024
            -1024 200 200 200 200 200 200 200 200 -1024
            200 200 200 200 200 200 200 200 200 200
            200 200 200 200 200 200 200 200 200 200
            200 200 200 200 200 200 200 200 200 200
            200 200 200 200 200 200 200 200 200 200
            -1024 200 200 200 200 200 200 200 200 -1024
            -1024 200 200 200 200 200 200 200 200 -1024
            -1024 -1024 -1024 200 200 200 200 -1024 -1024 -1024
        ]
        answer = [
            0 0 0 1 1 1 1 0 0 0
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            1 1 1 1 1 1 1 1 1 1
            1 1 1 1 1 1 1 1 1 1
            1 1 1 1 1 1 1 1 1 1
            1 1 1 1 1 1 1 1 1 1
            0 1 1 1 1 1 1 1 1 0
            0 1 1 1 1 1 1 1 1 0
            0 0 0 1 1 1 1 0 0 0
        ]
        phi, iter = ChanVese.classical_chan_vese(ellipse; max_iter=500)
        seg = phi .> 0
        @test seg == answer
    end
end


text:
using RCall
using Distributed
addprocs(2)
@everywhere include("ProblemGenerator.jl")

#Need a function that creates distributions of heterogenous states
@everywhere AddNoise2States(σ²) =
  [truncated(Normal(μ, σ²), 0, Inf) for μ in nonZeroSpeciesValues]

#Create a function that calculate the peak IFN distribution given:
 #The number of cells initially infected
 #The amount of noise desired for the initial protien concentrations
@everywhere function HeteroSimulation(primeCellPercent, varIC, prob)
    #Make a copy of the initial conditions
  u0 = copy(prob.u0)
    #Get number of primary cells wanted from percent
  primeCellNum = round(Int, nCells * primeCellPercent)
  secondaryCellNum = nCells - primeCellNum
    #Randomly place them in the ABM
  cellLocations = shuffle([ones(primeCellNum); zeros(secondaryCellNum)])
    #Update the number of initially infected cells (DNA conc.)
  u0[:, :, 2] = m2c.(1e3 .* reshape(cellLocations, N, N))

    @show varIC
    #Add noise to the nonzero initial conditions
  noiseDistributions = AddNoise2States(varIC)
  for (i, index) in enumerate(nonZeroSpeciesIdx)
    u0[:, :, index] = rand(noiseDistributions[i], N, N)
  end

    #Make a copy of the parameter container
  θ = deepcopy(prob.p)
    #Reset the number of cells infected initially
  θ.cellsInfected = fill(Inf, N, N)
  θ.cellsInfected[findall(u0[:, :, 2] .> 0.0), 1] .= 0.0
    #Reset the mass balances
  θ.mass = [u0[:, :, i] for i in nonZeroSpeciesIdx]

    #Define the new problem and solve
  probHetero = remake(prob; u0 = u0, p = θ)
  solHetero = solve(probHetero, CVODE_BDF(linear_solver = :GMRES))

    #Seperate by primary and secondary cells
  isPrimary = u0[:, :, 2] .> 0.0
  maxIFN = [maximum(solHetero[coord, 7, :]) for coord in cellIndicies]

  return [maxIFN[isPrimary], maxIFN[.!isPrimary]]
end

#Run the simulation as either Stochastic/Deterministic and Hetero
@everywhere probDet = ModelSetup(:ISD, :notStochastic, :Hetero)
@everywhere probStoch = ModelSetup(:ISD, :Stochastic, :Hetero)

varICRange = 0.0:0.1:1.0
varIClength = length(varICRange)
primeCellPercent = 0.63
primeCellNum = round(Int,nCells * primeCellPercent)

maxIFNDetHet = pmap(x -> HeteroSimulation(primeCellPercent, x, probDet), varICRange)
maxIFNStochHet = pmap(x -> HeteroSimulation(primeCellPercent, x, probStoch), varICRange)

function DataFrameMaker(sols)
  #Seperate out the primary and secondary cells
  primary = [sols[i][1] for i=1:varIClength]
  secondary = [sols[i][2] for i=1:varIClength]

  #Convert into a dataframe with 3 columns (percent,IFN, and celltype (primary or secondary))
  primary = DataFrame(percent=repeat(1:varIClength,inner=primeCellNum),
                      IFN = vcat(primary...),
                      CellState = repeat(["Primary"],length(vcat(primary...))) )
  secondary = DataFrame(percent=repeat(1:varIClength,inner= nCells-primeCellNum),
                        IFN = vcat(secondary...),
                        CellState = repeat(["Secondary"],length(vcat(secondary...))) )
  #Combine the two datadrames and return th result
  return vcat(primary,secondary)
end

ISDDetHet = DataFrameMaker(maxIFNDetHet)
ISDStochHet = DataFrameMaker(maxIFNStochHet)

@rput ISDDetHet
@rput ISDStochHet
R"""
library(ggplot2)
library(ggpubr)
GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin,
                           draw_group = function(self, data, ..., draw_quantiles = NULL) {
                             data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
                             grp <- data[1, "group"]
                             newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
                             newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
                             newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])

                             if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
                               stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
                                                                         1))
                               quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
                               aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
                               aesthetics$alpha <- rep(1, nrow(quantiles))
                               both <- cbind(quantiles, aesthetics)
                               quantile_grob <- GeomPath$draw_panel(both, ...)
                               ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
                             }
                             else {
                               ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
                             }
                           })

geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ...,
                              draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE,
                              show.legend = NA, inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin,
        position = position, show.legend = show.legend, inherit.aes = inherit.aes,
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}

commonFigureOptions <- list(geom_split_violin(scale="width",trim=FALSE),
  theme_pubr(border=TRUE),
  xlab("Initial Condition Variance"),
  ylab("IFN (nM)"),
  scale_x_discrete(labels= unlist(lapply(seq(0.0,1.0,0.1),toString))),
  stat_summary(fun.y = median, fun.ymin = median, fun.ymax = median, geom = "crossbar", width = 0.25, position = position_dodge(width = .25)),
  theme(plot.title = element_text(hjust = 0.5),aspect.ratio = 0.5)
  )

p1 <- ggplot(ISDDetHet, aes(factor(percent), IFN, fill = CellState)) +
  ggtitle("ISD: Deterministic (63% of Cell Infected)") +
  commonFigureOptions

p2 <- ggplot(ISDStochHet, aes(factor(percent), IFN, fill = CellState)) +
  ggtitle("ISD: Stochastic (63% of Cell Infected)") +
  commonFigureOptions

figure <- ggarrange(p1,p2,
                    labels = "AUTO",
                    common.legend = TRUE, legend = "right",
                    align = "h",
                    ncol = 1, nrow = 2)

ggsave("./Figures/Figure5new.pdf")
"""


text:
## Aetherlang interpreter exceptions and errors

# NOTE THAT AetherlangException IS NOT A SUBTYPE OF Exception
abstract type AetherlangException end

struct AetherlangError <: AetherlangException
    message::String
end

function aetherlang_print_exception(e, line::Vector{String}, dimname::String)
    print("\n\n"*(typeof(e) <: AetherlangException ? e.message : repr(e))*"\nAETHERLANG ERROR AT LINE:\n\t"*join(line, ' ')*"\n(COMMENTS REMOVED)\nIN DIMENSION "*dimname*'\n')
    exit()
end

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
module CameraGeometryTest

# write your own tests here
using FactCheck, CameraGeometry, Base.Test, Images

isinteractive() || FactCheck.exitstatus()

end


text:
# Function to create a Count Time Series Model specification object
# Input similar to Mspec, but allows for different inputs

"""
    Model(model, distr, link, pastObs, pastMean, X,
          external, zi)
Wrapper function to define count data models.
Default setting is an an IID Poisson process without regressors or zero inflation.

Structs have entries:

* `distr`: "Poisson" or "Negative Binomial" (Vector for INARMA)
* `link`: Vector of length two, "Linear" or "Log" (Vector for INARMA)
* `pastObs`: Lags considered in autoregressive part
* `pastMean`: Lags considered in MA/past conditional mean part
* `X`: Matrix of regressors (row-wise)
* `external`: Indicator(s) if regressors enter externally
* `zi`: Indicator, zero inflation Y/N

# Examples
```julia-repl
Model(pastObs = 1:2, pastMean = 1) # INGARCH(2, 1)
Model(pastObs = [1, 2], distr = "NegativeBinomial") # NB-INARCH(2)
Model(model = "INARMA", pastMean = 1, zi = true) # Zero inflated INMA(1)
```

For further details, see [Documentation](https://github.com/ManuelStapper/CountTimeSeries.jl/blob/master/CountTimeSeries_documentation.pdf)
"""
function Model(;model = "INGARCH",
    distr = "Poisson",
    link = "Linear",
    pastObs = Array{Int64, 1}([]),
    pastMean = Array{Int64, 1}([]),
    X = Array{Float64, 2}(undef, 0, 0),
    external = Array{Bool, 1}([]),
    zi = false)

    if !(model in ["INARMA", "INGARCH"])
        error("Invalid model type.")
    end

    if typeof(distr) == String
        if !(distr in ["Poisson", "NegativeBinomial"])
            error("Invalid distribution.")
        else
            distr = fill(distr, 2)
        end
    elseif typeof(distr) == Array{String, 1}
        if length(distr) == 0
            error("Invalid distribution.")
        end
        if length(distr) == 1
            distr = fill(distr[1], 2)
        end
        if length(distr) > 2
            distr = distr[1:2]
        end

        if !(distr[1] in ["Poisson", "NegativeBinomial"])
            error("Invalid (first) distribution.")
        end

        if !(distr[2] in ["Poisson", "NegativeBinomial"])
            error("Invalid (second) distribution.")
        end
    else
        error("Invalid distribution")
    end

    if typeof(link) == String
        if !(link in ["Linear", "Log"])
            error("Invalid link.")
        else
            link = fill(link, 2)
        end
    elseif typeof(link) == Array{String, 1}
        if length(link) == 0
            error("Invalid link.")
        end
        if length(link) == 1
            link = fill(link[1], 2)
        end
        if length(link) > 2
            link = link[1:2]
        end

        if !(link[1] in ["Linear", "Log"])
            error("Invalid (first) link.")
        end

        if !(link[2] in ["Linear", "Log"])
            error("Invalid (second) link.")
        end
    else
        error("Invalid link")
    end

    if typeof(pastObs) == UnitRange{Int64}
        pastObs = collect(pastObs)
    end

    if length(pastObs) == 0
        pastObs = Array{Int64, 1}([])
    end

    if typeof(pastObs) <: Integer
        pastObs = [pastObs]
    end

    if typeof(pastObs) <: Matrix{T} where T<:Integer
        if size(pastObs)[1] == 1
            pastObs = pastObs[1, :]
        elseif size(pastObs)[2] == 1
            pastObs = pastObs[:, 1]
        else
            error("Invalid specification of pastObs.")
        end
    end

    if !(typeof(pastObs) <: Array{T, 1} where T<:Integer)
        error("Invalid specification of pastObs.")
    else
        if any(pastObs .<= 0)
            error("Entries of pastObs must be positive integers.")
        else
            pastObs = sort(unique(pastObs))
        end
    end

    if typeof(pastMean) == UnitRange{Int64}
        pastMean = collect(pastMean)
    end

    if length(pastMean) == 0
        pastMean = Array{Int64, 1}([])
    end

    if typeof(pastMean) <: Integer
        pastMean = [pastMean]
    end

    if typeof(pastMean) <: Matrix{T} where T<:Integer
        if size(pastMean)[1] == 1
            pastMean = pastMean[1, :]
        elseif size(pastMean)[2] == 1
            pastMean = pastMean[:, 1]
        else
            error("Invalid specification of pastMean.")
        end
    end

    if typeof(pastMean) != Array{Int64, 1}
        error("Invalid specification of pastMean.")
    else
        if any(pastMean .<= 0)
            error("Entries of pastMean must be positive integers.")
        else
            pastMean = sort(unique(pastMean))
        end
    end

    if length(external) == 0
        external = Array{Bool, 1}([])
    end

    if typeof(external) == Bool
        external = [external]
    end

    if typeof(external) == Array{Bool, 2}
        if size(external)[1] == 1
            external = external[1, :]
        elseif size(external)[2] == 1
            external = external[:, 1]
        else
            error("Invalid specification of external.")
        end
    end

    if typeof(external) != Array{Bool, 1}
        error("Invalid specification of external.")
    end

    r = length(external)

    X = X .+ 0.0

    if typeof(X) == UnitRange{Int64}
        X = reshape(collect(X), (1, length(X)))
    end

    if length(X) == 0
        X = zeros(0, 0)
        if r > 0
            error("Invalid specification of external and X.")
        end
    else
        if typeof(X) <: Vector{T} where T<:AbstractFloat
            if r > 1
                error("Invalid specification of external and X.")
            end
            if r == 0
                external = [true]
                r = 1
            end
            X = reshape(X, (1, length(X)))
        end

        if typeof(X) <: Matrix{T} where T<:AbstractFloat
            if r > 0
                if size(X)[1] != r
                    error("Invalid dimensions in X.")
                end
            else
                external = fill(true, size(X)[1])
                r = size(X)[1]
            end

        else
            error("Invalid specification of X.")
        end
    end

    if r > 0
        iI = findall(.!external)
        rI = length(iI)
        iE = findall(external)
        rE = length(iE)

        if model == "INGARCH"
            if link[1] == "Linear"
                if any(X .< 0)
                    error("Negative regressors invalid for linear link.")
                end
            end
        end

        if model == "INARMA"
            if link[1] == "Linear"
                if rI > 0
                    if any(X[iI, :] .< 0)
                        error("Negative regressors invalid for linear link.")
                    end
                end
            end

            if link[2] == "Linear"
                if rE > 0
                    if any(X[iE, :] .< 0)
                        error("Negative regressors invalid for linear link.")
                    end
                end
            end
        end
    end

    if sum(distr .== "NegativeBinomial") == 2
        if r == 0
            distr[2] = "Poisson"
        else
            if rE == 0
                distr[2] = "poisson"
            end
        end
    end

    if length(zi) != 1
        error("Invalid specification of zi.")
    end

    if typeof(zi) in [Array{Bool, 1}, Array{Bool, 2}]
        zi = zi[1]
    end
    if typeof(zi) != Bool
        error("Invalid specification of zi.")
    end

    if model == "INGARCH"
        if length(pastObs) + length(pastMean) == 0
            return IIDModel(distr[1], link[1], X, external, zi)
        end

        if length(pastMean) == 0
            return INARCHModel(distr[1], link[1], pastObs, X, external, zi)
        end

        if (length(pastObs) == 0) & (length(pastMean) > 0)
            error("INGARCH(0, q) not supported.")
        end
        return INGARCHModel(distr[1], link[1], pastObs, pastMean, X, external, zi)
    elseif model == "INARMA"
        if (length(pastMean) == 0) & (length(pastObs) > 0)
            return INARModel(distr, link, pastObs, X, external, zi)
        end

        if (length(pastMean) >= 0) & (length(pastObs) == 0)
            return INMAModel(distr, link, pastMean, X, external, zi)
        end

        return INARMAModel(distr, link, pastObs, pastMean, X, external, zi)
    else
        error("Invalid model framework chosen.")
    end
end


text:
export @parallel_repeat, @parallel_alternate

function pfor_ntimes(niters::Int, f, N::Int)
    for iter in 1:niters
        for c in Base.splitrange(N, nworkers())
            @spawn f(first(c), last(c))
        end
    end
    nothing
end

function pfor_ntimes(niters::Int, f1, N1::Int, f2, N2::Int)
    for iter in 1:niters
        for c in Base.splitrange(N1, nworkers())
            @spawn f1(first(c), last(c))
        end
        for c in Base.splitrange(N2, nworkers())
            @spawn f2(first(c), last(c))
        end
    end
    nothing
end

macro parallel_repeat(niters, loop::Expr)
    if !is(loop.head,:for)
        error("malformed @parallel_repeat loop")
    end
    var = loop.args[1].args[1]
    r = loop.args[1].args[2]
    body = loop.args[2]
    quote
        pfor_ntimes($(esc(niters)),
                    $(Base.make_pfor_body(var, body, r)), length($(esc(r))))
    end
end

macro parallel_alternate(niters, loop1::Expr, loop2::Expr)
    if !all([is(loop.head,:for) for loop in [loop1, loop2]])
        error("malformed @parallel_alternate loop")
    end
    var1 = loop1.args[1].args[1]
    r1 = loop1.args[1].args[2]
    body1 = loop1.args[2]
    var2 = loop2.args[1].args[1]
    r2 = loop2.args[1].args[2]
    body2 = loop2.args[2]
    quote
        pfor_ntimes($(esc(niters)),
                    $(Base.make_pfor_body(var1, body1, r1)), length($(esc(r1))),
                    $(Base.make_pfor_body(var2, body2, r2)), length($(esc(r2))))
    end
end


text:
using Espresso
using XGrad
using GradDescent
using Distributions
using NNlib
using Requires

# otherwise CuArray type isn't resolved during loading
@require CuArrays using CuArrays

include("utils.jl")
include("diffrules.jl")
include("inplacerules.jl")
include("losses.jl")
include("modelopt.jl")
include("mlutils.jl")


text:
using URIParser

struct Address
    scheme::String
    host::String
    port::Integer
end

function Base.show(io::IO, address::Address)
    print(io,  "Address(scheme: $(address.scheme), host: $(address.host) port: $(address.port))")
end

function from_uri(uri::String, default_port::Integer=0)
    parsed = URI(uri)
    return Address(parsed.scheme, parsed.host, if parsed.port == 0 default_port else parsed.port end)
end

Base.:(==)(x::Address, y::Address) = x.scheme == y.scheme && x.host == y.host && x.port == y.port


text:
@testset "TensorGP.jl" begin

    tt = MPT([rand(1,4,2),rand(2,5,3),rand(3,4,1)]);
    @test order(tt) == 3;

end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
using LazySets: default_cddlib_backend

@test default_cddlib_backend(Float64) == CDDLib.Library(:float)
@test default_cddlib_backend(Float32) == CDDLib.Library(:float)
@test default_cddlib_backend(Rational{Int}) == CDDLib.Library(:exact)


text:
module mfeit1
using StaticArrays
using Elfel
using Elfel.RefShapes: RefShapeTriangle, manifdim, RefShapeInterval
using Elfel.FElements: FE, refshape, FEH1_T3
using Elfel.FElements: bfun, bfungradpar
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!
using Elfel.FEIterators: FEIterator

using MeshCore
using MeshSteward: Mesh, load, nspacedims, baseincrel
using Test
function test()
    fe = FEH1_T3()
    mesh = load(Mesh(), "qmesh.mesh")
    fesp = FESpace(Float64, mesh, fe)

    @test doftype(fesp) == Float64

    refc = [[1, 2, 5],
    [1, 5, 4],
    [4, 5, 8],
    [4, 8, 7],
    [7, 8, 11],
    [7, 11, 10],
    [2, 3, 6],
    [2, 6, 5],
    [5, 6, 9],
    [5, 9, 8],
    [8, 9, 12],
    [8, 12, 11], ]
    it = FEIterator(fesp)
    k = 1
    for c in it
        @test isapprox(c._nodes,  refc[k])
        k = k + 1
    end
    @test length(it) == 12

    # @show summary(mesh)

    numberfreedofs!(fesp)
    numberdatadofs!(fesp)
    @test ndofs(fesp) == 12
    @test nunknowns(fesp) == 12

    for i in [1, 4, 7, 10]
        setebc!(fesp, 0, i, 1, 0.0)
    end
    numberfreedofs!(fesp)
    numberdatadofs!(fesp)
    @test ndofs(fesp) == 12
    @test nunknowns(fesp) == 8

    # @show fesp._irsfields

    refd =  [[9], [1], [2], [10], [3], [4], [11], [5], [6], [12], [7], [8]]
    it = FEIterator(fesp)
    for el in it
         # @show el._dofs
         # @show refd[el._nodes]
         @test isapprox(el._dofs, [refd[n][1] for n in el._nodes] )
    end
end
end
using .mfeit1
mfeit1.test()

module mfeit2
using StaticArrays
using LinearAlgebra
using MeshCore
using MeshSteward: Mesh, load, nspacedims, baseincrel
using Elfel
using Elfel.RefShapes: RefShapeTriangle, manifdim, RefShapeInterval
using Elfel.FElements: FE, refshape, FEH1_T3
using Elfel.FElements: bfun, bfungradpar
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!
using Elfel.FEIterators: FEIterator, eldofs
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
using Test
A = [0.6744582963441466 0.2853043149927861 0.27460710155821255;
0.3781923479141225 0.2838873430062512 0.6316949656630075;
0.19369805365903336 0.8926164783344779 0.07006905962860177]
function test()
    fe = FEH1_T3()
    mesh = load(Mesh(), "qmesh.mesh")
    fesp = FESpace(Float64, mesh, fe)

    for i in [1, 4, 7, 10]
        setebc!(fesp, 0, i, 1, 0.0)
    end
    numberfreedofs!(fesp)
    numberdatadofs!(fesp)

    ass = SysmatAssemblerSparse(0.0)
    start!(ass, 12, 12)
    ke = LocalMatrixAssembler(size(A, 1), size(A, 2), 0.0)
    it = FEIterator(fesp)
    for el in it
        init!(ke, eldofs(el), eldofs(el))
        for j in 1:size(A, 2), i in 1:size(A, 1)
            ke[i, j] += A[i, j]
        end
        assemble!(ass, ke)
    end
    S = finish!(ass)

    D = fill(0.0, 12, 12)
    it = FEIterator(fesp)
    for el in it
        for j in 1:size(A, 2), i in 1:size(A, 1)
            D[el._dofs[i], el._dofs[j]] += A[i, j]
        end
    end
    @test isapprox(D - S, 0 * I)
end
end
using .mfeit2
mfeit2.test()

module mfeit3
using StaticArrays
using LinearAlgebra
using MeshCore
using MeshSteward: Mesh, load, nspacedims, baseincrel
using Elfel
using Elfel.RefShapes: RefShapeTriangle, manifdim, RefShapeInterval
using Elfel.FElements: FE, refshape, FEH1_T3
using Elfel.FElements: bfun, bfungradpar
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!
using Elfel.FEIterators: FEIterator, eldofs, eldofs, eldofentmdims, eldofcomps
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Test
A = [0.6744582963441466 0.2853043149927861 0.27460710155821255;
0.3781923479141225 0.2838873430062512 0.6316949656630075;
0.19369805365903336 0.8926164783344779 0.07006905962860177]
function test()
    fe = FEH1_T3()
    mesh = load(Mesh(), "qmesh.mesh")
    fesp = FESpace(Float64, mesh, fe, 3)

    for i in [1, 4, 7, 10]
        setebc!(fesp, 0, i, 1, 0.0)
    end
    numberfreedofs!(fesp)
        numberdatadofs!(fesp)

    it = FEIterator(fesp)
    @test em = eldofentmdims(it) == [0, 0, 0, 0, 0, 0, 0, 0, 0]
    @test eldofcomps(it) == [1, 2, 3, 1, 2, 3, 1, 2, 3]
    ref = [ [33, 1, 2, 3, 4, 5, 11, 12, 13],
     [33, 1, 2, 11, 12, 13, 34, 9, 10],
     [34, 9, 10, 11, 12, 13, 19, 20, 21],
     [34, 9, 10, 19, 20, 21, 35, 17, 18],
     [35, 17, 18, 19, 20, 21, 27, 28, 29],
     [35, 17, 18, 27, 28, 29, 36, 25, 26],
     [3, 4, 5, 6, 7, 8, 14, 15, 16],
     [3, 4, 5, 14, 15, 16, 11, 12, 13],
     [11, 12, 13, 14, 15, 16, 22, 23, 24],
     [11, 12, 13, 22, 23, 24, 19, 20, 21],
     [19, 20, 21, 22, 23, 24, 30, 31, 32],
     [19, 20, 21, 30, 31, 32, 27, 28, 29]  ]
    i = 1
    for el in it
        @test isapprox(eldofs(it), ref[i])
        i = i + 1
    end

end
end
using .mfeit3
mfeit3.test()

module maggfeit1
using StaticArrays
using Elfel
using Elfel.RefShapes: RefShapeTriangle, manifdim, RefShapeInterval
using Elfel.FElements: FE, refshape, FEH1_T3
using Elfel.FElements: bfun, bfungradpar
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FEIterators: FEIterator, elnodes
using Elfel.Utilities: @_check

using MeshCore
using MeshSteward: Mesh, load, nspacedims, baseincrel
using Test
function test()
    mesh = load(Mesh(), "qmesh.mesh")
    fesp1 = FESpace(Float64, mesh, FEH1_T3())
    fesp2 = FESpace(Float64, mesh, FEH1_T3())
    fesp3 = FESpace(Float64, mesh, FEH1_T3())

    refc = [[1, 2, 5],
    [1, 5, 4],
    [4, 5, 8],
    [4, 8, 7],
    [7, 8, 11],
    [7, 11, 10],
    [2, 3, 6],
    [2, 6, 5],
    [5, 6, 9],
    [5, 9, 8],
    [8, 9, 12],
    [8, 12, 11], ]
    it1 = FEIterator(fesp1)
    it2 = FEIterator(fesp2)
    it3 = FEIterator(fesp3)
    it = zip(it1, it2, it3)
    @test length(it) == 12

    k = 1
    for i in it
        @test elnodes(i[1]) == elnodes(i[2]) == elnodes(i[3])
        k = k + 1
    end

end
end
using .maggfeit1
maggfeit1.test()


text:
ESC = "\033[" #0x1B
RESET = "\033[m"

RED = "31"
BRIGHTRED = "31;1"
GREEN = "32"
BRIGHTGREEN = "32;1"

#println(ESC*"3;42;31mTexting"*RESET)

println("#")
println(ESC*GREEN*"mTexting"*RESET)
println(ESC*RED*"mTexting"*RESET)
println("@")

max_stars_repo_name: aole/RogueJ

max_issues_repo_name: aole/RogueJ

max_forks_repo_name: aole/RogueJ


text:
using Libdl

if VERSION >= v"1.3" && !Sys.iswindows() && !haskey(ENV, "SCIPOPTDIR")
    # Skip build in favor of SCIP_jll
    exit()
end

depsfile = joinpath(dirname(@__FILE__), "deps.jl")
if isfile(depsfile)
    rm(depsfile)
end

function write_depsfile(path)
    open(depsfile, "w") do f
        print(f, "const libscip = ")
        show(f, path)
        println(f)
    end
end

libname = if Sys.islinux()
    "libscip.so"
elseif Sys.isapple()
    "libscip.dylib"
elseif Sys.iswindows()
    "libscip.dll"
else
    error("SCIP is currently not supported on \"$(Sys.KERNEL)\"")
end

paths_to_try = []

# prefer environment variable
if haskey(ENV, "SCIPOPTDIR")
    push!(paths_to_try, joinpath(ENV["SCIPOPTDIR"], "bin", libname))
    push!(paths_to_try, joinpath(ENV["SCIPOPTDIR"], "lib", libname))
end

# but also try library path
push!(paths_to_try, libname)

found = false
tried = String[]
for l in paths_to_try
    try
        d = Libdl.dlopen(l)
        global found = true
        write_depsfile(l)
        break
    catch e
        push!(tried, "$(l): $(e.msg)")
    end
end

if !found && !haskey(ENV, "SCIP_JL_SKIP_LIB_CHECK")
    error("""
Unable to locate SCIP installation. Tried:

$(join(tried, "\n\n"))

Note that this must be downloaded separately from scip.zib.de.
Please set the environment variable SCIPOPTDIR to SCIP's installation path.
""")
end


text:
module Slurm

# package code goes here
include("api.jl")
include("error.jl")
include("hostlist.jl")

end # module


text:
using TextMining
using FactCheck
using Compat

reload("feature_vector_tests.jl")
reload("cluster_tests.jl")
reload("data_set_tests.jl")
reload("clustering_tests.jl")
reload("distribution_tests.jl")
reload("naive_bayes_tests.jl")
FactCheck.exitstatus()


text:
# Rocket
using NLOptControl, PrettyPlots
n = define(numStates=3,numControls=1,X0=[0.0,0.0,1.0],XF=[10.0,0.0,NaN],XL=[NaN,-0.1,0.0],XU=[NaN,1.7,1.0],CL=[-1.1],CU=[1.1]);
n.s.ocp.tfMax = 15.; # tFmin = 5.
states!(n,[:h,:v,:m],descriptions=["height (t)","velocity (t)","mass (t)"]);
controls!(n,[:T],descriptions=["thrust (t)"]);
dx=Array{Expr}(3,);
dx[1] = :(v[j])
dx[2] = :((T[j]- 0.2*v[j]^2)/m[j])
dx[3] = :(-0.01*T[j]^2)
dynamics!(n,dx)
configure!(n;(:solverSettings=>(:name=>:KNITRO)),(:finalTimeDV=>true))
@NLobjective(n.ocp.mdl,Min, n.ocp.tf)
optimize!(n)
allPlots(n)


text:
using Markdown

WebIO.render(md::Markdown.MD) = blocknode(md)

blocknode(x) = wrapnode(x)
inlinenode(x) = wrapnode(x)

blocknode(md::Markdown.MD) = vbox(map(blocknode, md.content))(className="jl-markdown")
inlinenode(x::AbstractString) = x

blocknode(md::Markdown.Header{n}) where {n} =
     node("h$n", map(inlinenode, md.text)...)

inlinenode(md::Markdown.Code) = node("code", md.code)

blocknode(md::Markdown.Code)  =
    node("pre", md.code)

blocknode(md::Markdown.BlockQuote) =
    blockquote(map(blocknode, md.content))

function blocknode(md::Markdown.List)
    el = md.ordered == -1 ? "ul" : "ol"
    node(el, map(x->node("li", inlinenode.(x)...), md.items)...)
end
blocknode(md::Markdown.Paragraph) =
    node("p", map(inlinenode, md.content)...)

inlinenode(md::Markdown.Paragraph) =
    node("div", map(inlinenode, md.content)...)

inlinenode(md::Markdown.Italic) = node("em", map(inlinenode, md.text)...)
inlinenode(md::Markdown.Bold) = node("span", map(inlinenode, md.text)...)(fontweight("bold"))
inlinenode(md::Markdown.Link) = node("a", map(inlinenode, md.text)..., href=md.url)

inlinenode(md::Markdown.Image) = node("img", src=md.url, alt=md.alt)
blocknode(md::Markdown.Image) = node("img", src=md.url, alt=md.alt)

function katexscope(formula, display)
    imports = ["https://cdn.jsdelivr.net/npm/katex@0.10.1/dist/katex.min.js",
               "https://cdn.jsdelivr.net/npm/katex@0.10.1/dist/katex.min.css"]
    s = Scope(imports=imports)
    formula = strip(formula, '$')
    onimport(s, js"""function (katex) {
             katex.render($formula, this.dom, {displayMode: $display})
             }""")
    s
end

function inlinenode(md::Markdown.LaTeX)
    style(katexscope(md.formula, false)(node("span")), Dict("display"=>"inline"))
end
function blocknode(md::Markdown.LaTeX)
    katexscope(md.formula, true)(node("div"))
end


text:
module TerminalLoggers

using Logging:
    AbstractLogger,
    LogLevel, BelowMinLevel, Debug, Info, Warn, Error, AboveMaxLevel

import Logging:
    handle_message, shouldlog, min_enabled_level, catch_exceptions

export TerminalLogger

include("StickyMessages.jl")
include("TerminalLogger.jl")

end # module


text:
"""
fpatan is the derivative of atan
used in the demo above
"""
function fpatan(x)
    return 1 / (1 + x * x)
end


text:
"""
    TPVal{T,H,L}

A static parametrized counterpart to the `TwicePrecision` type. Any instance
of `TPVal` should be converted to `TwicePrecision` as soon as it is not serving
as a parametric storage type.
"""
struct TPVal{T,H,L} end

"""
    gethi(x::Union{TPVal{T}, TwicePrecision{T}}) -> T

Returns the `hi` component of a twice precision number. Works for both
statically set `TPVal` and `TwicePrecision`.
"""
@inline gethi(::T) where {T<:TPVal} = gethi(T)
Base.@pure gethi(::Type{<:TPVal{T,H,L}}) where {T,H,L} = H::T
gethi(x::TwicePrecision) = getfield(x, :hi)

"""
    getlo(x::Union{TPVal{T}, TwicePrecision{T}}) -> T

Returns the `lo` component of a twice precision number. Works for both
statically set `TPVal` and `TwicePrecision`.
"""
@inline getlo(::T) where {T<:TPVal} = getlo(T)
Base.@pure getlo(::Type{<:TPVal{T,H,L}}) where {T,H,L} = L::T
getlo(x::TwicePrecision) = getfield(x, :lo)

Base.eltype(::TPVal{T,H,L}) where {T,H,L} = T
Base.eltype(::Type{<:TPVal{T,H,L}}) where {T,H,L} = T

TPVal(x::TwicePrecision{T}) where {T} = TPVal{T}(x)
TPVal{T}(x::TwicePrecision) where {T} = TPVal{T,T(gethi(x)),T(getlo(x))}()
TPVal{T}(x::TwicePrecision{T}) where {T} = TPVal{T,gethi(x),getlo(x)}()

Base.TwicePrecision(x::TPVal) = TwicePrecision(gethi(x), getlo(x))

@inline (::Type{T})(x::TPVal) where {T<:Number} = T(gethi(x) + getlo(x))::T

function Base.convert(::Type{TwicePrecision{T}}, x::TPVal{T}) where {T}
    return TwicePrecision{T}(gethi(x), getlo(x))
end
function Base.convert(::Type{TwicePrecision{T}}, x::TPVal) where {T}
    return TwicePrecision{T}(convert(T, gethi(x)), convert(T, getlo(x)))
end
Base.convert(::Type{T}, x::TPVal) where {T<:Number} = T(x)

# hacky catch for converting TwicePrecision to TPVal for static parameters
tp2val(x::TwicePrecision{T}) where {T} = TPVal{T}(x)
tp2val(x) = x


text:
# Copyright (c) 2020 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export nlogML

"""
`nlogML(Nt,a,A,y)`

Computes the negative log Maximum Likelihood (ML) normalized by the sequence length (``\\log(P(\\mathbf{X}))/N_t``) of a Hidden Markov Model (HMM) with `Ns` states, transition matrix `A` (a `Ns` × `Ns` matrix), initial probabilities `a` (`Ns`-vector) and observation probabilities `y` (a `Ns`×`Nt` matrix) for `Nt` observations. All `Array`s must have the same element type.

`nlogML(Nt,t2tr,A,y)`

Returns the negative log Maximum Likelihood (ML) normalized by the sequence length (``\\log(P(\\mathbf{X}))/N_t``) with a constrained path (see [`forward`](@ref)).
"""
nlogML(Nt::Integer,c) = -sum(log.(c))/Nt
nlogML(Nt::Array,cs) = sum( [-sum(log.(c))/Nt[i] for (i,c) in enumerate(cs)])

function nlogML(Nt::Integer,a,A,y)
  alpha, c = forward(Nt,a,A,y)
  return nlogML(Nt,c)
end

function nlogML(Nt::Array,a,A,y)
  alphas, cs = forward(Nt,a,A,y)
  return nlogML(Nt,cs)
end

export nlogMLlog
"""
`nlogMLlog(Nt,A,a,y)`

Same as [`nlogML`](@ref) but expects `A`, `a` and `y` to be in the log-domain.
"""
function nlogMLlog(Nt,a,A,y)
  alpha, logML = logforward(Nt,a,A,y)
  return -sum(logML)
end


text:
"""
Probability functions for DPP and k-DPP.

Methods:
--------
    pmf(dpp::DeterminantalPointProcess, z::Array{Int})
    logpmf(dpp::DeterminantalPointProcess, z::Array{Int})
    pmf(dpp::DeterminantalPointProcess, z::Array{Int}, k::Int)
    logpmf(dpp::DeterminantalPointProcess, z::Array{Int}, k::Int)
"""

"""Compute the log probability of a sample `z` under the given DPP.
"""
function logpmf(dpp::DeterminantalPointProcess, z::Array{Int})
    L_z_eigvals = eigvals(dpp.L[z, z])
    return sum(log.(L_z_eigvals)) - sum(log.(dpp.Lfact.values .+ 1))
end

"""Compute the log probability of a sample `z` under the given k-DPP.
"""
function logpmf(dpp::DeterminantalPointProcess, z::Array{Int}, k::Int)
    L_z_eigvals = eigvals(dpp.L[z, z])
    return sum(log.(L_z_eigvals)) .- log(elem_symm_poly(dpp.Lfact.values, k)[end, end])
end

"""Compute the probability of a sample `z` under the given DPP.
"""
function pmf(dpp::DeterminantalPointProcess, z::Array{Int})
    exp(logpmf(dpp, z))
end

"""Compute the probability of a sample `z` under the given DPP.
"""
function pmf(dpp::DeterminantalPointProcess, z::Array{Int}, k::Int)
    exp(logpmf(dpp, z, k))
end


text:
metric = Euclidean()
dists = [metric(query, d) for d in data]
ss = searchstructure(BruteForce, data, metric)

skip3(i) = i % 3 == 0
skip3inv(i) = !skip3(i)
skip3bulk(i, j) = skip3(i + j)

@testset "Search structure attributes" begin
    @test datatype(ss) === datatype(typeof(ss)) === eltype(data)
    @test getmetric(ss) === ss.metric
end

@testset "NeighborNumber" begin
    t = NeighborNumber(k)

    # Using bruteforcesearch function
    idxs1, ds1 = results1 = bruteforcesearch(data, metric, query, t)

    @test idxs1 == sortperm(dists)[1:k]
    check_search_results(data, metric, results1, query, t)

    # Using BruteForce instance
    @test search_allfuncs(ss, query, t) == results1

    # Again with skip function
    idxs2, ds2 = results2 = bruteforcesearch(data, metric, query, t, skip3)

    @test idxs2 == filter(skip3inv, sortperm(dists))[1:k]
    check_search_results(data, metric, results2, query, t, skip3)

    @test search_allfuncs(ss, query, t, skip3) == results2
end

@testset "WithinRange" begin
    t = WithinRange(r)

    # Using bruteforcesearch function
    idxs1, ds1 = results1 = bruteforcesearch(data, metric, query, t)

    @test Set(ds1) == Set(filter(<=(r), dists))
    check_search_results(data, metric, results1, query, t)

    # Using BruteForce instance
    @test search_allfuncs(ss, query, t) == results1

    # Again with skip function
    idxs2, ds2 = results2 = bruteforcesearch(data, metric, query, t, skip3)

    @test idxs2 == filter(skip3inv, idxs1)
    check_search_results(data, metric, results2, query, t, skip3)

    @test search_allfuncs(ss, query, t, skip3) == results2
end

@testset "Bulk search" begin
    for t in [NeighborNumber(k), WithinRange(r)]
        for skip in [nothing, skip3bulk]
            test_bulksearch(ss, queries, t, skip)
        end
    end
end


text:
module ProgramTunnel_fs
using Formatting

export ProgramTunnelInfo, hello, recvText, sendText, reverseRole!

mutable struct ProgramTunnelInfo

    recv_fn      :: AbstractString
    send_fn      :: AbstractString
    done_recv_fn :: AbstractString
    done_send_fn :: AbstractString

    chk_freq     :: AbstractFloat
    timeout    :: AbstractFloat
    timeout_limit_cnt :: Integer
    buffer_cnt :: Integer

    recv_first_sleep_max :: AbstractFloat
    recv_first_sleep :: AbstractFloat
    recv_first_cnt   :: Integer

    rotate       :: Integer
    recv_trackno :: Integer
    send_trackno :: Integer

    paths        :: AbstractArray
    error_sleep  :: Float64

    history_length :: Integer

    function ProgramTunnelInfo(;
        recv          :: AbstractString     = "ProgramTunnel-Y2X.txt",
        send          :: AbstractString     = "ProgramTunnel-X2Y.txt",
        done_recv     :: AbstractString     = "ProgramTunnel-DONE-Y2X.txt",
        done_send     :: AbstractString     = "ProgramTunnel-DONE-X2Y.txt",
        chk_freq      :: AbstractFloat                  = 0.05,
        path          :: Union{AbstractString, Nothing} = nothing,
        timeout       :: AbstractFloat                  = 10.0,
        buffer        :: AbstractFloat                  = 0.1,
        recv_first_sleep_max :: AbstractFloat = 5.00,
        recv_first_sleep :: AbstractFloat = 0.0,
        reverseRole   :: Bool = false,
        rotate        :: Integer = 100,
        history_length:: Integer = 20,
    )

        if chk_freq <= 0.0
            ErrorException("chk_freq must be positive.") |> throw
        end

        paths = []
        for i = 1:rotate
            push!(paths, joinpath(MPI.path, format("{:03d}")))
        end

        PTI = new(
            recv,
            send,
            lock,
            chk_freq,
            timeout,
            ceil(timeout / chk_freq),
            ceil(buffer / chk_freq),
            recv_first_sleep_max,
            recv_first_sleep,
            ceil(recv_first_sleep / chk_freq),
            rotate, 1, 1, paths, 1.0, history_length,
        )

        if path != nothing
            appendPath(PTI, path)
        end

        if reverseRole
            reverseRole!(PTI)
        end

        makeDirs(PTI)

        return PTI
    end
end

function makeDirs(PTI::ProgramTunnelInfo)
    for path in PTI.paths
        makepath(path)
    end
end

function cleanDir(PTI::ProgramTunnelInfo, i::Integer)
    rm(PTI.paths[i], recursive=true, force=true)
    makepath(PTI.paths[i])
end

function cleanHistory(PTI::ProgramTunelInfo)
    cleanDir(PTI, ((PTI.send_trackno-1 - PTI.history_length -1) % PTI.rotate) + 1)
end

function appendPath(PTI::ProgramTunnelInfo, path::AbstractString)
    for i = 1:PTI.rotate
        PTI.paths[i] = joinpath(path, PTI.paths[i])
    end
end

function reverseRole!(PTI::ProgramTunnelInfo)
    PTI.recv_fn, PTI.send_fn = PTI.send_fn, PTI.recv_fn
    PTI.done_recv_fn, PTI.done_send_fn = PTI.done_send_fn, PTI.done_recv_fn
end

function sendText(PTI::ProgramTunnelInfo, msg::AbstractString)
#    local double_chk_msg

    if strip(msg) == ""
        throw(ErrorException("Cannot send empty msg."))
    end

    while true
        try
            open(joinpath(PTI.paths[PTI.send_trackno], PTI.send_fn), "w") do io
                write(io, msg)
            end
        catch ex
            println(string(ex))
            println("keep sending msg.")
            sleep(PTI.error_sleep)
        end
        break
    end

    while true
        try
            open(joinpath(PTI.paths[PTI.send_trackno], PTI.done_send_fn), "w") do io
                write(io, "DONE")
            end
        catch ex
            println(string(ex))
            println("keep sending done.")
            sleep(PTI.error_sleep)
        end
        break
    end

    PTI.send_trackno = (PTI.send_trackno % PTI.rotate) + 1

end

function recvText(PTI::ProgramTunnelInfo)

    recv_fn = joinpath(PTI.paths[PTI.recv_trackno], PTI.recv_fn)
    done_recv_fn = joinpath(PTI.paths[PTI.recv_trackno], PTI.done_recv_fn)

    get_through = false
    sleep(PTI.recv_first_sleep)

    if ! (isfile(recv_fn) && isfile(done_recv_fn))
        PTI.recv_first_sleep -= PTI.chk_freq
        PTI.recv_first_sleep = max(0.0, PTI.recv_first_sleep)
        get_through = true

        println("[recvText] Message is already there. Adjust recv_first_sleep to : ", PTI.recv_first_sleep)
    else
        for cnt in 1:(PTI.timeout_limit_cnt - PTI.recv_first_cnt)

            sleep(PTI.chk_freq)

            if isfile(recv_fn) && isfile(done_recv_fn)
                get_through = true

                if cnt <= PTI.buffer_cnt

                    println("[recvText] Good guess of the recv_first_sleep : ", PTI.recv_first_sleep)

                elseif PTI.recv_first_sleep < PTI.recv_first_sleep_max

                    # Out of buffer, need to adjust: increase PTI.recv_first_sleep
                    PTI.recv_first_sleep += PTI.chk_freq
                    PTI.recv_first_sleep = min(PTI.recv_first_sleep_max, PTI.chk_freq)
                    println("[recvText] Out of buffer. Adjust recv_first_sleep to : ", PTI.recv_first_sleep)

                else
                    println("[recvText] Out of buffer. But reach to recv_first_sleep_max : ", PTI.recv_first_sleep)
                end

                break

            end

        end
    end

    if ! get_through
        ErrorException("[recvText] No further incoming message within timeout.") |> throw
    end

    result = ""

    while true

        try
            open(PTI.recv_fn, "r") do io
                result = strip(read(io, String))
            end
        catch ex
            println(string(ex))
            println("Keep receiving...")
            sleep(PTI.error_sleep)
            continue
        end

        if result == ""
            println("Empty msg received. Maybe due to IO delay. Sleep and do it again.")
            sleep(PTI.error_sleep)
            continue
        end

        break
    end

    PTI.recv_trackno = (PTI.recv_trackno % PTI.rotate) + 1

    cleanHistory(PTI)

    return result
end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Omega
using Test

@test Bool(anyₛ([3.0 >ₛ 2.0, 1.0 ==ₛ 0.3]))
@test !Bool(!anyₛ([3.0 >ₛ 2.0, 1.0 ==ₛ 0.3]))
@test Bool(allₛ([3.0 >ₛ 2.0, 110.0 <ₛ 203.3]))
@test !Bool(allₛ([3.0 >ₛ 5.0, 110.0 <ₛ 203.3]))


text:
using FastSpike
using Test
# TODO: check etwork output with a calculated ground-truth.
N = 100
batch_size = 1

for network in [Network, DelayNetwork]
    @testset "$network of Homogeneous Izhikevich neurons" begin
        connection = randomConnection(N, 0.1)
        #Izhikevich neurons parameters
        a = 0.02
        b = 0.2
        c = -65
        d = 2
        # Define Network
        net = network(
            neurons=Izhikevich(dt=1, a=a, b=b, c=c, d=d, v_thresh=30.0),
            batch_size=batch_size,
        )
        neurons = add_group!(net, N; name="neurons")
        connect!(net, neurons, neurons, connection)
        run!(
            net;
            input_voltage=ones((batch_size, N)),
            input_spikes=ones(Bool, (batch_size, N))
        )

        @test true
    end

    @testset "$network of Heterogeneous Izhikevich neurons" begin
        connection = randomConnection(N, 0.1)
        size_n1 = Int(N / 2)
        size_n2 = N - size_n1
        n1 = Array{Bool}([ones(size_n1); zeros(size_n2)])
        n2 = Array{Bool}([zeros(size_n1); ones(size_n2)])
        #Izhikevich neurons parameters
        r = rand(N, 1)
        a = reshape(0.02 .* n1 + (0.02 .+ 0.08 .* r) .* n2, (1, N))
        b = reshape(0.2 .* n1 + (0.25 .- 0.05 .* r) .* n2, (1, N))
        c = reshape((-65 .+ 15 .* r .^ 2) .* n1 + (-65) .* n2, (1, N))
        d = reshape((8 .- 6 .* r .^ 2) .* n1 + 2 .* n2, (1, N))
        # Define Network
        net = network(
            neurons=Izhikevich(dt=1, a=a, b=b, c=c, d=d, v_thresh=30.0),
            batch_size=batch_size,
        )
        neurons = add_group!(net, N; name="neurons")
        connect!(net, neurons, neurons, connection)
        run!(
            net;
            input_voltage=ones((batch_size, N)),
            input_spikes=ones(Bool, (batch_size, N))
        )

        @test true
    end

    @testset "$network of LIF neurons" begin
        connection = randomConnection(N, 0.1)
        # Define Network
        net = network(
            neurons=LIF(),
            batch_size=batch_size,
        )
        neurons = add_group!(net, N; name="neurons")
        connect!(net, neurons, neurons, connection)
        run!(
            net;
            input_voltage=ones((batch_size, N)),
            input_spikes=ones(Bool, (batch_size, N))
        )

        @test true
    end
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: HairpinIncr.mso								$
#*------------------------------------------------------------------
type Hairpin_HeatTransfer
	Hairpin_HeatTransfer()=begin
		N=outers.N
		Npoints=outers.Npoints
		new(
			area(Dict{Symbol,Any}(
				:Brief=>"Cross Sectional Area for Flow",
				:Default=>0.05,
				:Lower=>1e-8
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Hydraulic Diameter of Pipe for Heat Transfer",
				:Lower=>1e-8
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of zones",
				:Default=>2
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of incremental points",
				:Default=>3
			)),
			fill(temperature(Dict{Symbol,Any}(
				:Brief=>"Incremental Local  Temperature",
				:Lower=>50
			)),(Npoints)),
			fill(positive(Dict{Symbol,Any}(
				:Brief=>"Incremental Reynolds Number",
				:Default=>100,
				:Lower=>1
			)),(N)),
			fill(heat_trans_coeff(Dict{Symbol,Any}(
				:Brief=>"Incremental Film Coefficient",
				:Default=>1,
				:Lower=>1e-12,
				:Upper=>1e6,
				:DisplayUnit=>"W/m^2/K"
			)),(N)),
			fill(fricfactor(Dict{Symbol,Any}(
				:Brief=>"Incremental Friction Factor",
				:Default=>0.05,
				:Lower=>1e-10,
				:Upper=>2000
			)),(N)),
			fill(positive(Dict{Symbol,Any}(
				:Brief=>"Incremental Nusselt Number",
				:Default=>0.5,
				:Lower=>1e-8
			)),(N)),
			fill(positive(Dict{Symbol,Any}(
				:Brief=>"Incremental Prandtl Number",
				:Default=>0.5,
				:Lower=>1e-8
			)),(N)),
			fill(positive(Dict{Symbol,Any}(
				:Brief=>"Incremental Phi Correction",
				:Default=>1,
				:Lower=>1e-3
			)),(N)),
			fill(velocity(Dict{Symbol,Any}(
				:Brief=>"Incremental Tube Velocity",
				:Lower=>1e-8
			)),(N)),
			fill(enth_mol(Dict{Symbol,Any}(
				:Brief=>"Incremental Stream Enthalpy"
			)),(Npoints)),
			[:As,:Dh,:N,:Npoints,],
			[:Tlocal,:Re,:hcoeff,:fi,:Nu,:PR,:Phi,:Vmean,:Enth,]
		)
	end
	As::area
	Dh::length
	N::DanaInteger
	Npoints::DanaInteger
	Tlocal::Array{temperature}
	Re::Array{positive}
	hcoeff::Array{heat_trans_coeff}
	fi::Array{fricfactor}
	Nu::Array{positive}
	PR::Array{positive}
	Phi::Array{positive}
	Vmean::Array{velocity}
	Enth::Array{enth_mol}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Hairpin_HeatTransfer
function atributes(in::Hairpin_HeatTransfer,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
Hairpin_HeatTransfer(_::Dict{Symbol,Any})=begin
	newModel=Hairpin_HeatTransfer()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
#------------------------------------------------------------------------------#
#                                    Module                                    #
#------------------------------------------------------------------------------#

"""
# Description

    module seHTC

SymEngine backend (se) module to the Hermite Polynomial Tensor package.

# Usage

This module is `include`d by the `HermitePolyTensor` package and is not  usually
used in isolation, or directly; rather through the  interfaces  defined  on  the
`HPT` submodule, which is loaded automatically by the package:

```julia-repl
julia> using HermitePolyTensor

julia> typeof(HPT.seHTC)
Module

```
"""
module seHTC

#------------------------------------------------------------------------------#
#                                   Imports                                    #
#------------------------------------------------------------------------------#

import SymEngine

#------------------------------------------------------------------------------#
#                                   Exports                                    #
#------------------------------------------------------------------------------#

# Auxiliary Functions
export SVHP

# Tensor Component Functions
export HTC

#------------------------------------------------------------------------------#
#                             Auxiliary Functions                              #
#------------------------------------------------------------------------------#

"""
# Description

    SVHP(
        n::Int; ax::Int = 1, st::Bool = true
    )::SymEngine.Basic

Returns a symbolic `n`-th order single-variable Hermite polynomial in  the  axis
`ax`,  axis  1,  2,  and  3  corresponding  to  the  x,  y,  and  z  directions,
respectively. A "probabilist's" polynomial is returned if `st`  is  `true`  (the
default), otherwise, a "physicist's".

This function is not used directly; rather through the interfaces defined on the
`HPT` submodule, which is loaded automatically by the package, by specifying the
SymEngine ("se") backend:

# Usage

```julia-repl
julia> using HermitePolyTensor

julia> HPT.SVHP(3, be = "se")
-3*x + x^3

julia> HPT.SVHP(4, ax = 2, be = "se")
3 - 6*y^2 + y^4

julia> HPT.seHTC.SVHP(5, ax = 3, st = false) # Direct access
120*z - 160*z^3 + 32*z^5

```
"""
function SVHP(
    n::Int; ax::Int = 1, st::Bool = true
)::SymEngine.Basic
    # Validations
    if n < 0
        throw(DomainError("n = $n outside the valid domain [0, ∞)"))
    elseif ax <= 0 || ax >= 4
        throw(DomainError("ax = $ax outside the valid domain [1, 3]"))
    end
    # Execution
    if n == 0
        return SymEngine.Basic(1)
    else
        v = SymEngine.symbols("xyz"[ax:ax])
        ξ = sqrt(v^2) # ξ: U3be
        if st
            ω = (2*SymEngine.PI)^(-1//2) * SymEngine.exp(-ξ^2/2)
        else
            ω = (2*SymEngine.PI)^(-1//2) * SymEngine.exp(-ξ^2)
        end
        r = SymEngine.diff(ω, v, n)
        r *= (-1)^n / ω
        return SymEngine.expand(r)
    end
end

#------------------------------------------------------------------------------#
#                          Tensor Component Functions                          #
#------------------------------------------------------------------------------#

"""
# Description

    HTC(
        id::String; D::Int = 2, st::Bool = true
    )::SymEngine.Basic

Returns a single symbolic *component* of a Hermite polynomial tensor.

The indices are specified by `id`, a `String` containing x's, y's,  and/or  z's;
invalid characters on `id` are ignored, while `id`  without  x's,  y's,  or  z's
represent a zero-rank  tensor  (whose  only  component  is  `1`  for  all  space
dimensions).

The Euclidean space dimensionality is specified  by  `D`.  The  tensor  rank  is
implied by the provided `id`, and is always the total  count  of  *valid*  index
characters, considering `D`.

A "probabilist's" polynomial is  returned  if  `st`  is  `true`  (the  default),
otherwise, a "physicist's".

# Usage

This function is not used directly; rather through the interfaces defined on the
`HPT` submodule, which is loaded automatically by the package, by specifying the
SymEngine ("se") backend:

```julia-repl
julia> using HermitePolyTensor

julia> HPT.seHTC.HTC("xxy")
-y + x^2*y

```
"""
function HTC(
    id::String; D::Int = 2, st::Bool = true
)::SymEngine.Basic
    # Validations
    if D <= 0 || D >= 4
        throw(DomainError("dimension D = $D outside the valid domain [1, 3]"))
    end
    # Execution
    xc = count(i -> (i == 'x'), id)
    yc = D > 1 ? count(i -> (i == 'y'), id) : 0
    zc = D > 2 ? count(i -> (i == 'z'), id) : 0
    ret = SVHP(xc, ax = 1, st = st)
    if yc != 0
        ret *= SVHP(yc, ax = 2, st = st)
    end
    if zc != 0
        ret *= SVHP(zc, ax = 3, st = st)
    end
    return SymEngine.expand(ret)
end

end # module